Introduction

Glaciological modelling of Pleistocene ice sheets has hitherto been concerned to model equilibrium conditions at glacial maxima, either as an end in itself, or as an aid to comprehension of the geological products of glaciation (Reference RobinRobin, 1964; Reference Boulton, Boulton, Jones, Clayton, Kenning and ShottonBoulton and others, 1977; Reference SugdenSugden, 1977; Reference Denton and HughesDenton and Hughes, [^c1981]); to model the pattern of build-up of ice sheets (Reference Andrews and MahaffyAndrews and Mahaffy, 1976); or to model the response of ice sheets to externally imposed climatic changes (Reference WeertmanWeertman, 1976; Reference PollardPollard, 1978; Reference OerlemansOerlemans, 1980; Reference Budd and SmithBudd and Smith, [1981]).

Our intention is different. We argue that glacier systems respond to climatic change in an inherently more predictable fashion than most other major parts of the Earth’s surface system, and thus that the former extents of glaciers, identified by geological evidence, should be a valuable, and as yet largely unused indicator of palaeoclimate. In theory they should be more reliable indicators of regional palaeoclimate than the plant and animal species which are conventionally used. These palaeobiological techniques utilize empirical matching of modern species distributions with climate, though it is clear that complex series of variables may govern those distributions.

This article represents the application of a recently developed isothermal viscous flow theory (Reference Morland and JohnsonMorland and Johnson, 1980, Reference Morland and Johnson1982; Reference Morland, Morland, Smith and BoultonMorland and others, 1984) to the inference of mass-balance conditions from geological evidence of steady-state stages of former ice sheets. The model applied here is one in which the fundamental physical balances are satisfied in a properly formulated flow problem which is the basis for further development of a temperature-dependent, time-dependent model which will yield direct climatic inferences rather than indirect climatic indicators such as mass balance.

Theoretical Model

The equilibrium profile of an ice sheet for prescribed distributions of surface accumulation and ablation and of basal drainage, with a given basal sliding law, has so far been determined only under the restrictions of steady plane or axi-symmetric flow and uniform temperature (Reference Morland and JohnsonMorland and Johnson, 1980, Reference Morland and Johnson1982; Reference JohnsonJohnson, 1981). We adopt the steady plane isothermal flow solution to construct ice-sheet profiles for various accumulation and ablation patterns, recognizing that significant temperature variation through the sheet and the strong temperature dependence of ice response to stress must influence the profile detail. Our main qualitative conclusions are based only on the predicted span of the ice sheet, but nevertheless await confirmation by solution of the full thermomechanically coupled flow.

In rectangular axes Oxyz with Oxz horizontal, origin O at the ice sheet margin and flow parallel to the plane Oxy, the ice-sheet cross-section has surface y = h(x) and bed y = f(x), illustrated in Figure 1a. The bed is defined as a smooth contour representing a basal boundary for the global ice-sheet flow and ignores roughness on a scale small compared to the depth. The Reference Morland and JohnsonMorland and Johnson (1982) solution requires that the bed slope y = f′ (x) is everywhere small, so that the effects of any finite slope topography cannot be taken into account. The surface y = h(x) is stress-free relative to an assumed uniform atmospheric pressure, and subject to surface accumulation q per unit horizontal cross-section per unit time. Negative q defines ablation. We suppose that q = q(h), depending only on altitude, and that basal drainage is negligible.

Fig. 1. Ice sheet coordinate systems.

• (a) The Cartesian coordinate system Oxyz has its origin at the ice sheet margin. System Ox_Dyz has its origin at the extreme edge of the bed form f(x) where x = x_D − x_m. t_n and t_s are the normal and tangential tractions on the surface y = h(x), τ_b, u_b, and v_b are the tangential traction, tangential velocity, and normal velocity on the bed y = f(x).

• (b) Description of isostatic depression.

The other basal condition is the relation between the tangential traction τ_b and tangential velocity u_b, which define conditions on the smooth lower boundary y = f(x) of the global flow. These represent the mean effects of more complex ice and/or bed material movement on the roughness scale of the interface ignored in the global flow. We will use the conventional terms sliding velocity for u_b, and sliding law for the relation between τ_b and u_b, with the understanding that actual sliding of cold ice over the bed material is not necessarily implied. If real sliding is zero and if there is no thermally enhanced shear strain-rate, then the sliding relation must predict low mean velocity u_b. Since temperature variation is not incorporated in the present model, possible basal thermal conditions must be ignored in the adopted sliding law. There is no established law, but on the basis that a solution with a small surface slope can exist up to a margin, Reference Morland and JohnsonMorland and Johnson (1980, Reference Morland and Johnson1982) show that there must be additional dependence on overburden (basal) pressure P_b and that τ_b is linear in P_b as P_b → 0. We adopt a particular (restricted) form

(1)

where μ(0) is finite and non-zero, for which the function μ(p_b) has been determined empirically for two very different ice sheets (Reference Morland, Morland, Smith and BoultonMorland and others, 1984).

The ice is assumed to be an incompressible non-linearly viscous fluid with uniform density ρ_i and shear response

(2)

where is the rate of strain, D₀ = 1 a⁻¹, is a dimensionless deviatoric stress with unit σ₀ = 10⁵ N m⁻², and a(T) is the temperature-dependent rate factor. We adopt the polynomial representation for ω(J₂) determined by Reference Smith and MorlandSmith and Morland (1981) from Reference GlenGlen’s (1955) uniaxial stress data at T = 273.13 K, defining a bounded viscosity at zero stress, and the exponential representation of a(T) determined from Reference Mellor and TestaMellor and Testa’s (1969) data to evaluate a constant a(T) at the uniform temperature chosen to represent that of the dominant flow region. These are the functions used in the determination of the basal sliding relation (Reference Morland, Morland, Smith and BoultonMorland and others, 1984) where details are presented, together with a reformulation of the Reference Morland and JohnsonMorland and Johnson (1980 Reference Morland and Johnson1982) solution in fixed normalized variables suitable for any description of a large ice sheet.

