1. Introduction
The idea of a thermally regulated mechanism of glacier surges stems from Robin (1955). The concept is very simple. If an ice sheet is frozen at the base, then the ice flux is very small and it will thicken in time due to the surface accumulation. As il does so, the insulation afforded by the ice allows the base to warm until it reaches the melting point. At this point, the ice can slide and the resultant increased ice flux, if large enough, will draw-down the ice. In principle, this can re-freeze the base, so that a cyclic oscillation ensues. Clarke and others (1977) and Yuen and Schubert (1979) extended this idea by invoking thermal run-away the thermomechanical coupling between ice flow and temperature through the temperature-dependence of the flow law causes a multiple-valued relation between ice flux and thickness.
To understand the import of this, consider the (lumped) mass-conservation equation
where 
is accumulation rate, Q is proportional to ice flux and h is mean ice thickness. If Q is a multi-valued (e.g. S-shaped) function of h, such that h′(Q) < 0 for Q
c < Q < Q
m then, if Q
c < a < Q
m, periodic oscillations will occur which are manifested as surges, with quiescent phases on the lower (slow) branch Q < Q
c being punctuated by rapid surges on the upper (fast) branch Q > Q
m This mechanism has been elaborated by Fowler (1987 b).
However, thermal run-away may not be realistically possible in ice sheets and glaciers. Most simply, the multivaluedness occurs when the basal temperature is above the melting point and so is unphysical. Indeed, Fowler and Larson (1980) showed that, with realistic thermal-boundary conditions, the (unique) steady state was linearly stable, in contrast to the situation for Equation (1.1), if Q c < a < Q m.
Nevertheless, Robin’s basic concept retains its conceptual validity and, indeed, is likely to be relevant to the surges of Trapridge Glacier (Clarke and others, 1984). More recently, a similar concept has been advanced by MacAyeal (1993a) to explain hypothetical surges of the Hudson Strait ice stream which are thought to have occurred during the Last Ice Age (Andrewsand Tedesco, 1992; Clark, 1994), and for which indirect evidence exists in the deep-sea sediment record (Heinrich, 1988; Bond and others, 1992; Grousset and others, 1993). In a simple model, MacAyeal (1993b) described the oscillation by postulating that, when the base of the ice is molten, it relaxes to a dynamic state corresponding to modern ice-stream dynamics, as evidenced in West Antarctica.
As a description, this is fine but it involves an assumption, namely that a temperate base is sufficient to cause fast glacier flow in the Hudson Strait ice stream; clearly, this may not be the case, as is attested by the existence of non-streaming parts of the West Antarctic ice sheet which appear to be wet-based (e.g. Ice Stream C) and the plentiful glaciers with temperate bases which do not surge (e.g. Trapridge Glacier at the present time). The purpose of this paper is therefore to examine the interplay between thermal regime and basal sliding in greater detail, with a view to substantiating MacAyeal’s vision of the Hudson Strait surges.
2. Ice-Sheet Model
In seeking a simple but realistic model for ice-sheet motion, we make use of a lumped parameter version of a boundary-layer theory developed by Fowler (1992), based on an original idea due to Nye (1959), and subsequently developed by Lliboutry (1981). This is that large ice sheets have dynamics which can be characterized by “shear layers” near the base where the shear is largest due to the high stresses and high temperatures there. In addition, relatively high velocities cause thermal gradients to be elevated in a thermal boundary layer near base, with the shear layer lying inside this. Because of the boundary-layer nature of the flow, it is possible to produce a parameterized model which encapsulates the dynamics.
Physics of the model
We begin by describing the physical processes and their effects; the mathematical detail follows this. The horizontal motion of ice sheets is driven by the shear stresses at depth generated by the surface slope of the ice (i.e. ice flows down the “hill” generated by its surface). This flow is thermally activated by the dependence of viscosity on temperature, and this dependence is sufficiently strong that one can effectively consider that all the shear takes place where the ice is warmest, i.e. near the bed. We call the region near the bed where the shear is concentrated a shear layer and, in this layer, the horizontal velocity u rises from its value u b (the basal sliding velocity) at the bed to a far field (i.e. far from the bed) value u∞, which is esssentially independent of depth. Also in this layer, frictional heat causes an imbalance between the basal heat flux from the bed to the ice q b, and the heat flux delivered from the shear layer to the ice above, q T. Approximate relationships between these variables can be derived.
An important constituent of the model is the determination of the heat flux into the plug flow of ice above the shear layer. If u∞, is very small, this should be conductive (and thus q T ∝ h −1 where h is ice thickness) but, if u∞ is very large, we expect a thin thermal boundary layer to exist near the base, though we can show that the shear layer lies within this thermal boundary layer (see Fig. 1). In this case, we find and in our simple model we simply add the conductive to the convectively induced heat transfer.
We thus propose relations between q T, q b, u ∞, u b and h based on the dynamics of ice, which involve the basal shear stress τ and the basal temperature T b, We have four such relations and, in order to find an expression for the horizontal ice flux Q « hu„ we require two further relations amongst the variables. These follow from our description of the sliding process at the base. If the base is frozen, then u b = 0 and q b is the prescribed geothermal heat flux. However, if the base is temperate, then we prescribe T b, = T m, the melting temperature, and a suitable sliding law relating τ to u b. This relation also involves the effective pressure N (ice-overburden pressure minus subglacial drainage water pressure), which must be determined by a suitable drainage theory. In particular, Ν will depend on the water supply Q w, which is itself related to the ice flow, being proportional to the excess heat generated at the bed (geothermal plus friction minus heat lost to the cold ice above)
Fig. 1. Schematic representation of an ice sheet, based on conditions in the Laurentide ice sheet. The location of basal shear and thermal boundary layers is also indicated.
The sliding law is crucial to our conclusion that a multi-valued ice flux vs depth relation is possible and indeed likely. When the bed is temperate, sliding occurs but, if the velocity is low, then little frictional heat is generated and there is little water flow. Crucial to our results is the notion that the effective pressure N decreases as water supply increases. This says that as water pressure increases, the hydraulic flow increases, which makes sense, but is also contrary to the drainage characteristics of Röthlisberger channels, the hypothesis is based on recent theories of drainage over deformable sediment beds. With this assumption, low values of Q w correspond to high N, or low water pressures, and for sliding of ice over a deformable till, we expect that the basal velocity u b increases as water pressure increases, i.e. N decreases (the till becomes more mobile). Thus, low u b, implies low Q w implies high N, which is consistent with the assumption of low u b.
However, as the ice thickness increases, the conductive heat flux into the ice decreases and the water supply increases. An increase of Q w leads to a decrease of N and therefore enhances sliding. A positive feed-back exists, because the increased sliding velocity causes enhanced water production via frictional heat production, and in fact we find that the result is a run-away, in the sense that, for sufficiently large ice thickness, the “slow” mode (low u b, low Q w, high N) is not viable and a fast branch (high u b, high Q w low N) takes over. As a consequence, we find a multiple-valued relation for ice flux Q as a function of ice thickness h.
Mathematical model
The geometry we consider is shown in Figure 1. To fix ideas, we consider specifically a test section cf the Laurentide ice sheet which extends from an ice divide down the outlet of Hudson Strait, but the concept is applicable elsewhere such as for flowlines extending south into the ice lobes which terminated the ice sheet in the Great Plains area of North America. Τo the south of the divide, the Canadian shield provides a hard base and we conceptualize the ice here as being relatively stable. North of the divide, the test section we consider lies over Hudson Bay and is channelized into the ice stream in Hudson Strait. Let the horizontal ice velocity be u and the vertical coordinate be z (with z = 0 as the base). Glen’s flow law may be written in the approximate form appropriate for a shearing flow
where we use the Frank-Kamenetskii (1955) approximation (for) to simplify the Arrhenius term exp(–E/RT) Here, A m is the rate factor at the pressure-melting temperature T m, Ε is the activation energy, R is the gas constant and τ is the basal shear stress, which can be taken as constant through the shear layer.
In the thermal boundary layer, heat advection balances conduction but, within the thinner basal shear layer, we can ignore heat advection, so that the temperature is given approximately by
where subscripts denote partial derivatives, thus uz = ∂u/∂z, etc.; hence, integrating with respect to z,
where q T is the heat flux delivered from the shear layer to the thermal boundary layer of the ice sheet, q b is the basal heat flux at z = 0, and u∞, and u b are the far-field and basal velocities.
An integral relation between heat flux and τ is obtained by using Equation (2.1) in Equation (2.2). Multiplying Equation (2.2) by Tz and integrating, we obtain
where we have used the fact that the exponential term tends to asymptotically small values when z extends beyond the shear layer (Fowler, 1992).
The basal shear stress is given by τ = pgh sin α, where ρ is density, g is gravity, h is ice depth and sin α is the surface slope. In our parametric approach, we take h as divide height, l as the length of the test section (the diameter of Hudson Bay) and we take sin α ≈ h/l, so that
In reality, u∞ must change from zero at the divide to a large value at the outlet In our lumped approach, it is taken to be an average value. Then the mean depth is h/2, the area of the test section is ½lh, while the ice-drainage flux is ½hu∞. It follows that mass conservation can be expressed as , whence
where a is the average accumulation rate over Hudson Bay. In general, it may depend on the atmospheric adiabatic lapse rate, in which case a would be a function of the mean ice depth h/2.
In the thermal boundary layer, shear heating is negligible (it is concentrated in the warmer shear layer), so that (since the boundary-layer temperature relaxes quickly to a steady state)
where x is the horizontal coordinate (see Fowler (1992) for further details). The temperature profile necessarily evolves with x and it is indeed possible (with u = u∞, w = 0, T b = constant, for example) to find a similarity solution for Τ (indeed, this can even he done if u∞ = u∞(x)). With Τ = T b on z = 0 and Τ = T∞, as z → ∞ (i e. outside the boundary layer), the similarity solution is
(and f(ξ) = erfc(ξ)), whence the average heat flux delivered to the thermal boundary layer (from the shear layer) is . Since we thus have
where we simply add the second term to include the heat flux due to conduction when h or u∞ is small (T a is the surface temperature). Notice that a simpler approach to obtain qT would be to take a linear temperature profile through the thermal boundary layer thickness δT, thus the first term would q T = k(Tb —T∞)δT, where from Equation (2.7), we would dimensionally estimate . This gives the same result without the pre-multiplicative factor and is probably just as good.
The Equations (2.3), (2.4), (2.5), (2.6) and (2.8) provide five equations for the unknowns u∞, qb, τ, h and qT. However, the variables u b, T b, and T∞, remain to be prescribed. The temperature T∞ at the edge of the thermal boundary layer is derived from the divide temperature at an earlier time. Since heat conduction may be small (as measured by the smallness of δT compared to the ice depth), the surface temperature is simply advected downwards. However, we will assume for simplicity that the surface (divide) temperature is constant, and tbus we prescribe
Basal conditions
Finally, we constitute the basal values of T and u. There are four possible states of the basal system which we will consider.
(i) Frozen base
We prescribe
where G is the prescribed geothermal heat flux.
(ii) Sub-temperate base
When the base reaches the pressure-melting point, sliding is initiated, but at first there is no net water production (we ignore any possible flux from the surface). In this case, we have
where u b* is the “fully developed” sliding velocity detailed next.
(iii) Temperate base
When 0 < q b < G ub, then “fully developed” sliding occurs. This is the sliding commonly prescribed using a “sliding law” of the form τ - f(ub , N), where Ν a the effective basal water pressure (= pi – pw , where pi is ice-overburden pressure and p w is basal water pressure). Sliding theories for hard beds developed in the 1960s and 1970s led to that enunciated by Lliboutry (1979). In its simplest expression, this gives τ/N = f(u b/Nn ) (Lliboutry, 1987), where n is the exponent in Glen’s law. When f is a power law, we have the relation
(consistent with Lliboutry’s result nr + s = 1), derived theoretically by Fowler (1987a), with typical values r ≈ 1/4, s ≈ 1/4. Laboratory experimental results of Budd and others (1979) are also consistent with Equation (2.12), with (at high N) r = s = 1/3, (1 bar = 105 Pa, 1 year = 3 × 107 s); at low N, more appropriate values were r = s = 1. Bindschadler (1983) also found some consistency with Equation (2.12) using measurements on Variegated Glacier and indeed found a similar value of c.
However, these results apply for ice sliding over a lumpy bed, where the principal resistance is due to the flow past the bumps. For the situation of relevance here, where the principal resistance is due to deformation of a layer of saturated subglacial sediment, a simple sliding Law is just
where η is the viscosity of the till and h T is its depth. The effective water-pressure dependence arises through influence on the till viscosity. A law proposed by Boulton and Hindmarsh (1987) based on seven measurements at Breidamerkurjökull is in which case we could write Equation (2.13) in the form of Equation (2.12), with r = 1/a ∕ 0.75, s = b/a ∕ 1.4 and c = (Ah T)−r. Boulton’s law is controversial, however, and values of τ and N appropriate to Ice Stream B (τ = 0.15 bar, N = 0.5 bar) give a velocity of 9 m year−1 for till of thickness 8 m, as opposed to the 500 m year−1 observed. Alley and others (1987) estimated the viscosity of till under Ice Stream Β as 1010 Pa s.
Despite concerns with the Boulton-Hindmarsh rheology, it is the simplest type of rhealogy which mimics the which increasing water pressure has, namely to decrease till viscosity. However, it suffers one important conceptual limitation. Since the actual deformation of till involves particles moving round and past each other, dilation is necessary to accommodate the flow. Therefore, it is only realistic to suppose that the viscosity tends to zero as Ν → 0 for an unconfined flow (when it dilates to a slurry). For flow in a confined channel (such as a till layer), deformation causes dilation and consequently normal stresses, so that it is properly a non-Newtonian material. The simplest way in which a non-zero viscosity at N = 0 can be included is to replace N in Equation (2.12) by N + N 0, so that
There is little to constrain such an N 0 but a least-squares fit to Boulton and Hindmarsh’s limited data yields N 0 = 0.045 bar, comparable to the yield stress, and suggests that normally a relationship such as Equation (2.12) may although there is very little to constrain the choice of c, r and s.
The effective pressure Ν depends on the volume flux Q of water through the drainage system. Here we are concerned with the drainage from a till-based system, for which Walder and Fowler (1994) suggested that a system of “canals” scoured from the till will exist, within which Ν is essentially a (decreasing) function of Q; since Q increases downstream, Ν will decrease and so also will the viscosity, as is inferred for Ice Stream B (Alley and others, 1987). Taking into account basal topography, we then expect a reticular network of interwoven canals with spacing controlled by the wavelength of topographic variability.
The existence of a basal system of canals is by no means certain. A modified viewpoint is that the water simply ponds and drains as a patchy but connected water film (Alley, 1989). A sliding law based on the corresponding hydraulic system has been proposed and used by Alley (1990); he also found multiple dynamic states and indeed the crucial point is that the basal velocity increases with the water supply. Thus, while the precise quantitative nature of the basal hydrology is of importance, from the point of view of the sliding law, it is only a detail.
The excess heat production at the bed is G + τu b - q b and this produces a downward melt velocity of (G + τu b – q b)/ρ w L where L is latent heat and ρ w is the density of water. If the spacing of an assumed reticular drainage network is w d then the flux of water per channel is given (in parameterized form) by
and Ν is given by a relation such as that of Walder and Fowler (1994). For canals incised into relatively fine-grained sediments, they put (with sin α ≈ h/l)
where D s is a “characteristic” size of suspended grains in the water flow and . A representative value of D s might be D s = 0.1 mm = 10−4 m, in which case a slope of 10−3 and a water flux of 1 m3 s−1 gives N ≈ 0.3 bar. It must he emphasized that Equation (2.16) u possibly a crude representation and the relationship itself describes only one possible drainage scenario; others are arguably just as plausible.
In using Equation (2.12) or Equation (2.14), together with Equations (2.15) and (2.16), to prescribe τ in terms of u b, we have to consider what happens as Q w → 0. Equation (2.16) cannot strictly apply for small Q w since it implies that Ν → ∞ as Q w, → 0. In fact, for very small values of Q w, Ν will increase rapidly from zero (flotation) at the onset of channelized flow to a value , where is a critical effective pressure relating creep properties of till to those of ice (Walder and Fowler, 1994). Based on the Boulton-Hindmarsh rheology, , while for less viscous till, could be higher. At any rate, this suggests that the correct limiting value of Ν to use Q w → 0 is . A practical way to model this is to define via
and then replace Q w, in Equation (2.16) by Since we expect to be “small”, there is little difference made in doing this. We then use Equation (2.12) and the fully developed sliding velocity u b* is defined by
Now in fact, since is “small” and is “large”, we have that ub * is also small, and in practice we can ignore the sub-temperate region ((ii)above) altogether. We thus take the sliding law in the form
while Qw , > 0, and note that as Qw , → 0, then (formally) Ν → ∞ and ub → 0. The second of these allows for a continuous join to the frozen regime where ub = 0.
Molten base
A further regime is possible, which we simply ignore. If the heat flux to the ice reaches zero, i.e. qb = 0, then a basal layer of temperate ice forms (e.g. Hutter and others, 1988), and in this case the whole ice-dynamics model must change, since the ice viscosity depends on moisture content (Lliboutry, 1976). In such a situation, moisture is produced within the ice and drainage of this water to the basal system can only enhance the instability we envisage. However, it is not specifically included in our model and is in fact unlikely to occur with the present assumptions.
3. Analysis
The Equations (2.3), (2.4), (2.5), (2.6), (2.8) and (2.9), together with the boundary conditions (2.10) or (2.14), (2.15) and (2.16), are scaled by choosing
(that is to say, we define dimensionless variables u*, τ*, etc. by writing in the sequel, we drop the asterisks for convenience). We choose balances via
From these we derive
and for values
we find
These values are approximately representative of fast-flowing ice streams and are perhaps relevant to the conditions of a temperate bed. The corresponding dimensionless model equations are
together with
or
The parameters are defined by
and we use values
We then have
and μ3 is rather uncertain. Comparison of observed stresses ∼ 0.15 bar and velocities 500 m year−1 at Ice Stream Β would suggest μ3 ∼ 5, for example.
Cold approximation
We use the fact that μ1 ≪ 1 to approximate the relations in Equations (3.6) and (3.7) We have that , and then
where we take T a = −1 as constant (more generally it might he a function of h). Thus and the ice flux is
It increases monotonically with h until T b = 0, which is at h ≈ δ, and where . Thus, the ice is virtually stagnant if the base is frozen.
Temperate approximation
When Tb reaches zero at the base, we adopt the temperate conditions (3.8), which gives (with
These imply, successively, say, say, and thus we have to solve
These equations are valid providing (the dimensionless version of corresponding to a temperate base). From Equations (3.14), and . The cold-temperate transition is when ub = 0, so the temperate conditions are valid if ub ≥ 0, i.e. equivalently (since also q ≈ 1 when and, since the transition occurs when h ≈ δ, this is u > μ1δ6.
The temperate-molten transition is when q = 0. However, the third of Equations (3.14) implies that q > 0, so that in this model, basal molten zones do not occur.
We now consider Equation (3.15). Denote the left-and righthand sides as P(h, u),R(h, u) and also put , thus . For given u > 0, p and hence Ρ is monotone increasing with positive h, while R is monotone decreasing (we assume s < 12) Therefore, there is a unique positive h for each value of u.
Alternatively, fix a value of h > 0. The function p is a quadratic in u½ with no roots if two positive roots if δ < h < h c and one positive root if h < δ. On the other hand, R increases with u and, since then if r > s/3 (as we assume) R > Ρ for large enough (since as). By considering intersections of P and R, we then see that for large μ3 there will be a unique value of u (given h) for h > δ, while for small μ3,u(h) will be multi-valued. The two different possibilities are illustrated in Figure 2, using values given in Equation (3.11) for μ3 = 1.0 and μ3 = 20.0. For the values in Equation (3.11) we find that a multi-valued flux/depth relation occurs if μ3 < μ* ≈ 16.98. The lefthand nose of the curve, at h = h m, is to the left of the cold-temperate transition (i.e. h ≈ δ) if .
An analytic approximation to the relation can be obtained when δ3 is small. In that case, we can approximate the solution of Equations (3.15) as follows.
Fig. 2. Two flux curves for he parameters in Equation (3.11), together with μ3 = 1 (multi-valued) and μ3 = 20 (single-valued). With μ3 = 1, the slow nose is at h = h c = 2.86, and the fast nose is at h = h m = 1.44. These correspond to depths of 2145 and 1080 m, respectively. At hc, the velocity is 0.034 (∼ 17 m year−1), while at the same value on the fast branch it is 1.18 × 103 (∼ 590 km yeaf−1 or 2 cm s−1). At h m, the velocity is 5.24 (∼ 2.6 km year−1). One can show that the consequent surge duration is approximately 0.03, corresponding to 116 years (based on Equation (3.5)).
The fast branch is
where
for example, if r = s = 1/3, then , so we expect Equauons (3.17) to be accurate even for moderately low μ3. The lower and middle branches are given by
so the nose is at h = h c as defined before. At the upper (i.e. larger Q) nose, h = h m, u ≫ 1 and Equation (3.15) is approximated by
and it follows from this that at the nose
where
If r = s = 1/3, then υ ≈ 27, so that again we expect hc < δ for moderately low μ3. In particular, the the surging velocity depends sensitively on μ3.
We see that the sliding law can give rise to a phenomenon which we may term hydraulic run-away. If h reaches hc then there is an unstable transition to a fast-sliding mode by the enhanced production of meltwater, although the ice above remains cold. We thus have a mechanism for relaxation oscillations of the ice sheet.
4. Discussion
The parameterized model presented here is of course crude but nevertheless it encapsulates the basic physics of sliding over deformable sediments and shows that thermally regulated surging is a feasible phenomenon. Hydraulic run-away relies on two features. First, sliding is enhanced by increased water production. Thus, the increase in frictional heating with velocity is a positive feed-back. Countering this is the negative feed-back due to the increase in heat flux from the base with increased velocity. At low u b, the second of these dominates, since but at higher values of ub frictional heating predominates and a fast-sliding mode becomes possible.
As applications, we consider the Hudson Strait mega-surges to be a prime example (MacAyeal, 1993a, b). We should also like to consider whether thermally regulated surging of Trapridge Glacier can be looked at from the same point of view. In this case, the base is temperate and basal water is evacuated subglacially as ground-water flow. Currently, the glacier would reside on the slow branch. A difference in the valley-glacier case is that basal meltwater is likely to be mainly surface-derived and whether the positive feed-back of ice velocity on water pressure can operate is unclear.
Lastly, we consider sheet flow over a deformable bed, i.e. where the width of the flow is much greater than the depth and where the flow is fed from a stable catchment. we have in mind the Siple Coast of Antarctica. Suppose upstream that Q(h) is a monotone function (e.g. due to lack of a basal till cover) but that downstream we have Q = Q(h) as in Figure 2. With conditions of continuity of Q and h at the join, we conjecture that, if this Q value lies on the unstable part of the downstream Q(h) curve, a uniform flow is laterally unstable and will break up into fast- and slow-moving parts, and that these are precisely ice streams A study of this postulated instability is deferred to future work.
Acknowledgements
We are grateful to R. Alley and D MacAyeal, in particular, for sharing their insights with us. C. Johnson acknowledges the award of a U.K. NERC CASE (with the British Antarctic Survey).
