Motion of Sub-Freezing Ice Past Particles, with Applications to Wire Regelation and Frozen Soils

Abstract

Existence of a very thin layer of adsorbed water adjacent to particles embedded in ice allows relative motion between ice and particles even at sub-freezing temperatures if there are either applied stresses or macroscopic temperature gradients. Theoretical analysis of such motion involving a single sphere demonstrates that such motion is dominantly due either to “viscous” deformation in the ice or to mass transport in the liquid layer at temperatures below the nominal pressure melting-point, depending on the ratio of the sphere’s radius to a temperature-dependent “transition radius”. This result should also hold for motion of a cylinder (for which the creeping flow problem has no known analytical solution). Reviewing data on wire regelation at sub-freezing temperatures in the context of this analysis suggests that all “anomalous” data correspond to cases in which wire radii were greater than the transition radius, leading to dominance of ice-deformation effects. Ice motion past very small particles, on the other hand, is essentially accommodated entirely by mass transfer through the liquid layer. This result lends support to the “rigid-ice” model of frost heaving as proposed by R.D. Miller and co-workers, and permits approximate analysis of ice movement through a porous soil. In all cases involving relative motion between ice and particles at sub-freezing temperatures, the existence of macroscopic temperature gradients plays an important role.

Résumé

L’existence d’une très fine couche d’eau adsorbée sur des particules enchassée dans la glace permet un mouvement relatif entre la glace et les particules même à des températures inférieures au point de regel lorsqu’existe une contrainte ou un gradient macroscopique de température. Une analyse théorique d’un tel mouvement dans le cas d’une seule sphère montre que cela est principalement dû au transport de masse dans la couche liquide à des températures inférieures au point de fusion correspondant à la pression effective, en fonction du rapport du rayon de la sphère et d’une température associée au ”rayon de transition”. Ce résultat peut aussi être établi pour le mouvement d’un cylindre (pour lequel on ne connait pas de solution analytique). Une revue des données de regel sur les fils à des températures inférieures au point de congélation, dans le contexte de cette analyse suggère que tous les résultats ”anormaux” correspondent aux cas où les rayons étaient supérieurs au rayon de transition en privilégiant les effets de la déformation de la glace.

Zusammenfassung

Das Vorhandensein einer sehr dünnen Schicht adsorbierten Wassers an Partikeln, die in Eis eingelagert sind, lässt eine relative Bewegung zwischen dem Eis und den Partikeln auch bei Temperaturen unter dem Gefrierpunkt zu, wenn entweder Druckkräfte oder makroskopische Temperaturgradienten vorhanden sind. Die theoretische Analyse solcher Bewegung, beschränkt auf eine einzelne Kugel, zeigt, dass sie in erster Linie entweder auf “viskoser” Deformation im Eis oder auf Massentransport in der flüssigen Schicht bei Temperaturen unter dem nominellen Druckschmelzpunkt beruht, in Abhängigkeit vom Verhältnis des Kugelradius zu einem temperaturabhängigen “Übergangsradius”. Dieses Ergebnis sollte auch für die Bewegung eines Zylinders gelten (für die das Problem des Kriechfliessens keine analytische Lösung besitzt). Die Betrachtung von Daten zum Einfrieren von Drähten bei tiefen Temperaturen im Rahmen dieser Analyse lässt vermuten, dass alle “ungewöhnlichen” Daten solchen Fällen zugehören, bei denen die Drahtradien grösser als der Übergangsradius sind, was zum Überwiegen von Effekten der Eisdeformation führt. Andrerseits wird die Eisbewegung über sehr kleine Partikel im wesentlichen ganz durch Massentransport durch die flüssige Schicht bewirkt. Dieses Ergebnis stützt das Modell des “starren Eises” für Frosthebungen, vorgeschlagen von R.D. Miller und seinen Mitarbeitern, und gestattet eine Näherungsanalyse der Eisbewegung durch porösen Boden. In allen Fällen von relativer Bewegung zwischen Eis und Partikeln bei Temperaturen unter dem Gefrierpunkt spielt das Vorhandensein makroskopischer Temperaturgradienten eine wichtige Rolle.

1. Introduction

The long-standing notion that glaciers cannot slide over their beds at sub-freezing temperatures (i.e. at temperatures below the pressure melting-point) has recently been challenged. Shreve (1984) modified the Nye–Kamb version of glacier-sliding theory (Nye, 1969, 1970; Kamb, 1970) to account for the existence of a liquid-like layer at the ice–rock interface at sub-freezing temperatures, and showed that sliding rates, albeit very small, could lead to substantial total displacement over a period of many years. Field studies by Echelmeyer and Zhongxiang (in press) and by Hallet and others (in press) provide evidence of glacier sliding where basal ice is at least locally at sub-freezing temperatures, although the exact mechanisms of sliding in these cases remain poorly understood.