Starting with a fixed aspect ratio or magnitude of surface slope ε₀ = 0.005, and fixed accumulation, thickness, and semi-span magnitudes

(3)

define dimensionless variables by

(4)

where (u, v) are the longitudinal and normal velocity components, and σ_xy, p are the shear stress and mean pressure. The units for ū, ⊽, σ_xy, p are respectively 200 m a⁻¹, 1 m a⁻¹, 10⁵ N m⁻², 200 × 10⁵ N m⁻². Then

(5)

where

(6)

and the basal sliding relation (1) is written

(7)

The profile equation becomes

(8)

subject to the margin conditions

(9)

where

(10)

and Ǫ₀ = Ǫ(margin), μ₀ =μ0. Eliminating ū_b through Equations (7) and (5) gives a second-order non-linear ordinary differential equation for H(X), which can be integrated from the margin using the two initial conditions (9).

Functions μ(pb) in the basal sliding relation (7) have been determined from the profiles and accumulation/ablation distribution of two very different ice masses, namely the western flank of the Greenland ice sheet at lat. 70°N., and the north-west flank of the Devon Island ice cap (Reference Morland, Morland, Smith and BoultonMorland and others, _1984). That is, H(X) and Q(H) are prescribed, so that μ _b is determined by Equation (8), and μ(pb), τ_b by Equation (5). The Greenland data cover a much wider range of overburden pressure μ(pb)μ(pb)μ(pb)μ(pb), and are adopted for our initial sequence of profile constructions; they yield the curve labelled 1 in Figure 2. We wish also to determine the effects of significantly changing the sliding function μ(pb), spanning and extending beyond the limits represented by the Greenland and Devon Island ice cap data, so have used the wide variety of functions μ(pb) shown in Figure 2. For analytic convenience they are defined by varying the parameters of a linear-polynomial-linear construction (with continuous second derivatives) analogous to Equations (36) and (37) of Reference Morland, Morland, Smith and BoultonMorland and others (1984).

Fig. 2. Different function μ(pb) used to test effects of changing basal conditions. Curve labels are for identification only.

Finally, the bed profile Y = F(X) is either a given function of X, or is linear in H to represent isostatic sinking. If the bed undeformed by the ice sheet is flat, F = 0, then isostatic depression is represented by

(11)

where ρ_b is the bed density. Alternatively, if the ice sheet has advanced over a bed with initial profile Y = F₀(X) and local depression is proportional to the overburden pressure (see Fig. 1b), then

(12)

Equations (11) and (12) are identical for F₀ = 0 if c = ρ_i/ρ_b, and can then be expressed in a common form

(13)

which we use to model depression effects of the ice sheet.

Mass-Balance Characteristics

We use the conventional terminology net mass balance for the surface accumulation density q (frequently denoted b_n) which is a prescribed function of altitude h in the adopted model. We thus require a general statement about the possible range of patterns of q against h. Figure 3a shows plots of q(h) for a large number of modern glaciers. These accumulation/ablation distributions are a function of latitude and of the location of the glacier with respect to moisture sources in the major storm tracks. High rates of ablation in lower latitudes (e.g. Reference Budd and SmithBudd and Smith, [1981]) require high accumulation rates to balance them, so that net balance gradients (∂b_n/∂h = ∂q/∂h) tend to increase as latitude decreases. In high latitudes high rates of winter accumulation along the axes of major storm tracks are balanced by high summer ablation produced by the heat imported to high latitudes by these same systems, so that gradients also tend to increase towards storm tracks. A third effect is that of ice-sheet size whereby less moisture is capable of penetrating into the interior of a large ice sheet. This effect is shown in Figure 3a by the progressive diminution of net accumulation with altitude. For the purposes of modelling we have schematized these patterns into a sequence shown in Figure 3b ranging from low gradients (“cold, dry, continental” conditions) to high gradients (“warm, wet, maritime” conditions).

Fig. 3. (a). Net mass-balance distributions for a number of modern glaciers. (b) Model net mass-balance distributions used in glacier reconstructions. Labels defined in Table I.

• 1. Blue Glacier, Washington (1964–67)

• 2. South Cascade Glacier, Washington (1957–70)

• 3. Sentinel Glacier, B.C. (1965–70)

• 4. Place Glacier, B.C. (1965–70)

• 5. Woolsey Glacier, B.C. (1965–70)

• 6. Peyto Glacier, Alberta (1965–70)

• 7. Ram River Glacier, Alberta (1965–70)

• 8. Nigardsbreen, Norway (1969–73)

• 9. Hellstugubreen, Norway (1964–73)

• 10. Gråsubreen, Norway (1964–73)

• 11. Ålfotbreen, Norway (1964–73)

• 12. Storbreen, Norway

• 13. Storsteinfjellbreen, Norway

• 14. Blåisen, Norway

• 15. Cainhavarre, Norway

• 16. West Greenland, lat. 70°N (modelled)

• 17. Folgefonni, Norway. (1964–65 & 67–68)

• 18. Antarctica. Wilkes Station Profile

• 19. Thule Ramp, Greenland

• 20. Antarctica, Casey Profile

