Introduction

The origin of ice streams, which discharge 70% and 90% of ice from the Greenland and Antarctic ice sheets respectively, is widely recognized among glaciologists as a major unsolved problem. Another unsolved problem, not recognized by most glaciologists, is the disposition of a vast reservoir of gravitational potential energy in these ice sheets, resulting from a density inversion in ice near the bed that puts a thick layer of higher-density colder ice above a thinner layer of lower-density warmer ice. A solution to both problems is proposed in which the ‘ceiling’ of colder ice sinks into the ‘basement’ of warmer ice along ice flowbands by means of a unique form of polycrystalline thermal convection, thereby transforming these flowbands into ice streams and transforming the reservoir of gravitational potential energy into the thermal and kinetic energy associated with the fast motion of ice streams.

This investigation begins with a brief history of speculation surrounding the possibility of thermal convection in polar ice sheets, including implications for ice-sheet dynamics, abrupt climate change, and mantle dynamics associated with plate tectonics. Then the unique conditions for thermal convection in a polycrystalline polar ice sheet are reviewed. The investigation concludes with a proposal to understand the origin of ice streams in polar ice sheets as a necessary interaction between horizontal motion in advective gravitational sheet flow and vertical motion in convective gravitational streamflow.

Background

The possibility of thermal convection in polar ice sheets occurred to me while I attended lectures on ‘World Geology’ given by J. Tuzo Wilson at The Ohio State University during the spring semester of 1968. Wilson pioneered the geological paradigm shift that has resulted in a new synthesis now known as plate tectonics. He published three brief papers in those early years that were seminal. Reference WilsonWilson (1963) proposed deep-mantle thermal convection that sent plumes punching through a spreading Pacific lithosphere plate to account for the age of volcanoes in the Hawaiian island chain. Reference WilsonWilson (1965) proposed transform faults that connected oceanic rifts to oceanic trenches, thereby showing that Earth’s lithosphere consisted of interacting moving tectonic plates. Reference WilsonWilson (1966) presented evidence that the Atlantic Ocean had closed and reopened since the Cambrian, and proposed a theory in which lithosphere rifts opened oceans and lithosphere trenches closed oceans, producing mountain ranges along suture zones, a theory now accepted and called the ‘Wilson cycle of orogenesis’. Not one mathematical equation cluttered Wilson’s work.

Wilson believed thermal convection in Earth’s mantle caused lithosphere plates to move across Earth’s surface. In fluids heated from below, thermal convection begins when heat can be transported upwards more efficiently by mass transport than by thermal conduction within a static fluid. This transformation takes place when a ‘Rayleigh number’ reaches a critical value (Reference StruttStrutt, 1916). A 2164 m hole to bedrock had just been cored through the West Antarctic ice sheet in 1968 at old Byrd Station (80° S, 120° W) by the United States Army’s Cold Regions Research and Engineering Laboratory (CRREL). Temperatures and densities had been logged (Reference Gow, Ueda and GarfieldGow and others, 1968). This and earlier work on the creep (Reference GlenGlen, 1955) and conductive (Reference PounderPounder, 1965) properties of ice allowed calculation of the Rayleigh number at Byrd Station. It was just past the critical value for fluids (Reference HughesHughes, 1970). I was so bold as to suggest that, if thermal convection in the Antarctic ice sheet were to initiate and propagate like a chain reaction through the ice, it could lead to a surge of the ice sheet and perhaps to a new global Ice Age as Antarctic ice flooded the World Ocean. Early radio-echo records recorded vertical reflections that I interpreted as plumes of convective flow, an interpretation that had initial support from Reference Bentley and CraryBentley (1971).

This led to a series of papers that investigated the nature of the Rayleigh criterion, as applied to polycrystalline rocks instead of viscous fluids, and to the Antarctic ice sheet in particular (Reference HughesHughes, 1971, Reference Hughes1972a,Reference Hughesb). The idea of the Antarctic ice sheet as a simplified miniature mantle accessible at all depths by deep drilling and open to three-dimensional seismic and radar mapping over time led to other papers proposing a migrating tetrahedral pattern of thermal convection in Earth’s mantle that seemed to fit the changing pattern of lithosphere plate boundaries and continental drift over geologic time (Reference HughesHughes, 1973a,Reference Hughesb, Reference Hughes1975). Flood tide for this work was an invited paper based on the Reference WeertmanWeertman (1967) theory of mantle convection, but applied to the Antarctic ice sheet (Reference HughesHughes, 1976).

Ebb tide set in after Reference HarrisonHarrison (1971a,Reference Harrisonb) showed that radio-echo ‘convection plumes’ were actually specular echo patterns unrelated to ice flow. Reference PatersonPaterson (1981, p. 171) stated, ‘the most convincing evidence against convection comes from radio-echo sounding records; these show layers within the ice that are continuous over distances of hundreds of kilometers’. I countered by pointing to the absence of these layers in the ‘echo-free zone’ of ice near the bed (Fig. 1). This zone is below the density inversion and therefore is the ice most vulnerable to thermal convection, with the undisturbed layered ice riding passively above it (Reference HughesHughes, 1985). The boundary between the upper layered ice and the echo-free zone was essentially horizontal, even above a rugged bed, as would be expected for convection below a thermally caused density inversion. Reference Robin and MillarRobin and Millar (1982) described ‘near-vertical cusps and fingers’ that disrupted layering just above the echo-free zone, which I interpreted as the lowermost layers being broken up and drawn down by convecting systems that thoroughly scrambled deeper layers and thereby produced the echo-free zone (Reference HughesHughes, 1985; Reference Hughes1998, p. 140–144).

Fig. 1. The ‘echo-free zone’ in the Antarctic ice sheet. Top: A radio-echo record along flow at site T on map. Disturbed horizontal boundary is consistent with density-based limit on transient convection plumes in echo-free ice and inconsistent with advective ice flow over rugged bed topography. Bottom: A radio-echo record across flow at site B on map. Downwarped upper reflection horizons are consistent with steady-state convection rolls aligned with flow in echo-free ice, and inconsistent with advective ice flow over a smooth bed. From Reference Robin and MillarRobin and Millar (1982) and Reference HughesHughes (1998).

