1. Model Description

This paper presents results on an analysis of the Greenland ice sheet to the climate scenario inferred from the GRIP data from 25 0000 years BP to present. The results are based on a moidified and extended three-dimensional ice-sheet model. Its precursor was originally developed by Herterieh (1988) and improved substantially by Calov (1994), Rheologically, the ice is described as a very viscous power-law fluid with ice-flux temperature coupling (Hutter, 1983; Morland, 1984). The upper surface is free to evolve during all calculations, being affected by surface-snow addition and melting processes. The lithosphere is modelled as a heat-conducting rock layer which may respond isostatically to the overburden load as is made more speciflc below. At its lower boundary, the rock layer is subjected to a spatially and temporally constant geothermal heat flow of 42 × 10⁻³ W m⁻² and at the ice-bedrock interface the temperature is continuous, while the transverse heat flow and the sliding frictional heat balance the melting rate, if present. To this purely thermal lithosphere response the following simple bedrock-sinking model is added: it is envisaged that the rock layer is composed of vertical thin rigid cylinders, which can freely move relative to one another in the vertical direction; in this motion they are subject to gravity and buoyancy forces alone that develop when the respective pile is submerged into the asthenosphere. A Maxwellian model for the submergence depth of the lithosphere into the asthenosphere with a prescribed relaxation time accounts for the slow viscous deformation of the asthenosphere. The thermomechanically coupled viscous equations for the Stokes’ flow of the deforming ice are simplified for a flat small aspect-ratio ice sheet by employing the shallow-ice approximation of Morland (1984) and Hutter (1983). The emerging equations are implemented by use of centred finite differences on a staggered regular Arakawa grid. The finite-difference model for both the heat conduction and the bedrock sinking is transformed to a constant interval by implementing the so-called σ transformation, and is non-uniform. i.e. variable grids are used with decreasing mesh size in the vertical direction as the ice—bedrock interface is approached from below; tins achieves better satisfaction of the thermal-transition conditions at the bed. The atmospheric driving mechanisms are incorporated as spatio-temporal prescriptions of the atmospheric temperature, the accumulation as well as the melting rate. The driving temperature is additively decomposed into a spatial distribution of modern surface temperature as presented by Ohmura (1987) plus a purely time-dependent part T _D(t), that mimics the climate changes through the ice ages; i.e. the difference between the present mean temperature of the globe and its temperature at some earlier time t. Precipitation is based on Ohmura and Reeh’s (1991) map of present-day precipitation. This distribution is pre-multiplied with a function that is linearly correlated with T _D(t) such that it is reduced to half at the climate minimum, as shown by Galov (1994). Melting is implemented through a “positive degree-day function” (Braithwaite and Olesen, 1989) but it is improved by accounting for the remelting/refreezing processes of superimposed ice (Reeh, 1991), also as implemented by Calov (1994).

Data for present surficial and basal topographies are based on data asembled by Letréguilly and others (1990). The exact model description has been given by Calov and Hutler (1996). We only mention that the horizontal. △_x, △_y and vertical, △ _z , grid sizes are 40 km and 100 m, respectively, that Glen’s flow law is used with an exponent n = 3 and Paterson’s (1994) temperature-dependence of the rate factor is used; for Pleistocene ice, this rate factor is enhanced by an enhancement factor of E = 3 to account for the larger mobility of ice-age ice.

This paper gives a selection of results compiled in a video tape that shows the evolution of the computed topography of the Greenland ice sheet through the last two ice ages. Computations start 250 kyears ago with an ice-sheet geometry and temperature distribution as well as bedrock topography corresponding to those climate conditions in equilibrium, i.e. as if they had persisted for ever prior to that time. The climate driving from 250 kyears BP to the present is derived by the mean global atmospheric temperature as inferred from the δ¹⁸O temperatures deduced from the ice core drilled during the European Greenland Ice Core Project (GRIP) at 72°34ʹ N, 37°37ʹ W from 1990 to 1992 (see Johnson and others, 1992; Dansgaard and others, 1993). The correlation between T _D(t) and the ¹⁸O isotope ratio is

(1)

The time series of T _D as obtained through this correlation during the last 250 kyears is shown in Figure 1.

Fig. 1.

T_D(t) from GRIP ice core.

Similar calculations performed for one glacial interglacial cycle using another thermomechanical ice-sheet model (Huybrechts and Oerlemans, 1988) and different atmospheric forcing derived from ´¹⁸O measurements at the Greenland ice margin has been published by Letréguilly and others (1991).

2. Steady State as Initial Set-up