These patterns are normalized according to Equations (4) for use in the differential Equations (8), but our final results are expressed in physical variables. In Figure 3b the normalized representations Q(H) for each pattern q(h) are shown. In each case we have assumed linear dependence of Q on H in the ablation zone, i.e. for H ⩽ H_e, where H_e is the normalized equilibrium-line altitude (ELA), and also in the higher altitude zone, i.e. for H > H_s. A cubic polynomial is used for H_e ⩽ H ⩽ H_s to obtain continuity of Q(H) and Q′(H) at H = H_e and H = H_s. If we let g_a be the ablation gradient, g_s the gradient of the linear accumulation segment, and Q_s the value of Q at H = H_s, the representations Q(H) are given by

(14)

where

(15)

For each of the accumulation/ablation patterns Q(H) in Figure 3b, g_a, g_s, and Q_s are fixed. Furthermore, although the value of the ELA H_e varies according to the position of the origin on the H-axis in Figure 3b, the difference H_s-H_e remains unchanged. Thus a₁, a₂, and a₃ are constant for a given accumulation pattern. The chosen value for the difference H_s-H_e is 0.25, equivalent to a difference of 500 m in physical altitude. This value is assumed to hold for all choices of accumulation patterns and is a round estimate made from a study of the data in Figure 3a. We have listed in Table 1 the values of the parameters Q_s, g_s, g_a = a₁, a₂, and a₃, assuming H_s-H_e = 0.25, which determine the six different labelled patterns Q(H) appearing in Figure 3b. For some patterns we have considered two values of the gradient g_s.

Equation (8) is now solved for H(X), both for a flat bed F = 0 and an isostatically depressed bed F = vH, with v = −1/3 (equivalent to a density ratio ρ_b/ρ_i = 4). Initially we only solve the half-problem, i.e. from the left-hand margin to the divide. The full profile from margin to margin is assumed to be symmetric about the ice divide, equivalent to choosing the same mass-balance pattern Q(H) on both sides of the divide and with the same equilibrium-line altitude H_e. We use the net balance gradients shown in Figure 3b and the ELA as independent variables and investigate their separate influence on glacier size and shape.

Influence of Equilibrium-Line Altitude (ELA) and Ablation Gradient on Glacier Size and Shape for a Horizontal Bed

We first use the basal boundary condition inferred for the West Greenland ice sheet at lat. 70°N (Reference Morland, Morland, Smith and BoultonMorland and others, 1984), and solve Equations (8) for a given net accumulation/ablation gradient and ELA for a horizontal bed. Figure 4a shows the glacier spans and central heights which would be in equilibrium with the patterns defined by Equations (14) and (15), and Table I. The diagram indicates the range of accumulation/ablation conditions which could be inferred from the given span of a glacier in a steady state resting on a horizontal plane. It does not imply that a given climatic condition will necessarily produce a glacier of a certain size. The conditions which might lead to the evolution of an ice sheet of given size require a time-dependent analysis.

Fig. 4. (a) Glacier span versus ELA for different net balance distributions for a horizontal bed in the absence of isostasy. The bold lines connect points on the curves of similar central height. Note that glaciers with higher net balance gradients have steeper profiles.

(b) As 4(a) but with full isostatic compensation.

Reference WeertmanWeertman (1961) has demonstrated how the ELA will be related to the snow-line during glacier expansion or decay. (We assume that there is a horizontal or gently dipping plane within the atmosphere above which perennial snow will remain on mountains or glaciers penetrating through it. The line of intersection of this plane with a physical feature or any surface is the snow-line.) In the Northern Hemisphere the snow-line dips towards the north. Consider the situation in Figure 5a where a snow-line S₁S′₁ cuts the curve of ELA versus span for a glacier of span A at point X. If the overall elevation of the snow-line is reduced to S₂S′₂, the accumulation area is increased and the ablation area reduced, a positive balance will ensue, and the glacier will advance to B for which the curve of ELA versus span cuts the new snow-line at Y. An increase in snow-line height will produce the reverse trend.

Fig. 5. The effect of snow-line elevation and slope on glacier span.

(a) The effect of changing snow-line elevation on a glacier hitherto in equilibrium, when the snow-line slope is steeper than that of the curve of ELA against span (pattern 4(1/2)).

(b) The effect produced when snow-line slope is less steep than the curve of ELA against span.

(c) Modern and reconstructed glacial stage snow-lines from Europe and North America. Note that apart from the Labrador example, they are less steep than the curve of ELA against span of 4a and 4b.

TABLE I. Parameters for Normalized Accumulation/Ablation Patterns Q(H) Defined by Equations (14) and (15)

In Figure 6a–d we show the computed glacier profiles which would be in equilibrium with different equilibrium-line altitudes and net mass-balance patterns taken from Figure 3b. We find that glaciers with high net balance gradients, typical of maritime climates, have relatively high surface slopes reflecting the higher shear stresses required to generate flow velocities large enough to balance high accumulation rates. An important attribute of an advancing or retreating glacier is the ELA trajectory, the successive positions of the ELA on the glacier surface. For a given mass-balance pattern the ELA trajectory for successive steady-state profiles is a critical trajectory. In an advancing glacier, the ELA trajectory will lie below the critical trajectory and cannot intersect it, and vice versa for a retreating glacier. This critical ELA trajectory is thus a fundamental parameter relating climate and glacier behaviour in a given area.

Fig. 6. (a)-(d) The sequence of steady-state profiles for different net balance distributions associated with an increasing snow-line elevation. For a given mass-balance pattern the ELA for a given span must lie beneath the ELA trajectory for the glacier to advance, or above it to retreat. If it coincides with the ELA trajectory, the glacier front will be stationary.

(e) Comparison of ELA trajectories appropriate to different net balance distributions.