The echo-free zone was also ice in which oxygen isotope stratigraphy down core holes through the Greenland and Antarctic ice sheets displayed a characteristic ‘spikiness’ that I interpreted as sills in a dike–sill pattern of thermal convection (Reference HughesHughes, 1977), referring to Reference HughesHughes (1976). These spikes had also been interpreted as recording rapid changes in climate (Reference Johnsen, Dansgaard, Clausen and LangwayJohnsen and others, 1972), the accepted explanation (e.g. Reference MayewskiMayewski and others, 1997), now called Dansgaard–Oeschger events.

In summary, over nearly four decades I have argued without success for thermal convection in polar ice sheets, the Antarctic ice sheet in particular, as a key to unlocking mysteries in the dynamics of ice sheets, Earth’s mantle, and climate change. In this my latest attempt, I now propose that thermal convection provides an explanation for the origin of ice streams.

Thermal Buoyancy Stress

The basic argument for thermal convection is as follows. Advective gravitational motion along flow lines in the horizontal x direction is resisted by a deviator shear stress σ_xz that increases with depth in the vertical z direction:

(1)

Where ρ ₁ is mean ice density, g is gravity acceleration and θ is ice surface slope; whereas, apart from advection motion, vertical convective gravitational motion along z is resisted by a deviator normal stress given by

(2)

Where Δρ ₁ is the density decrease caused by thermal expansion of ice below the density inversion (Reference HughesHughes, 1976). Reference GowGow (1970) showed that Δρ ₁ = 4.0 kg m⁻³ below the density inversion at old Byrd Station (Fig. 2). The upper third of ice has lower densities as firn is compressed into ice in the first 300 m and as air bubbles are dissolved in ice in the next 700 m. Thermal convection is possible only below the density inversion at 1000 m depth, where the decrease in density is entirely explained by the temperature increase (Reference HughesHughes, 1970). Equation (1) gives σ_xz = τ ₀ = 45.6 kPa as the advective gravitational driving stress using maximum surface slope and total ice thickness z = h ₁ = 2164 m (Reference DrewryDrewry, 1983). Equation (2) gives as the convective gravitational driving stress due to thermal buoyancy using Δρ ₁ = 4.0 kg m⁻³ through ice thickness z = d _C = 1164 m below the thermal density inversion.

Fig. 2. Ice densities at in situ temperatures and surface pressures from the Byrd Station core hole in West Antarctica (Reference GowGow, 1970). Inset: similar data from the GISP2 core hole in Greenland (Reference GowGow and others, 1997).

Since τ ₀ drives horizontal advective flow, what prevents virtually identical σ _C from driving vertical convective flow? An answer to this question was based on the fact that Reference GowGow (1970) measured ice densities from ice cores removed from the core hole at old Byrd Station, so ice densities were measured at the surface atmospheric pressure, not at in situ pressures down the core hole. Reference Marion and JakubowskiMarion and Jakubowski (2004) measured the compressibility of ice to 2.0 kbar (0.2 ×10⁹ Pa). Their results give Δρ ₁ = 1.2 kg m⁻³ for the in situ density decrease in the lower 1164 m of ice down the Byrd Station core hole, where ice temperature increases from −28.7°C to −1.6°C (Fig. 3). This gives a convective gravitational driving stress of 13.7 kPa, compared to the advective gravitational driving stress of 45.6 kPa. Thermal convection could remain a possibility below the density inversion. If so, thermal convection provides a mechanism for releasing the vast reservoir of gravitational potential energy associated with the density inversion, a reservoir untapped in any model of ice-sheet dynamics. Advective gravitational flow has little effect on σ _C because Δρ ₁ increases as d _C decreases. Figure 4 shows an example for Law Dome ice cap (67° S, 113° E) on the Wilkes Land coast of East Antarctica (Reference BuddBudd, 1969, fig. 4.16, p. 92).

Fig. 3. Ice densities at in situ temperatures and at both surface and in situ pressures from the Byrd Station core hole (Reference Marion and JakubowskiMarion and Jakubowski, 2004).

Fig. 4. Calculated steady-state temperature profiles along an ice flowline originating at the summit of the Law Dome ice cap, Wilkes Land, East Antarctica (Reference BuddBudd, 1969).

Similar conditions exist in central Greenland. In ice 3 km thick, Reference GowGow and others (1997) located the density inversion 1 km above a frozen bed for ice cores extracted from the Greenland Ice Sheet Project 2 (GISP2) core hole. The reservoir of gravitational potential energy associated with the density inversion is comparable for both polar ice sheets. Converting this potential energy into kinetic energy by ‘dropping’ the cold ice ceiling into the warm ice basement is hindered because the warm ice is not easily displaced. So what happens, if anything? Reference Chappellaz, Brook, Blunier and MalaizéChappellaz and others (1997) compared the oxygen isotope stratigraphy in the Greenland Icecore Project (GRIP) and GISP2 ice cores to bedrock, respectively cored at the summit and 28 km from the summit of the central high dome of the Greenland ice sheet. They were unable to correlate the oxygen isotope stratigraphy in the lowest 300 m of ice. This ice also contained silt and small rock clasts from the bed (Reference GowGow and others, 1997). Published explanations for these observations range from climate instability during the Eemian to spatial changes in bed topographic and thermal conditions. Within the realm of possibilities is that gravitational sinking of the cold ice ceiling has displaced all but the lowest 300 m of warm basement ice. Oxygen isotope stratigraphy will not be disturbed in ceiling ice if it sinks en masse, but basement ice might be disturbed as it is displaced by thermal convection.

