1.Notation

Symbol Relation Description

D₀ 3 Deformation-rate magnitude (1 year⁻¹)

D 3 Deformation-rate tensor

C 3 Dimensionless temperature

C_F 3 Prescribed C at H = 0

E 15 Dissipation term

F 1 Dimensionless bed profile

G₁ 13 “Advection” coefficient

G₂ 13 “Diffusion” coefficient

H 1 Dimensionless ice-sheet profile

J 3 Second deviator stress invariant

L 16 Dimensionless lapse rate

M 2 Dimensionless margin-translation speed

Q 1 Dimensionless accumulation

R 8 Snow line

Τ 3 Dimensionless time

U 1 Dimensionless horizontal velocity

W 1 Dimensionless vertical velocity

X 1 Dimensionless horizontal position

Z 1 Dimensionless vertical position

S 1 Deviatoric stress tensor

a 3 Temperature-dependent rate factor

c 3 Temperature

d₀ 1 Depth scale

f 1 Bed profile

g 7 Rheological function for ice

h 1 Ice-sheet profile

m 23 Constant in bed-roughness model

n 19 Constant in internal-energy model

p 1 Pressure

q 1 Accumulation

q_m 1 Accumulation magnitude

u 1 Horizontal velocity

w 1 Vertical velocity

x 1 Horizontal position

z 1 Vertical position

Δ 1 Temperature magnitude (20 K)

18 Internal energy

7

17 Geothermal heat flux

ψ 15 Specific heat capacity of ice

β 15 Dimensionless diffusivity

20 Latent-heat representation

γ 19 Constant in internal-energy model

ε 1 Aspect ratio

k 7 ε/s

λ 15 Thermal conductivity of ice

μ 8 Bed roughness

μ_max 8 Maximum temperature-dependent μ

μ_min 8 Minimum temperature-dependent μ

θ 7 k²

ρ 1 Density of ice

σ 1 Stress tensor

σ₀ 1 Stress magnitude = 10⁵ Pa

v 3

υ 6

ω 3 Rheological function for ice

2. Introduction

Reference GlenGlen’s (1955) demonstration that ice is a thermo-viscous fluid over glacial time-scales had the implication that the evolution of temperature within an ice sheet and the flow of ice sheets could be a strongly coupled problem. A contemporary analysis of Reference Robin deRobin (1955) showed how ice-sheet flow patterns might alter temperature fields in the presence of a significant lapse rate (vertical atmospheric temperature gradient). An excellent review of progress over intervening years may be found in Reference PatersonPaterson (1981); since then, advances have been made in the formulation and analysis of thermo-mechanically coupled ice sheets.

Morland (1984) and Reference HutterHutter (1983) investigated the properties of thermal boundary conditions necessary for a reduced model. By adopting suitable scalings they showed that O(l) variation of the temperature field did not imply incorporation of extra terms into the momentum balance, and they also showed clearly in formal terms how the ice sheet was divided into an upper zone, where heat transport was predominantly by advection, and a basal boundary layer, rather thin for thick ice sheets, where advection is balanced by conduction.

By varying a prescribed temperature field, Morland and Smith (1984) demonstrated a significant effect on the geometry of ice sheets. However, Reference Hutter, Yakowitz and SzidarovszkyHutter and others (1986, 1987) concluded that temperature does not significantly affect the geometry when the coupled problem is solved. One resolution of this inconsistency is to suggest that the ice-sheet thermo-mechanical system posesses a self-regulatory structure, which disallows those temperature distributions which produced these different results.

While the steady-state solution may have a self-regulatory structure, time-dependent solutions need not behave so; for example, Reference Veen van der and Oerlemansvan der Veen and Oerlemans (1984) have suggested that there may be a bifurcation in the evolution of thermo-mechanically coupled ice sheets. The purpose of this paper is to investigate the time-dependent thermo-mechanically coupled behaviour of an ice sheet, solving for the evolution of the temperature field with no approximations in order to illustrate the modes of operation of the system.

Many early approximate time-dependent models which deal with the evolution of models of ice sheets in response to various external inputs and boundary conditions have been presented. Only one of them, due to Jensen (1977), computed coupling of temperature with rheology. This algorithm uses an explicit time-marching scheme in which the maximum permissible time-step is often substantially less than 1 year.