As Reference WeertmanWeertman (1961) has shown, when the slope of the snow-line is greater than the slope of the curve of ELA against span (Fig. 5a), a fall in snow-line elevation S₁S′₁ to S₂S′₂ will cause the glacier to advance to a new equilibrium position, whereas if the snow-line slope is lower than the slope of the curve of ELA against span, a fall in snow-line will cause unstable advance of the glacier (Fig. 5b) which would advance until a physical barrier was reached or some other change stopped its expansion.

Most reconstructed modern snow-lines from the areas of Europe and America occupied by large ice sheets durig the last glacial period have a very low northerly slope (Fig. 5c). Reconstructed snow-lines from different phases of the last glacial period (Fig. 5c) have a similar slope to modern snow-lines in the same area with the altitude reduced between 700 and 900 m. Most of these snow-lines have slopes lower than the typical slope of the curves of ELA against span (Fig. 4). Thus it might be concluded that contemporary ice sheets were inherently unstable in Reference WeertmanWeertman’s (1961) sense. This does not necessarily follow. An assumption generally made in the discussion of ice-sheet growth and decay is that the net balance gradients will remain constant through time, and that only the elevation of the ELA will change (e.g. Reference Andrews and MahaffyAndrews and Mahaffy, 1976; Reference WeertmanWeertman, 1976; Reference OerlemansOerlemans, 1980). We suggest that this is unlikely to be true. A climatic change which results in reduction of snow-line elevation may also change the associated net balance gradient. We would expect, for instance, that during ice-sheet growth in northern Europe and North America, decreased evaporation from a cooling ocean surface will result in a drier, colder atmosphere. Thus we expect an initially large net balance gradient to be reduced. Also, extension of the ice sheets to great size as they expand towards their maxima will make it increasingly difficult for moisture-bearing storms to penetrate very far into their interiors, so that accumulation rates will be diminished. While we might also expect the deflection of well-defined storm tracks around the ice sheets’ southern margins to result in high winter accumulation rates there (e.g. Reference LiljequistLiljequist, 1974; Reference Williams, Barry, Weller and BowlingWilliams and Barry, [^c1975]), this would be offset by the penetration of the ice sheets to warmer latitudes by this stage producing high summer ablation rates, so that the “elevation-desert” effect in the ice sheet interior will be balanced by low net ablation rates near the margin. Thus, rather than using the conventional model in which the ablation gradient is constant during ice-sheet growth and only the ELA changes (Figs 5 and 6), we suggest a model in which not only will the ELA change, but the net balance gradients will be reduced as the ice sheet expands. Thus we analyse the response of a glacier to changes in a more complex climatic framework.

We first take the case of snow-lines steeper than the curve of ELA against span. Figure 7a–c shows a sequence of snow-lines S₁ to S₅ and contemporary curves of ELA against span E₁ to E₅, so that snow-line S₁ is contemporary with curve E₁ and so on, with lower curves of ELA against span reflecting colder, drier, more continental conditions (lower net balance gradient). Consider a glacier in equilibrium with a given accumulation/ablation distribution and a snow-line S₁S′₁ which intersects the appropriate curve of ELA against span E₁ at 1. If the snow-line is reduced in elevation from S₁S′₁ to S₂S′₂, and at the same time the net balance gradient is reduced so that at location 1, S₁ − S₂ > E₁ − E₂, the intersection between the new snow-line S₂S′₂ and the new curve of ELA against span E₂ indicates that an advance to 2 will occur. If, however, the change in snow-line for the same net balance gradient change is less, so that S₁ − S₂ = E₁ − E₂ at 1 (Fig. 7b), then the glacier will remain stationary. If S₁ − S₂ < E1 − E2 at 1, the glacier will retreat (Fig. 7c). The converse applies to snow-line rise. Thus, for a glacier in a steady state, snow-line lowering does not necessarily lead to advance, nor snow-line rise to retreat; the reverse effect may be produced. Lowering of snowline coupled with the increase in elevation of the curve of ELA against span due to an increase in the net balance gradient will always lead to glacier advance.

Fig. 7. The combined effect of snow-line elevation and mass-balance gradient which control whether a glacier will remain stationary, expand or contract. S₁ to S₅ represent a sequence of snow-lines and E₁ to E₅ a sequence of contemporary curves of ELA against span reflecting a progressively lower net balance gradient.

(a)-(c) Snow-line gradient steeper than ELA span curve. If at point 1: S₁ − S₂ > E₂ − E₂, glacier advances to 2; S₁ − S₂ = E₁ − E₂, glacier standstill; S₁ − S₂ < E₁ − E₂, glacier retreats to 2.

(d) Snow-line gradient lower than the curve of ELA against span. If the curve of ELA against span is reduced in elevation from E₁ to E₂, and snow-line altitude S₁ is reduced, reduction to S₂ produces unstable retreat, to S₃ a standstill, and to S₄ an advance.

(e) If the climate becomes drier leading to a reduction in net balance gradient, an increasingly large depression of snow-line is required to maintain a stationary ice front.

In the case of snow-lines which are less steep than the curve of ELA against span, which we believe may have been the usual case for the two great northern-hemisphere ice sheets of Europe and North America, slightly different considerations apply. Changing snow-lines may lead to unstable behaviour as Weertman suggested, or to the establishment of new steady states. If snow-line and net-balance gradients are reduced so that S₁ − S₂ < E₁ − E₂ (Fig. 7d) the reduction in snow-line is less than that required to maintain a stationary ice front position and the glacier will undergo unstable retreat. If S₁ − S₄ > E₁ − E₄, the reduction in snow-line is more than required to maintain a stationary ice front and the glacier will undergo unstable advance. Figure 7e shows the change in snow-line elevation required to maintain a steady-state position for two different changes of balance gradient. Large glaciers require much greater changes in snow-line elevation to sustain a stationary ice front, and as climate becomes drier, a greater snow-line depression is required.