This kind of convection in ice sheets has similarities with the view of thermal convection in Earth’s mantle advanced by Reference WilsonWilson (1963, Reference Wilson1965, Reference Wilson1966). His lithosphere plates advecting horizontally above his convecting mantle would correspond to a cold ice ceiling advecting above a warm ice basement. Lithosphere ‘hot spots’ correspond to basaltic convection plumes rising from the deep mantle that turn on and off, and when turned on may punch through the lithosphere to produce volcanoes, such as the Hawaiian volcanic chain (Reference WilsonWilson, 1963), when partial melting occurs at the top of a plume. The ceiling of cold ice in polar ice sheets is much too thick to be penetrated by rising plumes of warm basement ice, but the idea of thermal convection in a warm polycrystalline mantle, whether of basalt rock or ice, is the central concern. Plumes can rise as either isolated pipes, as with volcanic ‘hot spots’ within lithosphere plates, or as interconnected curtains (rippling sheets) along plate boundaries composed of rifts, transform faults, and trenches (Reference WilsonWilson, 1965, Reference Wilson1966).

For thermal convection in fluids, superimposed advective flow transforms polygonal platform ‘cells’ of sinking cooler fluid enclosed by interconnected ‘curtains’ of rising warmer fluid into elongated convection ‘rolls’ with the roll axis aligned in the direction of advective flow (Reference Low and BruntLow and Brunt, 1925; Reference JeffreysJeffreys, 1928; Reference DeardorffDeardorff, 1965; Reference Gallagher and MercerGallagher and Mercer, 1965; Reference Davies-JonesDavies-Jones, 1971). The viscosity contrast between rising and falling limbs of thermal convective flow is small in fluids, so these limbs have comparable scales. However, the contrast is great in polycrystalline rocks.

Rising convective flow confined to narrow pipes and curtains seems to be a feature of thermal convection in polycrystalline rocks, as Reference WilsonWilson (1963, Reference Wilson1965) proposed for Earth’s mantle. That should be the situation for thermal convection in polar ice sheets as well, as they are polycrystalline rocks. Creep deformation in ice and other rocks is closer to the plastic end of the visco-plastic creep spectrum, whereas fluid flow is at the viscous end, so a visco-plastic yield stress characterizes polycrystalline creep (Reference HughesHughes, 1998, p. 196–200). A yield stress tends to localize creep deformation in the parts of a polycrystalline system where warmer temperatures and an ‘easy-glide’ crystal orientation are found (Reference HughesHughes, 1998, p. 119–122). For thermal convection in polar ice sheets, as in Earth’s mantle, these conditions confine rising flow to narrow pipes or curtains within broad regions of slowly sinking flow. At or near ice domes, warm rising flow may be either isolated pipes in transient creep or interconnected curtains in steady-state creep.

Increasing advective flow downslope from ice domes would align curtains in the flow directions. Slowly sinking flow between faster rising flow confined to narrow curtains would occur en masse in broad bands that, like the rolls in advecting fluid convection, would be aligned in the directions of advective flow. These bands will be incipient ice streams because they are flanked by curtains of warm ice that have an easy-glide ice fabric favorable to both vertical convective flow and horizontal advective flow. This is the basic reason for linking thermal convection in polar ice sheets to the initiation of streamflow.

It is worth noting that permafrost displays a form of non-thermal convection in which annual freeze–thaw processes bring loose rocks to the surface in interconnected rising curtains that enclose polygons of more compact material. This occurs on level ground in the permafrost condition. When the ground begins to slope, these polygons become moving bands flanked by narrow walls of loose rising stones.

Thermal convection has been analyzed most rigorously when it occurs in purely viscous fluids heated from below (e.g. Reference StruttStrutt, 1916; Reference KnopoffKnopoff, 1964; Reference Somerscales and DropkinSomerscales and Dropkin, 1966; Reference Somerscales and GazdaSomerscales and Gazda, 1969). Heat transport by mass transport via thermal convection begins in fluids when a condition of neutral buoyancy is attained. This situation is characterized by a ‘critical Rayleigh number’ of 657 when the fluid layer has free top and bottom surfaces, of 1101 when one surface is free and one surface is rigid, and of 1708 when both surfaces are rigid. In polar ice sheets, the top and bottom surfaces are free when the density inversion is closest to the ice surface, as for low accumulation rates with no advective flow (Reference RobinRobin, 1955), and when basal ice floats on a subglacial lake (e.g. Vostok Subglacial Lake, East Antarctica). Advective flow lowers the density inversion, making it a rigid boundary topped by thick ice. A frozen bed makes the bottom surface a rigid boundary. A deep density inversion above a frozen bed makes both top and bottom surfaces rigid boundaries. Another concern is the reliability of applying a critical Rayleigh number derived for fluid convection to polycrystal-line convection. Reference WeertmanWeertman (1967) examined this problem for viscous thermal convection in Earth’s mantle, and obtained a critical Rayleigh number, Ra, of 300 for free top and bottom boundary conditions in a comparison with Ra = 657 obtained by Reference StruttStrutt (1916) for fluid convection. With these constraints in mind, the Rayleigh criterion for initiating and sustaining thermal convection can be examined. As will be shown, quite different results follow.

Initiating Thermal Convection in Ice Sheets

Adjective flow cannot remove the density inversion in the Antarctic ice sheet, but can change its height above the bed. Convective flow can remove a density inversion caused by thermal expansion as Antarctic ice warms toward the bed, provided that a critical Rayleigh number is attained. As originally derived to quantify the onset of thermal convection in a fluid heated from below (Reference StruttStrutt, 1916), Ra is the dimensionless ratio of the upward velocity of heat flow due to convection, u _C, to the upward rate of heat flow due to conduction, r _C, in a fluid of height h _C having creep viscosity η and thermal diffusivity κ. Taking as the vertical creep rate due to thermal convection, ΔT _C as the upward temperature decreases in a fluid layer of height h _C, volume V _C, mass M _C and mean density ρ ₁:

(3)

If Δρ ₁ is caused entirely by thermal expansion uncorrected for compressibility, the volume coefficient α_v of thermal expansion can be related to Δρ _{1 for temperature increase} ΔT _C as follows:

(4)

In this case, noting that Δρ ₁ is negative for positive ΔT _C,

(5)

Equation (5) gives Ra for thermal convection in incompressible fluids. Equation (3) gives Ra for polar ice sheets in which ice compressibility is an important correction.

In applying Equation (3), the first problem is how to apply to polycrystalline ice a viscosity η that applies to fluids. Figure 5 shows the visco-plastic creep spectrum based on the usual flow law for steady-state creep in ice and rock:

