Introduction

The Antarctic ice sheet seems to have undergone important changes during the Pleistocene glacial cycles (e.g. Reference Denton and HughesDenton and Hughes, 1981). In particular, the smaller marine-based West Antarctic ice sheet (WAIS) appears to be very susceptible to changes in environmental conditions. This ice sheet largely rests on a bedrock below sea level, even if isostatic depression were not accounted for, and may have completely filled the continental shelf during maximum glaciation. There seems to be a general agreement that this was indeed the case during the Wisconsin in the Weddell Sea (Elverhøi, Reference Elverhøi1981), although views differ whether the shallow Ross Sea was completely land-ice covered (Stuiver, Reference Stuiver, Denton, Hughes, Fastook, Denton and T.J1981), or there was only a minor expansion (Drewry, Reference Drewry1979). Also, due to specific bathymetric conditions (the bedrock slopes downward towards the interior), the WAIS may have collapsed altogether following an interglacial retreat (Reference Thomas and BentleyThomas and Bentley, 1978), or during periods of higher temperatures that weakened the ice shelves and reduced their buttressing effect (Mercer, Reference Mercer1978; Reference Thomas, Sanderson and RoseThomas and others, 1979). However, the effectiveness of the latter mechanism of grounding-line instability, remains subject to considerable debate (e.g. Van der Veen, Reference Van der Veen1985). The larger land-based East Antarctic ice sheet (EAIS), on the other hand, emerges as a comparably more stable feature. However, due to the huge amounts of ice involved here, even relatively small fluctuations have a large impact, for instance when global changes of sea level are considered.

Fluctuations of the Antarctic ice sheet may result from changes in envionmental conditions in several ways. Lower ice temperatures in the basal shear layers will reduce deformation rates, thereby resulting in a thickening. Changes in accumulation rates (ablation is probably negligible, especially when the interest is in colder climates than today) are another factor. Lower global sea-level stands reduce the thickness necessary for grounding in regions where the bed is below sea level. A similar effect is changes in bed elevation due to isostatic bedrock adjustment. These mechanisms raise a number of fundamental questions concerning the Antarctic ice sheet's basic behaviour: i) what is the relative importance of changes in these environmental conditions; ii) what are the response time scales connected with these changes; iii) what is the spatial distribution of ice-sheet fluctuations; and iv) how do all these effects interplay during a glacial-interglacial cycle? The answers to these questions are important in the light of the interpretation of climatic information in ice cores (Byrd, Vostok and Dome C), a better understanding of ice-sheet histories from proxy data, and the contribution of Antarctic ice-sheet volume changes in world-wide sea level.

A steady-state sensitivity study on the glacial-interglacial contrast, employing a time-dependent numerical model that covers the entire Antarctic ice sheet, forms the subject of another paper (Huybrechts, in preparation). In the present paper we concentrate on transient behaviour, and attempt to simulate the last glacial-interglacial cycle, starting at 160ka B.P. This experiment was, among other things, suggested by the observation that it took 30 000 to 40 000 years for the Antarctic ice sheet to grow to full glacial extent after sudden glacial conditions were imposed. Also, in view of the long time scales connected with thermo-mechanical evolution (Reference Huybrechts and OerlemansHuybrechts and Oerlemans, 1988), it appears unlikely that a steady-state description is really appropriate to reconstruct the ice-sheet history. As a second consideration, modelling a full glacial-interglacial cycle is thought to be a good validation for the model: the model should not only be able to simulate a glacial build-up, but be capable of reproducing a deglaciation as well. More speculatively, such an experiment may yield information on the present evolution of the ice sheet as well.

The Model

A full account of the mathematical equations and numerical techniques governing the model is to be published elsewhere (Huybrechts, in preparation). It is essentially a three-dimensional extension of the thermo-mechanical flowline model in Huybrechts and Oerlemans (1988) to include a numerical ice-shelf model, grounding-line dynamics and a more sophisticated treatment of the bedrock (mechanically as well as thermodynamically). In short, the model basically integrates a continuity equation for ice thickness forward in time, where the velocities are found by solving the stress equilibrium equations subject to appropriate assumptions in the respective ice flow regimes, and substituting for a temperature-dependent flow law, which is of Glen type with exponent n = 3.