Experience with earlier computations indicates that in order to minimize errors due to the missing climate forcing before 250 kyears, calculations through the ice ages are best started from a steady-state configuration at a climate minimum. The GRIP data at 250 kyears BP do not exactly correspond to the Illinoian minimum but this time corresponds reasonably to an ice-age minimum. Thus, the model was run for 100 kyears into a steady state with a global temperature corresponding to T _D from the GRIP data which at that time was 3.94°C below the present temperature (see Fig. 1). The calculations show that steady state is very accurately obtained after 100 kyears; indeed, changes in the total ice volume of the ice sheet during the last 23 kyears of model integration are less than 2%0. The topography is shown in Figure 3a—d. Similarly, the temperature distribution is very close to a steady state. Incidentally, since only steady-state conditions of the ice sheet are sought, the thermal response of the rock bed, that generally results in a prolongation of the response time to a certain climate does not need to be accounted for.

Fig. 3.

Modelled surface elevation views, a. Topography after 100 kyears of model integration (steady state). Initial set-up for calculations driven by GRIP data. b. Surface elevation 127 kyears ago. c. Surface elevation 21 kyears ago. d. Present-day modelled surface elevation (after 250 kyears of model integration).

Several typical features of the calculations of steady state are worth mentioning. Comparing the computed surface topography for the steady-state configuration under the present climatic conditions with the present topography indicates that the north dome generally lends to be too high and too extended, whereas the south dome is somewhat too small. Moreover, both are slightly shifted to the west. In the east and in the north, the margin positions are at the shoreline, i.e. farther out than they are today, whereas in the south and southwest the margin positions agree better with today’s observed margin positions.

Table 1 summarizes the heights above sea level of the summit of the north dome and the two peaks of the south dome from observation at the present time, for the steady-state configuration under today’s climatic conditions and for the steady initial slate described above. Evidently, the computed peaks of the surface topography 250 kyears ago and today differ only moderately from one another even though the atmospheric temperature was almost 4°C lower than at present. The reason is the precipitation model used; precipitation is proportional T _D(t) such that it is half of today’s precipitation at the Wisconsin minimum when T _D =-10°C. As the climate becomes colder, this implies two different processes influencing the snow balance: in regions where in general no surface melting occurs this results in a reduced accumulation rate. On the other hand, for regions having surface melting under present conditions, the snow balance is increased, because there is less or even no melting, and decreases owing to the reduced precipitation rate. The two counteracting processes might perhaps balance out. In summary, one would then expect that the maximum heights are reduced. However, dynamics cannot be ignored: because of the overall colder conditions 250 kyears BP, the ice is stiffer owing to the thermomechanical coupling and ice flux towards the margins is therefore reduced. The net effect of all these competing processes is a slight increase of the maximum heights of the summits of the north and south domes, in the latter somewhat even more pronounced than in the former. This tendency can also be corroborated in the three variables listed in Table 2. The ice-covered area and ice volume are larger for the steady conditions 250 kyears BP than for the present-day computed conditions, whereas the temperate areas are smaller, because the ice temperature is reduced. It should also be mentioned that the location of the north dome in the steady state of 250 kyears ago is two grid points or 80 km further north and one grid point or 40 km further west than its present position. The two peaks of the south dome do not move together; they rather move apart (relative to their present positions): the northern peak is 40 km further north and that in the south drifts 80 km further to the west. At present, it is not clear which combination of processes is responsible for these results but it is clear that the motion of the summits will affect the age structure of the ice at great depth.

Table 1.

Heights above sea level at the peaks of the north and south domes for the situations indicated and described in detail in the text

Table 2.

Ice volume, ice-covered area and temperate basal areas; for further explanation, see text

3. 250 000 Modelled Years of Climate Variation

While we have computed the geometry, i.e. topography and ice temperature, through the last 250 kyears and the video film shows the entire time series, we can only report a few key results here.

Figure 2 shows two time series of the ice-covered area and the total ice volume of the Greenland ice sheet as obtained using the climate forcing of the GRIP data indicated in Figure 1. Through the Illinoian, ice volume and basal area stay generally high with an accelerated decrease as the Eem Interglacial, reaches its climatic optimum at 127 kyears BP; the minimum volume and basal area occur with some retardation of perhaps 5 kyears, The increase in ice volume and basal area is rather rapid as the climate moves into the Wisconsinan Ice Age with new relative maxima at about 110 kyears BP. The volume through this last ice age more or less continuously decreases as does the ice-covered basal area, though less, with a rapid but relatively small increase as the Holocene is approached. The two time series appear somewhat counter-intuitive but they reflect to a large extent the response to the precipitation model; they are consistent with results from independent computations by Greve (1995). Interesting in this connection are also the results compiled in Table 3, which displays the elevation of the north and south domes 250 kyears BP, at the Eem climate maximum 127 kyears BP, at the glacial maximum, 21 kyears, and today. They do reflect a mixed behaviour and are not consistently higher at the climate minimum than at present.