However, the concept of ice “sliding” past obstacles at sub-freezing temperatures is not new in the literature on frozen soils, where it has been promoted since at least 1972 by R.D. Miller and co-workers. The most common mode of ice-lens formation, known as “secondary heaving” (Miller, 1972), is thought by Miller and co-workers to involve actual motion of sub-freezing pore ice towards the ice lens, this motion being rendered possible by the existence of very thin films of “adsorbed” water at ice–mineral grain interfaces. The exposition of the concept of pore-ice motion as given by O’Neill and Miller (1985) is so lucid that it is well worth quoting them at some length. Considering first a single grain embedded in ice,

“at temperatures somewhat below 0°C the grain ought to be surrounded by a film of unfrozen liquid in equilibrium with the ice. If a temperature gradient is imposed, thermal equilibrium of water and ice at the interface is inconsistent with mechanical equilibrium in the hydrostatic field induced by surface adsorption forces. Whereas the thermal gradient induces asymmetry of film thickness, the action of adsorption forces is to center the grain within its liquid shell. Thus the temperature field constantly acts to diminish the film thickness on the cold side, while surface forces seek to retain that film thickness by removing ice from the warm side and transporting the resulting unfrozen water to the cold side, where it refreezes. The grain ought to migrate up the temperature gradient and its velocity should increase as it moves into an ever warmer environment with a corresponding increase in average thickness of the film, expediting the flow of liquid by which the centering tendency is expressed” (O’Neill and Miller, 1985, p. 283).

The predicted migration of grains through ice was observed by Hoekstra and Miller (1967) and by Römkens and Miller (1973). As a corollary, O’Neill and Miller (1985, p. 283) proposed that:

“If individual grains migrate through stationary rigid ice, traveling up a temperature gradient, then rigid ice that largely fills interstices between stationary grains ought to migrate down a temperature gradient. If the ice is inherently rigid, this movement is not flow but continuous regelation. Crystalline ice, everywhere bounded by liquid in both adsorption and capillary space, is continuously melting and reforming in a manner consistent with the geometry imposed by the array of soil grains.”

This scenario of “thermally induced regelation”, which must be distinguished from the better–known phenomenon of pressure-induced regelation, is illustrated as well by Figure 1. It should be noted that the adsorbed liquid films are very much thinner (only a few nm at –l°C)than the liquid layer in the case of pressure-induced regelation (cf. Nye, 1967; Gilpin, 1979, 1980[c]).

Fig. 1. Thermally induced regelation of ice through the pore space of an idealized soil. Ice moves in the direction of lower temperature (after O’Neill and Miller, 1985).

O’Neill and Miller’s description of the pore ice as “rigid” deserves further examination. By “rigid”, O’Neill and Miller meant that plastic deformation of the pore ice is of negligible importance in secondary heaving, and indeed they drew on rheological data to develop a semi-quantitative argument (p. 285) showing that plastic deformation should be negligible for silty soils (grain-size c. a few tens of μ m), which are known to be highly susceptible to frost heaving. In problems of glaciological interest involving mixtures of mineral grains and sub-freezing ice – say, the mechanics of “cold-based” glaciers resting on permafrost – the assumption of “rigid” ice may fail. We may anticipate from Shreve’s (1984) results – and will later demonstrate – that “thermally induced regelation” through a porous material with a typical grain-size of more than c. 100 μm will involve an important element of plastic deformation as well.

In section 2, we develop a theoretical model of ice flow past a single sphere at sub-freezing temperatures. This development draws upon the analysis of the motion of temperate ice past a sphere (Watts, unpublished) and Gilpin’s (1979) model of the liquid layer at the interface between sub-freezing ice and an embedded particle. This study, besides complementing Shreve’s (1984) modified sliding theory, provides insight into the results of some experiments on wire regelation and also lays the basis for an approximate analysis of ice movement through porous materials, presented in section 3.

2. Analysis

Ice flow

The passage of a particle through sub-freezing ice is conceptually similar to the more familiar temperate-ice situation. Ice mass may be re-distributed either by plastic deformation or by a melting–refreezing process made possible by a thin, continuous liquid layer at the particle-ice interface. ^* . A major difference from the temperate-ice situation arises because there may be macroscopic temperature gradients within the ice. In the absence of forces applied to the particle, it is in fact such temperature gradients that cause particle motion. Furthermore, explicit analyses of flow in the liquid layer is required (cf. Shreve, 1984).

