1. Introduction

In polar ice sheets, ice flow and its thermodynamics are strongly related: temperature determines to a large extent the viscosity of ice, so that for a 10 Κ temperature change, strain-rates, for a given stress, change by an order of magnitude, thereby constituting a major factor in controlling the flow characteristics and shape of large polar ice sheets.

Moreover, since the interaction of temperature and flow comprises a positive feed-back mechanism with regard to (sudden) global climatic changes, a potential instability exists: the increasing ice temperature – increasing strain-rate – increasing dissipation feed-back loop may, under suitable circumstances, lead to a runaway situation in which a large amount of continental ice is discharged into the ocean in a relatively short time. Reference WilsonWilson (1964) has put forward the hypothesis that surges of the Antarctic ice sheet, due to this mechanism of creep instability, could be a cause of global glaciation by creating extensive ice cover over the oceans of the Southern Hemisphere. The resulting increase in reflectivity for solar radiation would then lead to the temperature drop necessary to initiate the formation of ice sheets on the continents of the Northern Hemisphere. This idea has been discussed by, amongst others, Reference HollinHollin (1980), who reviews evidence pointing to a major East Antarctic ice surge 95 000 years B.P. The critical point here is the interpretation of higher sea-level stands and, in the light of recent calculations with geophysical models, some inferences remain speculative (e.g. Mörner 1980).

In principle it is possible that strain heating crosses a certain threshold where runaway temperature increase ultimately leads to extensive basal melting, bringing a continental ice sheet into the surging mode. At that stage, sliding instability may come into play, i.e. a feed-back loop where frictional heating due to sliding increases the production of lubricating melt water. Although potentially important, this mechanism acts entirely independently of creep instability, which is a consequence of the deformation properties of ice. Here we concentrate on creep instability.

Reference Clarke, Nitsan and PatersonClarke and others (1977) give a thorough treatment of the temperature regime in a uniform slab of ice. They solve the non-linear heat equation for prescribed downward advection and applied driving stress, and show that under specific conditions multiple steady states occur (in fact, the elementary cusp catastrophe). They also carried out an analysis of linear stability, revealing that typical e-folding times for runaway temperature increase are in the 10³–10⁴ year range for large ice sheets. The major shortcoming of this analysis is that it is “local”: horizontal temperature advection and driving stress (ice-sheet geometry) are not allowed to react to the changing temperature and velocity fields. A basic damping mechanism is thus excluded. The same objection applies to the study by Reference Schubert and YuenSchubert and Yuen (1982), dealing with exactly the same problem, but pursuing the matter a little bit farther by stating that the creep-instability heating may lead to extensive basal melting, a surge and a new ice age. Their scenario calls on climatically enhanced accumulation to increase the thickness of the East Antarctic ice sheet above a critical value for the onset of explosive creep instability. In a recent paper (Reference Yuen, Saari and SchubertYuen and others 1986), time-scales for the onset of this explosion were investigated and shown to be extremely sensitive to the activation energy and the pre-exponential constant of the ice-creep law, with lower activation energies favouring shorter time-scales.

To study the interaction between the velocity and temperature fields and its influence on the dimensions of large ice sheets, it is necessary to solve the equations governing ice-flow mechanics, together with the thermodynamic equation. To see whether, in reality, conditions can be such that creep instability may occur, such a model should include full time-dependence as well.

To date, the most rigorous approach has been taken by Reference JenssenJenssen (1977), who tried to solve the thermo-mechanical equations for “shallow ice flow”. He introduced a scaled vertical coordinate, transformed the relevant continuity and thermodynamic equation (prognostic equations for ice thickness and temperature) and tried to solve the system numerically. In applying the scheme to the Greenland ice sheet, however, numerical instabilities occurred, forcing the calculations to be interrupted after 1000 years’ integration. Over the last few years, Reference Hutter, Yakowitz and SzidarovszkyHutter and others (1986) have presented numerical solutions for the steady-state two-dimensional (vertical plane) case, based on a treatment given by Reference MorlandMorland (1984). However, in the paper by Hutter it appears that numerical problems are still encountered, reducing the applicability of their solution to cases where the basal sliding velocity component has to be far larger than the movement due to internal deformation.