Fig. 2.

Ice-covered area (a) and ice volume (b) versus time.

In Figure 3b we have displayed the modelled surface at the maximum of the last Eem Interglacial 127 kyears BP when the mean atmospheric temperature was 4.73°C higher than today. Greenland loses nearly 23% of its ice-covered area; a large part of the base in the southwest is now ice-free and the margin reaches the shoreline only in isolated segments. Because of the increase in precipitation rate during warm periods, the height of the north dome is now slightly larger than 250 kyears BP when the computation is commenced. It appears, for the north dome, that changes in the snow balance are more significant than dynamical effects due to changes in ice temperature through the thermomechanical coupling. The adverted ice transport away from the summit cannot oppose the increase in accumulation. This behaviour differs from that of the south dome, where the ice loses more than 100 m of its height (see Table 3). Evidently here, the snow-balance-reducing effects (ice flux and melting) are stronger than is the increase in precipitation due to the higher temperature. Generally, the south dome seems to be more sensitive to changes in ice temperature than to changes in the accumulation/ablation pattern. This is so, because of the rather narrow topography and the change in the amount of basal area with temperate ice in the vicinity of the south dome (where the ice is at the pressure-corrected melting point and therefore enhanced sliding does occur). This implies a relatively large influence of the ice temperature on the ice flow away from the south dome, thus reducing its height. Altogether, the increasing accumulation in the north cannot equalize the large loss in the south, so that the ice sheet loses more than 20% of its initial volume. In Figure 3b, the separation of the two summits of the south dome can be located; relative to their initial position, they have moved cast by 40 km.

Table 3.

Modelled surface elevations of Greenland’s ice domes during the last 250 000 years

Figure 3c shows the ice sheet during the last glacial maximum 21 kyears BP when the driving atmospheric temperature was −9.04°C lower than it is today. The ice sheet now extends everywhere to the shoreline; ice shelves cannot form in this model, since they are inhibited by an automatic calving process that is implemented for all ice reaching the ocean. Both the north and south domes exhibit the expected behaviour; the north dome becomes shallower, whereas the south dome becomes higher. Also, the north dome moves 40 km further to the north while the south dome has not moved. More explicitly, after 229 kyears of integration, i.e. 21 kyears BP, the north dome is 120 km further north than it is today, and of the south dome, the southern summit is 40 km further north and the northern peak is 40 km further west than today.

Figure 3d displays the plan view of the computed ice sheet after 250 kyears computation corresponding to its present conditions. Both domes again show the expected changes: the north dome is still higher and the south dome lower than at the climate maximum. The summit moves 40 km to the west and south (relative to its position 21 kyears BP), i.e. the computed present position of the north dome is 40 km westward and 80 km northward of its present observed location. Similarly, the northern peak of the south dome is 40 km north and the southern peak is 40 km west of their present observed positions. We have also computed the steady-state configuration of the Greenland ice sheet to present climatic conditions (by integrating over 100 kyears under such constant driving conditions). Deviations of the surface topography of this computation from the present observed topography arc similar in magnitude to those between the present topography obtained from a 250 kyear integration with the climate driving of the GRIP data and the present observed topography. Generally, the north dome is too high and it extends too far to the north and to the west, whereas the south dome is not high enough. Furthermore, the ice volume and the ice-covered area are quite similar to those for the steady-state calculation. While modelled ice margins in the south and southwest are quite similar to those at present observed, the ice sheet extends too far to the shore in the middle, in the northeast and especially along the northern shore. However, the positions of the domes are very close in the two situations; only the southern peak of the south dome is 40 km apart in the two situations.

4. Conclusions

We describe here a limited excerpt of results for the evolution of the topography of the Greenland ice sheet through the last 250 000 years of its climate history. The complete results have been compiled in a video film that remains with the authors and contains a wealth of additional information that awaits scrutiny for its climatological implications. We have shown the evolution of the total ice volume and of the ice-covered basal area; their evolution does not necessarily go in parallel with the driving atmospheric temperature; the response of the ice sheet is a complex interplay between the precipitation pattern, i.e. snow balance and melting, and advective transport of ice as it is influenced by the thermomechanical coupling embedded in the model. For a detailed analysis of ice-core data, the climatic history and thermomechanical response of the Greenland ice sheet seems to be important.