Fig. 5. The visco-plastic spectrum for steady-state creep and two criteria for visco-plastic yielding. A sharp knee develops in creep curves when the visco-plastic exponent n increases, where is the strain rate at creep stress σ, is the strain rate at plastic yield stress σ ₀ and . Visco-plastic yielding can be defined at the knee or at the stress intercept of the tangent line at . For n = 3, σ _v = 0.386σ ₀ at the knee and σ _v = 0.667 σ ₀ at the stress intercept (Reference HughesHughes, 1998).

(6)

where is the strain rate produced by deviator stress σ, A is a hardness parameter dependent mainly on the temperature and fabric of polycrystalline ice or rock, n is a visco-plastic exponent in the range 1 ≤ n ≤ ∞, σ ₀ is a plastic yield stress, and is the strain rate for all values of n when σ = σ ₀. By convention, linear viscosity is defined as:

(7)

Applying Equation (7) to viscous flow, when n = 1. Applying Equation (7) to plastic flow when n = ∞, all values of η are possible, from η → ∞ as to η →0 as . By giving any viscosity, Equation (7) gives no viscosity. As seen in Figure 5, when n = 1, , and when n = ∞, only two viscosities exist, η = ∞ at and η = 0 at . These values for viscous flow and plastic flow are delivered by visco-plastic viscosity η _v defined as:

(8)

It seems reasonable to conclude that Equation (8) also gives the correct viscosity for 1 < n < ∞, including n = 3 for glacier ice. For rising convection plumes in ice sheets, for velocity u_z along direction z gives v. Originally viscosity was defined for fluids shearing between opposite-moving parallel plates such that for velocity u_x normal to direction z, which would apply to the viscosity of basal ice shearing over a frozen bed.

Visco-plastic viscosity η _v given by Equation (8) is a necessary but not sufficient condition for describing creep in polycrystalline ice or rock. A visco-plastic yield stress σ _v is also required (Reference HughesHughes, 1998, p. 196–200). It can be described as the stress σ in Figure 5 where the bend in the curvature of the σ versus curve is greatest, so that η _v is high below σ _v and low above σ _v. This occurs where radius of curvature R for the σ versus curve is a minimum. Alternatively, σ _v can be defined as the stress intercept of a straight-line tangent to the creep curve at σ = σ ₀. Both possibilities are shown in Figure 5 for n = 3. This feature of creep in ice has not been investigated experimentally to determine the proper yielding criterion. This is important because convective flow would be minor if σ _C < σ _v, giving a higher η _v, but could be major if σ _C > σ _v, giving a lower η _v. These considerations do not apply to convection in fluids, for which σ _v = 0.

Creep in polycrystalline ice and rock is described by η _v and σ _v, with η _v high below σ _v and low above σ _v. This change becomes less pronounced toward the viscous end of the visco-plastic creep spectrum and vanishes at n = 1 (see Fig. 5). In addition, η _v and σ _v change when recrystallization converts a grain fabric unfavorable for easy glide of dislocations from grain to grain to one favoring easy glide (Reference HughesHughes, 1998, p. 117–122). Figure 6 shows how strain ε varies with strain energy E at constant stress σ (top), stress σ at constant strain rate (middle), and time t at constant stress σ (bottom) in polycrystalline ice or rock when recrystallization converts a random grain fabric to grain fabrics allowing progressively more easy glide between grains. An upper unstable visco-plastic yield stress σ _v =σ _U occurs at strain ε _U in strain-hardened ice or rock just before recrystallization begins, and a lower stable visco-plastic yield stress σ _v = σ _S occurs at strain ε _S in strain-softened ice or rock when recrystallization is complete.

Fig. 6. Idealized deformation of polycrystalline ice or rock. Top: strain energy E versus strain ε at constant stress σ. Middle: stress σ versus strain ε at constant strain rate . Bottom: time t versus strain ε at constant stress σ. Recrystallization in specified stress fields for a random grain fabric R produces circle fabrics C, and fabrics with 4, 3, 2 and 1 poles as easy glide becomes more dominant. Subscripts U and S refer to unstable and stable conditions before and after recrystallization, respectively (Reference HughesHughes, 1998).

By analogy, a block on an inclined plane begins to slide at a basal shear stress τ ₀ that is higher than τ ₀ needed to sustain sliding after lowering the inclination. This gives rise to static and dynamics coefficients of friction, respectively.

Curves in Figure 6 are called strain energy curves for E versus ε at constant σ (top), flow curves for σ versus ε at constant (middle), and creep curves for t versus ε at constant σ (bottom). Elastic strain increases linearly with stress in the flow curves with no change in strain energy, and instantaneously in the respective energy and creep curves. For upward thermal conduction in ice sheets, this represents atoms that vibrate about their lattice sites with increasing amplitudes as temperature increases, so that heat is conducted upward as warm atoms in basal ice transmit their vibrational energy to cold atoms toward the ice surface. Upward thermal convection begins when vibrational amplitudes become great enough to free atoms from their lattice sites, allowing them to move upward so that heat transport is by mass transport. This upward motion of atoms occurs in chains of line defects that dislocate the regular spacing of atoms in the ice crystal lattice. For this reason, these line defects are called dislocations (Reference Weertman and WeertmanWeertman and Weertman, 1964, Reference Weertman, Weertman and Nabarro1980). Strain energy in the E versus ε plots in Figure 6 (top) increases as moving dislocations pile up at grain boundaries in polycrystalline ice and rock. This is why E remains constant as ε increases for vibrational elastic strain, and why elastic strain takes place in zero time for the creep curves in Figure 6 (bottom). When dislocations form and begin to move, dE/dε increases, dσ/dε decreases and dt/dε increases as dislocations pile up at grain boundaries to cause strain hardening. The increase of dt/dε from zero at t = 0 is the same as a decrease of strain rate from the infinite strain rate during elastic deformation. Since this occurs at constant stress in creep curves, Ra, defined by Equation (3), must initially be infinite even though the thermal buoyancy stress is small and finite.

The first major finding of this analysis is that initiating thermal convection in ice sheets, or other polycrystalline rock masses heated from below such as Earth’s mantle, cannot be prevented, but it will weaken over time. Can it be sustained over time?

