Introduction

In Greenland, the processes limiting the size of the ice sheet are the ablation (melting of ice) and calving of ice into the sea, in roughly equal proportions. Large ablation areas can presently be found along the ice-sheet edge. Because of this, it is generally believed that the Greenland ice sheet is vulnerable to a climatic warming. Furthermore, the present ice mass raises the surface to more than 3000 m above the bedrock, which creates a climate very different from that which would have existed on bare bedrock if the ice sheet were to have been removed. On these lower elevations without ice cover, it has been suggested (Reference Oerlemans and Van der VeenOerlemans and Van der Veen, 1984) that the present climate would be too warm to allow the ice sheet to re-form. Since the Greenland ice sheet survived the transition from the last ice age to the present warmer climate, it would then be a relict ice mass remaining from a colder climate, in the sense that the ice sheet owes its existence solely to its past mass-balance history. Then, how stable is the ice sheet, and how would the ice sheet react to a warmer climate than today? Is the present Greenland ice sheet indeed a relict ice mass? If the ice sheet did not already exist, would it re-form under present climatic conditions?

The following study attempts to answer these questions by modelling the evolution of the Greenland ice sheet in present and warmer climates. We start with two different initial situations: first, with the present ice sheet, to find out under what climatic conditions it will disappear and, secondly, with bare bedrock, to investigate whether the ice sheet would reappear. For such a study, a good map of the Greenland ice sheet is essential, as well as a map of the bedrock topography. The existing maps of the Greenland ice sheet do not agree well with each other. A smoothed bedrock map, based on radio-echo-soundings performed by the Electromagnetic Institute (EMI), Technical University of Denmark, has been published by Reference Bogorodskiy, Bentley and GudmandsenBogorodsky and others (1985, p. 155). This map, however, is not sufficiently detailed for our work. In the text that follows, we will first describe the data sets we used to produce the maps; then, we will give an overview of the ice sheet and mass-balance models and, finally, we will describe the experiments carried out to test the sensitivity of the ice sheet under various climatic conditions.

The Maps

We used data of different origins for our maps. In order to minimize the number of interpolation procedures, we chose data that were already in digitized format, rather than printed maps. The elevations of the ice-free areas come from the ETOP5 World Data Base. For the ice sheet and bedrock elevations, we used data obtained by radio-echo-sounding flights undertaken by the Electromagnetic Institute (EMI) of the Technical University of Denmark. There are some problems associated with the EMI measurements because, at some spots where the flight lines cross each other, differences of up to 200 m in the surface altitude can be found. However, by comparing the EMI measurements from the southern part of the ice sheet with the surface-elevation data from Seasat and Geosat, Eckholm (unpublished) found that the differences are randomly distributed. By applying a smoothing filter, the map produced agrees reasonably well with that published by Reference OhmuraOhmura (1987). We did not use the satellite altimetry that exists south of 72°N. However, the map we obtained reproduces rather well the main features of the satellite-altimetry map of Reference Bindschadler, Zwally, Major and BrennerBindschadler and Zwally (1989).

Fig. 1. Surface-elevation contours of Greenland: the ice-free areas and the bathymetry are described by the ETOP5 data; the ice-sheet surface is obtained from smoothing the data from the EMI radio-echo-sounding flights. The triangles indicate the three deep-drilling locations: Camp Century, Summit and Dye 3.

Fig. 2. Bedrock elevations, derived from the EMI radio-echo sounding of surface elevation and ice thickness.

Furthermore, the EMI data have the advantage of having the same origin as the ice-thickness measurements (which are the only extensive source of information on ice thickness), so that the two data sets can be expected to be coherent when used together.

