Introduction
It has long been recognized that the network of liquid veins and nodes that lines the grain intersections in polycrystalline ice plays an important role in material transport and the regulation of temperature within temperate glaciers (e.g. Reference LliboutryLliboutry, 1971; Reference ShreveShreve, 1972; Reference Nye and FrankNye and Frank, 1973; Reference FowlerFowler, 1984). More recently, these studies have been extended and applied towards the understanding of post-depositional diffusion of the trace constituents that are recovered from polar ice cores (e.g. Reference Rempel, Waddington, Wettlaufer and WorsterRempel and others, 2001, Reference Rempel, Wettlaufer and Waddington2002; Reference Barnes, Wolff, Mader, Udisti, Castellano and RöthlisbergerBarnes and others, 2003; Reference Rempel and WettlauferRempel and Wettlaufer, 2003). As the temperature rises towards the melting point, the permeability of the vein network increases and the rate of intergranular flow becomes significant. Current interest in the potential for material transfer between polar ice sheets and subglacial lakes, both under ambient conditions and during future exploration (Reference Souchez, Petit, Tison, Jouzel and VerbekeSouchez and others, 2000, Reference Souchez, Petit, Jouzel, de Angelis and Tison2003; Reference SiegertSiegert and others, 2001, Reference Siegert, Tranter, Ellis-Evans, Priscu and Lyons2003; Reference Carsey, Behar, Lane, Realmuto and EngelhardtBell and others, 2002; Reference PriscuPriscu and others, 2003), motivates a quantitative re-examination of the balances that control intergranular flow. Gradients in the fluid flux may be associated with englacial phase changes that have contributed to the large quantities of refrozen lake water found above Vostok lake, Antarctica, and the accreted ice layers seen in borehole images from Antarctic ice streams (Reference Carsey, Behar, Lane, Realmuto and EngelhardtCarsey and others, 2002). It has even been suggested that refrozen vein water may contribute to drilling problems such as those that hampered operations in the warm basal ice at NorthGRIP (North Greenland Icecore Project; personal communication from S. Johnsen, 2002).
Here, we present a simple conceptual model that predicts the rate of intergranular transport as a function of the local temperature and impurity content, focusing in particular on the lower regions of an ice sheet that is perched above a subglacial liquid reservoir. The central predictions are summarized in Figure 1, which shows a spatially varying flow rate that changes direction a short, but finite, distance above the glacier base. We discuss how spatial variations in fluid flux can be accommodated by melting and freezing along grain boundaries, some distance above the glacier base. Scaling arguments are presented to estimate the magnitude of these effects.
Fig. 1. Dimensionless transport rate uū ≡ U/Uc as a function of dimensionless height above the glacier base ξ ≡ z/zc from Equation (7), assuming uniform impurity loading, grain-size and temperature gradient so that cB = 1, b = 1 and T = ξ. Using nominal parameter values that are representative of conditions above subglacial lakes, Uc V 10-2 ma–1 and zc V 0.3m. The horizontal dotted line separates the regions of downward and upward transport; the transition occurs where z = zc, at dimensionless height ξ = 1. The vertical dotted line separates regions where englacial melting is expected, from regions where englacial freezing can accommodate the changing rate of intergranular flow.
Intergranular Flow
We consider the flow rate through the vein network of a glacier with temperature T(z) that increases towards the pressure-melting point Tm at its base. The glacier floats atop a subglacial liquid reservoir, such as a lake. We assume that the ice pressure is approximately hydrostatic and that the liquid pressure at the base of the ice Pl(z = 0) is equal to the ice overburden pressure Pi(z = 0) ≈ ρigH, where H is the glacier thickness, g ≈ 10ms–2 is the acceleration of gravity, and ρi ≈ 920 kgm–3 is the ice density. Above the base, the Gibbs–Thomson effect causes the liquid pressure to be reduced from the ice pressure by an amount that is proportional to the surface energy γ ≈ 3.5 × 10-2 Jm–2 and the curvature of the ice–liquid interface, so that Pl ≈ Pi – γ/R, where R is the vein radius. For an idealized vein network that runs through polycrystalline ice with liquid volume fraction ϕ and grain diameter b, Reference Nye and FrankNye and Frank (1973) found that the permeability to fluid flow is described by k = ϕ2b2/χ, where the geometrical factor χ ≈ 2 × 103. Darcy’s law gives the fluid transport rate (positive upwards) as
where z is distance above the glacier base, η ≈ 1.8 × 10-3 Pas is the liquid viscosity, and ρl ≈ 103 kgm–3 is its density. Equation (1) describes how the density difference between liquid water and ice leads to a buoyancy-driven flow down through the veins towards the glacier base, and how this gravitational effect is modulated by the effects of surface energy when there are significant gradients in vein radius. An essential feature of this system is that these gravitational and surface-energy effects are of opposite sign near the glacier base, and this leads to a flow reversal so that the net transport at the lake interface is upwards into the ice.
To quantify U, we need to first determine the liquid fraction ϕ and the vein radius R. At temperature T < Tm, the penalty in latent heat that is required for the vein liquid to be at equilibrium with the surrounding ice is offset by the effects of surface energy, which are inversely proportional to the vein radius R, and by the colligative effects of soluble impurities, which depend on the intergranular solute concentration c (quantity of solute molecules per volume of vein liquid). We assume that the glacier contains soluble impurities at bulk concentration cB (quantity of solute molecules per volume of polycrystalline ice–liquid mixture) and that these impurities are completely segregated to the vein liquid so that the local solute balance requires that cB = ϕc. For polycrystalline ice with veins of radius R running along the edges of grains with diameter b ≫ R, the liquid volume fraction is ϕ ≈ R2/(αb2), where α ≈ 1 is a geometrical factor. The vein system is in equilibrium when (e.g. Reference NyeNye, 1991)
where Γ is the slope of the liquidus curve and L ≈ 3.3 × 105 J kg–1 is the latent heat. Dilute solution theory predicts that the liquidus slope for monovalent impurities is Γ ≈ 2 × 10-3Km3 mol–1.
To gauge the relative importance of impurities and surface energy in controlling the vein radius, we define the two characteristic radii 
and Rl = 
as lower bounds on the actual vein radii in the polycrystalline ice. (To calculate Rl in cases where a significant fraction of the impurities are located on grain boundaries and elsewhere, rather than in the vein liquid, we reinterpret cB as the quantity of solute molecules in the vein liquid per volume of polycrystalline ice–liquid mixture.) If no impurities were present, Equation (2) would predict R = Rγ, whereas if the effects of surface energy were negligible, equilibrium would require R ≈ Rl. For the range of grain-sizes and bulk impurity concentrations that are typically found in glacier ice, Rl ≫ Rγ for undercoolings Tm – T of a few milli-Kelvin and greater. This implies that the second term on the right side of Equation (2) is negligible.
We take R ≈ Rl and rewrite Equation (1) as
Equation (3) gives a continuum description of the rate that fluid is transported through the vein network across a unit horizontal cross-section of the polycrystalline matrix. Treating the liquid as incompressible, the product ρlU is equal to the fluid mass flux. To gain a better understanding of the controlling physical effects, it is useful to introduce appropriate scales that group the coefficients so that transport can be expressed in non-dimensional form.
Scaling Arguments
The temperature decreases away from the glacier base so dT/dz < 0 and the dependence of vein radius on temperature produces a contribution to the hydraulic gradient that tends to counteract the effects of buoyancy. Examining the right side of Equation (3), if we assume that temperature variations are the predominant cause of changes in vein radius so that 
then the hydraulic gradient that drives flow (i.e. the term in ‘{}‘) disappears at the point where
There is a wide variation in the bulk impurity concentrations found in meteoric ice, and it is difficult to constrain the salinity of subglacial lakes and reservoirs that are the source of accretion ice. There is also substantial variability in grain-sizes, which range from a few millimeters in diameter for most of the meteoric ice sampled in typical ice cores to perhaps a meter or so in diameter at the bottom of the Vostok core (Reference Souchez, Petit, Tison, Jouzel and VerbekeSouchez and others, 2000). We adopt nominal values of c0 ≈ 10-3 molm–3 and b0 ≈ 4 × 10-3 m to define the dimensionless concentration and grain-size as
We favor a choice for b0 that is much smaller than the grain-sizes found in the top of the 200m accretion layer above Vostok lake, since the samples recovered at that location are expected to have grown considerably during the approximately 20 kyr since the ice was emplaced (Reference Benatov and WettlauferBell and others, 2002). We recognize that the values of c0 and b0 appropriate for a given physical setting may be somewhat different from the nominal values chosen here, but view these as useful for illustrating the basic controlling phenomena.
For heat fluxes typical of those detected in polar ice cores, the magnitude of the temperature gradient is G ≈ 2 × 10-2Km–1. This is a reasonable approximation for the gradient near the lake interface provided that the release of latent heat that accompanies englacial phase changes is much smaller than the conductive heat flux. With the nominal grain-size and impurity content, and with ∂T/∂z ≈ –G, Equation (4) is satisfied and the hydraulic gradient reaches zero when the undercooling Tm – T is approximately 6 × 10-3 K. Assuming that the bulk impurity concentration and grain-size are spatially uniform, the height above the glacier bed at which the transport velocity reaches zero is
which is about 0.3 m for the nominal parameter values. We identify zc as a natural length scale for this problem and non-dimensionalize position using ξ ≡ z/zc. Above ξ = 1 the pressure gradient in the liquid is smaller than hydrostatic and a flow is driven downwards. Below ξ = 1, the effects of surface energy and the large gradients in vein radius cause liquid to be drawn up from below.
A natural velocity scale for fluid transport is obtained by considering the rate of downward flow that would be present at height zc if the hydraulic gradient were dominated by buoyancy so that
We define the dimensionless temperature as T ≡ (Tm – T)/(Gzc) so that T = ξ when the temperature gradient is ∂T/∂z = –G. The dimensionless transport rate is
In the idealized case where the grain-size, impurity content and temperature gradient are constant so that cB = 1, b = 1 and 
Equation (7) predicts that 
and a maximum downward transport rate of Umax = –(44/55)Uc ≈ –0.08Uc occurs at z = (5/4)2zc ≈ 1.6 zc (see Fig. 1). With the nominal parameter values, Uc ≈ 10–2 ma–1 so the maximum rate of downward transport Umax ≈ –10–3 ma–1 is expected approximately 0.5m above the bed. For b0 ≈ 10–2 m, Umax ≈ –0.1ma–1 about 0.1 m above the bed. The maximum downward transport occurs higher up when the bulk impurity concentration is lower, but the rate itself is reduced in the manner described by Equations (5) and (6).
Implications
The dimensionless transport rate from Equation (7) is shown as a function of height in Figure 1 for the case where the composition, grain-size and temperature gradient are all constant, i.e. cB = c0, b = b0 and Tm = T = Gz. In the region above ξ = 1, the transport rate is negative, indicating that the flow is towards the glacier bed. In this idealized case, we predict a downward buoyancy-driven flow that accelerates downstream so that uū ≈ –ξ–2 when ξ ≫ 1. As the effects of surface energy gain importance in the lower regions, the rate of acceleration decreases until the maximum downward rate Umax is reached at ξ ≈ 1.6. Since the ice and water can both be regarded as incompressible, the increasing rate of downward flow above this level suggests that additional liquid is supplied by englacial melting. The rate of downward flow decreases in the lowest portions of the glacier, reaching zero at the level where ξ = 1. The flow is upwards beneath this, transporting water at a rate that decelerates away from the bulk interface as ξ-5/2. To conserve water mass, this decelerating flow must be accompanied by englacial freezing. For example, all of the downward transport crossing from the level ξ ≈ 1.6 where U = Umax is solidified between this height and the level where the flow rate reaches zero at ξ = 1. The magnitude of upward transport can be much higher than 
near the lake interface for ξ < 1. However, the permeability model breaks down as the temperature increases to the point where the predicted vein radius approaches the grain-size, which itself becomes comparable to the height above the glacier base. To have confidence in the predictions of this simple, continuum model, we require that z ≫ b ≫ R; it is the former condition, z ≫ b, that typically controls the lower limit to which this continuum treatment applies.
Further consideration is required to understand the manner in which englacial phase transformations can accommodate the changing rate of intergranular transport. The geometry of the vein network is thermodynamically controlled, as described by Equation (2). Any changes to the vein radius that accompany melting or freezing along the vein walls are modulated by the effects of heat flow as the system adjusts to maintain equilibrium. The presence of thin liquid films along boundaries between adjacent grains may provide alternative locations where englacial phase changes can occur without significant alteration to vein geometry. Evidence for such grain-boundary melting comes from the rapid growth rates observed during laboratory recrystallization experiments (Reference Ohtomo and WakahamaOhtomo and Wakahama, 1983), the inferred changes in grain-growth mechanisms derived from ice-core analyses of fabric and grain-size at temperatures greater than about –10C (Reference Thorsteinsson, Kipfstuhl and MillerThorsteinsson and others, 1997; Reference De La Chapelle, Castelnau, Lipenkov and Duvalde la Chapelle and others, 1998), the transition in ice rheological properties, notably the relatively high activation energy for ice deformation under similar temperature conditions (e.g. Reference PatersonPaterson, 1994, p. 86), and theoretical calculations that demonstrate how the formation of a pre-melted liquid film serves to minimize the total free energy along the boundaries between adjacent ice grains (Reference Benatov and WettlauferBenatov and Wettlaufer, 2004). The suggestion made here is that the changes in flow rate are accompanied by phase changes that involve water exchange between the relatively high-permeability vein network and the significantly less permeable grain boundaries. This would enable the divergence in fluid transport rate to produce a divergence in ice velocity that does not necessarily require a significant deformation to the ice grains themselves, so that the approximation of hydrostatic ice pressure remains valid.
Focusing attention on the region where the intergranular transport is downwards, above ξ = 1, for the scenario depicted in Figure 1 the increase in transport rate towards Umax is concentrated between z ≈ 1.6zc and z ≈ 10zc. For the nominal parameter values, we require that each year a volume of liquid equivalent to a horizontal layer 1 mm thick is melted from the ice in this region of a few meters thickness, and subsequently refrozen in the 17cm thick region between z ≈ 1.6zc and z = zc. This is comparable to the rate of basal accretion inferred for the ice above Vostok lake (Reference Bell, Studinger, Tikku, Clarke, Gutner and MeertensBell and others, 2002) and suggests that, if the modeled conditions fall within the appropriate parameter regime, the recrystallized lake ice sampled by the Vostok ice core should have undergone approximately one cycle of melting and refreezing by this mechanism after emplacement nearer to the glacier base. Such recycling behavior has potential as a purification process, wherein minor quantities of impurities that are incorporated within the grain interiors during the initial solidification of lake water are subsequently segregated to the liquid network. This may have a role in explaining the drop in electrical conductivity measurement signal near the base of the Vostok core at the transition from meteoric to accreted ice (Reference Souchez, Petit, Jouzel, de Angelis and TisonSouchez and others, 2003). More comprehensive modeling efforts designed to reproduce this and related features show promise for unraveling the dynamics in the basal regions of ice sheets.
It has been suggested that frazil ice may accumulate at the reservoir interface (Reference Souchez, Petit, Tison, Jouzel and VerbekeSouchez and others, 2000). The process of compaction and recrystallization that accompanies such an accretion process would certainly alter conditions at the reservoir interface and affect the range over which upward flow can be driven by surface energy effects. Similarly, if the underlying reservoir is highly saline, then a dendritic mushy layer would be expected to form, with much higher liquid fractions than the current model is formulated to describe (e.g.Reference Worster, Davis, Huppert and MullerWorster, 1992). Nevertheless, even with such complications, the flow model presented here is still valid further above the interface, where the liquid pathways take the form of the network of veins that we have assumed.
Conclusions
Transfer of material from the glacier into the reservoir is expected only by melting right at the glacier–reservoir interface. Intergranular flow through glaciers is predominantly downwards, driven by the hydraulic gradient that accompanies the density difference between liquid water and ice. However, when the glacier caps a subglacial reservoir that is at a fluid pressure equal to the ice overburden pressure, the effects of surface energy lead to a reversal in flow direction above the glacier base. Because of this, the intergranular flow does not empty into the reservoir directly. This implies that the transport of trace constituents, including nutrients from the overlying meteroric ice and potential contaminants from drilling operations, should not be carried below this level into the lake itself by intergranular flow. Gradients in the rate of fluid transport lead to englacial melting that is concentrated in a region of a few meters thickness above the glacier base. This is situated immediately above a layer that is typically just a few tens of centimeters thick wherein solidification occurs. Recycling of intergranular water between these two layers can occur at a rate comparable to the accretion of intergranular ice above Vostok lake.
Acknowledgements
This work has benefited from many useful discussions with J.S. Wettlaufer, thoughtful editorial comments by D.R. MacAyeal and careful reviews by K. Hutter, R. Souchez and an anonymous referee.