Sustaining Thermal Convection in Ice Sheets

A second major task is to determine whether thermal convection in ice sheets can be sustained, once it begins. An answer to this question depends on how thermal convection is visualized. Figure 7 is a cartoon showing how a convection plume may form. Ice at the bottom of the Byrd Station core hole in central West Antarctica and the GISP2 core hole in central Greenland consists of large inter-grown ice crystals undergoing dynamic recrystallization (Reference GowGow, 1970; Reference GowGow and others, 1997), probably because strain energy accumulates more rapidly in warm ice that exists in a stress field that changes continually as ice moves over a bed having variable roughness and hydrological conditions that control ice–bed coupling. Dynamic recrystallization allows very high strain rates, so this ice favors strain rates that begin with the infinite strain rate of elastic deformation. Above this ice layer, simple shear over the bed produces ‘easy glide’ represented by the ‘single-maximum’ ice fabric and creep curve 1 in Figure 6 (bottom). Optic axes of ice grains tend to be vertical because of easy-glide planes having a hexagonal symmetry line parallel to the regional bed slope, which is of order 10⁻² to 10⁻³.

Fig. 7. Development of a thermal convection plume in the easy-glide basal layers of an ice sheet. Easy-glide horizontal layers deforming in simple shear (top) are deformed in bending shear to initiate a thermal plume (middle) and become easy-glide vertical layers encasing the rising plume (bottom).

A rising plume of warm basal ice, whether the plume is a pipe or a curtain, must upwarp the hexagonal planes having an easy-glide ice fabric produced by advective flow of ice moving in simple shear above the bed. A vertical plume of ice in convective flow will begin as transient creep and undergo strain hardening as it produces dislocations that pile up against grain boundaries not oriented to favor easy glide from grain to grain. Hard glide results from the dislocation ‘traffic jam’ at grain boundaries where strain energy accumulates and eventually causes recrystallization. This relieves the traffic jam by producing a new grain fabric favoring easy glide between grains in the rising plume of convecting ice. Figure 7 shows how upwarping ice layers in horizontal easy-glide advective flow over the bed produce ice layers in vertical easy-glide convective flow above the bed. Easy glide occurs in a sheath of ice that encloses the rising plume. Colder ice that does not recrystallize sinks slowly en masse between the plumes. Ice stratigraphy is preserved in sinking ice.

Figure 8 shows how the creep curve begins with elastic strain and passes through transient creep to slow, unstable steady-state creep in strain-hardened ice before recrystallization to fast, stable steady-state creep in strain-softened ice after recrystallization. Strain hardening leads to visco-plastic instability at time t _v. This triggers recrystallization that is complete at time t _r, after which an easy-glide ice fabric exists for the stress field causing the deformation. Strain ε over time t for the creep curve in Figure 8 can be written for t < t _v:

Fig. 8. Components of the creep curve in polycrystalline ice. Deformation begins with elastic strain ε _e at time t = 0 and continues with visco-plastic ductile strain, beginning with transient strain and m = 1/3 for ice, to which is added slow unstable steady-state strain , which is gradually replaced by fast stable steady-state strain during time t _r– t _v of recrystallization. These stages are identified as 1–4. In dynamic recrystallization, t _v = t _r = 0, and the only ductile strains are and .

(9)

for t _v ≤ t ≤ t _r:

(10)

and for t > t _r:

(11)

where ε _e is elastic strain, is transient strain, is steady-state strain before recrystallization, and is steady-state strain after recrystallization. Reference GlenGlen (1955) found that m = 1/3 for ice.

A rising convection plume in the polar ice sheets covering Antarctica and Greenland begins in basal ice where dynamic recrystallization would move both t _v and t _r to t = 0 on the time axis of Figure 8 so that:

(12)

where occurs in strain-hardened ice and occurs in strain-softened ice. The rising ice plume upwarps ice deforming in simple shear above the basal ice layer. Instead of easy glide in simple shear, ice in the plume now deforms as easy glide in bending shear, such as would occur between cards when a deck of cards is bent. Bending produces a fold in the deck of cards. A similar fold would occur if the plume rose as a curtain, and folding is self-enclosed if the plume rose as a pipe. The flow law for ice, in terms of vertical strain rate in Equation (6), would be:

(13)

where R _C is a scalar of strain rate, which for vertical strain rate is, following Reference ThomasThomas (1973a,Reference Thomasb):

(14)

Rising convection pipes are most likely near interior ice domes where advective ice flow is minimal. Vertical shear on the sides of pipes gives . Minor downslope advective flow along x and y contributes to and , and can bend (fold) the pipes. For rising flow in pipes, and . Setting and taking n = 3 for ice (Reference GlenGlen, 1955), Equation (14) becomes R _C ≈ [1 − 1/2 + 1/4 + 1 + 1]/[2 − 1/2] = 1.5 for convection pipes. Rising convection curtains are most likely where advective ice flow is strong away from interior ice domes. Like convective rolls in an advecting fluid, curtains will align with direction x of advective flow, so and . Setting and n = 3 for ice, noting that shears the curtain in the downslope direction, Equation (14) becomes R _C = [1 – ­1 + 1 + 1 + 1] = 3.0 for convection curtains.

Setting and in Equation (12), with and both applied to in Equation (13):

(15)

where A _v and A _r are the respective ice-hardness parameters before and after recrystallization. Therefore A _v ^mn > A _r ⁿ .

The strain rate obtained from Equation (15) for constant R _C, A _v, A _r and σ _C is:

(16)

Rayleigh number Ra defined by Equation (3) becomes:

(17)

Taking m = 1/3 and n = 3 for ice (Reference GlenGlen, 1955):

(18)

Equation (18) shows that Ra decreases from infinity at t = 0 to a stable steady-state value at t = ∞ of:

(19)