Figure 8 shows a glacier advancing via a sequence of steady-state stages from an initial condition of growth in a relatively maritime climate with a high net balance gradient to its maximum in a “continental” climate with a low net balance gradient. Each stage represents an intersection between a snow-line and the appropriate curve of ELA against span in Figure 4a. The ELA trajectory (Fig. 4a) shows the maximum elevations of the sequence of snow-lines required to produce this advance, or the sequence of minimum elevations required to produce retreat. In general, as is shown in Figure 7e, a greater snow-line depression is required to sustain advance as the climate becomes progressively drier.

Fig. 8. (a) Ice-sheet growth on a horizontal surface in which both ELA and net balance gradient change through condition 1 to 5a. Rapid climatic changes produce a rapid change in ELA trajectory for a fixed span.

(b) Ice-sheet growth with full instantaneous isostatic compensation. Note the dramatic reduction in central height which can accompany a frontal advance associated with a change from mass balance pattern 4(1/2) to 5a, a reduction in net balance gradient associated with colder, drier conditions.

In order to produce a halt stage or short retreat during an advance, climate could either become drier (and presumably colder), so that the ELA trajectory coincides with or is less than the snow-line elevation, or it could become warmer (and presumably wetter), with an increased snow-line elevation which coincides with or exceeds the ELA trajectory.

The coupling of mass balance, flow, shear stress, and glacier profile in a viscous glacier model produces a reduction in the centre height:span ratio as the balance gradient is reduced during glacier expansion into an increasingly continental climate. In Figure 8a we contrast glacier expansion during phases of gradually changing climate with expansion during a phase of rapidly changing climate. In this latter case, the ice front advances substantially whilst the centre height is substantially reduced. The effect of this is to produce a relatively small volumetric change in the ice sheet even though there may be a large change in span. Thus we should not expect the ice-sheet volume changes implied by the deep ocean isotopic record to be precisely reflected by changes in ice-sheet margins. Indeed, there is evidence that at least in North America, the ice-sheet span maximum was achieved long after the global ice-sheet volume maximum (Reference Ruddiman and McIntyreRuddiman and McIntyre, 1981).

Role of Isostasy

Figure 6c–d shows the profile contrasts between glaciers resting on horizontal beds and those for which there is full isostatic compensation. A higher ELA is required for a glacier resting on a horizontal bed compared with a glacier of similar span resting on a bed which has suffered significant isostatic depression (Fig. 4). The ELA trajectories required for growth in the latter case are lower than in the former.

Consider a glacier which advances rapidly to a given steady-state position (before significant crustal depression can occur) under conditions where the slope of the snow-line is greater than the slope of the curve of ELA against span. As the crust begins to sink beneath the glacier load, a progressively lower ELA will be required for equilibrium and the glacier will retreat until a new equilibrium is established. The interaction between the response rate of the glacier and the response rate of the crust could produce complex behaviour, but eventually an equilibrium may be established, although in practice it is doubtful whether climate will ever remain stable long enough to allow this to occur. Where the snow-line slope is less than the slope of the curve of ELA against span, crustal sinking beneath a glacier in equilibrium will lead to unstable retreat of the glacier which will not attain a subsequent equilibrium unless the climate changes.

A glacier advancing via a series of equilibrium stages and always in isostatic equilibrium with the underlying crust has an ELA trajectory shown in Figure 8b. If we assume that full isostatic compensation is never attained, then the maximum ELA trajectory required for an advance must lie between the trajectory for zero compensation (Fig. 8a) and that for full compensation (Fig. 8b).

Not only will the effect of isostatic depression be to act as a brake on glacier advance, but also as a brake on glacier retreat. A glacier retreating because the snow-line lies above the curve of ELA against span appropriate to substantial isostatic depression may reach an equilibrium halt stage or re-advance if crustal uplift causes relative lowering of the snow-line so that it cuts the curve of ELA against span appropriate to a lesser degree of isostatic depression (Fig. 6c, d).

Ice Sheets with Asymmetric Net Balance Patterns

Modern ice sheets frequently show strong asymmetry across the ice divide with a steep flank of small span opposite a low-slope flank of greater span. We investigate the possibility that different mass-balance characteristics across the divide may produce such asymmetric ice sheets. Consider an ice sheet with a similar net mass-balance pattern on both its flanks but with a different ELA on each side. We recall that for a given basal boundary condition and mass-balance distribution, the differential Equations (8) determines a unique semi-profile or symmetric full-profile, with the central height and semi-span also part of the solution. If we reconstruct both flanks independently, starting from the respective margins, we find that the central heights do not match (independent solutions in Fig. 9a).

Fig. 9. Effect of asymmetric mass balance, where net balance gradients are constant (pattern 5a) but ELAs differ.

(a). Where mass-balance patterns are different on opposite flanks of an ice sheet, independent solutions for each flank do not match at the ice divide. A “coupled” solution effectively increases the span and centre height of the B-side (lower ELA), whereas the A-side profile is identical for both independent and coupled solutions.

(b). The solid line is the curve of ELA against span for mass-balance pattern 5a for an independent A-side. If the ELA on the B-side is lowered, the dashed lines show curves of ELA against span for the B-side. The dotted lines connect points of similar height and thus give the spans to the north and south of the ice divide.