Three distinct flow regimes are considered, namely grounded ice flow, where the normal stress tensor is isotropic, ice-shelf flow, where there is no shearing in horizontal planes and a balance exists between the pressure gradient force and gradients in normal stress deviatore, and a transition zone in between at the grounding line, where all stress components contribute to the effective stress in the flow law. For numerical reasons, the model ice shelf actually extends all the way to the edge of the square grid, where boundary conditions for freely floating and uniform spreading ice shelves are applied. This means that ice-shelf calving is not considered explicitly, but “cutting off” the ice shelf closer to the coast appeared not to influence the position of the grounding line. Following an analysis by Bindschadler (Reference Bindschadler1983), basal sliding in the model is of Weertman-type corrected for the effect of subglacial water and restricted to regions that are at the pressure melting point. The equations for the ice dynamics are quite similar to Herterich (Reference Herterich1988) for the grounded ice and MacAyeal and Thomas (Reference MacAyeal and Thomas1986) for the ice shelf, although these models do not account for changes in ice thickness or ice-sheet domain (which is our prime interest, in fact).

The temperature of the flow properties of ice is then accounted for by an Arrhenius-type equation (Paterson, Reference Paterson1981), where the temperature distribution is found by solving the thermodynamic equation in the ice and the bed below (down to 2 km rock depth, where a geothermal heat flux of 54.6 mW m⁻² is the boundary condition). In the ice shelf, the effects of basal melting and/or basal accretion and of spatially varying density and ice stiffness were ignored. Instead, a steady-state linear vertical temperature distribution, with a surface temperature of −18° C (representing present day ice-shelf conditions and allowed to vary with climatic change) and a fixed basal temperature of −2°C, was assumed to exist within the ice shelves.

The Earth deformation model consists of a viscous fluid asthenosphere enclosed by a uniform, thin and elastic lithospheric shell. For this rheology, the following values were used; lithospheric flexural rigidity: 10²⁵Nm; mantle density: 3300 kg m⁻³; asthenosphere diffusivity: 0.5 10⁸ m² yr⁻¹. Since flexure of a rigid plate is a linear process, isostatic displacements due to varying ice and water loading can be calculated at each gridpoint separately, and are consequently superimposed in a square with sides of 800 km.

Model inputs are bed topography, surface temperature, accumulation rate, sea-level stand and an initial state (that may be a thin slab of ice of equal thickness). The model then essentially outputs the time-dependent three-dimensional ice-sheet geometry. Flow and temperature calculations are made on a three-dimensional grid with a horizontal resolution of 40 km and 10 layers in the vertical (gradually decreasing in ice thickness towards the base, where the shear concentrates), employing an efficient finite-difference scheme based upon the Alternating-Direction-Implicit method. The complete model has been run on a CRAY-2 computer, and takes around 18 minutes CPU time for a 10 000 year integration. Time steps are 40 years.

Data and Reference Runs

Bedrock elevation, surface elevation and ice thickness data for the present ice sheet come from an original 20 km-digitalization by the Australian group (Reference Budd, Jenssen and SmithBudd and others, 1984) of the Drewry (Reference Drewry1983) map folio series, and have been locally complemented beneath ice-shelf areas with more recent data (Herrod, Reference Herrod and Konnen1986). The square grid, laid out over a stereographic projection with standardparallel 71°S, is centred at the pole and comprises a 141 × 141 grid point with Δx = Δy = 40 km. Distortions due to the sphericity of the Earth are disregarded. With respect to the basic climatic data sets (10 m ice temperature and surface accumulation rate), a quadratic interpolation and smoothing on original measurements (approximately 1000 data points as compiled at Scott Polar Research Institute (SPRI), supplied to us by D. Drewry) were performed, because a satisfactory parameterization could not be obtained. In particular, the accumulation rate appeared to be hard to parameterize in terms of geometric model output. The surface temperature data are then perturbed in response to background temperature changes (assumed uniform over the ice sheet) and changing ice-sheet geometries according to an atmospheric lapse rate of 5.1K/100m below the 1500 m contour and 14.3/1000 m above the inland plateau. These values are suggested by a multiple regression study on the SPRI data (P. Fortuin, personal communication). Accumulation rates in different climates are then calculated proportional to the derivative of the saturation vapour pressure at the temperature prevailing above the surface inversion layer, following a model suggested by Lorius and others (Reference Lorius1985). For a 10°C decrease in surface temperature, the resulting glacial accumulation rate is typically around 50–60% of its present value.

Fig. 1. The Antarctic ice sheet, as modelled in an interglacial (left) and glacial (right) steady state. Elevation contours are 0.5 km apart; the lowest contour is the grounding line.