A different approach to the solution of the fully time-dependent problem was taken by one of the authors in a number of studies of a global nature (Reference OerlemansOerlemans 1982, 1983). Here, the vertical temperature profile was prescribed as a second-order polynomial in height relative to bedrock, with time-dependent coefficients. Together with a prescribed scaled velocity profile, this allows a much more efficient scheme for time integration, while retaining the essential “physics”. However, these thermodynamic calculations are not accurate enough in the lower layers (i.e. where the shear concentrates).

To carry out a more refined analysis, we recently developed a (three-dimensional) fully time-dependent ice-sheet model with high resolution in the lower layers. Section 2 provides a brief description of a flow-line version of the model and of the procedure that was adopted to solve the system numerically. In a following section a schematic (East Antarctic) continental flow line is forced with changing climatic conditions (mass balance and surface temperature) and its transient and steady-state behaviour is investigated. It will turn out from this basic sensitivity analysis that the strain-rate dissipation feed-back s.s. (i.e. without the inclusion of basal melting on ice-mass discharge) only moderately enhances the sensitivity of the basal temperature and mass-flux transient response and is apparently too weak to give rise to large-scale melting and surge-like runaway behaviour with regard to climatic change. It appears that the rather long time-scale on which an instability might develop permits the temperature and stress fields to damp and react to disturbances in a complex interactive fashion. Horizontal heat advection plays a crucial role in this process. Nevertheless, even if the temperature dependence of the flow properties of ice may only change the basal boundary conditions significantly in marginal situations, it still exerts a major influence on the evolution of ice thickness over time in a large polar ice sheet. To conclude and assess the role of thermo-mechanical coupling in more realistic situations, representative East Antarctic flow lines are modelled and their thermal regimes and overall profiles are discussed in terms of dependence on changing environmental conditions.

2. Model Description

The numerical grid-point model used in this study describes ice flow along a flow line and computes the fully coupled velocity and temperature fields in a two-dimensional vertical plane, while accounting for any flow divergence or convergence in the continuity equation. Since it is rather inconvenient to integrate the thermodynamic equation in the computational stage on a grid fixed in space, as the upper and lower ice boundaries will generally not coincide with grid-points, a new vertical coordinate ζ, scaled to the local ice thickness, has been introduced. Following Jenssen (1977), this stretched dimensionless vertical coordinate is defined by

such that ζ = 0 at the surface and ζ = 1 at the base for all values of x and t. In this definition, Η is ice thickness (m) and h is bedrock elevation (zero at sea-level) in the usual Cartesian coordinate system (x,z,t), in which the x-axis is horizontal (x = 0 at the divide) and the z-axis is pointing upwards. The way equations are transformed into this non-orthogonal system (x,ζ,t), which will be used below, can be found (for example) in Reference HaltinerHaltiner (1971). In this approach, velocity and temperature fields are now computed at levels of constant ζ (the layer interfaces). For the experiments discussed in this paper, the vertical domain was subdivided into ten layers concentrated towards the base, with a lowermost grid spacing of Δζ = 0.02, which proved to be sufficient to capture the essential characteristics of the variables in the model.

Ice deformation is assumed to result from shear strain; the shear-stress distribution in the model is given by

where the usual assumptions (sufficiently small bedrock and surface slopes, grid-point spacings an order of magnitude greater than ice thickness) are applied. In this expression g is the acceleration of gravity and ρ is the ice density (910 kg m⁻³), taken to be constant.

Substituting this equation in a Glen-type flow law where exponent n = 3, and integrating the resulting equation with respect to the vertical and using the above coordinate transformation, yields an expression for the horizontal velocity u(ζ):

where the basal boundary condition u(1) = 0.

The flow-law coefficient A(T*) introduces the temperature dependence of ice deformation in the model. Apart from factors such as crystal-fabric and impurity content, that might be equally important, laboratory experiments suggest that A(T*) obeys an Arrhenius relationship:

In this expression a is specified below, R is the universal gas constant (8.314 J mol⁻¹ K⁻¹), Q is the activation energy for creep, and T* is the absolute temperature corrected for the dependence of the melting point on pressure (Τ* = Τ + 8.7 10⁻⁴ Ηζ, where Τ is expressed in K). m is a tuning parameter and proved to be necessary in order to adjust slightly the height-to-width ratio in some of the model runs. Its value generally comprised between 1 and 10 and it may be thought of as an implicit way to include the softening effect of crystal-fabric and impurity content. Unless indicated otherwise in the text, the following values were chosen for a and Q: T* < 263.15 K, a = 1.14 10⁻⁵ Pa⁻³ year⁻¹, Q = 60 kJ mol⁻¹; Τ* ≥ 263.15 Κ, a = 5.47 10¹⁰ Pa⁻³ year⁻¹, Q = 139 kJ mol⁻¹, so that A(T*) lies within the limits put forward by Reference Paterson and BuddPaterson and Budd (1982).

Vertical motion (in m/year, negative downwards), which occurs as a result of accumulation and vertical strain, is calculated from the incompressibility condition, yielding

with the kinematic boundary condition at the upper surface given by

and M is the local mass-balance (m/year), positive in the case of accumulation. Generally the calculated w at the ice-bedrock interface approached zero (the present case, if there is no sliding) within 0.01 m/year, revealing that the scheme employed conserves mass very well.

In order to adjust the flow parameter A during the calculation, the temperature distribution within the ice sheet is obtained from the thermodynamic equation:

Here Τ is absolute temperature (K), t is time (years), k is thermal conductivity (6.62 10⁷ J m⁻¹ K⁻¹ year⁻¹) and c _p is specific heat capacity (2009 J kg⁻¹ K⁻¹). In this equation, heat transfer is considered to result from vertical diffusion (first term), horizontal advection (second term), vertical advection (included in the third term, together with the various correction terms that account for time-dependent layer geometry) and internal deformational heating, due to horizontal shear strain-rates. Boundary conditions follow from the mean annual air temperature at the upper surface. At the base, ignoring heat conduction in the bedrock below, the geothermal heat flux is incorporated in the basal temperature gradient, taken to be 2 K/100 m. This corresponds to a heat flux of 4.2 × 10⁻²W m⁻², considered to be a representative value for East Antarctica. Whenever the pressure-melting point is reached, temperatures are set equal to that value. At the divide, all horizontal derivatives are zero.

In the model experiments discussed below, the bedrock response has also been taken into account, as the basal temperature distribution turns out to depend critically on total ice thickness. We considered a viscous asthenosphere, with lithosphere deflection that followed local hydrostatic equilibrium and an ice – mantle-rock density ratio of 0.3:

In this expression, h₀ is the undisturbed bedrock topography and D _a (10⁸ m² year⁻¹) is a diffusion coefficient. With a typical length scale of 1000 km, this equation leads to a relaxation time for bedrock adjustment of around 10 000 years (Reference Oerlemans and van der VeenOerlemans and van der Veen 1984).

Finally, time-dependent evolution of ice thickness is described by

where HU, the vertically integrated mass flux per unit width, is given by

and b(x) is the width distribution along a flow line. In this paper, only “half” an ice sheet is actually computed, where boundary conditions are taken as zero thickness gradient (divide) and there is prescribed ice thickness at the edge (zero, unless indicated otherwise in the text).

In order to solve (equ.9) numerically, this equation is written as a diffusion equation, with a “vertically integrated” diffusion coefficient that is given by the scalar component of (equ.10). We opted for a finite-difference approach, of the type implicit-in-time / central-in-space. A staggered grid in space was employed, in effect calculating mass fluxes between grid points with a mean diffusivity. Smoothing in this way turns out to keep the integration stable. The resulting finite-difference equations are usually quite readily solved by an explicit integration scheme. However, such a scheme has the important drawback that in order to preserve stability, time steps necessarily have to be small. Alternatively, a semi-implicit scheme (differing from a fully implicit method, in that diffusivity is evaluated at the old time step) turned out to work very satisfactorily. Although not unconditionally stable, partly due to strong non-linear coupling through the source terms, this scheme allows much larger time steps to be taken. The resulting tri-diagonal linear systems are then easily solved by a Gaussian elimination method.