In solving the asymmetric problem we solve Equations (8) to obtain the profile h(x) from margin to divide as before, thus determining the central height h_c. The remainder of the profile is then determined, with the new ELA, by starting from the divide with h = h_c, h′ = 0. This latter problem is over-specified, however, in that the initial value h = h_c does not in general coincide with the central height obtained by constructing the profile from the far margin. The outcome of this over-specification is that the slope h′ (x) becomes infinite as the other margin is approached, and thus the ordinary differential equation, valid only for small slopes, no longer applies. Although the physically realistic finite-slope solution near this second margin must then be determined by consideration of a (passive) vertical boundary layer, the singular solution (dashed curve in Fig. 9a) provides us with a good estimate of the span. In Figure 9a, we see that while the uncoupled profile constructed from the margin has a span of 314 km and a central height of 1673 m for the A-side and a span of 219 km and a central height of 1208 m for the B-side, the coupled solution from the left margin, retaining the A-side, has an increased span of 286 km on the B-side. The dotted lines in Figure 9b connect selected A-side spans with the appropriate B-side spans (similar central height) where the B-side ELAs are 250 m and 500 m lower than the A-side ELAs.

We have also examined the effect of changing net balance gradients whilst keeping ELAs constant on opposing flanks. The effects are very similar to those shown in Figure 9b, but not so evident. We conclude that major asymmetry in an ice sheet is unlikely to result directly from contrasting mass-balance conditions on opposite flanks, but more from the nature of the basal boundary on each side of the divide.

Effect of the Basal Boundary Condition

We now present a selection of results obtained by changing the basal sliding function μ(pb) over the wide range shown in Figure 2. Figure 10a shows the profiles over a horizontal bed with no depression for accumulation pattern 4(½) in Figure 3b with ELA at 1500 m. The labels correspond to those of the μ(pb) curve in Figure 2. If μ is low, basal movement is easy, and only a relatively low shear stress and associated surface slope are needed to discharge a given accumulation, which is reflected by the comparison of profiles 5^*, 10, 4, 9, 8 with profiles 2 and 3. However, because on most glaciers we find that the net balance does not change dramatically beyond a small height above the ELA (Fig. 3), changes in basal boundary condition in the accumulation area have little effect on glacier span, although the profile shape may be quite different. Thus, basal boundary condition 5^* (Fig. 2), with much easier basal movement in the accumulation area, has a similar span to 2 (Fig. 10a), although 5^* produces a much flatter profile. This is because similar μ(pb) values in the ablation area produce similar snout profiles, and similar ablation-area extents, which, for a steady state, must be balanced by similar accumulation-area extents. Glacier span is, however, sensitively dependent upon the basal boundary condition in the ablation area, which affects surface profile, ablation-area extent, and thus accumulation-area extent. Thus, boundary condition 11 (Fig. 2), with low μ(pb) values in the snout region, produces a very much larger span than the other conditions for a similar mass-balance distribution (Fig. 10a). The evidence from most modern glaciers, which tend to have similar snout profiles, is that μ(pb) values do not vary dramatically in the ablation area, and thus, for a given net balance gradient a predictable relationship between ELA and span exists. Thus, although the basal boundary condition used in this application is derived from the Greenland ice sheet, we conclude that span is not sensitively dependent on the precise values used in the basal sliding law, and therefore our reconstructions do not have to assume a similarity between the modern Greenland ice sheet and former ice sheets.

Fig. 10. (a) Different glacier profiles for different basal boundary conditions and a fixed ELA of 1500 m and fixed net balance gradient. The numbers on the profiles show the different basal boundary conditions used (see Fig. 2). Although central heights vary considerably, in most cases spans are very similar, except those (8,11) where the μ(pb) values are low in the terminal area (Fig. 9).

(b) Curves of ELA against span for different basal boundary conditions.

With a reasonable range of accumulation/ablation conditions and relatively high (“normal”) μ(pb) values in the snout region, the present solution fails to produce a large ice sheet with a span of the order of 1000 km or more. The large size of the Antarctic ice sheet reflects high rates of marginal calving which require a very broad accumulation area to balance it. The earlier Quaternary ice sheets in Europe and North America, 1500–2000 km in span, had no such advantage over much of their margin, yet we believe that they were near to equilibrium at their maximum extents (see later discussion). This may arise through low terminal μ(pb) values, which could have been achieved by deformation of subglacial sediments (Reference Boulton and JonesBoulton and Jones, 1979).

Effect of an Irregular Bed

Equation(8) for the dimensionless profile H(X) can be solved for an irregular bed where the dimensionless bed form F is a prescribed function of X. As an example we take part of the bed of the ice sheet which flowed over northern Britain during the last glacial period (Fig. 11b). We match this with an 11-degree Chebyshev polynomial, and use an alternative co-ordinate system Ox_Dyz (Fig. 11b). The dimensionless bed form F(X) is based on the representation f₀ in Figure 11b. In the absence of isostasy, the relation is simply

(16)

where x_m is the distance from the origin of the Ox_Dyz system to the origin of the Oxyz (OXYZ) system (Fig. 11b). Note that when the displacement x_m is zero, then the coordinate systems coalesce and X_D = x = ℓ₀X. In the isostatic case the dimension less bed form used in Equations (8) is

(17)

which, with ρ_b/ρ_i = 4, becomes

(18)

Fig. 11. (a) Devensian (last glacial period) flow patterns over northern England and the location of the York and Escrick moraines.

(b) The bed profile for the chosen northern England flow line (dashed) and the modelled bed form f₀ fitted by a Chebyshev polynomial (continuous line). The double coordinate system is also shown (see Fig. 1a).

The curves of ELA against span and comparative profiles in Figures 12 and 13 show that large-scale irregularities reduce the span and increase the profile steepness for a given net balance gradient, reflecting the additional resistance to glacier flow provided by major topographic features.

