Introduction

Palaeoclimate modelling is of primary interest because it allows us to gain a greater physical insight into the behaviour of the climate system during the Earth’s history and possibly to help in future climate prediction. For example, the simulation of the last glacial-interglacial cycle seems to be particularly interesting, because this cycle is well documented and consequently its simulation may help to validate climate models. More generally, the climatic state of the Earth is far better defined for the late Pleistocene than for any other geological period (Reference Imbrie, Berger, Imbrie, Hays, Kukla and SaltzmanImbrie and others, 1984). This period is marked by an alternation of ice-sheet advances and retreats (glacial-interglacial cycles). In particular, spectral analyses of the climatic records of the last 450 000 years (e.g. Reference Imbrie, Berger, Imbrie, Hays, Kukla and SaltzmanHays and others, 1976) have suggested a strong link between orbital forcing and the glacial-interglacial cycles.

Up to now, several theories have been presented to explain the link between orbital forcing and the glacial-interglacial cycles. The most popular one is the Milankovitch theory (Reference MilankovitchMilankovitch, 1941; Reference BergerBerger, 1988), which requires that, for a glacial age to occur, northern high-latitude summers must be cold enough to prevent winter snow cover from melting, in such a way to allow a positive value in the annual budget of snow and ice, and to initiate a positive feed-back cooling over the Earth through a further extension of the snow cover and a subsequent increase of the surface albedo.

In order to examine the Milankovitch hypothesis, a two-dimensional climate model coupled to a one-dimensional asthenosphere-ice-sheet model has been developed and used to simulate the last glacial-interglacial cycle (Reference Gallée, van Ypersele, Fichefet, Marsiat, Tricot and BergerGallée and others, in press). In this paper, after a short description of the model, the main results of the last glacial-interglacial cycle simulation are discussed and the mechanisms responsible for the glaciation and the déglaciation in the model are analysed with special attention to the feed-backs.

Present-Day Climate