A comparable approach was applied to the solution of the thermodynamic equation. Here, the terms involving ζ-derivatives were made implicit, leading once more to a set of tri-diagonal equations that had to be solved at every grid point. However, a note about the horizontal advective terms is in order here, because replacing the derivatives by central differences turned out to generate oscillations in the solution. This problem is usually circumvented in diffusion-convection equations with a high Peclet number (i.e. a constant proportional to the ratio of advective velocity and diffusivity) by introducing an artificial horizontal diffusion process (e.g. Reference Mitchell and GriffithsMitchell and Griffiths 1980). We use an “up-wind” differencing scheme that can be shown to introduce an artificial horizontal diffusivity equal to uΔx/2, where Δx is the grid spacing. Apart from stabilizing the integration, this influences results only marginally, as in most cases the associated artificial heat transfer turns out to be an order of magnitude smaller than the horizontal heat-advection term. In the experiments discussed here, a spatial resolution of 50 km was chosen, allowing time steps of 50–100 years.

For a more detailed description of the full three-dimensional model and the numerical procedures involved, the reader is referred to Reference HuybrechtsHuybrechts (1986).

3. Basic Sensitivity Versus Environmental Conditions

In order to investigate basic sensitivity of a polar ice sheet with respect to environmental conditions (accumulation and sea-level temperature) and to assess the role of thermo-mechanical coupling, a series of (in the first instance) steady-state experiments was set up. For this purpose we considered a “clean” flow line of uniform width, with a flat bedrock profile that extended from the ice divide at x = 0 to a fixed-edge position at x = 1500 km. First of all, a reference experiment was defined as follows: mean annual sea-level temperature and atmospheric lapse rate were set at 253 K and 12 K/1000 m respectively. Over the Antarctic ice sheet, the accumulation rate appears to be strongly correlated with surface temperature T _s[K] (e.g. Reference Robin.de, Johnsen and Robin.deRobin and Johnsen 1983), as the amount of precipitation is limited by the amount of water vapour that can be carried by the atmosphere at any temperature:

so that at a surface elevation of 3000 m, M = 0.05 m/year, T _s = 217 Κ and at sea-level M = 0.625 m/year respectively. With this parameterization the mass balance depends (in climatic-change experiments) both on elevation and on sea-level temperature.

Initially, computations started with zero ice thickness and bedrock elevation. Depending on the mass balance, a steady state was usually established after 200 000–300 000 years’ integration. With Δx = 50 km, Δt = 50 years and ten layers in the vertical, this required (including graphical output) about 300 CP s for every 100 000 years and 130 000 octal execution field length on a CDC-Cyber 175–750 computer. A realistic height-to-width ratio was obtained by setting m = 10 in (equ.4).

Basic sensitivity of the model to accumulation and/or sea-level temperature changes is shown in Figure 1. Keeping the accumulation distribution fixed with respect to the 253 Κ run (plusses) and increasing/decreasing sea-level temperature only results, not surprisingly, in decreased/ increased mean ice thicknesses, whereas the reverse is true for mean basal temperature. The sensitivities found here reflect the combined effect of many competing feed-back mechanisms, such as the positive feed-back of the changes in elevation on surface temperature, and the opposing effects of increased ice thickness that lead to basal warming, together with the influence due to differences in deformational heat release (although relatively unimportant in this case). However, since the differences in mean basal temperatures more or less follow the applied forcing, these additional feed-backs more or less cancel each other out in this particular experiment. In the 263 Κ model run, the whole base was actually near or at the melting point, explaining in part the lower sensitivity found at that side. Circles in Figure 1 refer to a sensitivity experiment in which the temperature distribution of the ice sheet was kept constant, so in effect showing the effects of changes in accumulation on ice thickness due to ice dynamics only. In this case there exists a negative feed-back of accumulation on surface elevation: higher elevations tend to result in decreased precipitation rates.