Recently, an algorithm using an implicit time-marching scheme has been developed by Reference Hindmarsh, Morland, Boulton and HutterHindmarsh and others (1987) and Reference Hindmarsh and HutterHindmarsh and Hutter (1988), where it is possible to take time-steps of 500 year. This method is between one and two orders of magnitude faster than explicit methods, and it makes possible the computation of the thermo-mechanical evolution of ice sheets through cycles of 100 000 year feasible without the use of a super-computer.

We use the Morland-Hutter (Reference HutterHutter, 1983; Morland, 1984) reduced model to describe the mechanical behaviour of the ice sheet, and to compute the evolution of temperature using the complete dissipation-advection-diffusion equation. The surface-boundary condition is prescribed temperature, while the basal boundary condition is geothermal heat flux and basal frictional heat.

Hindmarsh and Hutter found that in every case part of the base of the ice sheet reached the melting point. In their algorithm, no account is taken of the change in the physical basal boundary condition for the energy balance which occurs at this temperature; heat flowing in is used to melt ice, either to produce basal melting or else water between ice crystals. This logical branch is easy to represent when (inefficient) explicit methods are being used, but would almost certainly have a deleterious effect on the Newton-Raphson iteration scheme used by Hindmarsh and Hutter. Instead, so as to maintain sub-zero temperatures, we have adopted a method similar in intent to the enthalpy method (e.g. Reference Elliot, Ockendon and MathElliot and Ockendon, 1982), though a more complete resolution of the problem will require the development of new models (e.g. Hutter, 1982; Reference FowlerFowler, 1984; Hutter and others, 1988).

The aim of this paper is to investigate the modes of operation of thermo-mechanically coupled ice sheets. Section 3 presents the governing equations. Section 4 discusses the boundary conditions. In section 5, we use the Morland-Hutter model to predict scales of operation of the various processes, and we illustrate our predictions in section 6. A discussion and conclusion are presented in sections 7 and 8.

3. Governing Equations and Algorithm

Consider the plane, gravity-driven flow of an ice sheet (Fig. 1). Let x be the horizontal coordinate, z the vertical coordinate in a plane Cartesian system, z pointing upwards.

Fig. 1. Illustrating the coordinate system used.

The lines z = h(x,t) and z = f(x,t) are the equations describing the bounds of the domain, the domain being confined to the region h(x,t) − f(x,t) ≥ 0, with equality defining the margins. Both surfaces can freely deform in time, though we only treat the case where the deformation of f is prescribed. The problem is to determine the profile z = h(x,t) as a function of position and time, the velocity distribution (u,w) in the directions (x,z), respectively, and the evolution of the temperature c.

We follow Morland (1984) and use a scaling

where d₀ is a typical thickness for the ice sheet, q_m is the accumulation magnitude, and ε ‹‹ 1 is a characteristic ice-sheet aspect ratio. Thus, X,Z,H,F are the dimensionless equivalents of x,z,h,f and similarly U,W are the dimensionless counterparts of the velocities u,w. Time t is scaled by d₀/q_m into T, which implies a dimensionless strain-rate (D) with unit q_m/d₀. Stresses (σ), deviatone stresses S, and pressure p are scaled by pgd₀ where ρ is the density of ice and g the acceleration due to gravity. Henceforth, we will always refer to the normalized form of the stresses and strains. Temperature c is normalized by C = c/Δ, Δ = 20 Κ, where c is in °C.

The margin translation speed M is found to be

where the superscript A = L,R indicates evaluation at the left and right margins.

We use a viscous relation with a temperature-dependent rate factor

where J = trS²/s², s = σ₀/ρgd₀, v = σ₀q_m/D₀ρgd^{2 0} with σ₀ = 10⁵ Pa, D₀ = 1 a⁻¹, and adopt the Smith-Morland (1981) model functions

α₁ = 0.7242, β₁ = 11.9567, α₂ = 0.3438, β₁ = 2.9494.

Satisfying mechanical balance to lead order in ε yields

where the subscripts denote tensor components. From Equations (3) and (5) we obtain

where

The basal boundary condition is provided by a sliding law (see section 4) for a discussion of the type proposed by Morland and others (1984), viz.

where p =k(h-f) and μ( p ) is a function representing bed roughness; using this with Equations (6) and (8) we obtain

the second equality only being valid to leading order in ε. By integrating Equations (6) with (7) subject to the basal condition (9), we obtain U/(Z),

Then by differentiating U with respect to X, using the incompressibility condition