No matter what Rayleigh numbers are calculated from Equations (18) and (19), including Ra = ∞ at t = 0, there exists a critical Rayleigh number Ra_C at neutral buoyancy. Take Ra = Ra_C ≈ 10³ as the critical Rayleigh number for steady-state thermal convection when t →∞, with Δρ ₁ = 1.2 kg m⁻³ below the density inversion in polar ice sheets, g = 9.8 m s⁻² and κ = 37.2 m² a⁻¹. Reference Denton and HughesDenton and Hughes (1981, p.246–249) determined that A = 140kPa a^1/3 in warm basal Antarctic ice. Bending shear in glacial ice produces an easy-glide ice fabric that was reproduced in the laboratory at a shear stress of σ = 117 kPa at −3°C, and for which recrystallization increased the strain rate from to , a 12.8-fold increase (Fig. 9; Reference Hughes and NakagawaHughes and Nakagawa, 1989). The corresponding decrease in A is A _r = A/(12.8)^1/3 = 59.8 kPa a^1/3. Solving Equation (19) for h _C gives h _C = 1266 m using R _C = 1.5 for convection pipes and h _C = 1103m using R _C = 3.0 for convection curtains. These values compare with h _C ≈ 1000 m below the density inversion at Byrd Station, central West Antarctica, (Reference GowGow, 1970) and at GISP2, central Greenland (Reference GowGow and others, 1997). It seems that steady-state thermal convection may be marginally stable, especially for rising curtains, and could therefore scramble ice below the density inversion and produce the echo-free zone in Figure 1. Strong radar systems now produce weak, discontinuous radar-reflecting horizons below the density inversion in Antarctica and Greenland. This would be consistent with the fragile nature of steady-state thermal convection in these ice sheets. Reference Alley, Gow, Meese, Fitzpatrick, Waddington and BolzanAlley and others (1997) reported a fold in Greenland ice near the bottom of the GISP2 core hole. They attributed it to possible variations in bed topography for a frozen bed, but a convection pipe or curtain also produces a fold according to the cartoon in Figure 7 and advective flow could tilt the fold. Ice at this depth is also where oxygen isotope stratigraphy from the GISP2 and GRIP cores, only 28 km apart on the high central dome of the Greenland ice sheet, fails to match (Reference Alley, Gow and MeeseAlley and others, 1995; Reference Chappellaz, Brook, Blunier and MalaizéChappellaz and others, 1997). This is consistent with thermal convection in ice near the bed, at least intermittently in transient creep.

Fig. 9. Creep curves for simple shear in polycrystalline ice at −3°C and shear stresses of 117 kPa (curve A) and 55 kPa (curve B).

Transient creep occurs throughout creep deformation but it is eventually overtaken by steady-state creep, as seen in Figure 8. Strain hardening causes a reduction in strain rate during transient creep. In dynamic recrystallization, steady-state creep is by easy glide rather than by hard glide. Transient creep in hard glide dominates until time t _C obtained by setting t = t _C in Equation (16) and equating the two right-hand terms. Taking m = 1/3 and n = 3 for ice (Reference GlenGlen, 1955):

(20)

Take σ _C = 13.7 kPa due to the density inversion, A _r = 59.8kPa a^1/3 for easy glide after recrystallization, and Av as a time-dependent viscosity before recrystallization (Reference PatersonPaterson, 1981, fig. 3.7, p.41) that gives A = 8 × 10¹³ Pa s = 2.5 × 10³ kPa a as the viscosity of ice in the viscous part of the deformation map for creep in ice (Reference Duval, Ashby and AndermanDuval and others, 1983; Reference HookeHooke, 2005, fig. 4.16, p. 64). For time-dependent viscous creep in Equation (15) this would be A _v = 2.5 × 10³ kPa a a ^{m − 1} = 2.5 × 10³ kPa a ^m . Taking m = 1/3 for ice, Equation (20) then gives t _C = 0.20 years for convection pipes with R _C = 1.5 and t _C = 0.10 years for convection curtains with R _C = 3.0. Thermal convection dominated by transient creep in hard glide lasts about 1 month. Then it is overtaken by steady-state creep in easy glide.

Figure 9 shows creep data obtained by Reference Hughes and NakagawaHughes and Nakagawa (1989) for polycrystalline ice having a random grain fabric and deformed in simple shear at −3°C. Recrystallization had not taken place after 1100 hours when the applied shear stress was 55 kPa, but began after 150 hours and was complete by 550 hours when the stress was 117 kPa. The single-maximum easy-glide ice fabric was as strong as those produced by easy glide in bending shear at a similar temperature behind the face of a calving ice wall on Deception Island (60° S, 60.6°W) off the Antarctic Peninsula. The ice wall and shear bands date from the 12 August 1970 volcanic eruption. Ice was sampled in December 1970. The creep curve for 55 kPa in Figure 9 shows that transient creep remained an important strain component after 1100 hours (46 days) of creep. It therefore seems possible that transient creep could be important >1 month after a convection plume is initiated by a thermal buoyancy stress of 13.7 kPa for warm ice in dynamic recrystallization that continually produces an easy-glide ice fabric from strain-hardened ice. In that case, t _v = 0 and gives in Figure 8 for creep strain by easy glide in recrystallized ice.

Figure 10 plots the creep curve given by Equation (15) as the sum of transient creep in hard glide and steady-state creep in easy glide for dynamic recrystallization in a rising curtain of convecting ice, using values of m, n, R _C, σ _C, A _v and A _r determined here. Note that the two contributions to creep are equal at t = 0.10 years. For dynamic recrystallization at −3°C and a buoyancy stress of 13.7 kPa, creep strains from transient creep in hard glide and from steady-state creep in easy glide are equal at a strain of about 0.004. This is two orders of magnitude less than the strain of about 0.4 at 117 kPa (Fig. 9). Dynamic recrystallization in warm ice allows strain softening to occur simultaneously with strain hardening linked to transient creep at the beginning of creep deformation, not at some later time when more strain energy has accumulated.

Fig. 10. A creep curve for dynamic recrystallization at −3°C given by Equation (18). The creep curve is the sum of ductile strains for hard glide in transient creep before recrystallization and for easy glide in steady-state creep after recrystallization.