Before embarking on the time-dependent model run, an interglacial reference experiment was first set up in order to assess the pre-exponential proportionality factor in the flow law, that serves a tuning purpose (Reference Huybrechts and OerlemansHuybrechts and Oerlemans, 1988). Calculations started initially with the present ice-sheet configuration. Assuming that the observed bed topography is in steady state with present day ice thickness and still keeping the grounding line fixed to its present position, the coupled temperature and velocity fields were then run 100 000 years forward in time. The resulting ice sheet after the model has relaxed to a stationary state for another 100 000 years, and the grounding line is now allowed to vary according to a flotation criterion, is shown in Figure 1 (left panel). Regarding the fact that the model contains only one tuning parameter and has a maximum number of degrees of freedom, the result appears quite acceptable. Somewhat larger deviations from the present observed state include thicker ice in the Antarctic Peninsula and in the Pine Island and Thwaites Glacier catchment areas (WAIS) and a slight recession of the grounding line at the seaward edge of the most overdeepened EAIS outlet glaciers (Totten, Ninnis and Mertz glaciers). However, this does not necessarily mean that the model is in error, as the present ice sheet may not be in steady state. Also shown in Figure 1 (right panel) is the glacial steady-state run, in which sea level was stepwise lowered by −130 m, background temperatures were decreased by 10°C and concomitant reduced accumulation rates were applied. This ice sheet now extends almost all the way to the continental shelf edge, corresponding with a volume gain of 7.4 10⁶ km³ of ice (32.42 against 25.02), or alternatively, an Antarctic contribution to global sea-level lowering with 15.9m (including the correction made for ice merely displacing ocean water). More details of these and other steady-state experiments will however be presented elsewhere (Huybrechts, in preparation).

Simulation of the Last Glacial-Interglacial Cycle: Results and Discussion

In order to simulate the last glacial-interglacial cycle, two forcing functions were used to drive the model. The Vostok temperature signal (Reference JouzelJouzel and others, 1987) was chosen to perturb the 10 m temperature and accumulation fields and is assumed to be independent of latitude. The sea-level record presents some ambiguities, because benthonic oxygen isotope records from deep sea cores record temperature change as well as global ice volume (Shackleton, Reference Shackleton1987). According to Chappell and Shackleton (Reference Chapell and Shackleton1986), this temperature effect may partially explain the discrepancy between the ice volume record that these records imply, and that derived from the altitude of dated coral terraces around the world, an effect that would even be enforced when the average isotope composition of the former ice sheets is taken into account (Reference Mix and RuddimanMix and Ruddiman, 1984). For these reasons, the more familiar “saw-tooth” sea-level function was chosen, piecewise linearly connecting the points (−140 ka, −130 m; −125 ka, 0 m; −16 ka, −130 m; 6 ka, 0 m) respectively. The maximum deviation from either curve extreme (New Guinea terraces and oxygen isotope ice-volume record) is then generally within 20–25 m (Figure 2, upper panels).

Fig. 2. Forcing (upper panels) and evolution of some model variables (lower panels) during the last glacial–interglacial cycle. The vertical bar in the panel below refers to the glacial (G) and interglacial (I) steady states at Byrd (80°S, 120° W; WAIS). The corresponding bars for Vostok (78.5°S, 106.8°E; EAIS) should be at −9 m and −1 m respectively. They almost coincide due to the counteracting effects of accumulation and temperature changes.

The model calculations start at 160 ka B P., when the ice sheet has been run initially to a stationary state with a sea-level depression of 130 m and a temperature perturbation of −5°C. The evolution of some large-scale variables is shown in Figure 2 (middle panels). As a rough picture, the evolution of Antarctic ice volume appears to be mainly controlled by changes in global sea level, and further modulated by fluctuations in accumulation rate and ice temperature. An apparent feature here is the remarkable low during and shortly after the Sangamon (Eem) interglacial, when the modelled WAIS is substantially smaller than today (although there is no complete collapse). This modelled reaction seems to be primarily a consequence of initial conditions (a fully depressed steady-state “glacial” bed, still enforced by the phase lead of the temperature and accumulation rise) together with specific bathymetric conditions beneath the WAIS (bedrock sloping downwards towards the interior), allowing the grounding line to retreat substantially more than during the last deglaciation. However, it also suggests that a larger disintegration of the WAIS is possible as long as the bed is depressed deeply enough. How the evolution of surface elevation in the East and West Antarctic ice sheets looks in this experiment is demonstrated in Figure 2 (lower panel). Typical elevation changes at Vostok station (EAIS, full line) are generally within 100 m, whereas these modelled fluctuations at Byrd (WAIS, dashed line) are substantially larger, due to important changes in the grounded ice domain.