Fig. 1. Steady-state mean ice thickness and mean basal temperature versus environmental conditions. + : fixed accumulation depending on x only; x: accumulation also depending on surface temperature; o; fixed englacial temperature distribution; 253 Κ : reference experiment.

However, in the full-model run (crosses), where sea-level temperatures and (therefore) accumulation rates vary, and full thermo-mechanical coupling is considered, warmer conditions actually lead to slightly increased ice thicknesses, in spite of increased dissipation and the rather significant dependence of ice thickness on changes in surface temperature alone. As the concomitant rise in mean basal temperature appears to be quite small, this means that on the whole, and in this particular situation, increased advection of colder ice towards the basal shear layers following a climatic warming largely compensates for the combined warming effects of: (i) increased surface temperature, (ii) increased ice thickness, and (iii) increased mass flux and dissipation. So, generally speaking, the mutual coupling of temperature and velocity seems to have a stabilizing effect on the steady-state configuration of a large polar ice sheet like East Antarctica that is exposed to different climatic environments within realistic ranges: the integrated effects on mean steady-state ice thickness of roughly doubling accumulation rates and raising sea-level temperature by 10 Κ (a typical glacial–interglacial climatic shift) appear to be of similar magnitude, but of opposite sign.

Fig. 2. Differences in ice thickness and basal temperature following a linear 10 Κ glacial–interglacial climatic shift over 10 000 years starting at t = 0. (1): x = 0 km; (2): x = 40 km; (3): x = 800 km; (4): x = 1200 km; x = 0 at the divide.

Up to now only steady and mean states have been investigated. How these ice sheets actually evolve from one state to another at different locations is a different matter, realizing that the various terms in the heat equation all operate on different time-scales. Figure 2 shows the evolution of ice thickness and basal temperature at four locations (labelled in increasing order outwards to the margin) following a linear 10 000 year glacial-interglacial temperature increase of 10 K, starting from a steady-state 243 Κ “glacial” ice sheet (this is the state marked by “x” in Fig. 1). Basal temperatures (lower panel) initially increase, with higher values down-stream, but later they are counteracted by advection of cooler ice, resulting in a slight basal cooling. After that, the temperature rises slightly again, influenced by the arrival of a strongly diffused “warm wave”, as a result of increased surface temperatures conducted towards the base. Note also that near the divide (where there is hardly any dissipation) basal layers are actually cooling, as a result of increased vertical advection of cold ice. To demonstrate the specific role of the various terms in the heat equation in this experiment, Figure 3 is included, showing a typical situation in the basal layers (here for ζ = 0.90) at x = 800 km (corresponding to curve (3) in Fig. 2). Note that because the basal layers slope slightly upward at this particular location (as a consequence of the initial flat bedrock profile and isostatic adjustment), the vertical velocity component, and hence vertical advection, are of minor importance here.

The evolution of ice thickness (upper panel in Fig. 2) follows the opposite trend to some extent. Increased accumulation rates initially lead to increased ice thicknesses. This trend is then reversed as a result of enhanced shear strain due to a warmer base. Eventually, this “lowering wave” travels up-stream and reaches the divide in a much weakened form about 20 000 years later. In the light of these results, it may thus very well be possible that, at present, (mainly outer) parts of the Antarctic ice sheet are actually lowering as a delayed response to the increased accumulation rates associated with the last interglacial warming. Apparently, dependence of ice deformation on temperature adds a very long time-scale to the system. Furthermore, in these experiments no sign of runaway behaviour was found. The dissipation – strain-rate feed-back seems to be far too weak to give rise to massive basal warming, at least not after a typical glacial–interglacial climatic shift.

Fig. 3. Evolution of the heat transfer terms for x = 800 km and ζ = 0.90, corresponding to the experiment in Figure 2.

