Introduction
An important application of analyses of melt-water infiltration in cold snow is to question of what the run-off response of cold glaciers and ice caps will be in the event of future climatic warming. This is believed to be an important source of possible future sea-level rise, but is currently poorly constrained (National Research Council, 1985). The delay between the onset of increased melt-water input and the appearance of increased run-off to the ocean depends on the timing of processes which alter the hydro-logical structure of cold permeable snow and firn to allow substantial lateral transport of water to an established run-off zone at some lower elevation. Run-off paths may develop along near-surface impermeable horizons which form over uninfiltrated snow. Nearly impermeable horizons are observed to develop in cold, permeable firn under present-day conditions, where water contents are very low, so accurate modeling of water flow at low water contents may be important for accurate predictions of future run-off. This work is focused upon the basic heat and mass-transfer relationships at the interface between infiltrating water and cold snow.
Much work has been devoted to the analysis of water flow in snow, from a theoretical perspective (Reference ColbeckColbeck, 1972, Reference Colbeck1973a, b, Reference Colbeck 1974, Reference Colbeck1975a,Reference Colbeck b,Reference Colbeck 1978), through experiments (e.g.Reference Colbeck Colbeck, 1974), and through observations (e.g.Reference Marsh and Woo Marsh and Woo, 1984a, b).Reference Illangasekare, Walter, Meier and PfefferIllangasekare and others (in press) present a fully two-dimensional finite-element model of water flow in sub-freezing heterogeneous snow; all previous models are either one-dimensional (e.g.Reference ColbeckColbeck, 1972) or quasi-two-dimensional (Reference ColbeckColbeck, 1973a; Reference Marsh and Woo Marsh and Woo, 1984b), in which information about phenomena influencing flow in the horizontal direction is included parametrically in equations explicitly describing flow in the vertical direction.
Fully two– and three-dimensional models of water flow in snow are important tools for understanding how water flows through a heterogeneous snow-pack. Snow-packs (including glacier accumulation) typically have a layered structure which is hydrologically variable over horizontal distances somewhat larger than the magnitude of spacing of vertical variations. This heterogenous layered structure influences water flow in the snow-pack and is, in turn, altered by the water flow, through equilibrium–metamorphism process when the snow is at 0 ° C, and through refreezing when the snow is initially below 0 ° C. Refreezing processes are particularly significant in the infiltration of summer melt water on polar glaciers.
Observations by the authors at Mer de Glace Agassiz, Ellesmere Island, Canada (paper in preparation by W.T. Pfeffer), indicate that melt-water infiltration to depths of about 30 cm occurs annually, with resulting changes in the structure of the snow-pack including the formation of ice layers. Annual accumulation at the latitude and elevation of the Mer de Glace Agassiz study site is approximately 40 cm/year (20 g/cm2 year water equivalent), and annual melt is approximately 2.3 g/cm2 year (Reference KoernerKoerner, 1979). The liquid-water content in the snow appears to be very low (i.e. near residual values) during infiltration, and significant structural alteration probably occurs with the first arrival of melt water. Accordingly, we believe that an accurate understanding of events at the leading edge of melt-water infiltration is vital to a realistic model of the infiltration process in cold snow. This paper addresses the details of melt-water movement at the boundary between wet and dry (cold) snow, with the goal of developing the proper numerical (e.g. finite-element) formulation for melt-water advance in the presence of freezing.
To analyze water flow in an initially sub-freezing snow matrix, one must adapt the conventional porous media formulation to include refreezing phenomena. Refreezing changes both the liquid-water content (through loss of liquid water) and the porous medium structure (through addition of mass to the ice matrix). The basic character of refreezing has been analyzed by Reference ColbeckColbeck (1975a), and his analysis has been adapted in Reference Marsh and WooMarsh and Woo’s (1984b) quasi-two-dimensional model. Briefly, Colbeck’s analysis requires that water entering sub-freezing snow must first provide enough latent heat through refreezing to raise the temperature of the snow to 0 ° C, and then satisfy the residual saturation (the minimum water content required to allow water movement over capillary forces), before water is free to move further. Colbeck constructed an expression for the rate of advance of water into sub-freezing snow (the advance of the “wetting front”), given the water influx rate, the initial temperature and porosity of the snow, and the characteristic residual saturation. Marsh and Woo used this relationship to model the formation of ice layers in cold snow.
A fully two-dimensional numerical model of water flow into cold snow must be able to determine how much water is introduced into a spatial element over any particular time step, and how much of the water introduced will freeze in that time step. However, Colbeck’s analysis includes an implicit assumption, which makes his model difficult to use in this application. The rate of wetting-front advance is proportional to the water-influx rate regardless of its magnitude; implicit is the assumption that the local heat transfer is sufficiently rapid that the freezing process can warm the snow to 0°C at any rate dictated by the water influx. At low influx rates typical of snow-pack ripening, this is a good assumption. In cases of non-uniform flux rates, however, such as in fingering of the wetting front, the local water flux could be very high for short periods. In any case, some relationship must be established describing the rate of freezing of water entering cold snow to analyze accurately changes in saturation and changes in void geometry due to accretion of ice. From a numerical point of view, the rate of freezing must be known to avoid discontinuities in the vertical distribution of water content which instantaneous freezing would effectively produce.
The purpose of the present analysis is to provide a physical basis for estimating the rate of freezing of liquid water introduced into a cold porous matrix. The analysis can be approached at any of several levels, ranging from order-of-magnitude estimates to the more detailed capillary suction and freezing model presented here. Ultimately, the estimate used in a model must be a simple one, relating the rate of liquid-water loss due to freezing (with dimensions of ML −3 T −1) to local grain-size, porosity, and initial temperature. An elaborate analysis of the grain-scale phenomena controlling freezing cannot appear in an efficient macroscopic model of water flow, but it can provide reliable constraints for a simpler relationship to be used in such a model. The analysis will also assist in designing laboratory experiments which will be conducted to study this problem.
Model of Flow at the Wetting Front
In the previous analyses of water flow in snow where initially sub-freezing conditions are considered, the snow-pack is divided into two regions: a lower, “dry” region where the snow temperature is below 0°C and no liquid water is present, and an upper “wet” region, where the snow and water temperature is at 0°C and liquid water is present. Reference Marsh and WooMarsh and Woo (1984a, b) defined two fronts related to a fingering-front geometry, but snow at any point is either “wet” or “dry” in the sense defined above. In this analysis, we define a third, intermediate, zone which consists of snow at sub-freezing temperatures surrounded by water at 0 ° C. This zone is in thermodynamic disequilibrium: liquid water is freezing on to cold snow grains, and the temperature of the grains is rising to 0°C. This intermediate (“dynamic”) zone tends to be eliminated by the thermal equilibration process between water and ice, but is sustained by the addition of new liquid water from the wet zone above. The dynamic zone must propagate into new cold snow in the dry zone as newly equilibrated snow becomes an extension of the wet zone.
The pertinent processes in each of the three zones thus defined can each be described with a separate set of equations. In the wet zone, the pertinent process is water flow in a porous medium. Darcy’s equation governing water flow in an unsaturated medium applies, modified appropriately for grain growth caused by wet-snow metamorphism:
where Κ is the hydraulic conductivity of the snow, k rw is the relative permeability of snow with respect to water, φ is the total head, θ is the volumetric water content, and τ is time. The temperature equation does not apply since the snow/water mixture is thermodynamically constrained to be at 0 ° C. We note, in agreement with Colbeck (1986), that wet snow is not truly isothermal (grain-scale temperature gradients are present which are responsible for the wet-snow metamorphism process). In the context of heat and mass transfer on scales larger than the grain scale, however, a wet snow-pack is effectively isothermal, since the small temperature gradients average to zero over the larger scales and do not provide a mode for large-scale heat transfer.
In the dry zone, no water is present, and the temperature equation alone applies, with the thermal conductivity and heat capacity having values appropriate for snow of a given bulk density:
where k s ρ s, and C s are the effective thermal conductivity, density and heat capacity, respectively, for dry snow, Τ is the temperature, and τ is time. Since no water is present, no water-flow equation applies. We neglect convective transport of water vapor. This is an important component of the total energy transfer during periods of kinetic grain growth, but can be neglected during periods of melting, when macroscopic temperature gradients are small.
In the dynamic zone, both water flow and heat conduction are acting, but not in ways identical to the other zones. Water flow is again described by Darcy’s equation, but with the addition of a sink term representing loss of liquid water due to freezing, whose value as a function of time is given independently by an analysis of the freeze-on process:
where the variables are as defined in Equation (1), with the addition of m i which is the mass of water lost to freezing per unit volume of snow. The freeze-on process also enters Darcy’s equation through the change in grain-size due to freezing, which influences the hydraulic conductivity.
Since the dynamic zone has two distinct temperature (T waterand T ice), the temperature equation does not have a conventional form like Equation (2). However, it is clear that the water temperature does not need an equation to describe it, since it is thermodynamically constrained to be at 0 ° C. The ice temperature, on the other hand, does change, in response to two energy fluxes: inter-grain conduction, and latent heat from freezing of liquid water at grain surfaces. When the temperature gradients associated with these two processes acting on a single grain are compared, it may be seen that the inter-grain conduction is relatively negligible: snow-pack temperature profiles in the Arctic are typically 0.1–1.0°C/cm (Reference Marsh and WooMarsh and Woo, 1984a; paper in preparation by W.T. Pfeffer), while the intra-granular temperature gradient in a 1 mm grain at –15°C in contact with water at 0°C is 300°C/cm. Accordingly, the term describing inter-grain conduction (∇ · k s ∇T) may be dropped from the temperature equation, leaving an equation expressing the fact that virtually all of the temperature change in a cold ice grain in contact with water is due to latent-heat release from freezing on that grain:
where the variables are defined as above, with the addition of the latent heat of fusion L. The physical picture represented by this equation is of a snow-pack composed of individual grains with negligible thermal intercommunication. The temperature change in any given grain is determined solely by the rate of freezing of liquid water on to that grain. The rate of freezing is determined by the temperature gradient in the grain and by the rate at which water is supplied.
To summarize the processes acting in the three zones: in the wet, or isothermal zone, mass (and latent heat) is transferred by water flow according to Darcy’s law and no sensible heat is transferred, since the system is thermo-dynamically constrained to be at 0 C. In the dry zone ahead of infiltrating water, energy transfer is by conduction along temperature gradients only. No water is present for mass transfer, and vapor transport is neglected. In the dynamic zone, energy transfer is by mass transfer and release of latent heat upon refreezing. Inter-grain conduction is neglected but this is small in comparison with the latent heat.
We re-emphasize here that the dynamic zone is in thermodynamic disequilibrium and that there is no supercooling of the liquid water. This is in contrast to a similar situation in frozen soils, where water and ice can be in equilibrium at sub-freezing temperatures due to pressure melting-point depression in small soil pores. Pores in snow are much too large (a representative minimum is 0.05 mm (Reference SommerfeldSommerfeld, 1983)) for pressure melting-point depression to be significant.
Equations (3) and (4) describing water flow and temperature in the dynamic zone are coupled through the expression (∂m i/∂t), the rate at which water is frozen on to snow grains. This can be evaluated a priori, given the initial temperature and some choice of snow-grain geometry. The appropriate choice of snow-grain geometry is difficult to determine, since it requires a knowledge of the three-dimensional grain geometry as well as knowledge of what the typical water pathways are through the geometry, and a way to reduce the real geometry to a tractably simple idealization. In the following section, we present models based on two simple geometries which can be constrained by existing measurements. More elaborate models can be constructed on these foundations, but require more elaborate information from stereological studies of snow. The models presented here give estimates of the rate of the freezing process which acts in the dynamic zone, and can be used to estimate the width of the dynamic zone under various circumstances. We re-emphasize that the dynamic zone is generally very small (on the order of millimeters), and probably negligible in seasonal snow-pack situations. It is important to understand this small zone correctly to model wetting fronts in numerical models, however, especially in models of polar snow-packs, where significant structural alteration in the snow may occur with the passage of the wetting front.
Numerical Models of Dynamic Zone and Results
Two idealized snow-pack geometries are investigated in this section: snow grains represented as isolated spheres and void space represented as capillary tubes in ice. The time for the ice to reach 0°C is calculated for various temperatures and various grain and void sizes in each of these models. Also shown is the special case of long continuous capillaries, where capillary suction determines the water flux through the snow. Saturation profiles can be calculated for chosen distributions of capillary tube sizes.
Snow grains as isolated spheres
The snow in the dry zone of the snow-pack is regarded as a collection of spheres of some initial temperature and size, or size distribution (chosen in these calculations from steorological data of Reference SommerfeldSommerfeld (1983)). In the dry zone, heat is conducted by inter-grain contacts (and conduction and convection in the pore space) according to Equation (2), where the effective thermal conductivity k S of the snow accounts for the combined modes of heat transfer in the ice and air (Reference Gray and MaleGray and Male, 1981). Once water enters the dry snow, however, the conduction term is dropped from the temperature equation, as described above, and heat transfer is limited to the release of latent heat of fusion. The rate at which this occurs can be calculated by solving the problem of a sphere initially of radius r 0 and temperature T 0, subjected to a Stefan boundary condition with Τ = 0°C at the ice–water interface and in the water, maintained for all time after τ = 0. This problem cannot be solved analytically, however, and either a numerical solution or a simplification of the problem is required.
The final radius of the ice sphere can be calculated simply by considering the difference in sensible heat between the initial volume of ice at T 0 and at 0 ° C (the “cold content”); the relationship between initial radius r 0 and final radius r f for initial temperature T 0 is:
Using this equation, it is found that for T 0 = –15°C (a representative near-surface temperature for High Arctic snow-packs in early spring (paper in preparation by W.T. Pfeffer), the increase in radius is only 3%. For T 0 = –30 ° C (a representative minimum temperature), the increase is approximately 6%. These are small changes in the water–ice interface position, and accordingly the Stefan boundary condition is dropped from the problem formulation, and the heat transfer is calculated by considering only the problem of a sphere of given intial radius and temperature held at 0°C at the boundary for times after τ = 0. The radius change as a function of time can then be calculated with little loss of accuracy by converting the heat flux at the sphere surface to the equivalent rate of mass accretion through the latent heat. An analytic solution exists to this problem; it is (Reference ArpaciArpaci, 1966, p. 287):
The associated rate of radius change is:
Figure 1 shows the time to equilibration for ice grains of various initial radius and temperature as calculated from Equation (6). The times are somewhat ill-defined since the growth rate diminishes as e−t , and negligibly small radius and temperature changes persist indefinitely. In calculating the equilibration times in Figure 1, the grain was considered to be equilibrated when (dr/dt)/r is less than 10−8/s. The equilibration time is substantially more sensitive to initial radius than initial temperature. This is expected, since a decrease in initial temperature acts to increase the cold content, but also increases the temperature gradient, while a decrease in initial radius decreases the cold content and increases the temperature gradient. As a general rule, grains of typical size in a snow-pack (0.01–0.2 cm in diameter) will equilibrate in 0.5 s or less.
Fig. 1. Time to thermal equilibrium for ice spheres of various initial radii and sub-freezing temperatures in contact with water at 0 ° C.
The effect of neglecting the Stefan boundary condition can be checked by comparing the calculated final radius with the precise final radius of Equation (5). The calculated final radius differs from the precise value by –0.02% for r 0 = 0.1 cm to –0.3% for r 0 = 0.02 cm.
Figure 2 shows the rate of freezing of water on to sphere populations for a given initial radius, temperature, and porosity. All five curves are for an initial temperature of –15 ° C, and an initial porosity of 0.6 (corresponding to a density of 0.366 g cm −3); the individual curves are for initial radii of 0.025, 0.0375, 0.05, 0.0625, and 0.0875 cm. For a given initial radius, the density is matched by choosing the appropriate number density of spheres; the freeze-on rate (g s−1 cm−3) is calculated by multiplying the rate for a single sphere by the number density. The initial freeze-on rate for all initial radii lies between 0.25 and 0.6 g s−1 cm−3, and decreases to less than 0.1 g s−1 cm−3 after 0.1 s. The total freeze-on change corresponds to a density change of 0.0343 g cm−3, or 0.094% of the initial density.
Fig. 2. Rate of freezing of water on to sub-freezing sphere populations of various initial radii. The initial temperature is –15 ° C in all cases. The rate for an individual sphere (MT−1) is converted to a rate per unit volume of snow (MT−1L−3) by multiplying by the number density of grains. The number density for each grain-size is chosen to maintain a fixed bulk porosity of 0.6.
Capillary tube model
Water pathways in the snow-pack are treated as capillary tubes of specified radius in cylinders of ice of specified wall thickness and initial temperature. A representative capillary radius (or distribution of radii) can be taken from pore-size measurements, while the wall thickness (or distribution of thicknesses) represents a typical grain radius. A more particular interpretation of the wall thickness can be made by determining a typical coordination number between grains and water paths, and assigning a wall thickness of equal area (in cross-section) to the fractional parts of grains contacting a single water path. This was not done here, however, and the wall thickness is related directly to an average grain radius.
Treating water pathways in the snow-pack as capillaries of constant radius is a substantial idealization of the true situation, for the porous part of a real snow-pack is composed of voids of many sizes linked together in sections whose lengths are roughly equal to the grain-size. By giving the water pathways in the snow a definite and simple geometry, however, the water flow can be calculated as well as the rate of freeze–on. The standard expression for Poiseuille flow in a capillary is used, with a simple modification to account for the freezing-on process on the inner wall of the capillary. (Poiseuille flow does not strictly apply to flow in this situation, since the freezing-on process with density change between phases violates the boundary condition that there be no radial component of fluid velocity at the wall. The radial component of flow at the wall is very small compared to the axial component, however, and this complication is ignored in this analysis.) Neglecting the gravity component of the total head, the standard expression is (Reference CoreyCorey, 1977);
where u w is the water velocity (and the capillary front velocity), R is the capillary tube diameter, r is the diameter at the capillary front, x w is the distance from the tube entrance to the capillary front, h is the standing head at the tube entrance, and σ and μ are the specific surface energy (relative to air) and viscosity of water at 0°C. In all the calculations which follow, h is set to 0. The water pressure at the top of the capillary is thus at 0 (atmospheric) pressure, and whatever input flux is demanded by the capillary-front advance is supplied at this pressure. Ordinarily, R and r are the same but, in the case of freezing-on, r is the initial radius of the capillary (at the leading edge of the water), while R has been reduced by freezing-on, and varies with distance along the wetted part of the capillary. The change in R is accommodated here by taking the average radius along the wetted part of the capillary at any given time. Since the freeze-on process reduced the capillary radius rapidly at first and quickly stabilizes to a steady value, the average radius R diminishes rapidly when water first enters the capillary, and stabilizes asymptotically to a steady value as the equilibrated part of the capillary occupies progressively greater fractions of the total wetted length.
Gravity forces have been neglected in the formulation of Equation (8), although these can be added simply by including a gravity-head term along with (2σ/r) and h. Gravity is neglected in this analysis, however, to isolate the effects of capillarity, and also since unbroken capillaries long enough for gravity to exceed capillarity are unlikely to occur in a real snow-pack (gravitational forces exceed capillary forces in a tube of radius r = 0.025 cm only for capillaries longer than 6.2 cm).
Equation (8) is integrated to give the position of the wetting front as a function of time:
where 〈R〉 is the average radius along the wetted length of the capillary.
Heat transfer due to refreezing is calculated in the radial direction; heat transfer in the axial direction is neglected. As discussed in the description of the dynamic zone, temperature gradients arising from conduction in the ice are much smaller than those between ice and freezing water, and are consequently neglected. As for the spherical problem, no analytic solution exists for the problem of freezing in a capillary when the Stefan boundary condition exists. The problem can be solved by neglecting the Stefan boundary condition, as is done for the sphere problem, but solutions (e.g. Reference Carslaw and JaegerCarslaw and Jaeger, 1959, sect. 13.4) require extensive numerical work to obtain roots to Bessel functions. Since a numerical treatment is required in any case, a finite-element solution was made here which includes the Stefan boundary condition. The finite-element analysis is for radial heat flow in an isotropic medium, bounded at small r by the Stefan condition at the interface with a water-filled hole of non-zero radius, and bounded at large r by the condition that ∂T/∂r = 0. The solution procedure is an adaptation by N. Humphrey (personal communication) of the method described by Reference Jarvis and ClarkeJarvis and Clarke (1974). The numerical procedure solves the problem of an infinite thick-walled cylinder of given initial temperature, radius, and wall thickness, subject to boundary conditions of fixed temperature (0°C) at the inner surface, and no heat flux across the outer surface. The outer boundary condition is imposed to approximate the symmetry of heat flow in a snow grain in contact with water on all sides. The outer wall of the capillary represents the centers of several snow grains surrounding a void, and the flow of water in the capillary represents water flowing in the void.
The interaction between the freeze-on process and the capillary flow is simulated by advancing the capillary front a distance (δx = u wδt), where u w is determined by Equation (8), and (δt) is 0.002 s in these calculations. After advancing the capillary front to its new position, each previous (δx) has its radius updated according to the heat-transfer solution, which provides the current radius and wall-temperature distribution for times after the onset of freezing. The average radius 〈R〉 is then recalculated and the front advanced again. The simulation gives, at each time step, the position and velocity of the capillary front, the capillary radius as a function of position behind the front, and the wall-temperature distribution at positions behind the capillary front. At any given position behind the capillary (or wetting) front, the rise in temperature at the insulated outer wall to 0 ° C corresponds to the thermal equilibration of the capillary at that position, or equivalently, the arrival of the “melting front” — the boundary between the dynamic zone and the wet zone, as defined above.
Advance of wetting and melting fronts
Figure 3 shows wetting- and melting-front positions, and average capillary radius, as a function of time after introduction of water into the capillary. Figure 3a is for an initial capillary radius of 0.025 cm, and Figure 3b is for an initial radius of 0.05 cm. For both capillaries, the wall thickness is 0.05 cm and the initial temperature is –15 ° C. A comparison of the wetting-front advance in the two capillary sizes shows that the larger capillary moves water further in during any given interval of time than the smaller capillary, and at a rate which diminishes more slowly. This is expected, since the water velocity (Equation (8)) increases linearly with r (for h = 0), but also occurs because the freeze-on process has a greater proportional effect on the radius at smaller initial radii, and thus reduces 〈R〉 in Equation (9) more significantly for smaller radii. This is especially true at times shortly after the introduction of water into the capillary, when the rapidly changing zone near the capillary front occupies a large part of the total wetted length. At longer times, when most of the wetted length of the capillary is equilibrated, the velocity diminishes due to the increasing length of the capillary only.
Fig. 3. Position (depth in capillary) of welting front and melting front as a function of time following introduction of water, and accompanying change in average radius in wetter length of capillary. Initial température = –15° C; wall thickness = 0.05 cm. (a) for r0 = 0.025 cm; (b) for r0 = 0.05 cm.
Fig. 4. Lag time between introduction of water into capillary, and thermal equilibration at top of capillary (arrival of melting front), for various values of initial radius (r0) to wall thickness (D). Initial temperature = –15°C.
A comparison of the melting-front arrival times for the different capillary sizes shows that the larger radius equilibrates more rapidly (the melting front is initiated at the entrance to the capillary sooner). This can also be seen in the radius change versus time, where the radius change flattens out earlier for larger initial radius. This is because the heat transfer occurs more rapidly for larger radii, where the surface area of the inner boundary is larger. Once the melting front is initiated in the capillary, it propagates at the same speed as the wetting front did at the same position since the time lag between the two is a fixed interval for a given initial radius and temperature. This lag time is also a function of initial radius, wall thickness, and temperature. The lag time is shown in Figure 4 as a function of r 0 for T 0 = –15°C.
Front advance as a function of capillary size
Figure 5 shows the position of the wetting and melting fronts after 2.5 s for ten different initial radii between 0.005 and 0.10 cm. The relationship between front positions and initial radius is shown here over a range of capillary sizes, and gives a measure of the gap between the wetting and melting fronts at a fixed time for the range of initial radii, as well as an average velocity of front propagation over 2.5 s. Also shown in Figure 5 is the minimum initial radius (0.0207 cm) which will equilibrate without freezing shut at T 0 = –15 ° C and an initial wall thickness of 0.05 cm. The two capillaries with initial radii below this limit (0.005 and 0.0125 cm) froze shut after 0.03 and 0.21s, respectively. This result is significant, since it shows that there is a limiting ratio of initial radius to wall thickness ro/D for any given sub-freezing initial temperature below which capillaries will not contribute to the transmission of water. This limit can be calculated easily from Equation (5), and is shown in Figure 6 for a range of initial temperatures between 0 and –40 ° C. Using Figure 6 to find that the minimum r 0/D is 0.414 at Τ = –15 ° C, one can calculate that for cold, coarse-grained snow with D (interpreted as half the average grain-size) equal to 0.1 cm, capillaries below 0.04 cm will freeze shut, while for moderately fine-grained snow (D = 0.025 cm) capillaries below 0.01 cm will freeze shut.
Fig. 5. Position of wetting and melting fronts at 2.5 s after introduction of water into capillary, for values of initial radius between 0.0 and 0.1 cm. Initial temperature = —15 ° C. Crosses at upper left denote furthest extent of water in capillaries which froze shut before 2.5 s.
Fig. 6. Minimum value of the ratio of initial radius (r0) to wall thickness (D) which will not freeze shut before reaching thermal equilibrium, as a function of initial temperature.
Water flow and saturation profiles in a collection of capillaries
When water is introduced into a collection of capillaries of different sizes, the water will move at different rates in the different capillaries and, if the position in the capillaries can be related to depth in the medium, saturation profiles leading up to the wetting front can be calculated as a function of time. In the following calculation, we assume the porous snow-pack to consist of a bundle of parallel
capillaries. Since gravity forces are not included, the length along the capillary bundle should not be interpreted as depth more than any other direction. Straight capillaries can be expected to produce anomalously high wetting-front speeds, since the capillary paths in a real snow-pack are tortuous and interrupted by voids. Nonetheless, the results give insight into the behavior of flow through distributions of capillaries, and the approach can be adapted to a more realistic model of capillaries linking voids by using more detailed information on the geometry of void space in a snow-pack.
Reference SommerfeldSommerfeld (1983, fig. 7) gave frequency distributions of grain lengths and void lengths in new snow. Both the grain- and void-length distributions are very strongly peaked (33%) at a minimum size of 0.005 cm, and fall rapidly and monotonically to minimum frequency at maximum sizes of 0.053 cm for grains and 0.13 cm for voids. His void-length distribution was used here to synthesize a population of capillaries with initial radii in eight size categories between 0.025 and 0.10 cm with the same distribution statistics (Table I). Even though capillary sizes below 0.02 cm are the most numerous in the frequency distribution, they were not included in the synthesized distribution since they freeze shut quickly (as seen in Figure 5). The wall size was kept at 0.05 cm, despite the fact that the grain-size corresponding to this wall thickness is too large for the measured grain-size distribution. This was done to allow comparison with the other calculations presented here, and to restrict the number of parameters being varied. It is difficult in any case to choose a single wall thickness which will be representative of the entire range of observed grain-sizes. A more complete model will combine a range of capillary sizes with a range of wall thicknesses.
Fig. 7. Rate of freezing of water on to inner wall of capillaries of four initial radii (as marked), per unit length of capillary, as a function of time following introduction of water. Initial temperature is –15 ° C.
Table. I. Frequency Distribution Of Voids In Snow, Adapted From Sommerfeld (1983)
The freeze-on rate per unit length of capillary for four of the ten capillary sizes considered is shown in Figure 7. The freeze-on rate for a collection of capillaries is found by weighting the individual capillary results by their respective relative frequencies, and multiplying the result by the total number density of capillaries in snow of a chosen bulk porosity (0.6, as in the corresponding example for spheres). The result for the collection of capillaries is shown in Figure 8. The rate of freezing-on is less overall in capillaries than in spheres of the same initial porosity and temperature, but persists for slightly longer. This is consistent with the fact that the ratio of volume-to-contact area with water is greater for spheres than for capillaries, so an equivalent cold content in the ice must be transferred to the water at a slower rate in the capillary geometry. The difference in total freeze-on rate between the sphere model and capillary model also gives information on the extent to which one needs to characterize the grain and void geometry accurately from measurements. The initial rates of freezing-on (e.g. at τ = 0.05 s) differ by a factor of roughly 20, depending upon which initial radius for the sphere model is compared to the capillary model (note that a distribution of sphere radii could be chosen). The total time for the freeze-on process, on the other hand, varies at most by a factor of six between all the geometries considered, with the bulk of the process occurring over times varying by a factor of two. In cases where the time variation of the freeze-on rate needs to be known accurately, the grain/water-channel geometry must be described accurately. If only the total time for freeze-on and the average rate need to be known, an accurate characterization of the geometry is not so critical.
Fig. 8. Rate of freezing of water per unit volume in a collection of capillaries as a function of time following introduction of water. The size distribution and total volume of capillaries is described in the text. Initial temperature is —15 ° C.
By weighting the individual capillary results by the product of their respective relative frequencies and total void-area fractions, the total flux through the capillary collection can be calculated. Figure 9 shows saturation profiles in the hypothetical snow-pack described above whose void space consists of parallel capillaries. The size distribution of capillaries is the same as that used in the earlier calculations. The profile is diffusive, with the distance between the wetting front and the fully saturated front increasing with time, while the profile overall moves more slowly with increasing time. The latter fact is expected, since the pressure gradient driving the flow diminishes as the wetted capillary length increases. The diffusive character of the profile is due to the fact that the larger capillaries in the population move water faster than the smaller capillaries, and hence faster than the aggregate of all capillaries, which carry the fully saturated flow. The propagation speed of the profile is anomalously fast compared to observed rates in snow-packs, owing to the assumption of parallel straight capillaries. This can be corrected by using a model of linked capillaries and voids, but requires more extensive data on snow-pack geometry.
Fig. 9. Saturalion profiles in collection of capillaries as a function of time following introduction of melt waler at top. The size distribution and total volume of capillaries is described in the text. Initial temperature is —15°C.
Conclusions
We have developed the notion of a dynamic zone separating dry, sub-freezing snow from wet “isothermal” snow, characterized by thermal disequilibrium between liquid water and sub-freezing ice, and by a rate of freezing of water controlled by the ice/water geometry and initial temperature. Rates of freezing of water on to sub-freezing ice grains are generally very high, with equilibration times for typically sized grains being of the order of 0.5 s. The time to equilibration, and especially the rate of freezing-on during equilibration, is sensitive to the assumed ice/water geometry.
Two ice/water geometries were investigated: isolated spheres of ice surrounded by water, and capillary tubes of ice into which water is pulled by capillary suction. Isolated spheres are relatively easy to parameterize from observational data: initial temperature, a distribution of radii, and bulk number density of spheres must be given. To what extent spheres adequately represent snow grains will depend on the situation being modeled. The rate of freezing on to isolated spheres can be calculated accurately, and the rate included as liquid-water loss in the Darcy flow equation. The Darcy equation requires knowledge of the relative permeability k rw, however, which is in general poorly established for snow and is especially uncertain at the low and rapidly changing water content in the dynamic zone. Furthermore, it is questionable whether Darcy’s equation can be used to describe water flow in a situation where it may happen that water content is so low that not all the potential water pathways are being used.
The capillary model provides a local mechanism for moving water (capillary suction) so that Darcy flow need not be used. However, this characterization requires more snow-pack data. Initial temperature, grain-size, and void-size distributions must be known, as well as information on whether all sizes of grains are uniformly in contact with all sizes of voids, and how long and tortuous the capillary passages are.
The sphere model reaches thermal equilibrium at shorter times than the capillary model (cf. Fig. 1 to Fig. 4 and Fig. 2 to Figs 7 and 8). It is difficult to compare the two models precisely, however, since the parameterizations are different, and precisely the same geometrical situation cannot be put into both models. Translating equilibrium times into a width for the dynamic zone requires knowledge of the water-input velocity from the wet zone; this was partially addressed by the capillary-front model, where the water velocity is an intrinsic property of the model, and which produces a separation as a function of time between the capillary front (wetting front) and the melting front, where thermal equilibrium has developed in the capillary (Fig. 5). The capillary model can be carried further to calculate saturation profiles as a function of time (Fig. 9). The calculations based on the capillary model are very hypothetical, however (especially the saturation profiles), and serve primarily to give a qualitative picture of the behavior of the capillary geometry under freezing conditions. More snow-pack information is needed to constrain properly the several parameters required for the capillary model, and indeed to guide the choice of a capillary model as opposed to a sphere model, or some other model not considered here. Estimates have been provided here for the timing of the freeze-on process in the dynamic zone, but more detailed modeling waits experimental information on the actual geometry of water flow at the wetting front.
Acknowledgements
We wish to thank N. Humphrey for supplying the finite-element program used in calculating the freezing in a capillary, and R.A. Sommerfeld for supplying the void-frequency distribution data. The work was improved by discussions with S. Colbeck, and with W. Briggs (Applied Mathematics Program, University of Colorado at Denver). This work was supported by grant no. DE-FGO2–87ER60570 from the CO2 Research Division of the U.S. Department of Energy.