Figure 11 plots the strain rate and the Rayleigh number Ra for convective creep using these values and h _C = 1103 m for a rising curtain. Note that the Ra curve tracks the curve, as required by Equation (3), both being infinite at time t = 0. Figure 11 shows that Ra exceeds any critical Rayleigh number Ra_C for thermal convection in fluids until t = 0.10 years, when transient creep in hard glide equals steady-state creep in easy glide for the rising convection curtain. At t = 0.10 years, Ra = 1732 and is nearly identical to the critical Rayleigh number Ra_C = 1708 for rigid top and bottom boundaries when initiating thermal convection in a fluid heated from below (Reference KnopoffKnopoff, 1964). After 10 years, Ra = 1297, compared to Ra_C = 1101 when the fluid has one rigid and one free boundary. There is little change in Ra for the rising curtain of ice after 10 years, with Ra = 1296 at infinite time. These values will be somewhat lower for rising pipes of ice, so convecting pipes should be rare compared to convecting curtains. Similarly, crustal hot spots within lithosphere plates and linked by Reference WilsonWilson (1963) to pipes of mantle convection are rare compared to crustal rifts along plate boundaries linked to curtains of mantle convection (Reference HughesHughes, 1973a,Reference Hughesb, Reference Hughes1975). Hot spots produce sporadic volcanic eruptions whereas rifts discharge lava continuously. This would be the case if pipes operate as intermittent transient convection and curtains operate as continuous steady-state convection. The same can be expected for thermal convection in polar ice sheets, given the lower Ra value for pipes due to the lower R _C in Equation (14). Surface expressions of these features are suppressed on ice sheets because convective activity occurs below the density inversion, but convection in the icy mantle of Jupiter’s moon Europa may produce similar surface features (Reference PappalardoPappalardo and others, 1998; Reference SullivanSullivan and others, 1998; Reference Pappalardo and CarrPappalardo and Carr, 2004; Reference Freeman, Moresi and MayFreeman and others, 2006).

Fig. 11. The time dependence of convective strain rate, , from Equation (16) and Rayleigh number, Ra, from Equation (17) for h _C = 1103 m. Here m = 1=3 and n = 3 for ice.

Convection Curtains Facilitate Streamflow

Sheet flow becomes streamflow in the tributaries of ice streams. How this is accomplished remains unknown, but a thawed bed, wet basal sediments and channelized bed topography are helpful (Reference Bindschadler, Bamber, Anandakrishnan, Alley and BindschadlerBindschadler and others, 2001). These are the conditions that favor convection curtains aligned with advective ice flow. Vaughan and others (2001, fig. 5) show the dendritic nature of tributaries that converge on Pine Island and Thwaites Glaciers, the two fastest ice streams in Antarctica. Lateral shear zones separate fast streamflow from slow sheet flow across the sides of ice streams. Strain heating in lateral shear zones would tend to elevate ice temperature, and hence the density inversion, in an ice stream by warming ice flowing laterally across the shear zones, according to modeling experiments by Reference Raymond, Echelmeyer, Whillans, Doake, Alley and BindschadlerRaymond and others (2001) and confirmed for Rutford Ice Stream, West Antarctica, by Reference Doake, Frolich, Mantripp, Smith and VaughanDoake and others (1987). Lateral flow assists return convective flow.

Figure 12 shows temperature profiles on opposite sides of a lateral shear zone that separates Whillans Ice Stream, West Antarctica, from the Unicorn, an ice ridge between the two main tributaries of Whillans Ice Stream (Reference Kamb, Alley and BindschadlerKamb, 2001, fig. 4). The temperature of Whillans Ice Stream increases with depth virtually from the surface in ice 1080 m thick. That and a wet bed of low-strength till give Whillans Ice Stream close to free top and bottom boundary conditions that allow the lowest critical Rayleigh number, Ra = 657, for initiating thermal convection in fluids heated from below. In contrast, the temperature inversion in the Unicorn is 300 m below the surface and the bed is frozen, giving rigid top and bottom boundary conditions, which require Ra_C = 1708 to initiate thermal convection in fluids heated from below. Rising warm convective flow in lateral shear zones cools and sinks in ice streams, not ice ridges.

Fig. 12. Temperature profiles in Whillans Ice Stream (Ice Stream B) and Unicorn Ice Ridge (Ridge B1–B2), West Antarctica (Reference Kamb, Alley and BindschadlerKamb, 2001).

Figure 13 is a cartoon showing a transverse cross-section of an ice stream at three downslope locations. Convection curtains aligned with advective ice flow are shown as transverse cross-sections in lateral shear zones, with a soft thawed bed on the lowering ice-stream side and a hard frozen bed on the higher ice-ridge side. Along with the higher density inversion on the ice-stream side, based on the temperature profiles in Figure 12, this gives close to free–free and rigid–rigid top and bottom boundary conditions across the lateral shear zones. This is why the sinking limbs of convection would occur in ice streams rather than in ice ridges between ice streams. Sinking ice empties the reservoir of gravitational potential energy above the density inversion and converts it into kinetic energy of motion that facilitates, and may even initiate, streamflow. Since ice streams discharge 90% of Antarctic ice and 70% of Greenland ice, most of the potential energy reservoir is drained if this form of thermal convection acts to lower and thin the reservoir as shown in Figure 13.

Fig. 13. A transverse cross-section of an ice stream in which convection curtains rise in the lateral shear zones. Horizontal lines represent stratigraphic horizons in the ice for the reservoir of gravitational potential energy (P) above the density inversion (D) and for thermal convection in ice below the density inversion (C), as manifested by warm rising curtains of ice (R) in lateral shear zones (L) and en masse sinking of ice between shear zones. Sinking lowers layer P into layer C, converting gravitational potential energy into kinetic energy of motion as the ice stream moves downslope (top to middle to bottom).

It can be argued that ice streams empty the potential energy reservoir in any case, regardless of whether they are flanked by rising convection curtains in their lateral shear zones. Nevertheless, the easy-glide ice fabric that favors horizontal advective ice flow equally favors vertical convective ice flow. A vertical buoyancy deviator stress is produced by the density inversion, so a vertical strain rate will exist according to Equation (16), and the initial Rayleigh number given by Equation (17) will be infinite, thereby exceeding any possible critical Rayleigh number needed to initiate thermal convection. When the stress field is allowed to warp free top and bottom surface, Reference WeertmanWeertman (1967) found Ra_C = 300, less than half of Ra_C = 657 with no warpage. Figure 12 shows that ice in Whillans Ice Stream is warmer and lower than ice in the Unicorn. If this is the sinking limb warmed by ice in the rising convection curtains, then we have a demonstration of a complete convective system in ice streams and their lateral shear zones.