In order to illustrate in more detail the behaviour of the model ice sheet with respect to sudden climatic changes and the role of the dissipation – strain-rate feed-back, another experiment was carried out, in which at time = 0 sea-level temperature was suddenly increased from 243 to 253 Κ (including the associated doubling of the accumulation rate). In this (drastic) experiment, some changes were made in the formulation of the original model: the pressure-melting-point limit on ice temperature was removed and the activation energy in the creep law was set at 60 kJ mol⁻¹, irrespective of temperature. This then clearly isolates the effect of horizontal heat advection and facilitates comparison with respect to previously mentioned studies on creep instability (Clarke and others 1977, Yuen and others 1986).

The resulting evolution of ice thickness and basal temperature at location x = 1000 km are shown in Figure 4.Curves labelled (2) refer to an experiment in which feed-back on strain-rates was excluded, in effect using the 243 Κ “glacial” ice-sheet temperature distribution to calculate the flow parameter A. Comparing (2) with the full-model run (1) clearly shows the enhanced basal warming that results from the positive feed-back of strain-rates on dissipation. The resulting basal warming of about 4 Κ corresponds to an approximate tripling of the horizontal velocity, from 17 to 50 m/year (not shown here). Farther down-stream, basal temperatures slightly exceed the pressure-melting point, though this does not have much influence on the behaviour of the model. There appears to be no exponential time-scale connected with this feed-back loop, and after about 8000 years of rising temperature and lowering elevation the trend is reversed in much the same way as depicted in Figure 2. Although difficult to assess from this experiment, part of the stabilizing effect is no doubt due to a reduction of shear stress following decreasing ice thickness. What turns out to be more important, however, is the damping effect connected with horizontal heat advection, as is immediately clear from the model runs labelled (3), in which the horizontal heat-advection term was “switched off. In this case, creep instability did develop, forcing calculations to be halted after 5000 years’ integration. Decreasing the time step confirmed that this instability is not a mere numerical artefact. This runaway situation compares well with the (super-)Exponential model behaviour as found by Yuen and others (1986). Although the qualitative behaviour did not change to a great extent, increased activation energies and higher initial velocities were found to shorten significantly the characteristic time-scale for explosion. For instance, setting Q = 139 kJ mol⁻¹ for temperatures above 263.15 Κ (corrected for pressure melting) led to an explosion within 1000 years at this particular location.

Fig. 4. Basal temperature

and ice-thickness

evolution following a sudden 10 Κ climatic warming at x = 1000 km. (1): full model; (2): flow parameter independent of temperature; (3) : no horizontal heat advection.

Thus it may be concluded from this experiment that the development of a creep instability is essentially a consequence of a basic shortcoming in model formulation, namely the neglect of horizontal heat advection. Moreover, if widespread basal melting were to result from enhanced warming due to the velocity-temperature feed-back, only marginal areas in the outer part, where temperatures are already close to melting, would be affected. This point will be taken up in the discussion of the East Antarctic ice sheet.

4. Calculation on Selected East Antarctic Flow Lines

To investigate implications for the East Antarctic ice sheet, with emphasis on the glacial–interglacial contrast, four flow lines thought to be representative of different glaciological environments were modelled (Fig. 5). Bedrock and surface elevation were digitized on a 50 km grid along a flow line, taken from the Antarctic glaciological and geophysical folio (Reference DrewryDrewry 1983). Present-day accumulation rates (not shown) were taken from a preliminary map produced at the Scott Polar Research Institute, that also provided a tape with the original accumulation and surface-temperature point measurements (personal communication from D.J. Drewry). Other parameters for these flow lines are listed in Table I:

Fig. 5. Location of modelled flow lines.

Bed elevation was assumed to be in isostatic equilibrium with present-day ice thickness. Figure 6 shows the geometry of the modelled flow lines, representing steady states under present interglacial conditions. Although the experiments conducted so far clearly indicate that the Antarctic ice sheet is probably never in a steady state, the tuning parameter m could be chosen (independently of x) such that the modelled surface elevations approached measured surface elevations within 100 m or less, which we feel proves the validity of the model equations. Of these flow lines, A (Vostok flow line) and D (Dome Argus – Lambert Glacier) are converging and debouching into an ice shelf, where present-day ice thickness at the grounding line is a boundary condition. Flow lines Β (down-stream of Dome C – Dumont d’Urville) and C (ridge Β – near Casey Station) are completely land-based and of a diverging nature (with zero ice thickness at the margin as a boundary condition). As our present interest is focused around thermo-mechanical response behaviour, changes in the length of flow lines due to changes in sea-level or possible grounding-line movement were isolated from the experiments.