The radio-echo-sounding data (along flight tracks) and the ETOP5 data (every 5 min of arc) were interpolated on to a regular grid and combined, and subsequently smoothed to produce a map of Greenland (Fig. 1). A map for the bedrock was obtained by subtracting the ice-thickness data from the surface elevation (Fig. 2). The central part of the ice-sheet base is at present below sea level and is much lower than the surrounding mountains near the coast. However, since the ice thickness in the area approaches values of up to 3000 m, the ice-sheet base would rise by up to 1000 m if the ice load were to be removed. We computed the ice-free, unloaded bedrock after isostatic rebound, by attributing a density of 0.910Mg m⁻³ to the ice, and 3.370 Mg m⁻³ to the mantle (assuming that the bedrock is now in equilibrium). The result is plotted in Figure 3. This changes the mean altitude of the Greenland bedrock considerably (Table 1).

In order to evaluate the accuracy of our data, we compared their values with measurements at a few points. Stations of known elevation are easy to find on land as well as on ice. Ice-thickness measurements, however, are not so abundant. As a matter of fact, they are limited to the sites of Dye 3 and Camp Century, where deep drillings have reached the bedrock. The result of the comparison is given in Table 2.

Fig. 3. Uplifted bedrock elevations after isostatic rebound.

At Camp Century, the altitude of the surface appears to be overestimated by 107 m. Since the surface slope of this area is of the order of 100 m per 20 km, this difference corresponds to the uncertainty in the location of the station in the gridded data. 100 m also seems to be the order of magnitude of the errors in other points, also on land. Many of the ETO elevations for the coastal stations are slightly below sea level, but this is probably because the topography is averaged inside each grid cell. The ice thickness, on the other hand, is underestimated at Camp Century by 124 m, and at Dye 3 by 40 m. From this, it can be concluded that both surface altitudes and ice thicknesses are probably accurate within about 100 m.

The Ice-Sheet Model

The ice-sheet model is three-dimensional, time-dependent and calculates the fully coupled temperature and velocity fields. Ice deformation is assumed to result from shear strain and longitudinal stresses are disregarded. Calving of ice into the ocean is simulated by letting the ice disappear at the coast, i.e. the ice thickness is set to 0 at the coast edge. Other than that, there are no constraints on the ice-sheet geometry, which is freely generated by the model in response to the model inputs. Model inputs are bed topography, surface temperature, mass balance, thermal parameters and an initial state. The model also includes the response of the underlying bedrock to changing ice load. The Earth deformation model is based on a viscous asthenosphere (Reference Oerlemans and Van der VeenOerlemans and Van der Veen, 1984) and lithosphere deflection is given by considering local isostatic equilibrium. The model was originally developed for the Antarctic ice sheet, but has been adapted to simulate conditions on Greenland. A complete description of the model has been given by (Reference HuybrechtsHuybrechts 1986,in press), and Reference Huybrechts and OerlemansHuybrechts and Oerlemans (1988).

The equations are solved numerically, using finite differences, on a horizontal grid with squared cells of 20 km. With 14 layers in the vertical for the ice flow and thermodynamic calculations, this adds up to a total of nearly 200 000 grid points. We used a scaled vertical coordinate in order to avoid boundary problems in the temperature calculations. The upper layer has a thickness of 10% of the ice thickness, and the lowermost layer is 2%. This allows a much more refined description of ice deformation in the basal layers, where the shear concentrates. The diffusion equations are solved by using an Alternating-Direction-Implicit scheme. This approach has the advantage that larger time steps can be taken than the more conventional explicit integration schemes. For the fine grid that was chosen, the time step is set at 2 years, which means that only a powerful computer can handle the number of operations involved. We used the CRAY-2 computer of the University of Stuttgart. On this machine, the model needs 40 min CPU time for a 10 000 years integration.

The Mass Balance

The mass balance is the most important factor determining the state of the ice sheet. In the model, the components of the mass balance (accumulation and ablation) are calculated separately and its perturbations are parameterized in terms of temperature. Temperature is thus the principal forcing variable. Although run-off and accumulation result from quite complex processes, involving the general circulation pattern in the atmosphere and the energy balance at the ice-sheet surface, such a simplification is necessary to keep calculation times on the fine grid in use within acceptable bounds.