Fig. 12. Curves of ELA against span for glaciers with zero and full isostatic compensation for the selected northern England flow line in Figure 11b and for a horizontal bed f₀ = 0

Fig. 13. (a) A series of equilibrium profiles and the ELA trajectory for an ice sheet of constant net balance gradient (pattern 4(1/2)) advancing from south-west Scotland to the Vale of York. (Zero isostatic compensation).

(b) As for (a) but with full isostatic compensation.

(c) A comparison of different possible profiles of an equilibrium ice sheet terminating in the Vale of York, for several patterns of net balance distribution and their appropriate ELAs.

(d) Different spans produced for different net balance patterns in northern England for an ELA of 1500 m.

Reconstructions of Quaternary Ice Sheets

Several attempts have been made to use glaciological models to reconstruct Quaternary ice sheets. Reference Andrews and MahaffyAndrews and Mahaffy (1976) made a series of assumptions about the evolution of climate in North America during the last glacial period to estimate likely patterns of growth of the Laurentide ice sheet. Reference Boulton, Boulton, Jones, Clayton, Kenning and ShottonBoulton and others (1977) also assumed a climatic regime to reconstruct the thermal and dynamic characteristics of the Late Devensian ice sheet over Britain. Reference SugdenSugden (1977) modelled the Laurentide ice sheet by assuming a longitudinal profile which was a mean of profiles from the Antarctic and Greenland ice sheets and with accumulation patterns derived from the same regions. Reference Denton and HughesDenton and Hughes ([^c1981]) modelled all the major Quaternary ice sheets at their presumed maximum extent by making assumptions about basal thermal conditions and the appropriate basal flow laws. Reference WeertmanWeertman (1976), Reference OerlemansOerlemans (1980), and Reference PollardPollard (1978) have attempted to model the growth and decay of ice sheets in order to test if their major cycles of growth and decay could primarily be produced by Milankovitch solar radiation cycles.

Our aim is to use geological evidence of glacier form and extent to reconstruct climate. We assert that for a glacier of known surface profile flowing over a known bed there is a unique combination of ELA and net balance gradient. However, we do not know of geological evidence which might enable us to reconstruct the surface profile of a former glacier to which our two-dimensional flow analysis would be appropriate. Such evidence is only available for valley glaciers, which must be analysed three-dimensionally. However, we have abundant evidence of spans of former glaciers for which a two-dimensional approximation would be appropriate. Here we claim that for such a glacier in equilibrium, a predictable set of combinations of ELA and balance gradient are appropriate. We thus use the known extents of inferred steady-state glaciers as indicators of palaeoclimate.

The geological evidence for the former presence of a steady-state glacier margin has been reviewed by Reference Boulton and RobinBoulton (1983). It is:

• (a) major accumulation of melt-out till and/or flow till in a well-defined zone;

• (b) series of closely spaced push moraines;

• (c) laterally extensive evidence (hundreds of kilometres) in the form of moraines and sediments of an ice-marginal position.

Using such evidence, a sequence of equilibrium stages were identified during the retreat of the European ice sheet, to which a glacio-climatic model could be applied.

As an example we choose the ice sheet which covered much of northern England during the maximum of the last (Devensian) glacial period. In north-eastern England a major laterally continuous moraine, the Escrick moraine, is succeeded, 10 km to the north at York, by another similar moraine. Internally these moraines show evidence of glacier pushing, and they contain flow tills interbedded with glacial outwash (Reference GauntGaunt, 1970). They clearly represent substantial halt stages in the retreat of the late Devensian ice sheet from this area, which is known to have occurred after about 18000 B.P. (Reference Penny, Penny, Coope and CattPenny and others, 1969). From this we assume that a steady-state equilibrium was attained here. Erratic transport paths, and drumlin and striae orientations all indicate that the ice sheet flowed into this area from the north-west (Reference CharlesworthCharlesworth, 1957; see also Fig. 11a) and that it had its principal source in south-west Scotland where Reference CornishCornish (in press) has recently argued that a major relatively static ice centre was located during the Late Devensian glacial maximum. Further, reconstructed flow lines suggest parallel flow in this sector of the ice sheet, and thus two-dimensional flow can be assumed. The bed over which the ice sheet flowed was highly irregular (Fig. 11).

We determine the conditions under which a steady-state ice sheet centred over south-western Scotland could have extended into the Vale of York, and also the conditions which would have enabled equilibrium to be attained for ice sheets terminating north of the maximum extent. This was done both for zero isostatic depression and full isostatic compensation. It is interesting to compare these results with those obtained for a horizontal bed and a horizontal bed plus isostasy (Fig. 12). For a given snow-line elevation and net balance gradient, the extent of the ice sheet on the irregular bed is less than it would be on a horizontal bed. The lowering of snow-line required to sustain a given span on isostatic depression of the crust is substantially reduced for an irregular bed.

Figure 13a–b shows a series of equilibrium stages of a glacier expanding along the northern England flow line with and without isostatic depression of the crust. Their ELA trajectories for a single mass balance pattern are compared. Figure 13c shows the different profiles which would be produced by the combinations of ELA and net balance gradients required to maintain a steady-state ice sheet extending to the Escrick moraine. Figure 13d shows the extents of steady-state ice sheets with different accumulation/ablation patterns but with a constant ELA of 1500 m.