A two-dimensional climate model which links the Northern Hemisphere atmosphere, ocean mixed layer, sea ice and continents has been validated over the present-day climate (Reference Gallée, van Ypersele, Fichefet, Marsiat, Tricot and BergerGallée and others, 1991). It is a latitude—altitude model. In each latitudinal belt, the surface is divided into at most five oceanic or continental surface types, each of which interacts separately with the sub-surface and the atmosphere. The oceanic surfaces are free of ice or covered by ice, while the continental surfaces are covered or not by snow and there is a Greenland ice sheet. The atmospheric dynamics are represented by a zonally averaged quasi-geostrophic model, which includes a new parameterization of the meridional transfer of the quasi-geostrophic potential vorticity and a parameterization of the Hadley sensible-heat transport. The atmosphere interacts with the other components of the climate system through vertical (luxes of momentum, heat and water vapour. The model explicitly incorporates detailed radiative transfer, surface-energy balances, and snow and sea-ice budgets. In particular, a parameterization of the effects of snow metamorphism on its albedo is included. In short, snow metamorphism is due to water-vapour diffusion in the snow layer and is responsible for an increase of the snow grain-size (e.g. Reference ColbeckColbeck, 1983), which leads to an albedo decrease (e.g. Reference WarrenWarren, 1982). This process is referred to hereafter as the “snow-aging” process. The shifts of the taiga/tundra boundary modify the surface albedo, because of a larger snow albedo over tundra than over taiga. In the model, these shifts are parameterized as a function of the July continental temperature. The vertical profile of the upper-ocean temperature is computed by a mixed-layer model which takes into account the meridional convergence of heat. Sea ice is represented by a thermodynamic model including leads and a new parameterization of lateral accretion. With this parameterization, sea ice is allowed to form in latitude bands where the zonally averaged water temperature is above the freezing point. This allows us to take into account the heterogeneous distribution of the sea-surface temperature within a zonal belt.

Simulation of the present climate shows that the model is able to reproduce the main characteristics of the present-day general circulation. The seasonal cycles of the oceanic mixed-layer depth, the sea-ice extent and the snow cover are also well reproduced. Sensitivity experiments show the importance of the meridional sensible-heat transport by the Hadley circulation in the tropics, of the seasonal cycle of the oceanic mixed-layer depth and of the sea-ice formation in latitude bands where the average water temperature is above the freezing point.

The Last Glacial-Interglacial Cycle

The two-dimensional climate model was asynchronously coupled to a model of the three main Northern Hemisphere ice sheets and their underlying bedrock in order to assess the influence of several factors, including snow-surface albedo over the ice sheets, which contribute to ice ages when the model is forced by astronomically driven insolations and by CO₂ concentration taken from the Vostok ice-core data. The model was first run from 122 kyear BP to the present, taking into account the ice-sheet feed-back but keeping the CO₂ concentration at the Eemian level ≃270 ppmv. This experiment is later referred to as EXP1 and the simulated total-ice volume deviation from the present-day value is displayed in Figure 1. Large variations of ice volume are simulated between 122 and 55 kyear BP, with a fast latitudinal extension of the North American and Eurasian ice sheets starting at 120 kyear BP. The model simulates a maximum ice-volume deviation from the present-day value of 42 × 10⁶ km³ at 19 kyear BP and a total deglaciation as well, this simulation of both the glaciation and deglaciation being a crucial test for evaluating the potentiality of a model to simulate glacial-interglacial cycles. Moreover, the simulated evolution of the individual three northern ice sheets is generally in phase with geological reconstructions. However, the major discrepancy between the simulated and the observed climate lies in the underestimation of the temperature variations and of the southward extent of the ice sheet.

Fig. 1. Variations over the last glacial-interglacial cycle of (i) the deviation from the observed present-day values of the simulated total continental ice volume over the Earth, assuming the Antarctica ice-volume reconstruction adapted from Denton (full line), and (ii) the variation of the global sea-water oxygen ratio given by Reference Labeyrie, Duplessy and BlancLabeyrie and others (1987) and Reference Duplessy, Labeyrie, Blanc, Wanner and SiegenthalerDuplessy and others (1988). The forcing used for this simulation is only the insolation variation at the top of the atmosphere (Reference BergerBerger, 1978).

A second experiment (referred to as EXP2) was undertaken to address the climatic effects of CO₂ variations over the last glacial-interglacial cycle. Including the CO₂ variations over the last 122 kyear improves significantly the temperature simulation when compared to the results of EXP1, while the ice-volume reconstruction is less significantly improved (Reference Gallée, van Ypersele, Fichefet, Marsiat, Tricot and BergerGallée and others, in press). At the Last Glacial Maximum (LGM, at 19 kyear BP) the simulated temperature of the Northern Hemisphere was 3.4°C colder than the present-day simulated value, and the simulated total ice-volume deviation was 53 × 10⁶ km³ (see Fig. 2). A comparison of the results of EXP1 and EXP2 at LGM allows us to attribute 50% of the temperature change and 30% of the ice-volume change in the Northern Hemisphere to the long-term CO₂ change. In addition, a feed-back analysis of the LGM simulated in EXP2 shows that 67% of the computed cooling is due to the surface-albedo increase (with 27% due to the water-vapour feed-back) and 33% is related to the CO2 amplification (with 13% due to the water-vapour feed-back).

Fig. 2. Variations over the last glacial-interglacial cycle of (i) the deviation from the observed present-day values of the simulated total continental ice volume over the Earth, assuming the Antarctica ice-volume reconstruction adaptedfrom Denton (full line), and (ii) the variation of the global sea-water oxygen ratio given by Reference Labeyrie, Duplessy and BlancLabeyrie and others (1987) and Reference Duplessy, Labeyrie, Blanc, Wanner and SiegenthalerDuplessy and others (1988) (dashed line). The forcing used is the insolation variation at the top of the atmosphere (Reference BergerBerger, 1978) and the C02 variation reconstructed from the Vostok ice core (Reference Barnola, Raynaud, Korotkevich and LoriusBarnola and others, 1987). (From Reference Gallée, van Ypersele, Fichefet, Marsiat, Tricot and BergerGallée and others (in press).)

The Feed-Back Mechanisms

Entering the glaciation

Summer insolation is peaked around 123 kyear BP at high northern latitudes. The minimum was reached 11–12 kyear later. For June, it decreased from 545 to 440 W m⁻², a difference which amounts to near 20%. At these latitudes and months, the insolation is characterized by a strong precession signal, but in December obliquity dominates the spectrum at 60° Ν (maxima and minima alternatively occur at 147, 129, 115 and 88 kyear BP). Each time that this type of decrease in insolation occurs, positive feed-backs amplify the response of the climate system to such changes in the external forcing and lead to the build-up of ice sheets not only up to the minimum of insolation but a few thousands of years later (4–5 kyear) due to the inertia of the ice sheets. During the same time, negative feed-backs progressively slow the formation of the ice sheets up to the beginning of their retreat when the ice maximum is reached.

The mechanisms playing a significant role during the ice-sheet initiation phase in EXP2 (where insolation and CO₂ vary with time) are illustrated by comparing, in Table 1, their individual impacts between 122 and 120 kyear BP in the 65–70° Ν latitude band. For example, the variation in the snow-field continental fraction increases from 0 to 1 in July between 122 and 120 kyear BP, revealing a slower melting of the summer snow at 120 kyear BP, The annual maximum of the snowfield albedo at 120 kyear BP is also larger than at 122 kyear BP. This is due to the fact that the snow over the tundra has a greater albedo than that over the taiga and that the tundra extends at the expense of taiga between 122 and 120 kyear BP, because of the decrease of the mean July continental temperature.

Table 1. Comparison between 122 and 120 kyear BP of the insolation, of the CO₂ concentration and of the simulated continental and atmospheric climate characteristics illustrating the ice-sheets initiation in the 65–70° N latitude band. Insolation is given for June (Reference BergerBerger, 1978), CO₂ concentration is an annual mean value (Reference Barnola, Raynaud, Korotkevich and LoriusBarnola and others, 1987), maximum snowfield albedo occurs in winter and other climate characteristics are given for July

These continental temperature variations also impact the heat fluxes above the continental surface. When summer insolation decreases, the continental temperature (T) decreases and the melting of the snowfield is delayed. At the same time, the taiga is replaced by the tundra which increases the albedo of the vegetated surface covered by snow. Both the snowfield and the tundra therefore lead to an increase in the surface albedo, creating a positive feed-back. This cooling is also responsible for a decrease in the upward infrared, of the latent-heat and of the sensible-heat fluxes at the surface (cf. lines 7–9 in Table 1) which feed back negatively on T, but reduces the amount of heat absorbed by the atmosphere. Consequently, the temperature and the water-vapour content of the atmosphere decrease, leading to a decrease in the downward infrared radiation at the surface. In summary, the radiative budget there becomes negative, which induces a net accumulation at the surface, the ablation being suppressed. Moreover, the atmospheric cooling directly increases snowfall in this region. As a consequence of the ablation suppression and of the snow-precipitation increase, an ice sheet starts to form. In the 65–70° Ν latitude band on the North American ice sheet, snow precipitation increases from 0.38 to 0.44 m (ice equivalent) year⁻¹ between 122 and 120 kyear BP while ice ablation decreases from 1.32 to 0 m (ice equivalent) year⁻. These values show that the ablation variations play a major role in the net ice-accumulation variations. This behaviour of the model has been confirmed by sensitivity experiments (for more details, see Reference Gallée, van Ypersele, Fichefet, Marsiat, Tricot and BergerGallée and others, in press). The initiation phase lasts roughly 2 kyear.

With the altitude of the ice sheets increasing, the temperature at their surface decreases (a positive feed-back) which decreases progressively as the snowfall at their top (a negative feed-back). At the same time, the ice sheets extend over the continent. Because of the continentality effect, snowfall at the top of the ice sheets in the interior of the continents decreases even more (a negative feed-back). For example, in the 65–70° Ν latitude band over the North American ice sheet, snow precipitation decreases from 0.44 to 0.25 m (ice equivalent) year between 120 and 110 kyear BP. Consequently, the averaged snow precipitation over the ice sheets decreases and this slows down their growth. The maximum volume of ice is reached, at 110 kyear BP, 4 kyear later than the minimum of insolation.

One should also mention (i) the important amplification of the temperature decrease since 114 kyear BP, due to the decrease in the CO2 concentration in the air, and (ii) the increasing depression of the bedrock below the ice sheets because of their growth.

The deglaciation process

As the insolation already started to increase 4 kyear before the maximum of the ice is reached, this increase begins to become important particularly at the southern edge of the ice sheets which begin to melt. At the same time, the absence of snowfall, particularly over the ice sheets, triggers the “snow-aging” process which is efficient in regions of the ice sheets where the temperature is above −10°C. This decreases significantly the albedo of both the non-melting and melting snow cover, particularly in the south. In the meantime, snow accumulation has remained positive over the northern parts of the ice sheets which definitively creates a north-south flux of ice within the ice sheets. At the southern edge, melting is faster than ice flow and the isostatic rebound is not fast enough to level up the ice sheet, which could prevent future ablation there. This is not the case north of the ice sheets where the isostatic rebound maintains the ice-sheet thickness, allowing the north-to-south flux of ice to continue to be efficient. As a consequence, the ice-sheets’ volume decreases, the height of the ice sheets decreases, and temperature at their surface increases which increases the global ablation.

As an example, the North American ice-sheet evolution between 110 and 98 kyear BP is displayed in Figure 3, for EXP2. Between 110 and 104 kyear BP, the ice-sheet retreat consists essentially in a northward motion of the southern front of the ice sheet; after 104 kyear BP we observed a decrease in both the ice-sheet height and extent. The ice-sheet retreat leads to a decrease in the zonal surface albedo and a further replacement of the tundra by the taiga which feed back positively on the albedo decrease, directly and by reducing the albedo of the continental surfaces covered by snow. This leads finally to a minimum of ice volume reached about 5 kyear later than the insolation maximum.

Fig. 3. The latitudinal extent of the North American ice-sheet crest elevation above present sea level, between 110 and 98 kyear BP, for the experiment EXP2.

Conclusions

It is found in our model that the variations of the ablation are more important for the ice-sheet response than the snow-precipitation variations. A key mechanism in the déglaciation after the last glacial maximum appears to be the “aging” of snow which significantly decreases its albedo. The other factors which play an important role are the ice-sheet altitude, the insolation, the taiga cover, the ice-albedo feed-back, the ice-sheet configuration (“continentality” and “desert” effect), the isostatic rebound, CO2 changes and the temperature-water-vapour feed-back.

These processes and feed-backs are only those which play a major role in the model. Other important ones are missing, like the deep-ocean circulation, clouds, particularly over the ice sheets and dust in the atmosphere. Further experiments therefore remain to be done in order to determine the relative importance of all these processes.