In particular, precipitation is a difficult process to model. It depends not only on a cyclonic activity, depression paths and moisture content, but accumulation rates over the Greenland ice sheet are also determined by such factors as temperature (since colder air can carry less precipitable moisture), continentality and orientation of the ice sheet with respect to prevailing winds (orographic effect). These factors are much too complicated to be parameterized. Alternatively, one could think of using a general circulation model. However, those models only describe the climate of the whole planet, and for a period limited to a few years only. Also, they are not sufficiently detailed for a specific area such as Greenland, and cannot handle evolution over many thousands of years. So, a more practical approach is to use observed data and to perturb the resulting accumulation distribution by a prescribed change in the temperature. We used a map by Reference Ohmura and ReehOhmura and Reeh (1991) for the mean annual precipitation rate. The accumulation rate is then prescribed to vary by changes in surface temperature DT according to

The fact that part of the precipitation falls as rain is neglected. This is however not considered a serious constraint, since only a small fraction of the annual precipitation is involved, during the summer months, and a large part of the rainfall may be expected to refreeze into superimposed ice, before it can run off to the coast.

A 5.3% change in accumulation rate for every 1 deg change in mean annual temperature is suggested by correlating annual accumulation rates measured in shallow ice cores in central Greenland with the corresponding δ¹⁸Ο values (Reference Clausen, Gundestrup, Johnsen, Bindschadler and ZwallyClausen and others, 1988), and converting δ¹⁸Ο to temperature by means by means of the factor 0.62‰ × δ¹⁸O/deg valid for present-day central Greenland conditions (Reference DansgaardDansgaard, 1961). For a 5 deg temperature shift, the corresponding accumulation rates change by around 30%. The highest degree of uncertainty in the present approach is connected with the evolution of precipitation in warmer climates and/or different ice-sheet geometries; if the ice sheet changes substantially, or even disappears, this may well affect not only the intensity but also the pattern of the accumulation distribution. However, very little is known in that respect. In order to evaluate the sensitivity of the model results to the above accumulation parameterization, a different relation was also used, where the increase in accumulation is limited to the present value:

The ablation model (Reeh, unpublished) is based on the degree-day method. As shown by Reference Braithwaite, Olesen and OerlemansBraithwaite and Olesen (1989), there is a high correlation between positive degree days and melt rates at West Greenland ice-margin locations. The annual number of positive degree days that represents a melt potential is then calculated from the July temperature and the mean annual temperature, that is parameterized in terms of altitude and latitude using data compiled by Reference OhmuraOhmura (1987). The positive degree-day model also accounts for the daily cycle and for random temperature variations from the regular, long-term annual cycle. The calculated number of positive degree days (PDD) is then used to melt, in the following order, snow, superimposed ice and ice, with degree-day factors of 0.003 m ice melt/PDD for snow, and 0.008 for ice (Reference Braithwaite and ThomsenBraithwaite and Thomsen, 1984). The first 60% of melted snow forms superimposed ice, and the rest runs off. In the event that superimposed ice remains at the end of the melting season, the warming effect due to the release of latent heat during its formation is included in the ice-sheet surface-temperature calculations.

Fig. 4. Modelled ice-sheet elevations for the present climate.

In order to test the model, as well as to have a reference run against which to compare the ice-sheet evolution for different climates, Figure 4 shows the present ice sheet simulated in steady state. The modelled ice sheet then appears to be slightly larger and thicker than the actual ice sheet, occupying 1.78 × 10⁶ km² and containing 3.21 × 10⁶ km³ of ice as compared to 1.67 × 10⁶km²and 2.83 × 10⁶ km³ for the present ice sheet. However, considering that the complete ice-sheet geometry is internally generated and that response of the ice sheet to increased accumulation rates since the beginning of the Holocene is not precisely known, the similarity with the present ice sheet is actually quite good.

The Model Experiments