For the period of the last glacial maximum or about 18000 B.P., Reference BüdelBüdel (1949) and Reference Kaiser and WrightKaiser (1969) have reconstructed European snow-lines, suggesting a snow-line depression of about 1000 m compared with the present day. This would give a snow-line elevation between about 600 and 800 m a.s.1. on an ice sheet in steady state terminating in the vicinity of York. An ELA of 600–800 m a.s.l. would require (from Fig. 12) a low ablation gradient to sustain this steady state. Modern glaciers with a similarly low ablation gradient are located in areas of relatively cold, dry, continental climate. However, conditions similar to those at the southern margins of the European and North American ice sheets may not occur at the present day. A low ablation gradient could also be achieved in a region of high summer temperature and cold, snowy winters with a rapid fall off of precipitation over the ice sheet.

We suggest that climatic conditions during the Devensian late glacial re-advance between 11000 and 10000 B.P. are likely to have been similar to conditions during the build-up stage of the Late Devensian ice sheet over Britain. Studies of corrie moraines in south-west Scotland, formed during the late glacial re-advance, suggest a contemporary snow-line of 500–600 m. Figure 12 suggests that a small glacier centred over the mountains of south-west Scotland with an ELA of this value would require a relatively high ablation gradient, typical today of a maritime climate.

Thus we suggest that the growth of an ice sheet from centres in the highlands of northern Britain took place in a relatively maritime climate, but that as it expanded the climate became increasingly continental. We model the initial stage by a steady-state glacier with its margin some 50 km from the ice divide, an ELA of 500–600 m and an appropriate net balance gradient (Fig. 12), and assume that intermediate stages of spans 100 km, 150 km and 175 km have net balance gradients intermediate between the initial value and the final value for a glacier of span 200 km and ELA 1010 m (Fig. 14). This value has been used, as the geologically inferred value of 600–800 m a.s.l. would require an even lower net balance gradient than pattern 5a, and we have not computed such a pattern. We have assumed full isostatic compensation throughout. The glacier surface slope becomes lower as it expands, since the force required to drive decreasing accumulation through the system becomes less. The basal velocity distribution is strongly affected by the topography.

Fig. 14. Series of equilibrium stages for a glacier advancing along the northern England flow line when climate changes from relatively “maritime” to strongly “continental”.

(a) Zero isostatic compensation.

(b) Full isostatic compensation.

Conclusions

• (1) For an isothermal ice sheet in a steady state the surface profile, including its span and central height, is a function of the accumulation/ablation pattern, the basal boundary condition, and the form of the glacier bed.

• (2) For a given glacier bed form and accumulation/ablation function a change in the basal boundary condition beneath the accumulation area will not affect the ice-sheet span, although the central height may be affected. Thus, provided that there is not a considerable between-glacier variation in the basal boundary condition in the ice-sheet snout region, we conclude:

  • (a) An ice sheet of given span will be in equilibrium with a predictable series of accumulation/ablation patterns.

    If in addition the ELA is known from geological evidence, we conclude:

  • (b) An ice sheet of given span will be in equilibrium with a unique accumulation/ablation pattern.

    If we also know the ice sheet terminal profile from geological evidence, we conclude further that:

  • (c) An ice sheet of given span will be in equilibrium with a unique ELA and net balance gradient.

• (3) The state of stability or instability of an ice sheet, and whether it will advance, retreat, or remain stationary as a result of a climatic change, depends not only on whether snow-line is raised or lowered, but also on the associated changes of net balance gradient. Lowering of snow-line does not necessarily lead to glacier expansion, nor raising to contraction. A small lowering of snow-line coupled with reduction of the net balance gradient may lead to glacier contraction. Although the climate has become cooler, the glacier has been starved of accumulation.

  The concept of an ELA trajectory is introduced, which is the path which would be followed by the equilibrium-line altitude of a sequence of steady-state glaciers during advance or retreat. It marks the maximum possible elevation of the snow-line during glacier advance, or its minimum during retreat. The ELA trajectory will depend on the sequence of climates which accompanies glacier expansion or decay.

• (4) We expect the net balance gradient on an ice sheet to diminish as the ice sheet expands (increasingly “continental” climate). As the volume to span ratio is smaller for a glacier of given size when its net balance gradient is less, we expect the rate of volume increase as span increases to become less. A sudden change of climate from maritime to continental may thus lead to advance of the ice-sheet front with little volumetric change. There is evidence that the volume maximum of some Pleistocene ice sheets preceded their greatest areal spread. This may reflect an important climatic change.

• (5) This latter effect results from the assumption, supported by observations on ice sheets and ice shelves, that in the long term ice behaves as an incompressible non-linearly viscous fluid. If the net balance gradient increases without an increase in the glacier’s net balance, a greater shear stress is required to drive a greater discharge through the system, and thus the longitudinal surface profile becomes steeper.

  Much previous modelling is derived from or based on (Reference WeertmanWeertman's 1961 Reference Weertman1976) analyses based on an assumption of perfect plasticity. In such models the speed of ice flow adjusts itself to changes in net balance gradient without a change in the shear stress. It effectively de-couples a highly coupled system which is at the heart of the interaction between atmospheric and ice-sheet dynamics.

• (6) The snow-lines required for glacier growth are lower the stronger the isostatic component of crustal sinking.

• (7) Application of the model to the flow of the ice sheet which covered part of northern Britain during the last glacial period suggests that its initial growth was favoured by a relatively high net balance gradient (a relatively “maritime” climate) whilst at its maximum extent the net balance gradient was much smaller suggesting a relatively “continental” climate.

• (8) Without iceberg calving at the margin, such as occurs around the margin of the modern Antarctic ice sheet, and using a basal boundary condition appropriate to modern ice sheets, we are unable to model those extensive parts of the large last glaciation ice sheets in the northern hemisphere which terminated on land, and which geological evidence suggests were in a steady state. We conclude that basal movement must have been much easier for these ice sheets.