The first part of this analysis directly parallels that of Watts (unpublished, p. 25-27). Figure 2 shows the geometry to be considered. For convenience, the sphere of radius r is considered fixed at the origin of coordinates, the usual spherical coordinates r, θ, and ψ being used. The ice flows steadily past the sphere with velocity V, the far-field velocity being of magnitude U in the z direction. Stresses σ _ij are measured relative to the no-flow configuration, with the convention that tensile stresses are positive; pressure p in the ice is defined by –p = σ _kk /3 (using the summation Figure 2 convention for subscripts). The ice is assumed to be an incompressible fluid with Newtonian viscosity η _i. lnertial effects are so small that the equation of motion reduces to

(1)

which combined with the incompressibility condition

(2)

implies that ∇² p = 0. Symmetry conditions, along with the constraint that p → 0 as r → _∞, lead to a solution for the pressure

(3)

Fig. 2. Flow of sub-freezing ice past a sphere. Standard spherical coordinates r (radial distance), θ (polar angle), and ψ (azimuthal angle) are used. The liquid-layer thickness, greatly exaggerated here, would actually be non-uniform.

where A is a constant to be determined. Equation (3) may now be substituted into Equation (1) to solve for V, and standard relationships from fluid mechanics (e.g.Bird and others, 1960, p. 90) used to find the stresses. The aforementioned far-field condition on flow is equivalent to the boundary condition

(4)

Furthermore, the liquid layer at the ice–sphere interface lubricates that boundary, hence

(5)

Equation (5) is sufficiently accurate as long as

(6)

where h is the liquid-layer thickness. With these boundary conditions, the velocities and stresses are readily found to

(7a)

(7b)

(7c)

(7d)

with all other σ _ij vanishing.

If the ice were temperate, the constant A would be determined (cf. Nye, 1967; Watts, unpublished) by solving for the temperature distribution and imposing the constraint that the temperature and normal stress σ _rr everywhere on r = R must satisfy the pressure-melting relationship ∝σ _rr ,where is measured relative to 0°C. In the present case of sub-freezing ice, and σ _rr are not so simply related. Determination of A requires that we solve not only for the temperature field but also for the liquid-layer thickness h everywhere at the sphere–ice interface, h and being functionally related (cf. Gilpin, 1979, p. 239). We will find two expressions for h and fix the constant A by requiring that these two expressions be equivalent.

Flow within the liquid layer

We will now adopt the model of the liquid layer at sub-freezing temperatures proposed by Gilpin (1979) and applied by him to problems of wire regelation, particle migration, and ice lensing (Gilpin, 1979, 1980[a], [b], [c]). Gilpin’s fundamental assumption – one motivated by a variety of experimental phenomena – is that the chemical potential of water, but not of ice, is lowered in close proximity to a solid surface, the thermodynamic effect being given by (Gilpin, 1979, p. 238)

(8)

where the chemical potential μ _w of water is given by the difference between the chemical potential μ _wB of bulk water and the change in chemical potential μ _wσ due to the surface. For mathematical convenience, Gilpin assumed μ _wσ to be given by

(9)

where a and α are constants and y is the distance measured normal to the surface and towards the ice. Equation (9) is valid only for y > y ₀ where “y ₀ is of the order of a few molecular dimensions” (Gilpin, 1979, p. 238).

The chemical attraction of water to the surface will be manifested in part by an increase in water pressure near the surface (Gilpin, 1979, p. 239; Equation (5)):

(10)

where v _w is the specific volume of water and _w is the water pressure relative to some datum, here taken as the no-motion, no-temperature-gradient state. Equilibrium at the ice–water phase boundary leads to the additional condition (Gilpin, 1979, p. 239; Equation (11)):

(11)

where v _i = specific volume of ice, Δv = v _i –v _w , σ_iw = surface tension of ice–water interface, K = mean curvature of ice–water interface, L = latent heat of fusion, T _a = absolute temperature (K), and _wh = water pressure at ice–water interface.

Equation (11) shows that the “melting” temperature depends not only on pressure but also on liquid-layer thickness. Variations in liquid-layer thickness, hence water pressure, lead to flow along the liquid layer. Mass conservation requires that such flow be balanced by melting or freezing at the phase boundary. This conservation relationship may be readily put into mathematical form.

Let q _W (θ) be the "upward" mass-flow rate through the liquid layer at the polar angle θ (cf. Fig. 2). Assuming that the water behaves in a linearly viscous fashion and that Equations (6) are valid, we may write, in analogy to Gilpin (1979, p. 240; Equation (17)):

(12)

As long as the liquid-layer thickness is constant in time, q _w (θ) must be balanced by ice influx q _i (θ), given by the surface integral

(13)

where n is the unit outward normal from the phase boundary and the integral is taken over that part of the surface spanned by polar angles θ to π. Substituting Equation (7a) into Equation (13) and integrating, then equating the fluxes q _i(θ) and q _w(θ), we find after some rearrangement

(14)