The purpose of the model experiments is to evaluate the reaction of the Greenland ice sheet to a climatic warming. More specifically, we want to study (1) what would be left of the ice sheet in a warmer climate, (2) what increase of temperature would make the ice sheet disappear, and (3) conditions in which the ice sheet would re-form if it had to start from bare, ice-free ground.

In the first experiment, the steady state of the ice sheet was computed for a range of temperature rises between 1 and 8 deg relative to the present climate. The ice sheet appears sensitive to a climate warming: even a small increase of temperature will cause a decrease of the ice volume. However, on the map scale, for an increase up to 2 deg, those changes are hardly detectable. For warmings of up to 3 deg, most of the changes occur at the ice margin, particularly in the southwestern part, which retreats by up to 60 km (Fig. 5a). The surface elevations in the central part of the ice sheet, on the other hand, appear to be relatively unaffected. It takes a temperature increase of 4 deg for changes in the ice sheet to become more apparent: the ice sheet then splits up into two parts, one large part covering central and northern Greenland, and a much smaller ice cap over the southern mountains. It is interesting to note that the deep-drilling site of Dye 3 becomes ice-free. With a 5 deg temperature rise, the main ice sheet is shown to shrink even more, and northern Greenland, including the area of the Camp Century deep-drilling site, then also becomes ice-free. For a temperature increase of 6 deg, the remaining ice sheet disappears completely, leaving behind only four small ice caps on the southern and eastern mountains. It then takes a temperature increase of 2 deg more for the remaining ice masses to disintegrate completely.

Fig. 5. Modelled steady-state ice sheet for temperature increases of 3–6 deg as compared to the present. The ice-sheet boundaries are not explicitly drawn, but the steep front can be recognized where three or four elevation contours ( every 200 m) run close together.

Fig. 6. The ice sheet after 50 000 years of modelled evolution for temperature increases of 0–4 deg relative to the present. Here, the ice-free topography after isostatic rebound served as an initial configuration.

In the second experiment (Fig. 6), we investigated the conditions under which the ice sheet would re-form. These runs started from the bare bedrock (after isostatic rebound) for various temperature variations between −10 and +4 deg. Surprisingly, under present conditions, the ice sheet re-forms into an ice sheet very similar to the one simulated by the model when starting with the present ice sheet. The ice sheet will even re-form for climates warmer than today: it re-forms completely for a climate up to 2 deg warmer, and possibly also for a climate 3 deg warmer (in this case, equilibrium is still not reached after 50 000 years). The Greenland ice sheet is thus clearly not a relict ice mass left over from a previously colder climate.

Fig. 7. Dependence of ice-sheet volume on a temperature rise relative to the present. Open triangles: present ice sheet used as initial configuration. Solid circles: ice-free Greenland used as initial configuration.

The volume changes are summarized in Figure 7. As seen already in Figures 5 and 6, the temperature-ice-volume curve is not a one-to-one relationship, but presents a hysteresis. At all temperature increases between 0 and 6 deg, we have two possible steady states of the ice sheet, a situation that has been theoretically suggested by Reference Oerlemans and Van der VeenOerlemans and Van der Veen (1984), although not for the high temperature range found in our study. At a 4 deg temperature increase, for example, one of the steady states consists of an almost ice-free Greenland, while the other consists of a large central-northern ice sheet and a small southern ice cap. Which situation will occur (large ice sheet or no ice sheet) is then not only determined by climatic conditions but also by the preceding mass-balance history.

As mentioned before, the accumulation parameterization is probably the most delicate part of the mass-balance model. The present distribution on the ice sheet is based on existing measurements in Greenland, and is well described. In a landscape with little or no ice, however, it is very possible that the precipitation pattern would be different. The accumulation is parameterized only in terms of temperature changes, which means that changes in the accumulation distribution due to orographic changes are not accounted for. As long as the ice sheet exists in Greenland, these changes are not important and can be safely neglected. In situations without an ice sheet, this may introduce some errors. Such errors are likely to be found near the western coast, where the presence or absence of ice means a large difference in the surface topography. In southern Greenland and along the eastern coast are high mountains. Whether they are covered with ice or not does not affect the altitude distribution much, and changes of accumulation due to orographic changes are expected to be minor there. By neglecting accumulation changes due to orographic changes in the simulations with a bare bedrock as an initial condition, we may then be overestimating the rate of development of the ice sheet. However, the effect on the final steady state is probably minor.