Fig. 6. Fig. 6. Steady-state “interglacial” model geometry of selected flow lines (see Fig. 4).

Fig. 7. Steady-state “interglacial” basal temperature distribution of selected flow lines (see Fig. 4). Thick lines represent pressure melting.

Figure 7 shows the corresponding steady-state basal temperature distribution during present interglacial conditions. In general, basal temperatures tend to rise towards the edge, and appear to be controlled to a great extent by local ice thickness, due to the insulating effect of ice. Obviously, if widespread basal melting were to result from climatic warming and enhanced accumulation, flow line C would turn out to be the most vulnerable. This flow line passes through Wilkes Land, where the thickest ice on Earth (up to 4.5 km) is found. Time-dependent modelling confirms the reality of this statement. Least vulnerable appears to be the flow sector debouching into Lambert Glacier and Amery Ice Shelf, where the base is well below melting. Note the warm base close to the divide of flow line B.

Fig. 8. Steady-state interglacial–glacial surface elevation and basal temperature differences of selected flow lines (see Fig. 4).

Figure 8 shows results of a comparison between a glacial and interglacial steady state of surface elevation and basal temperatures. For the glacial climate run, sea-level temperature was lowered by 10 K. Accumulation rates along the respective flow lines were then calculated in a way similar to that described by Reference Lorius, Jouzel, Ritz, Merlivat, Barkov, Korotkevich and KotlyakovLorius and others (1985). First, the temperature of formation of precipitation T _f was estimated, following indications by Robin and Johnsen (1983) and Reference Lorius, Jouzel, Ritz, Merlivat, Barkov, Korotkevich and KotlyakovLorius and others (1984, and references there): T _f close to the temperature T _i that prevails over the surface inversion layer, inversion strength linearly related to T _s and typically 15 Κ over central Antarctica, leading to: Tf[K] = 0.67 T _s[K] + 88.9. Then accumulation rates were calculated from the product of its present value, times the ratio of the derivatives of the saturation vapour pressure, with respect to T _f for the time of precipitation, and with respect to T _f for present conditions. For surface temperatures prevailing over central Antarctica, this glacial–interglacial temperature shift roughly halves present-day accumulation rates. As a general picture, surface elevations tend actually to be somewhat lower during steady-state interglacial conditions, whereas the rise in basal temperatures remains below the applied sea-level forcing (10 K). Ice thickness, and consequently bedrock topography, appears to be the main factor in explaining the differential response to the imposed glacial–interglacial contrast. In short, elevation increase as a result of imposed warming is most pronounced in the parts of the East Antarctic ice sheet where on the whole the ice is thickest and, not surprisingly, basal temperature changes the least. Advection due to increased accumulation seems to counteract (proportionally to initial ice thickness) the warming effects of increased ice thickness and higher surface temperatures, as a result of vertical diffusive heat conduction and increased deformational heating.

The changes in surface elevation (amounting to a 0.7 fraction of total ice thickness in the steady state) all show a steeper edge and a smaller slope immediately up-stream. This results from the warmer base and higher velocities in the outermost areas and may also be due in part to a fixed outer boundary condition (where the surface-elevation shift is zero, not shown on the graph). Although strongly modulated by local ice thickness, the curves of glacial–interglacial temperature shifts all display a somewhat similar shape, with a strong maximum towards the edge and a smaller one close to the divide. The minimum at the divide results from the combined effects of increased ice thickness and surface temperatures due to vertical diffusive heat conduction and increased vertical advection of cold ice. This minimum in absolute terms appears to be controlled by ice thickness, in the sense that the thicker the ice, the more important the cooling effect, which follows increased accumulation relative to warming due to vertical diffusive-conduction effects. A second minimum about half-way along the modelled flow line appears to be connected with the development of a reversed vertical temperature gradient in the upper layers, due to enhanced horizontal heat advection. Farther down-stream the overall impact of dissipation (with a contribution due to the strain-rate – dissipation feed-back) becomes apparent.