We may evaluate the integral in Equation (14) by recognizing that Equation (10) may be re-cast as

(15)

and using Equation (11) to eliminate _wh, For h >> y ₀ , we find

(16)

where the contribution of the term dK /dθ may be shown to be negligible because of the assumption that h varies little with θ (Equation (6)).

Temperature–stress relationship at the phase boundary

We now derive the expression that replaces the pressure-melting condition used in Watts’ (unpublished) analysis of a sphere moving through temperate ice. Equation (11) may be re-written as

(17)

where _i, the “pressure” in the ice normal to the phase boundary, is measured relative to the no-motion, no-temperature-gradient state. When Equation (6) holds, _i = – σ _rr (R); hence, using Equation (7d), Equation (17) becomes after differentiation and re-arrangement

(18)

Equations (16) and (18) constitute two equations for the three unknowns h, >, and A. We now use energy-conservation considerations to find a third expression.

Temperatures at the ice–sphere interface

The temperature distribution at the ice–sphere interface is readily found by adapting the analysis of Gilpin (1979, p. 241-42) to the problem at hand. The essential difference between our formulation and Gilpin’s is that he implicitly used the rigid-ice model; with the present model the normal component of ice velocity at the interface is of magnitude (cf. Equation (7a)) (U + A/R) cos θ instead of simply Ucos θ. The temperature is therefore (cf. Gilpin, 1979, p. 242)

(19)

where _u= “undisturbed” temperature at center of sphere, k _i ,k _s = thermal conductivities of ice and sphere, respectively, and G _T = imposed temperature gradient dT/dz in ice far from the sphere (cf. Fig. 2). Equation (19) is valid if convected heat, as well as work done by the sphere’s motion, are of negligible importance in comparison to conduction, as is almost always true (Philip, 1980, p.195-98), and if the temperature drop across the liquid layer is negligible. This latter condition requires (Nye, 1967) that h/R « k _w /k _s , where k _w is the thermal conductivity of water; this condition is nearly always met as long as we restrict our attention, and justifiably so, to geological materials. A further assumption in deriving Equation (19) is that there is no net heat flow either towards or away from the sphere. This last assumption requires some further discussion.

Drake and Shreve (1973, p. 66) found in their wire-regelation experiments with ice at 0°C that, at sufficiently large driving forces, the wires left behind themselves a trace of water and vapor, the volume of which could be explained only if there had been a net flow of heat toward the wires. Formation of the trace was associated with a transition from a “slow” to a “fast” mode of wire motion. A similar sort of transition was observed in wire-regelation experiments using ice at temperatures below 0°C (Telford and Turner, 1963; Gilpin, 1980[c]), although these authors did not discuss any water trace behind the wires. The exact reasons for the transition remain uncertain (Drake and Shreve, 1973, p. 69; Gilpin, 1980[c], p. 446-47). It seems plausible that such a transition might also occur for a sphere moving through ice, although there are no experimental data to test this hypothesis properly. Our analysis should therefore be considered restricted to the “slow” mode of motion.

We now substitute Equation (19) into Equations (16) and (18) to eliminate , finding after some re-arrangement

(20)

(21)

Liquid-layer thickness

Before proceeding to the solution of Equations (20) and (21), it is useful to recast these equations into dimensionless form. Following Gilpin (1979, p. 292), the liquid-layer thickness will be re-expressed as

(22)

where h _c is the equilibrium thickness for a stationary sphere, given by

(23)

and the dimensionless deviation h’ from this thickness must be much less than one, as implied by Equation (6b). The characteristic temperature will be taken as (Gilpin, 1979 n 242)

(24)

We may now write three characteristic velocities that arise out of Equations (20) and (21), two of them identical to those defined by Gilpin (1979, p. 242), viz.:

(25)

and the third, new to our analysis:

(26)

Gilpin (1979) pointed out that the velocities V _cη and V _ck are characteristic of the regimes in which particle motion is limited by, respectively, flow through the liquid layer or heat conduction through and around the particle. The “new” velocity V _cc , which characterizes the regime in which ice deformation limits particle motion, arises because we now explicitly consider the ice to be deformable.

We finally introduce a characteristic temperature gradient, again following Gilpin (1979, p. 242):

(27)

We may substitute Equations (22) – (27) into Equations (20) and (21), and integrate to find two expressions for h’ (θ), the (dimensionless) variation in liquid-layer thickness. Because we are looking only at “slow-flow” cases for which h’ « 1, certain approximations may be made in performing the integrals. The lowest-order solutions, derived in Appendix A, are

(28)

and

(29)

Equating these last two expressions, we find that

(30)