and the kinematic basal condition

where B* = B*(X,H,T) is the rate at which ice is melting at the base, an expression for W(Z) is obtained of the form

By substituting Ζ = Η in relations (10) and (11) and substituting these expressions in the kinematic surface condition

where Q* = Q*(X,H,T) is the rate of accumulation on the surface, an evolution equation for Η

is obtained, where

dHdF ) (dHdFdC

and it may be shown that G₂(·,·,0,·,·) = 0, G₁,(0,·,·,·,·,·) = 0.

The temperature C evolves according to the spatially two-dimensional advection-dissipation-diffusion equation

Here, β is a constant diffusion coefficient, E a reaction term, and ε a typical aspect ratio introduced during the scaling. With the above scaling and with Equation (3) we have (see Morland (1984) for the derivation):

where ψ = 2 × 10³ J kg⁻¹ Κ⁻¹ is the specific heat capacity of ice and λ = 7 × 10⁷ J Κ⁻¹ m⁻¹ a⁻¹ is the thermal conductivity of ice. We consider the case symmetric about the divide, and select ε = 0.005 (thus, horizontal conduction may be ignored, see relation (14)), d₀ = 2000 m, q_m = 1m a⁻¹, which gives β = 0.021 and a unit time-step of 2000 year. The reaction term E arising from internal dissipation is as given in Equation (15). Given U, H, and

as being all O(1), then its magnitude at most is about 0.5.

Boundary conditions consist of (i) a Dirichlet condition applied to the upper surface according to

where C_F is a constant and L is the lapse rate (with the scaling we use, the dimensionless lapse rate is equal to the lapse rate in K/100m); and (ii) a Neumann condition applied to the base representing a constant geothermal heat flux plus the heating due to sliding, given by the product of the shear stress at the base y(Z = F) and the sliding velocity U(Z = F). In scaled form, we find that the basal Neumann condition is given by

where

is the geothermal heat flux = 1.576 × 10⁶ J m⁻² a⁻¹= 0.05 W m⁻².

It will be seen that the dimensionless basal power generation Ω = U(H - F)∂H/∂X is O(1) at the margin, though longitudinal gradients ∂Ω/∂X are a bit above this. With the scaling used in Equation (17), the quantity ∂² C/∂Z∂X may reach O(1/ε), implying similar magnitudes for longitudinal derivatives of the temperature and the rate factor a(C) (relation (4)). Inspection of relationship (14) shows that horizontal conduction may still be ignored; it may also be shown that the reduced model momentum balance is still respected, though with error O(ε) rather than Ο(ε²). We must take this as a warning that there are likely to be problems at the margin.

The ice-sheet inception is from a snow-pack where the snow line Ζ = R(X,t) intercepts with the topography (Fig. 2). Above the snow line, ice accumulates while below it any ice will melt. For our calculations we introduce the accumulation function Q(H*), where Η* = Η − R(X,t) defined by

Fig. 2. Contours of the accumulation Q*(X,Z), and the snow line Ζ = R(X). The ice sheet is symmetric about X = 1 and is initiated from the area X > 0.95. In areas above the Q = 0.5 contour the accumulation is 0.5 m/a.

with the function Q(H*) for 0 < H* < 0.25 being defined by the cubic fitted through these points with slopes defined above. As an initial temperature at Τ = 0, we use C₀(X,Z) = C_F + LZ.

Hindmarsh and Hutter devised a finite-difference fixed domain method in which the ice sheet is mapped on to a rectangle; thus, stretched coordinates are used in both directions in contrast to Jensen (1977), who used a fixed domain method in the vertical domain only. A consistent second-order accurate Crank-Nicholson scheme is used to step the calculation through time. Multi-point up-winding is used to eliminate the spurious modes associated with hyperbolic equations. The resulting equations are non-linear and are solved by a Newton-Raphson scheme which includes expansions with respect to the mapping terms. The resulting matrix has a nested bordered form and is inverted by a pre-conditioned conjugate-gradient scheme.

We use a discretization of 21 points in the vertical and 31 in the horizontal. The time-step is not constant, as at points where the variables are changing rapidly and convergence of the Newton-Raphson iteration becomes difficult, the algorithm automatically cuts the time-step. In all calculations the maximum time-step was set to ΔT = 0.25, corresponding to a time-step of 500 year. These discretizations are not fine enough to prevent cumulative error over glacial time-scales from becoming very large; however, the purpose of the calculations is to illustrate the modes of operation of ice sheets which do not require extremely accurate calculations.

Hindmarsh and Hutter found that in most instances the ice reached the melting point locally. Thus, the preceding model breaks down unless amendments are introduced. Physically, two situations may arise (Reference FowlerFowler, 1980): either a temperate volume will form with ice and water co-existing, or the ice will only reach melting at the very base. In the former situation, the enthalpy density of the ice is increased by increasing the water content. This may be redistributed (Hutter, 1982; Reference FowlerFowler, 1984) by movement of the interstitial water. In the latter case, the enthalpy content of the ice does not increase, but the water is supposed to drain in some unspecified way.

While the latter situation may arise beneath an ice sheet, we do not have a formal guarantee that it will do so under all the situations considered. However, modelling two-phase flow requires computation of the water content using constitutive relations which are not yet well established (Hutter, 1982; Hutter and others, 1988). It also introduces a new moving-boundary problem within the (moving) ice-sheet domain, and thus changes the structure of the problem, making it even more complex.

The simplest remedy is to adopt the “enthalpy” methods (Reference ElliotElliot and Ockendon, 1982) where the effect of the latent heat is represented by an increased specific heat over a narrow temperature range. Normally, the specific heat is set to a value appropriate to the molten phase at temperatures above the melting point, but in our case this is not meaningful, and we let the specific heat capacity tend to infinity as the temperature approaches the melting point. This is equivalent to letting the water drain from the ice sheet without any adjustment to the continuity equation, and provides an infinite sink for heat entering the system.

The formal specification is as follows. Let ∧ represent the internal energy of the ice. Then, for simplicity considering vertical diffusion and dissipation only, the variation of ∧(C) with time following the motion is given, using Equations (14) and (15), by

We select for the variation of enthalpy with temperature

where γ and n are arbitrary constants and 13.65 is the scaled melting temperature of ice. The first term in Equation (19) represents normal linear variation, while the last is a hyperbolic term which tends to infinity as C tends to zero; this prevents temperatures from rising above zero.

Differentiating Equation (19) with respect to temperature (remembering we are differentiating the modulus of a negative number) gives d∧/dC

and, by defining a new function by

substitution into Equation (18) gives

This relation replaces Equation (14).

In our calculations we have used two models (Fig. 3, graphs A and B) for the dependence of

, one requiring

Fig. 3. Pseudo-latent heat

(A) Model 1, (B) Model 2, and (C) basal roughness μ(C)/μ_max plotted against scaled temperature.

rather less heat to raise temperatures near melting point than the other. The validity of the approximations depends upon the context and is discussed further below.

4. The Basal Boundary Condition and Movement of The Margins

In this paper we have used the form (i.e. relation (8)) of the sliding law introduced by Morland and others (1984) which permits sliding motion at the margin despite the presence of zero shear stress, and requires the bed resistance to disappear with overburden. A mathematically similar form has been used by Reference Budd and SmithBudd and Smith (1981), while a (linearized) Reference WeertmanWeertman (1957) form,

for example, does not permit sliding at the margin. The functional form of these sliding laws has experimental verification (Morland and others (1984) from the Greenland ice sheet; Reference Budd, Keage and BlundyBudd and others (1979) from laboratory data), but we do not wish to present it as based upon any greater physical significance than that of well-posedness, though it can be readily changed so that effective pressure replaces pressure. A sliding law does not necessarily imply a discrete interface at which slip is occurring, but merely implies some enhanced basal motion, possibly on sub-grid dimensions. The theoretical complexity of constructing sliding laws for ice sheets has been illustrated by Reference LliboutryLliboutry (1987a), who postulated the necessity for three scales of re-normalization.

Ice-sheet or glacier advance at the margin may be due to sliding or by ice over-riding. The latter requires an infinite slope at the margin and a significant role for longitudinal stresses in the momentum balance. Reduced ice-sheet models which do not permit sliding at the margin do, nonetheless, exhibit marginal advance. This may be justifiable (though it has not been rigorously proved or even attempted for glacier models) by appeal to a “weak” principle (Elliot and Ockendon, 1982) which says that in some average sense the margin motion is computed correctly. It should be noted, however, that the truncation error of such schemes, where it can be derived, is often the spatial discretization interval to a power less than one (e.g. Elliot, 1987), an error rather greater than that normally associated with discretizations of parabolic equations. This matter will be considered in a future paper.

By explicitly permitting margin sliding, it is possible to derive a formula for the translation speed of the margin (e.g. Hindmarsh and others, 1987; Hindmarsh and Hutter, 1988). This permits fixed domain methods to be used (Hindmarsh and others) which retain second-order accuracy. These algorithms, when run with specification of a very resistant bed, predict very little margin motion and very steep slopes at the margin, which is contrary to the result predicted by “weak” formulations. This poses the question as to whether the “weak” formulation is a convergent approximation to the real problem in the case of no marginal sliding. We shall address whether this is important through the results of our calculations.

Probably all glaciologists would agree that temperature affects “sliding” to some degree, even though there is little agreement as to the mechanisms. We have included calculations in which the roughness coefficient μ is temperature-dependent. The dependence is purely heuristic, and is designed to see the influence of longitudinal variations in the bed roughness that are as sharp as the Morland-Hutter model permits. We adopt the model equation

This yields maximum sliding resistance at low temperatures and gives a sliding resistance of μ_min at melting. The higher m is, the shorter the transition zone between cold sliding and warm sliding is. We use μ_min/μ_max = 0.05 and m = 7.5 (Fig. 3, graph C).

5. Rates of Operation and Magnitudes of The Components

The formulation of the model above shows that there are a large number of parameters and boundary conditions which may be varied in order to alter the response of the ice sheet. These include the distribution of atmospheric temperature, geothermal heat, ablation and accumulation, topography, isostatic response, and the sliding law. We have not yet analysed the response of the ice sheet to variations in all these parameters, nor is it obvious that this is the best way to analyse this system or to present its operation. Instead, we present what we believe to be the key features of the ice-sheet system, identify the time-scales over which they operate and illustrate the results with some of the calculations we have carried out.

These time-scales are deduced from relations (14) and (13). The coefficients of advection, conduction, and reaction are 1, β, and 1, respectively, each of these having the meaning of inverse dimensionless time. Thus, 1/β and 1 are dimensionless conduction or dissipation times.

The temperature-evolution equation has three types of contribution: conduction, advection, and dissipation. The conduction coefficient β = 0.021, which implies a conduction time of 50 or 100000 years with the scaling we have used here. Near the margin, the dissipation time is O(1), and, as will be seen, the contribution from basal friction is even greater. The result of this is that dissipation time-scales are 2000 year or less. An advection time may be estimated by computing the time it takes for heat to be advected into the basal layer where conduction is important (Reference MorlandMorland, 1984). For an ice sheet 2 km thick, we obtain time-scales of the order of 10000 years.

Reference Lliboutry, Martin and SandersonLliboutry (1987b, p. 217) postulated two thermal time-scales for ice sheets, one belonging to the basal boundary layer and one for the ice sheet as a whole. We argue that, since the boundary layer arises as a response to (a) dissipation terms or (b) advection (when advection and conduction are of the same magnitude; Reference MorlandMorland, 1984), it is not useful to attribute it a time-scale in its own right, but rather to let this time-scale arise naturally from the time-scales of the processes which cause it to be.

It is evident that the basal temperature of the ice will couple with the temperature field in the bedrock so as to alter the geothermal heat flux into the ice. We shall show that the problem is one-dimensional and then estimated likely deviations in the heat flux. The maximum depth of interference of a 100000 year temperature cycle into the Earth is no more than 3 km (Reference BirchBirch, 1948). Similar results have been obtained by Reference RitzRitz (1987) and Reference WaddingtonWaddington (1987). This cold wave would never reach the base of the lithosphere, and in consequence it is not correct to consider the problem as being steady. Since the depth of penetration is of the same magnitude as the ice-sheet thickness, we use the length scales of the ice-sheet problem to normalize the bed heat-transfer problem. The heat-transfer equation is then exactly as Equation (14) with no advection and any dissipation is due to radiative decay. The aspect ratio ε of the problem is the same and we may infer that horizontal conduction is negligible provided that the boundary condition (i.e. the ice temperature and the geothermal heat flux at the depth of maximum interference) has O(1) longitudinal second derivative. As will be seen, this is respected everywhere except perhaps at the margin. The coupled ice-bed problem may therefore be solved in the bed by integrations of the diffusion equation over the vertical dimension through time. If we ignore heat production in the lithosphere, we can make estimates of the gradients induced by possible temperature variations.

It may be shown (e.g. Reference Carslaw and JaegerCarslaw and Jaeger, 1959) that, if a sinusoidal surface-temperature variation is applied to an infinite medium, the maximum gradient of temperature [∂C/∂Z]_max is

where C_A is the amplitude of the applied temperature variation and T_L is the period of the variation (in dimensionless time) and we use the fact that the thermal conductivities of ice and rock are very similar. Using the Morland scaling and letting C_A = 0.5, i.e. an amplitude of 10 Κ. and assuming that a reasonable lower bound to time of the temperature change is 100 year, gives [∂C/∂Z]_max − 20 (cf. the scaled steady geothermal temperature gradient, 1.0). For time periods of 10 000 years, the maximum temperature gradient is 2. Thus, transient perturbations to the bed-temperature field will significantly affect the geothermal heat flux over short time periods. These transient values are comparable in magnitude with the maximum basal heating.

The height-evolution Equation (13) has three terms: a diffusive term, an advective term, and a source term. The diffusion term reaches its maximum at the divide, where it is O(1). The second derivative of height ∂²H/∂X² is usually at a minimum at the divide and is rarely above 0.2, implying a kinematical time-scale of 10 000 year or greater at the divide. At the margin, the diffusion disappears while the advection coefficient and slope are minimally O(1) and often approach 10, implying kinematical time-scales here of the order of hundreds or maximally 2000 year. The accumulation is O(1) at maximum, implying an accumulation time-scale of thousands of years, while marginal ablation is commonly O(10), implying time-scales of hundreds of years.

6. Calculations Illustrating The Modes of Operation

We may now see how these processes combine in determining the evolution of the temperature and surface profiles of an ice sheet. We initiate all the calculations from the snow-pack (Fig. 2). Unless stated otherwise, we use model 1 (Fig. 3, graph A) to represent the effect of latent heat, essentially because this was the smoother of the models and produced better performance from the algorithm.

Fig. 4a-d. Cross-sections of an ice sheet with μ = 1, ∂R/∂X = −0.25, C(Z) = −1.1 − Z. The characteristic pattern of a molten zone inside the margin evolves very early. Length unit = 400 km, height unit = 2 km.

7. Discussion

We summarize the results:

• (i) The temperature field just in from the margin is dissipation-driven. The heat production in most cases is sufficient to bring the zone just inside the margin to the melting point. Conduction is more efficient over these small thicknesses and advection of cold into the area is not important.

• (ii) The temperature field at the base of the divide area can adopt two modes. One is a cold-based mode, where the ice sheet is thicker and the low temperatures are due to the advection in cold ice from altitude. The other mode is warm-based and is associated with thin ice sheets and warm conditions. This is not a bifurcation in the strict sense (cf. Reference Veen van der and Oerlemansvan der Veen and Oerlemans, 1984), but there is a sensitive range of parameter values which stimulates a bifurcation.

We have not demonstrated a multiple steady-state solution and, if there were one, then the fastest agent of heat transfer to produce a change from cold-based to warm-based is conduction, which operates upon a far slower time-scale than the kinematical processes of adjustment at the divide, meaning that a collapse is extremely unlikely (cf. the initiation waves of Fowler (1987) which operate on time-scales of hours and days to start valley-glacier surges).

Ice sheets will grow until they have sufficient area in the ablation zone to balance the accumulation zone. The much higher velocities and the fact that ablation is an order of magnitude greater than accumulation means the areal extent of the accumulation is usually an order of magnitude greater than the ablation zone (see e.g. Fig. 2). The total ablation is, in consequence, quite sensitively dependent upon the geometry of the marginal areas, and this effect is amplified by the fact that ablation is greatest at the very margin.

The equilibrium margin slope is dependent upon bed roughness and marginal ablation (e.g. Hutter, 1983; Morland and Smith, 1984), while the margin-slope evolution is dependent upon spatial derivatives of these quantities. In our model, the marginal temperature is defined by a Dirichlet condition. In cases where the bed roughness is temperature-dependent, this therefore defines the marginal roughness (see the great variation caused by temperature on margin slope in Figure 9, for example).

There is a very large amount of evidence which indicates that processes of basal motion in marginal areas are rather different from those in core areas. We believe that the principal difference is the existence of large amounts of debris underneath the ice sheet (Reference Boulton and HindmarshBoulton and Hindmarsh, 1987) which can deform. There is evidence from many areas (Reference BoultonBoulton, 1987) that this was a significant process underneath at least the marginal areas of the Quaternary mid-latitude ice sheets. There is also evidence of foreland pushing by these ice sheets. We believe that this could only have occurred when the sediment was unfrozen; furthermore, the fact of foreland pushing implies that the stress distribution in the coupled ice-sheet/substrate system is not simply that of the shear-stress gradient balancing the pressure gradient; stresses are also being transmitted longitudinally.

The transmission of these stresses and the consequent strain-rates are dependent on the temperature of the foreland and in particular whether it is frozen or unfrozen. Geological evidence suggests that the foreland areas were permafrost areas. Exactly what this implies about the temperature of the ice at the margin is not clear as, according to our model, warm ice would be advected from the zones of high dissipation. (The fact that the ice cools down in our model is really a consequence of applying a forced boundary condition to a hyperbolic outlet.) The exact distribution of marginal temperature then becomes a very complex problem. Furthermore, the mid-latitude ice sheets were often flanked by very large pro-glacial lakes. These would have insulated the margins from the annual temperature variations and kept the area at the melting point.

The thrust of these arguments is that the processes of enhanced basal deformation at the margin is dependent, perhaps quite significantly in view of the differing strengths of frozen and unfrozen sediments, upon details of the margin stress and temperature fields. The governing equations for ice-sheet-profile evolution (Equation (13)) at the margin incorporate no diffusion, and are in principle sensitive to very small perturbations. The kinematical sensitivity of the margin is capable of amplifying these perturbations. It is likely that a more sophisticated model of margin mechanics would include diffusive or dispersive terms (as in Hutter, 1981) which might damp out sensitivity. We suggest that a more detailed model of the margin, including coupling with the bed, is necessary.

Fowler’s (1980) conjecture that only the base of an ice sheet will ever reach the melting point seems to be correct when vertical velocities are directed downwards. The bulging associated with compression will, however, advect warm ice upwards, meaning that the margins may become temperate. There is ample evidence to suggest that ice sheets do compress at the margin (Reference BoultonBoulton, 1987). This means that some kind of representation of temperate ice will be necessary in an ice-sheet model (see also the computations of Hutter and others (1986)). In dissipation-driven areas, it does not seem that the details of our enthalpy formulation of the model are important, but are rather more so in the divide area.

Comparison of our results of basal temperature with those of Oerlemans and van der Veen (1984, fig. 5.8) show that our computed warm patches lie somewhat closer to the margin than do those of Oerlemans and van der Veen. It is simple to explain this observation in terms of the different basal sliding laws used. In the cases where the basal slip is not temperature-dependent, using the Morland and others sliding law (relation (8)) predicts basal power generation (relation (17)) to vary as H(∂H/∂X)², while Oerlemans and van der Veen, who use a Weertman functional form (relation (22)), have the basal power generation varying as (H∂H/∂X)^m, where m ≥ 2. The geometry of ice sheets dictates that the Oerlemans and van der Veen position of maximum basal power generation must lie further in than ours. Since the margin is kinematically very sensitive, this implies that not only the magnitudes of the coefficients inserted into sliding laws are important but also the functional form, since if a Weertman form implies cold margins we should expect steeper ice-sheet surface gradients and larger ice sheets.

A further inference that may be made from the basal temperature profiles is that regional refreezing is likely to be limited to a marginal zone. These show that a moving basal particle of ice is warmed nearly as far as the margin. It is only when heat is being removed from the ice that it can use the latent heat available in water. This is not a statement that regelation will not occur but that zones of refreezing should be limited. This statement must be qualified in a number of ways. First, in our model we have taken no account of the fact that the rheology of ice is implicitly pressure-dependent, as the correct temperature to chose (Reference RigsbyRigsby, 1958) in the constitutive relation is the homologous temperature. Wet ice being moved into a region of lower pressure at a constant absolute temperature is having its homologous temperature lowered, and will therefore become stiffer, and susceptible to refreezing. Owing to the dissipation-dominated thermal behaviour near the margin, we doubt that the ice could cool down, but it may be that depressurizing in areas of thick ice from ice being forced upwards by topography could occur faster than heating by conduction, and refreezing would occur here. This may be the case at Byrd Station. However, if the marginal zone of refreezing is larger than our model suggests, the implication is that the marginal model needs improvement. We doubt that adopting a Weertman sliding law is sufficient.

The difference in profiles obtained by the prescribed temperature fields of Reference Morland, Smith and BoultonMorland and Smith (1984) and of Hutter and others (1986) may be explained in terms of the self-regulatory structures of ice-sheet temperature profiles. Morland and Smith produced very flat profiles when they had cold margins and warm divides. As we have shown, this temperature structure does not appear to be admissible because of the nature of the processes which determine the structure of temperature fields within ice sheets.

These results have several implications for geological interpretations:

• (i) The sensitivity of the marginal area increases the importance of the evidence from this area. The tectonic deformation of soft sediments beneath the ice sheet and the structures thereby formed contain information which is amplified by the kinematics of the margin to have implications on the scale of the ice sheet itself.

• (ii) Zones of refreezing are predicted to be limited except when topographic effects occur. If refreezing is important, this contains information about the functional form of the ice-sheet profile. The areal extent of warm ice, interpretable by geological evidence, is dependent quite sensitively upon the lapse rate, which may therefore be constrained by geological evidence.

Finally, we have an alternative position on whether there is a shear-heating instability in ice sheets (see, for example, Reference Clarke, Nitsan and PatersonClarke and others, 1977; Reference FowlerFowler, 1980; Reference Yuen, Saari and SchubertYuen and others, 1986). Since the marginal areas reach the melting point, there is always an approach to high temperature fueled by dissipation, i.e. the phenomenology of a shear-heating instability. However, because these high temperatures are attained very early, they are part of the ice sheet. We have shown that ice sheets can, nonetheless, grow with very warm margins. The divide areas produce very low shear heat and we have shown that kinematical adjustment is at least as fast as any warming process. Thus, the thermal phenomenology of a shear-heating instability will usually (always?) exist, but this ever-presence implies no spectacular kinematical consequences.

8. Conclusions

• (i) Continental ice-sheet temperature fields evolve according to processes occurring at three time-scales: a dissipative time-scale of the order of thousands of years, an advective time-scale of the order of tens of thousands of years, and a conductive time-scale of the order of 100 000 years. Dissipation is the most important process in the marginal areas, advection in the upper regions of the divide area, and conduction in the basal divide area.

• (ii) Ice-sheet kinematics have two time-scales: an advective one of thousands of years and a diffusive one of tens of thousands of years. They are most dominant in the marginal and central areas, respectively.

• (iii) Dissipation drives marginal temperatures up to the melting point very rapidly; however, the form of the sliding law we use probably exaggerates this effect.

• (iv) Temperatures at the base of the divide are able to adopt two modes: a cooling one, associated with thickening of the ice sheet and a warming mode associated with low ice sheets. The latter is produced by low lapse rates and low thermal inertia, but we have not conducted a full search of the possible parameter variations. We do not believe that the warming mode leads to an instability, as the warming is due to conduction which operates on a very much slower time-scale than does mechanical adjustment.

• (v) The marginal area shows extreme kinematical sensitivity and, with the model we use, acts to amplify the effects of the marginal mechanics. This fact is also true of other one-dimensional models. Longitudinal gradients are not negligible at the margin, and this problem is made more acute by the existence of deformable materials underneath these areas of ice sheets. The thermal regime at the marginal area could be of extreme importance, and blanketing by pro-glacial lakes could even be of importance in determining marginal motion. The kinematical sensitivity of the margin increases the importance of the deforming-bed problem.

• (vi) The small-aspect ratio of the bed-temperature interference problem means that horizontal conduction of heat is not likely to occur except in marginal areas. The thermal inertia of the bed is important and may alter the geothermal heat gradient significantly. We do not believe that this will alter the modes of operation of ice sheets, because conduction is the slowest process operating and the ice sheet can therefore adjust to these changes.

• (vii) Continental ice sheets may have had temperate margins.

• (viii) The shear-heated approach to melting of the margins always seems to be present and seems to rule out the possibility of a shear-heated collapse of large ice sheets.

• (ix) The study has highlighted the inadequacies of numerical and mechanical models of the marginal areas. A more refined model of the mechanics and the coupled bed/ice-sheet thermal problem is needed and will dictate the numerical strategy needed.

Acknowledgements

This work was carried out as part of U.K. Natural Environment Research Council grant GR3/6737, “The temporal and spatial dynamics in relation to their geological products”. We had useful conversations with A. Fowler about approximate representations of the effects of phase transitions and with L. Morland about the reduced model.