Fig. 8. Steady-state ice sheet for a temperature increase of 3 deg above the present, with the accumulation-temperature relationship given by Equation (2).

Changes in the general circulation pattern may also influence the precipitation pattern. To include those effects, however, is beyond the scope of this work. However, some experiments with a different precipitation-temperature relationship (constant precipitation in warmer climates; relation 2) were performed, in order to estimate the sensitivity of the accumulation variation on the results (Fig. 8). It appears that the resulting ice-sheet configurations are only moderately affected: in the warming experiment, the splitting of the ice sheet now occurs for a temperature rise of 3 deg instead of 4 deg. In the experiment with initially bare bedrock, the ice sheet will re-form for temperatures up to 1 deg warmer than at present instead of 2 deg. Both of these accumulation-temperature relations probably represent extreme situations, and the different results give an idea of the uncertainty caused by the lack of knowledge of how the accumulation rate will change in a warmer climate. It is certainly comforting that the results appear to be only moderately influenced.

The isostatic adjustment rate of the bedrock due to changing ice load is on the order of 20 000 years for Greenland. A different value will influence the rate at which equilibrium is reached in the simulation. However, when equilibrium is reached, the deflection of the bedrock is independent of this value.

In the experiments, we tried to run the model ice sheet to a steady-state equilibrium, in order to reach a reference state for each temperature. To do that, we computed the evolution of the ice sheet during 50 000 years for each temperature change. For the warming experiments (starting from the present ice sheet) shown in Figure 7, most of the ice-sheet volume variations occur during the first 10 000 years, and equilibrium is approximately reached after 30 000 years. However, for the initially bare bedrock experiments (see Fig. 9), 50 000 years was not always sufficient to reach equilibrium: for temperatures of 2 and 3 deg higher than present, the equilibrium is not reached even after 50 000 years. The building of an ice sheet is clearly a much slower process than its destruction. What mostly influences the evolution of the ice sheet is the mass balance, which depends on temperature. It is limited upwards, because the accumulation decreases with decreasing temperature, while it is not limited downwards, because for a climatic warming, the ablation increases faster than the accumulation. There is then a limit to the rate of formation while, in principle, there is no limit to the rate of destruction.

The rate of formation of the ice sheet also depends on the initial extent of the accumulation area on the ice-free bedrock, and hence on the topography. For DT = 0 and −10deg the ice sheet builds up much faster than for DT = 2 and 3 deg. In the latter cases, the accumulation area at the start is restricted to the high elevations of the eastern mountain tops, so that the initial ice caps are very small. They can grow only by advecting ice formed at higher elevations, and this in turn creates an altitude-mass-balance feed-back by raising the overall elevation and thus by enlarging the accumulation area. This process is slow and may well last more than 100 000 years. For DT = 0, the larger part of Greenland is in the accumulation zone right from the start, not only the eastern mountains but also the lower western mountains. The total volume of ice that can accumulate during the first 1000 years is much larger for DT = 0 than it is for a higher DT. As for DΤ = −10 deg, all of Greenland has a positive mass balance. This illustrates the effect of the topography. For climatic conditions similar to today or colder, a model with a lower resolution that does not reproduce the details of the mountains would still be able to initiate an ice sheet on the ice-free bedrock. However, for climatic conditions 2 or 3 deg warmer than at present, the mountains are responsible for the initiation of the ice sheet. Although the lower reaches of Greenland may not have a climate suitable for the development of an ice sheet, the mountains act like a source of ice that can cover most of Greenland with ice, given enough time. With less detailed topography, it would be impossible to model such an effect.