where R _v = V _c η/V _ck (Gilpin’s definition) and we introduce the “new” ratio R _c = V _c η/V _cc .R _v is almost always much less than unity (Gilpin, 1979, p. 245-46), so we may neglect terms in R _v in Equation (32). R _c = (l/2)(v _i /v _w)² (η _i/η _w)(h _c/R)³ may be larger or smaller than unity depending on the value of sphere radius R.

We finally need an expression relating the ice-flow rate to the force F exerted on the sphere. Any applied force must be balanced by a “drag” force, given by

(31)

The drag per unit cross-sectional area, Pd, is equal to F/πR ² ; using Equations (7d) and (30), neglecting terms involving R _v , and integratging, we find

(32)

Equation (32) is then our fundamental relationship expressing the ice-flow rate (or equivalently, the sphere’s velocity) as a function of the macroscopic temperature gradient, the applied force (or equivalently drag), and material parameters.

Relative efficacy of regelation and viscous deformation: the “transition radius”

Pressure-induced flow. It is useful to consider the special case G _T = 0 because it most clearly illustrates the way in which the sphere’s rate of motion depends on its size. When G _T = 0, Equation (32) becomes, after some re-arrangement:

(33)

For a given U, the drag is a maximum when R = R*, where

(34)

For R «R _*, motion of the sphere is accommodated primarily by regelation at sub-freezing temperatures and flow through the liquid layer (which, we emphasize, is only a few nm thick); for R >> R _* , the regelative process is very inefficient and the sphere’s motion is accommodated primarily by deformation of the ice. At R = R _* , neither the regelative nor the deformational process is particularly efficient and the resistance to motion is the greatest.

It should be noted that, not surprisingly, our expression (Equation (34)) for the transition radius (with G _T = 0) is identical, to within a small numerical constant, to the “transition wavelength” of Shreve’s (1984, p. 343; Equation (11)) sliding theory for cold-based glaciers. An analogous situation exists with respect to the Nye–Kamb glacier-sliding theory and Watts’ (unpublished) analysis for sphere motion through temperate ice.

Flow with macroscopic temperature gradient. The essential physics of the sphere migration are not altered by the existence of a non-zero value of G _T , although the transition radius is altered. Equation (32) may be re-written in dimensional form as

(35)

where all symbols, including R _* , are as defined previously. The transition radius, now denoted by R_tranS, is still defined by the condition ∂P _d /∂R = 0. This leads to a quartic equation for R _trans;

(36)

Explicit evaluation of the quartic is extremely tedious and unnecessary, as we may easily place bounds on the value of Rtrans. Only one positive, finite value of R _trans can satisfy Equation (36). Clearly, if G _T → 0, we must find Rtrans = R _*. Furthermore, if one examines the functional behavior of ∂R _trans/∂G _T and ∂² R _trans /∂G ² _T , it is easy to show that for G _T > 0, the extreme value of R _travs is 4^–1/³R_*, and that such an extremum is a minimum. (This extremum in fact occurs as G _T → ∞.) Hence, 4^–1/3 R _* ⩽ R _trans ⩽ R _* , which is a rather narrow range (4^–1/3 = 0.62).

In Figure 3, we show the upper bound on R _trans as a function of (–T) in the temperature below 0 C. The temperature dependence of R _trans is contained implicitly (cf. Equation (37)) in the temperature dependence of ice viscosity η_i, water viscosity η_w, and liquid-layer thickness h _c(The calculation of η_i, which is also stress-dependent, is described in Appendix B.) The characteristic liquid-layer thickness h _c is assumed to be given by (Gilpin, 1980[c]) h _c = h ₁(–T _c)^–1/2.4, where h ₁, = 3.5 nm K^1/2.4. This expression for h arises from Gilpin’s best fit of his theoretical predictions to his data on wire regelation, discussed next.

Fig. 3. Transition radius between regelation-dominated and creep-dominated flow. for various driving stresses P_d’ Effective viscosity of ice is dependent on both temperature and P _d

Transition radius: relevance for wire-regelation studies

It is difficult to compare the theoretical predictions above to experimental data, because the only set of experiments on motion of a sphere through ice (Townsend and Vickery, 1967) were conducted with the ice at 0°C and atmospheric pressure. Our predictions of transition radius for motion of sub-freezing ice past a sphere may, however, provide guidance for interpreting results of “wire-regelation” experiments (e.g. Telford and Turner, 1963; Drake and Shreve, 1973; Gilpin, 1980[c]; Tozuka and Wakahama, 1983). Intuitively, it seems likely that the motion of “very small” wires through ice will involved essentially no “viscous” deformation of the ice, whereas the regelation process will be very inefficient for “large” wires (cf. Nye, 1967; Tozuka and Wakahama, 1983). ^* Let us therefore assume that the predictions of transition radius given in Figure 3 are relevant to the wire-motion problem, and compare these predictions with experimental data.