Figure 3 displays the Antarctic ice sheet at several stages during the last glacial-interglacial cycle. A moderate lowering of global sea level with 35–40 m appears to be enough to initiate grounding in the Weddell Sea, when the ice shelf runs aground on a number of high banks near the present shelf front. The threshold for grounding in the Ross Sea, on the other hand, appears to be much larger and significant grounding here is only produced in later stages of a glacial period. Since the free water depth below the ice shelf generally increases outwards, grounding is also of a more gradual nature as compared to the situation in the Weddell Sea. As a consequence of the long time scales involved, the Antarctic ice sheet has not fully reached its maximum glacial extent until 16 ka B.P., in particular in the Ross Sea. At this time, the Antarctic ice-sheet volume has grown to 31.0 × 10⁶ km³, corresponding with a global sea-level lowering of 12.3m relative to the interglacial reference run. It is important to note that this figure is substantially less than the 35 m put forward by Nakada and Lambeck (Reference Nakada and Lambeck1988), indicating that the northern hemisphere ice sheets may have contained more ice than assumed in their analysis.

The deglaciation of the Antarctic ice sheet, that essentially appears as a partial disintegration of the WAIS, begins in the model around 8000 ka B.P. and is nearly completed by the present time. Although it is clear that this deglaciation occurs at least 3000 years later than what is generally believed, a recent discussion on the matter by Peltier (Reference Peltier1988) also seems to favour an important delay in West Antarctic melting by as much as 7000 years after the onset of northern hemisphere deglaciation. This time delay is also the reason why the maximum ice volume is reached in the early stages of the glacial-interglacial climatic transition: accumulation rates have increased, while the ice-sheet domain has not started to shrink and the surface warming signal has not reached basal shear layers yet. Grounding-line retreat in the model is thus triggered by a rise in world-wide sea level, thereby providing the link between glacial events on both hemispheres. Another important feature here is that this disintegration seems to develop its own dynamics (environmental boundary conditions do not change after 6 ka B.P.) and takes around 6000 years to complete, once the threshold for grounding-line retreat is surpassed. The ice-sheet retreat now follows in a reverse order and occurs earlier in the Ross Sea than in the Weddell Sea. A crucial role is also played by the bed adjustment process: it allows the grounding line to reach depressed inland areas at first, when later on isostatic rebound slows the retreat down and eventually results in a small advance. During this evolution, ice rises in the Ronne-Filchner ice shelf (Korff Ice Rise and Berkner Island) never completely disappear. In the EAIS, changes in ice thickness are less pronounced and reflect the combined effects of accumulation and temperature changes much like the flowline studies discussed in Huybrechts and Oerlemans (Reference Huybrechts and Oerlemans1988), on which a thinning wave due to postglacial grounding-line retreat is superimposed. If the model is run further for another 40 000 years, an ice sheet close to the present interglacial reference state shows up again.

Fig. 3. Stages in the modelled evolution of the Antarctic ice sheet during the last glacial-interglacial cycle. Contour interval is 0.5 km; the outermost line is the 200 m contour. Hsl: world-wide sea-level stand; ΔTsl: temperature change.

Conclusion

According to the results discussed in this paper, the model appears to be able to successfully account for major features of the Antarctic ice sheet during a complete glacial-interglacial cycle. The most important fluctuations are then concentrated in the West Antarctic ice sheet and are to a large extent controlled by changes in sea level. Grounding occurs more readily in the Weddell region than in the Ross region. These results support the hypothesis that the Antarctic ice sheet basically follows glacial events on the northern hemisphere through changes in eustatic sea level. Also, no sign of a catastrophic collapse of the entire WAIS following a glacial-interglacial transition was found, although on the basis of our model results, in the initializing phase, an interglacial disintegration of larger amplitude than the last one cannot be excluded, as long as the bed under the WAIS is depressed deeply enough. Changes in ice volume of the East Antarctic are also of some importance, but, on the longer time scales (>10 000 years), the effects of changing accumulation rates and ice temperature are approximately counteracting each other.

The results presented in this study could, however, still be improved in several ways. For instance, a better quality record for global sea-level depression, especially during the final stages of the Wisconsin, may yield a somewhat different maximal glacier stand, in particular in the Ross Sea. Model output at the present time could serve to remove the steady-state assumption on present bed heights. Also, the model may be used in an iterative fashion to date the Vostok signal more precisely and correct for shifts in surface elevation. In this respect, particle trajectories and ice- age depth profiles could easily be constructed, even in a transient situation. However, it is believed that these refinements are unlikely to influence model outcome in the sense that results would become entirely different.

Acknowledgements

This research would not have been possible without the support of the Alfred Wegener Institute, Bremerhaven for the use of their CRAY-2 facilities and for providing optimal working conditions during an 8-months stay. Special thanks are also due to H. Oerlemans for his encouragement and stimulating discussions. Useful comments by J. Determann, S. Kipfstuhl and H. Decleir and critical remarks on the manuscript by D.R. MacAyeal and K. Herterich were much appreciated. This work was supported by the Belgian National Fund for Scientific Research (NFWO) and in part sponsored by the Belgian Office of Science Policy under contract ANTAR/04.