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. Reference ShreveShreve (1984) modified the Nye–Kamb version of glacier-sliding theory (Reference NyeNye, 1969, 1970; Reference KambKamb, 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 Reference Echelmeyer and ZhongxiangEchelmeyer and Zhongxiang (in press) and by Reference Hoekstra and MillerHallet 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” (Reference MillerMiller, 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 Reference O’Neill and MillerO’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” (Reference O’Neill and MillerO’Neill and Miller, 1985, p. 283).
The predicted migration of grains through ice was observed by Reference Horiguchi and MillerHoekstra and Miller (1967) and by Reference Seeburger and NurRömkens and Miller (1973). As a corollary, Reference O’Neill and MillerO’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. Reference NyeNye, 1967; Reference GilpinGilpin, 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 Reference O’Neill and MillerO’Neill and Miller, 1985).
O’Neill and Miller’s description of the pore ice as “rigid” deserves further examination. By “rigid”, Reference O’Neill and MillerO’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 (Reference WattsWatts, unpublished) and Reference GilpinGilpin’s (1979) model of the liquid layer at the interface between sub-freezing ice and an embedded particle. This study, besides complementing Reference ShreveShreve’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.Footnote * . 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. Reference ShreveShreve, 1984).
The first part of this analysis directly parallels that of Reference WattsWatts (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
which combined with the incompressibility condition
implies that ∇2 p = 0. Symmetry conditions, along with the constraint that p → 0 as r → ∞, lead to a solution for the pressure
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.Reference Bird, Bird, Stewart and LightfootBird and others, 1960, p. 90) used to find the stresses. The aforementioned far-field condition on flow is equivalent to the boundary condition
Furthermore, the liquid layer at the ice–sphere interface lubricates that boundary, hence
Equation (5) is sufficiently accurate as long as
where h is the liquid-layer thickness. With these boundary conditions, the velocities and stresses are readily found to
with all other σ ij vanishing.
If the ice were temperate, the constant A would be determined (cf. Reference NyeNye, 1967; Reference WattsWatts, 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. Reference GilpinGilpin, 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 Reference GilpinGilpin (1979) and applied by him to problems of wire regelation, particle migration, and ice lensing (Reference GilpinGilpin, 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 (Reference GilpinGilpin, 1979, p. 238)
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
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 0 where “y 0 is of the order of a few molecular dimensions” (Reference GilpinGilpin, 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 (Reference GilpinGilpin, 1979, p. 239; Equation (5)):
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 (Reference GilpinGilpin, 1979, p. 239; Equation (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 Reference GilpinGilpin (1979, p. 240; Equation (17)):
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
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
We may evaluate the integral in Equation (14) by recognizing that Equation (10) may be re-cast as
and using Equation (11) to eliminate 
wh, For h >> y
0
, we find
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 Reference WattsWatts’ (unpublished) analysis of a sphere moving through temperate ice. Equation (11) may be re-written as
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
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 Reference GilpinGilpin (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. Reference GilpinGilpin, 1979, p. 242)
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 (Reference PhilipPhilip, 1980, p.195-98), and if the temperature drop across the liquid layer is negligible. This latter condition requires (Reference NyeNye, 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.
Reference Drake and ShreveDrake 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 (Reference Townsend and VickeryTelford and Turner, 1963; Reference Hallet, Hallet, Gregory, Stubbs and AndersonGilpin, 1980[c]), although these authors did not discuss any water trace behind the wires. The exact reasons for the transition remain uncertain (Reference Drake and ShreveDrake and Shreve, 1973, p. 69; Reference Hallet, Hallet, Gregory, Stubbs and AndersonGilpin, 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
Liquid-layer thickness
Before proceeding to the solution of Equations (20) and (21), it is useful to recast these equations into dimensionless form. Following Reference GilpinGilpin (1979, p. 292), the liquid-layer thickness will be re-expressed as
where h c is the equilibrium thickness for a stationary sphere, given by
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 (Reference GilpinGilpin, 1979 n 242)
We may now write three characteristic velocities that arise out of Equations (20) and (21), two of them identical to those defined by Reference GilpinGilpin (1979, p. 242), viz.:
and the third, new to our analysis:
Reference GilpinGilpin (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 Reference GilpinGilpin (1979, p. 242):
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
and
Equating these last two expressions, we find that
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 (Reference GilpinGilpin, 1979, p. 245-46), so we may neglect terms in R v in Equation (32). R c = (l/2)(v i /v w)2 (η i/η w)(h c/R)3 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
The drag per unit cross-sectional area, Pd, is equal to F/πR 2 ; using Equations (7d) and (30), neglecting terms involving R v , and integratging, we find
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:
For a given U, the drag is a maximum when R = R*, where
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 Reference ShreveShreve’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 Reference WattsWatts’ (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
where all symbols, including R * , are as defined previously. The transition radius, now denoted by RtranS, is still defined by the condition ∂P d /∂R = 0. This leads to a quartic equation for R trans;
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 ∂2 R trans /∂G 2 T , it is easy to show that for G T > 0, the extreme value of R travs is 4–1/3R*, 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 (Reference Hallet, Hallet, Gregory, Stubbs and AndersonGilpin, 1980[c]) h c = h 1(–T c)–1/2.4, where h 1, = 3.5 nm K1/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 Pd’ 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 (Reference Townsend and VickeryTownsend 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. Reference Telford and TurnerTelford and Turner, 1963; Reference Drake and ShreveDrake and Shreve, 1973; Reference Hallet, Hallet, Gregory, Stubbs and AndersonGilpin, 1980[c]; Reference Tozuka and WakahamaTozuka 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. Reference NyeNye, 1967; Reference Tozuka and WakahamaTozuka and Wakahama, 1983).Footnote * 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 Reference Telford and TurnerTelford and Turner (1963), Gilpin (1980[c], and Reference Tozuka and WakahamaTozuka and Wakahama (1983). Reference Telford and TurnerTelford 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.
Reference Tozuka and WakahamaTozuka 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 (Reference O’Neill and MillerO’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. Reference O’Neill and MillerO’Neill and Miller, 1985). Force-balance considerations analogous to those presented in the next section (cf. Reference PhilipPhilip, 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. Reference O’Neill and MillerO’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. Reference DullienDullien [c1979]; Reference ShreveSeeburger and Nur, 1984; Reference YaleYale, 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 Reference PhilipPhilip (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 Reference WattsWatts’ (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.Footnote *
The basic result we adapt from Reference PhilipPhilip’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 Reference GilpinGilpin’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 (Reference PhilipPhilip, 1980).Footnote † 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. Reference GilpinGilpin, 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
where p i = 1/vi is the density of ice.Footnote * 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
where β = (πR/H)2. 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. Reference GilpinGilpin, 1979, p. 241-42; and our Equation (16)):
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
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
where z 0 is the z-coordinate of the “row” of cylinders of interest and x 0 is an arbitrary point along such a row. The mean gradient of ice pressure G p is conveniently defined as (cf. Reference PhilipPhilip, 1980, p. 203)
Combining Equations (41) and (42), and our definition of P d for a cylinder, gives G p = –2rp d /h 2 , Equation (40) may therefore be written as
It is clear from Equation (43) that the ice-flow rate U is proportional to the gradient of a “generalized potential” Ф where
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 Reference GilpinGilpin’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
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. Reference Radd and OertleRadd and Oertle, 1973)
which is identical to the thermodynamic condition for cessation of flow in the liquid layer (Reference GilpinGilpin, 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 Reference O’Neill and MillerO’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 H2O 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 H20 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
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,
or relating h/r to porosity Φ by straightforward geometrical considerations,
For more general, non-cubical packings of cylinders, a reasonable functional form for K R may be
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. Reference Burt and WilliamsBurt and Williams, 1976; Reference Horiguchi and MillerHoriguchi 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. Reference MillerMiller, 1973). The real pore space would contain a geometrically complex mixture of ice and water. H2O 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 (Reference O’Neill and MillerO’Neill and Miller, 1985);
-
(ii) Through-flow in the adsorbed films (Reference GilpinGilpin, 1980[a]), as opposed to the local flow associated with the regelative mechanism as analyzed here.Footnote *
In detail, the various modes of H20 mass transfer will probably interact. Furthermore, the present analysis, when considered in the context of the work by Reference GilpinGilpin (1980[a] and Reference O’Neill and MillerO’Neill and Miller (1985), suggests that the gradient of the generalized potential ψ = 
i + (L
To the best of our knowledge, the above analysis leads to the first prediction of rates of pore-ice motion. Reference O’Neill and MillerO’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. Reference O’Neill and MillerO’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 (Reference GilpinGilpin, 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 Reference O’Neill and MillerMiller (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.
Acknowledgements
Discussions with B. Hallet on the physics of ice-lensing motivated the analyses presented above. B. Hallet, K. Hutter, and an anonymous reviewer carefully critiqued an earlier version of this paper. Correspondence with R.D. Miller helped clarify certain aspects of the pore-ice regelation problem. F. Bardsley drafted the figures and Quaternary Research Center staff helped prepare the typescript. Financial support was provided by U.S. National Science Foundation grant EAR83-19119.
Appendix A. Liquid-Layer Thickness
Equations (21) and (22) of the main text may be written in dimensionless form as:
where all symbols are as defined in the main text. The very lowest-order approximation to use in integrating these equations, and one valid only for h’« 1, is to take 1 + h’ ≈ 1, which makes the right-hand sides of Equations A-1 and A-2 independent of h’. This leads to the expressions given by Equations (30) and (31) of the main text.
The next level of approximation for h’ « 1 would be a linearized expansion with, for example, (1 + h’ ) 3 ≈ 1 + 3h’. Integration of Equations (A-l) and (A-2), in this case, obviously leads to complicated exponentials (if we neglect the non-linear terms) and makes determination of the constant A a complicated exercise. The resulting expression must then be expanded in terms of trigonometric functions in order to calculate the drag (cf. Equation (33) of main text). This “refined” expression for the drag does not differ greatly from the expression found from the simpler analysis. The simplicity of Equation (35) for the drag is considered adequate reason to stay with the lowest-order analysis.
APPENDIX B Calculation of the Apparent Ice Viscosity
We have assumed ice to have a Newtonian-viscous rheology, in spite of the proliferation of experimental evidence to the contrary, in order to facilitate our analysis, much as Reference NyeNye (1969, 1970) did in his glacier-sliding theory. Nonetheless, we may roughly account for the actual rheology of ice by treating the ice viscosity as stress-dependent. In particular, we will assume (cf. Reference ShreveShreve, 1984, p. 344, table I)
where η0 = constant, Q = activation energy for creep, R = gas constant, and T a = absolute temperature, and the meaning of <τ2> will be explained shortly.
The “effective shear stress” τ is defined by (e.g. Reference PatersonPaterson, 1981, p. 85)
where the τ’ ij are the deviatoric stresses in the ice, and the summation convention for repeated subscripts is implied. For the case of ice flow past a lubricated sphere, we found that the only non-zero stress component was σ rr , hence
where p is the mean stress, equal to σ rr /3 in this case.
We now assume that, in order to characterize the effective viscosity for ice deformation adjacent to the sphere, we may use « τ 2 >, the mean-square value of τ over the sphere’s surface:
However, it is easy to show that
Using Equation (B-6) in (B-5) and integrating, we find
and finally using Equation (B-7) in (B-1),
To fix the constant η0 , we follow Shreve and assume a viscosity of 1 bar a (3.16 – 1012 Pa s) at 0°C and an effective shear stress of 1 bar (0.1 MPa). Taking Q = 6 – 104 J mol–1 (cf. Reference Telford and TurnerShreve, 1984) and R = 8.314 J mol–1 K–1, we therefore find η0 = 1.06 – 1011 Pa3 s. Calculated values of R * are then found using Equation (B-8) for ice viscosity.