Conclusions

Even when the density decrease with depth due to rising temperature with depth in an ice sheet is corrected for the density increase due to increasing pressure with depth, a density inversion still occurs and gives rise to a thermal buoyancy stress that causes a significant vertical strain rate. Owing to the crystal anisotropy of ice, with creep deformation confined largely to ‘basal’ planes having hexagonal symmetry, rising convection plumes should be confined to pipes and curtains with slow sinking flow between. Pipes may occur beneath ice domes where surface accumulation rates and advective ice flow are minimal, but curtains aligned with advective ice flow are favored as ice spreads away from domes. Curtains are most likely along the sides of ice streams as slow sheet flow converges to become fast streamflow. This transition begins in the tributaries of ice streams, so each tributary may be flanked by rising convection curtains in its lateral shear zones. As tributaries merge to produce an ice stream, lateral shear becomes confined to the sides of the ice stream. Since lateral shear produces an easy-glide ice fabric that equally favors horizontal advective flow and vertical convective flow, rising convecting curtains will be confined to the outer lateral shear zones that remain when tributaries merge to become the main trunk of an ice stream. Sinking convective flow occurs en masse in the ice stream (see Fig. 13).

A chicken-or-egg question emerges as to whether thermal convection produces ice streams or merely becomes localized in ice streams. Both possibilities can be imagined. For thermal convection to produce ice streams it must precede streamflow. This is possible if sporadic isolated convecting pipes in transient creep under interior ice domes become nodes connecting a polygonal network of convecting curtains in easy-glide steady-state creep as advective flow moves ice downslope from ice domes (see Reference HughesHughes, 1976, fig. 8). Transient pipes would be analogous to those causing volcanic ‘hot spots’ within lithosphere tectonic plates (Reference WilsonWilson, 1963), and steady-state pipes connecting curtains would be analogous to hot ‘triple junctions’ at tectonic plate boundaries (Reference WilsonWilson, 1965). Advective flow away from ice domes would align convecting curtains in the downslope direction. Parallel curtains separated by a basal channel having weak ice-bed coupling would provide the weakened lateral coupling that might initiate streamflow. If ice streams form without convecting curtains rising in converging tributaries of streamflow, the thermal buoyancy stress below the density inversion might produce rising curtains of convecting ice in the lateral shear zones of ice streams and their tributaries. Vertical shear would be much less than longitudinal shear because the buoyancy stress would be much less than the downslope gravitational driving stress, of order 14 kPa compared to 45 kPa. The fact that ice streams are lower and warmer than their flanking ice ridges is consistent with descending convection flow between lateral shear zones.

Thermal convection would be favored by dynamic recrystallization in basal ice having an easy-glide fabric produced by the basal shear stress that resists advective flow of ice over the bed. A thermal buoyancy stress due to the density inversion in overlying ice is greatest at the bed and may locally upwarp basal ice where ice–bed coupling is weakest, such as in thawed patches on a frozen bed or water pockets on a thawed bed. Basal ice deforming in simple shear continues to deform in bending shear where upwarping occurs, thereby preserving the easy-glide ice fabric for both horizontal advecting flow and vertical convecting flow. Progressive upwarping allows easy glide in simple shear within a boundary layer of ice between rising pipes or curtains and the slowly sinking ice that surrounds them. Curtains should replace pipes as advective flow increases downslope from interior ice domes of an ice sheet, with curtains being aligned with advective flow. Curtains then become zones of ice in easy glide for both convective and advective ice flow, so it seems natural for ice streams to develop between curtains and above a bed weakly coupled to basal ice. Ice streams lower as the return limb of thermal convection. This would occur along bedrock channels blanketed by water-saturated till, the very places where Reference Bindschadler, Bamber, Anandakrishnan, Alley and BindschadlerBindschadler and others (2001) found tributaries of ice streams tended to develop.

Figure 14 is a cartoon model showing how an ice stream might be generated from the interaction between horizontal advective ice flow in downslope directions and vertical convective ice flow in curtains aligned with downslope flow. Warm basal ice rises in curtains, allowing ice to sink slowly between curtains. Ice advecting laterally across the curtains delivers warmed ice to replace slowly sinking cold ice, while basal water flows across the curtains because slowly sinking ice between the curtains produces a transverse pressure gradient at the base of overlying ice. Inflow of basal water reduces ice–bed coupling between curtains and allows the slowly sinking ice to be advected downslope faster than flanking ice that becomes more coupled to the bed as it loses basal water. An easy-glide ice fabric in the vertically rising curtains reduces ice–ice coupling across the curtains, thereby accelerating downslope advection of ice sinking between the curtains. This would be the beginning of streamflow. The lowering and warming cold ice empties the reservoir of gravitational potential energy above the density inversion, with ice streams being the ‘sluices’ that are opened by thermal convection.

Fig. 14. A cartoon showing how downslope motion of an ice stream is facilitated by rising convection curtains in its lateral shear zones and by basal water crossing the shear zones and uncoupling the ice stream from its bed. Slow-sinking convection flow occurs between the rising curtains and lowers the surface of the ice stream.

Whether this process can or should be incorporated into existing ice-sheet models is an open question. My personal interest in thermal convection and ice streams dates from my entry into glaciology in 1968, and from the start contained a suspicion that the two flow processes are linked. I have a doppelgänger screaming at me that they are. Many years ago I asked Barclay Kamb to investigate this possibility. He did and told me he could not compose a field experiment that would definitively answer the question. Yet my doppelgänger kept screaming at me that convection exists and does produce ice streams. I hope by this paper my doppelgänger will shut up and spare me the rage. Figures 12 and 14 illustrate a hypothesis that can be tested in the field. Faced with my imminent retirement, I invite others to take up this task.