Most experiments on wire movement through ice have been conducted in the pressure-melting regime; the best-known data on wire movement through sub-freezing ice would appear to be those of Telford and Turner (1963), Gilpin (1980[c], and Tozuka and Wakahama (1983). Telford and Turner (1963) used a 225 μm radius wire and ice in the temperature range –0.5° to –4° C, with an applied pressure of 46 bar. Figure 3 therefore suggests that the chosen wire radius was at or above the transition radius for the entire temperature range of the experiments, and that “viscous” deformation of the ice should have been very

Gilpin’s (1980[c]) own experiments on wire movement through sub-freezing ice involved wires with a variety of radii and pressures < c. 10 bar (cf. Gilpin, 1980[c], figs 6, 7, and 8). Even for those tests at very low temperatures (down to –35°C), the wire radii were significantly less than R _trans. Ice-deformational effects in Gilpin’s tests were therefore probably negligible.

Tozuka and Wakahama (1983, p. 4153) reported measurements using a 150 μn radius piano wire at pressures up to 50 bar and temperatures down to –1°C. These measurements were intended explicitly to examine the role of ice deformation in the overall wire motion. Tozuka and Wakahama estimated that the transition from regelation-dominated motion to creep-dominated motion occurred in the pressure range 15–30 bar at –0.7°C. A transition radius of 150 μm at these conditions is quite consistent with our estimates (cf. Fig. 3),

The transition radius and the “rigid-ice” approximation for frozen soils

The theoretical predictions of R _trans illustrated in Figure 3 are also directly relevant to the question raised in the introduction about the correctness of the rigid-ice formulation in the theory of frost-heaving (O’Neill and Miller, 1985). Although a soil is composed of a multitude of grains, with none likely to be spherical, the results of our analysis of ice flow past a single sphere should be illustrative of the basic physics involved in ice movement through a soil.

Physically plausible conditions under which substantial amounts of ground-freezing and heave occur almost never involve overburden loads p _0B of more than a few bars (cf. O’Neill and Miller, 1985). Force-balance considerations analogous to those presented in the next section (cf. Philip, 1980, p. 203) indicate that Pd/p _OB should be of order R/L _f, where L _f is the thickness of the zone through which ice moves. R/L _f is necessarily less than one. We therefore conclude that typical values p _d of of interest in frozen-ground phenomena are usually no more than a few bars. Figure 3 then shows that the transition size R _trans, should be no less than c. 100 μm, and at least several times that value at temperatures above about –2°C. Because soils that exhibit significant heave are characterized by grain-sizes of less than several tens of microns, it seems clear that viscous deformation of the ice should play a negligible role in the overall ice motion. O’Neill and Miller’s (1985) treatment of the pore ice as “rigid” therefore seems quite reasonable. Neglect of viscous deformation also permits us to construct an approximate theory for pore-ice motion.

3. Ice Movement Through a Porous Medium: An Approximate Analysis

Real porous media have exceedingly complex micro-structures that defy description. A number of workers have attempted to explain macroscopic properties of porous materials, such as hydraulic permeability, electrical conductivity, and elastic moduli, on the basis of simplified microphysical models (e.g. Dullien [^c1979]; Seeburger and Nur, 1984; Yale, unpublished). Such models always involve highly idealized descriptions of the pore-space geometry but, nonetheless, seem to describe the essential physical phenomena. Similarly, we cannot hope to describe in all detail the microstructure of a soil and rigorously model pore-ice motion through such a material; we can, however, construct an idealized model of a soil and thereby elucidate the basic physics of pore-ice motion.

The lack of importance of viscous deformation in pore-ice motion, as discussed above, allows us to use the “rigid-ice” approximation. This is exceedingly convenient, because Philip (1980) has presented results for ice motion past arrays of cylinders for the pressure-melting regime. We may directly adapt his analysis – in much the same way as we adapted Watts’ (unpublished) work in the preceding section – to solve the problem of thermally induced regelation of ice through an array of cylinders, which we will take as our highly idealized model of a soil. ^*

The basic result we adapt from Philip’s (1980) analysis is his expression for the temperature field at the ice–cylinder interface. Obviously, when there is a macroscopic temperature gradient (G _T ≠ 0), the mean temperature of individual cylinders may differ, but the results for variation of temperature about any cylinder remain valid. We then proceed as in Gilpin’s (1979) analysis to solve for the liquid-layer thickness and, finally, for the drag force.

For an infinite square array of cylinders, each of radius R and with array spacing H (Fig, 4), and assuming that the thermal conductivity of the cylinder k _s equals that of the ice k _i , the temperature field at any point may be represented mathematically by an infinite series, each term representing the contribution of an individual cylinder which behaves formally as a “thermal dipole.” The series may formally be summed by contour-integral methods but the solution would be in terms of elliptic functions (Philip, 1980). ^† Philip (1980) has presented what he terms “a good approximation” for the temperature field in terms of simpler functions. We have modified his solution to include the case of G _T ≠ 0 (cf. Gilpin, 1979, p. 242) and changed to a cylindrical coordinate system, where the origin may be taken at the center line of any cylinder of interest so long as we look only at temperature variations instead of absolute values.

Fig. 4. Idealized porous medium, composed of a cubic packing of cylinders

The temperature variation at the cylinder–ice interface may be expressed as

(37)

where p _i = 1/v_i is the density of ice. ^* Later calculations (e.g. of the drag force) would essentially involve integrating this expression, a task that appears quite daunting. One may obtain useful approximations to Equation (37) by expanding in terms of trigonometric functions. The larger the quantity r/h, that is, the closer the cylinders are to each other, the more terms must be kept to yield a sufficiently close approximation. For the problem at hand, we may elucidate the basic physics even with a fairly “low order” expression. This restricts us to fairly large cylinder spacings but makes the mathematical manipulations much easier.

The lowest-order expansion of Equation (37) that still includes the effect of multiple cylinders may be shown to be

(38)

where β = (πR/H)². This may be shown to be a very good approximation as long as h /2R (i.e. cylinder spacing/cylinder diameter) is greater than about 2.87 (whereas a cubic close-packed array would have h /2R = 1).

Equation (38) is combined with an expression for the variation in liquid-layer thickness (cf. Gilpin, 1979, p. 241-42; and our Equation (16)):

(39)

The procedure for calculating the drag proceeds exactly as in Gilpin (1979) and in the foregoing analysis, hence it need not be repeated here. We will again restrict our attention to cases in which the liquid-layer thickness does not vary greatly around the cylinder, i.e. the dimensionless thickness perturbation h’ « 1. The drag force per unit cross-sectional area of the cylinder, P _d , must now be understood as P _d = F/2lR, where F/l is the drag force per unit length of the cylinder. The relationship between U, P _d , nd G _T may be written as

(40)

In the absence of externally applied forces, the existence of a non-zero drag P _d for each cylinder in an infinite array requires a macroscopic gradient of ice pressure. This becomes clear from simple force-balance considerations. The drag per unit length F/l must be balanced by a gradient in ice pressure p _i such that

(41)

where z ₀ is the z-coordinate of the “row” of cylinders of interest and x ₀ is an arbitrary point along such a row. The mean gradient of ice pressure G _p is conveniently defined as (cf. Philip, 1980, p. 203)

(42)

Combining Equations (41) and (42), and our definition of P _d for a cylinder, gives G _p = –2rp _d /h ² , Equation (40) may therefore be written as

(43)

It is clear from Equation (43) that the ice-flow rate U is proportional to the gradient of a “generalized potential” Ф where

(44)

with _i , the ice pressure, and , the ice temperature, measured relative to datums that are conveniently chosen as atmospheric pressure and 0°C. The form of the potential Ф bears a curious resemblance to the potential that appears in Gilpin’s (1980[a]) frost-heave model. Gilpin assumed that the flow rate in the unfrozen liquid layer would be proportional to the gradient of water pressure in that layer. Using his thermodynamic description of the unfrozen water, he then found that the flow rate q is proportional to the gradient of the potential Ψ (our notation), where

(45)

independent of the details of grain packing in the porous material. Because the regelative ice flux is inextricably tied up with water flow in the unfrozen liquid layers, this suggests to us that an “exact” solution to the pore-ice regelation problem, rather than the approximate analysis presented above, would lead us to a driving potential Ф equal to Ψ. Strictly, this must be considered speculation; however, we note that there is a well-known experimental criterion that describes the conditions necessary for cessation of frost-heave (and presumably pore-ice motion); this criterion is (e.g. Radd and Oertle, 1973)

(46)

which is identical to the thermodynamic condition for cessation of flow in the liquid layer (Gilpin, 1980[a], p. 919).

It should be noted from Equation (43) that steady flow necessarily requires U to be uniform throughout the spacefilling array of cylinders. (Recall that we have followed O’Neill and Miller (1985) in assuming that the pore ice forms a connected network.) Because the liquid-layer thickness h _c varies with temperature, hence position, steady flow is only possible if the gradient of Ф (Equation (44)) varies spatially in such a way as to make U uniform. B. Hallet (personal communication) has suggested that the requirement of spatially uniform U holds only if “through–flow” of H₂O is restricted to the solid ice, and that in the pore space of a real frozen soil, with unfrozen “capillary” water in addition to adsorbed water films, steady flow of H₂0 could occur even if U varied spatially. K. Hutter (personal communication) has also pointed out that U need not be spatially uniform if particle spacings are not everywhere the same. However, without prescribing U in some fashion, it is difficult to see how an analysis of the pore-ice regelation process could have been developed. It is unlikely that the assumption of spatially uniform U is grossly in error.

We propose to define an apparent hydraulic conductivity K _R on the basis of Equation (43). This hydraulic conductivity is defined by a Darcy’s law type of expression

(47)

where the factor of v _w /v _i on the left-hand side “corrects” for the density difference between water and ice. For the geometry considered,

(48)

or relating h/r to porosity Φ by straightforward geometrical considerations,

(49)

For more general, non-cubical packings of cylinders, a reasonable functional form for K _R may be

(50)

where is a numerical constant of 0(1) and the function f depends on porosity and the geometry of packing.

The apparent hydraulic conductivity defined by Equation (48) is shown as a function of temperature in Figure 5 for several choices of r. The array was assumed to be in cubic close packing (h/r = 2). The hydraulic conductivity values in Figure 5 are rather low in comparison with most measured values for frozen soils (e.g. Burt and Williams, 1976; Horiguchi and Miller, 1980, 1983). This result is actually not particularly surprising. In the present model, all of the pore space – aside from the extremely thin unfrozen films of water at ice–grain interfaces – is assumed to be ice-filled. In a real porous medium, the great variability in pore sizes and shapes leads to a continuous variation in ice content as a function of sub-freezing temperature (e.g. Miller, 1973). The real pore space would contain a geometrically complex mixture of ice and water. H₂O mass transfer in a real soil should occur by two processes not considered in the present model:

Fig. 5. Apparent hydraulic conductivity of frozen porous medium, for various particle radii R. Cubic close packing is assumed.

• (i) Flow in the pore water outside of the very thin “adsorbed” films on grain surfaces (O’Neill and Miller, 1985);

• (ii) Through-flow in the adsorbed films (Gilpin, 1980[a]), as opposed to the local flow associated with the regelative mechanism as analyzed here. ^*

In detail, the various modes of H₂0 mass transfer will probably interact. Furthermore, the present analysis, when considered in the context of the work by Gilpin (1980[a] and O’Neill and Miller (1985), suggests that the gradient of the generalized potential ψ = _i + (L

/v _i T _a) should be taken as the “driving force” for all H₂O mass transport in a frozen porous medium. ^†

To the best of our knowledge, the above analysis leads to the first prediction of rates of pore-ice motion. O’Neill and Miller (1985, p. 286-87), who explicitly considered pore-ice motion in their numerical frost-heave simulations, avoided the need for an explicit physical model predicting the value of U (V _I in their notation) by using certain mass-balance considerations. This procedure, although not addressing the basic physics of the pore-ice regelation process, was very convenient for computational purposes. Our formulation does point out the fundamental physical controls on the regelation process but is unfortunately not convenient for computational simulations because of the idealizations involved in describing the packing of “grains”. Further work along these lines may lead to a more useful theoretical formulation. Such work would need to include a more realistic model of the pore space, which in a real porous material would contain both ice and water (e.g. O’Neill and Miller, 1985).

Discussion

The foregoing analysis has two important consequences. Prediction of the transition radius R _* as a function of temperature and applied load (or drag) provides clear guidance to future investigators who may seek to extend earlier experimental work on wire regelation at sub-freezing temperatures. The rates of wire movement are so exceedingly small at temperatures much below 0°C (Gilpin, 1980[c]) that experimenters might be tempted to use large loads to increase wire speeds. Our results make it clear that wire radii would have to be quite small, perhaps c. 10 μm, to avoid appreciable deformation in the ice. This may lead to restrictions on the types of wires used.

The formal analysis also points to the reasonableness of the rigid-ice approximation used by R.D. Miller and co-workers in their studies of frost-heaving. The notion that pore ice should be “rigid”, yet mobile, in a freezing soil has not gained general acceptance among workers in that field, in spite of what we view as highly persuasive arguments in its favor: arguments summarized in the recent paper by O’Neill and Miller (1985). Our analysis essentially predicts that pore-ice motion should occur in an undeformable porous medium, and specifies the functional dependence of the rate of motion on parameters such as temperature, temperature gradient, ice-pressure gradient, and grain-size. The model should be essentially valid even for a soil, which clearly is not undeformable (i.e. the soil grains might move relative to each other by processes other than ice-lens formation), as long as the rate of any relative grain motion is small compared to the rate of pore-ice movement.

The analysis developed here for ice movement through a porous material should also be useful for examining theoretically the way in which the basal ice of a glacier might “invade” the pores of an underlying layer of glacial till. We defer such discussion to a separate paper, in which we examine this issue for both temperate and cold-based glaciers.