Fig. 9. Growth of the ice sheet for various climatic conditions, starting with an ice-free Greenland as an initial condition.

Modelling equilibrium states of an ice sheet are useful to assess its sensitivity to climate change. However, they are unlikely to be encountered often in the glacial history. A temperature history derived from an ice-margin oxygen-18 record from central West Greenland (Reference Reeh, Oerter, Letréguilly, Miller and HubbertenReeh and others, in press) indicates that the temperature in the last interglacial was about 3–4 deg above the present value for a duration of about 10 000 years. The resulting ice sheet would then have looked like something in between those shown in Figure 5a and b, and would not have been far from the equilibrium state. Furthermore, Figure 10 shows that a temperature increase of at least 6 deg during a period of 20 000 years is required for the ice sheet to disappear completely. This is a much higher temperature increase as well as a much longer time period than the Eemian interglacial. This suggests that the Greenland ice sheet must have survived the Eemian interglacial. However, large areas may have been ice-free, and it is questionable whether the sites of Camp Century in northwest Greenland and particularly Dye 3 in southern Greenland were still ice-covered during that period, since they are located in areas of large margin fluctuations. This could explain the difficulties in interpreting the basal part of those ice cores. Central Greenland, on the other hand, seems to be the most favourable location for retrieving pre-Eemian ice, either by future deep drilling in the Summit area (GISP and GRIP), or by ice-margin studies (Reference Reeh, Thomsen and ClausenReeh and others, 1987).

We also conducted a short-term experiment to evaluate the response of the ice sheet to the anticipated greenhouse warming (Reference Huybrechts, Letréguilly and ReehHuybrechts and others, in press). Forcing the model with a temperature scenario rising exponentially to 8 deg by the year AD 2100, the ice volume would decrease by 68 500 km³. This is only a small fraction of the total ice volume in Greenland, however such wastage would be clearly observed; it corresponds to a mean ice-margin retreat of 3.2 km, and a worldwide sea-level rise of 17 cm.

Conclusions

The ice sheet is only moderately vulnerable to a climate warming: it would take a temperature increase of 6 deg, sustained over 20 000 years, or a temperature increase of 8 deg, sustained over 5 000 years, to allow the ice sheet to disappear completely. This increase of 6 deg is higher than the 3ߝ4 deg for the previous interglacial (Eemian). The model results then suggest that the ice sheet probably survived this period. This means that the Greenland ice sheet still contains pre-Eemian ice, which is available for palaeoclimatic studies. Furthermore, under present-day climatic conditions, if the ice sheet did not already exist, it would re-form on the ice-free Greenland topography. This indicates that the ice sheet is not a mere relict from a previously colder climate, but an active ice sheet that sustains itself. Although a climatic warming would immediately result in increased melting of the Greenland ice sheet and retreat of the ice margin, it seems that, from the more general point of view of the existence of the ice sheet, it is somewhat resistant to a climatic warming.

Fig. 10. Ice-volume evolution of the Greenland ice sheet for warmer than present conditions, starting with the present ice sheet as an initial configuration.

Acknowledgements

This paper is published with the approval of the Director, Alfred-Wegener-Institut. P. Huybrechts is supported by the Belgian National Fund for Scientific Research (NFWO) and is sponsored in part by the Belgian Office of Science Policy under contract ANTAR/04. This study would not have been possible without the excellent computing facilities (CRAY 2, University of Stuttgart) which were made available to us at the Alfred-Wegener-Institut. We thank A. Ohmura (ΕΤΗ Zürich, Switzerland) for communicating data that were unpublished at the time of the study, and the Electromagnetic Institute, Technical University of Denmark, Copenhagen, for providing the radio-echo-sounding data of the Greenland ice sheet. The Geodetic Section of the Alfred-Wegener-Institut assisted with the ETOP5 data. Useful comments from C. Lorius on an earlier draft of this paper were also much appreciated.