As shown above, thermo-mechanical coupling adds a very long time-scale to the system, implying that a steady-state description may be quite inappropriate for deducing the present state of the ice sheet with respect to “glacial” conditions. This is demonstrated in a last experiment, in which the four flow lines were forced with a linear 10 Κ glacial-interglacial shift over 6000 years, starting from a steady glacial state 16 000 years ago, as a general climatic trend deduced from Antarctic ice-core data (Robin 1983). A final “interglacial” steady state was then established only after 200 000–300 000 years’ integration. The resulting shifts in surface elevation at present (time = 0) are shown in Figure 9. They amount to around 100 m in central areas, whereas more coastal areas appear to have thinned slightly under the influence of an increase in basal temperature. It may be mentioned here that this situation compares well with the findings of Lorius and others (1984), on the basis of the total gas content of ice cores.

Fig. 9. Modelled difference in surface elevation of selected flow lines since the beginning of the present interglacial (see Fig. 4).

Of particular interest here is Vostok Station, as changes in elevation constitute a complicating factor in interpreting the palaeoclimatic δ¹⁸O signal of the ice core obtained there (Lorius and others 1985). According to our calculations, the change in elevation since the beginning of the interglacial 16 000 years ago amounts at this time to about +94 m. The model predicts that the surface will rise another 6 m during the next 6000 years; after that, if the climate were unaltered, the surface would be lowered by 150 m before settling down to a new steady state. According to this calculation the inferred 16 000–0 B.P. climatic contrast in the core would be underestimated by about 0.8 Κ, due to changes in elevation.

5. Conclusion

In this study, the response of the East Antarctic ice sheet to changing climatic conditions has been investigated by using an efficient numerical model with full thermo-mechanical coupling. With emphasis on the glacial-interglacial climatic contrast, our model calculations have shown that the mutual coupling of temperature and deformation has a stabilizing effect on the steady-state configuration of the East Antarctic ice sheet. In time-dependent experiments, the model shows a relatively fast response of ice thickness to increased accumulation rates, followed by a lowering trend (due to warmer bases) that travels slowly up-stream. Time-scales for this reversal are of the order of 10³ years in the margins and 10⁴ years in central zones.

According to the model results, creep instability seems a rather unlikely candidate to initiate surging of parts of the East Antarctic ice sheet, mainly because there appears to be no fast time-scale connected with the temperature-dependent creep behaviour of ice and the damping effect associated with, in particular, horizontal heat advection. As pointed out, only Wilkes Land (where at present the thickest ice on Earth is found) may be vulnerable to widespread basal melting, due to enhanced dissipation following a climatic warming. Here, a sophisticated treatment of basal hydraulics and sliding is required in order to show how rapid ice-mass discharge could increase.

Our model calculations indicate that surface elevations in central East Antarctica have risen by around 100 m since the beginning of the present interglacial, whereas coastal areas may have thinned slightly. Ice thickness, and consequently bed topography, appears to be a main factor in controlling the response of the basal temperature field. In addition, thermo-mechanical coupling appears to add a long time-scale to the system, indicating that the Antarctic ice sheet is probably never in a stationary state. This suggests the use of a model of the present type for interpretation of the climatic signal in ice cores. Particle trajectories and ice age–depth profiles can easily be calculated, even in a transient mode.

The strong dependence of basal temperatures and general climate sensitivity on geometry and topography necessitates the extension of this study to a fully three-dimensional calculation. With present-day supercomputers it is certainly possible to repeat the sensitivity experiments with a 100 km grid for the entire Antarctic ice sheet. Such a study is in progress.

Acknowledgement

Philippe Huybrechts is supported by the Belgian National Fund for Scientific Research (N.F.W.O.) and is sponsored in part by the Belgian Ministry of Science Policy under contract ANTAR/04.