Introduction

Glaciers and ice sheets respond sensitively to climate change, but their response is difficult to interpret. In this study, we try to improve our understanding of glacierclimate interactions by means of a three-dimensional numerical model of the Greenland ice sheet. Such models were previously developed by Reference MahaffyMahaffy (1976), Reference JenssenJenssen (1977). Verbitskii (1981), Reeh (1982), Reference Oerlemans and van der VeenOerlemans and Van der Veen(1984), Reference Grigoryan, Buyanov, Krass and ShumskiyGrigoryan (1985), Hall (1987), Reference Huybrechts and Oerlemans.Huybrechts and Oerlemans (1988), Reference Fastook and ChapmanFastook and Chapman (1989), Reference Lindstrom and MacAyealLindstrom and MacAycal (1989) and Reference CalovCalov (1994). The Jenssen, Huybrechts, Lindstrom and Calov models take the coupling of the temperature and velocity fields into account, by solving the model equations on a three-dimensional network. The model we present here belongs to that group.

The Ice-Sheet Model

The model we developed is a three-dimensional ice-sheet model, including the coupling between the ice temperature and the velocity. The evolution of the ice-sheet geometry results from a mass-balance and a bedrock-adjustment model. The basic equations governing the model are similar to those of other three-dimensional models such as that of Reference HuybrechtsHuybrechts (1990). In the shallow-ice approximation, ice deformation results from shear strain only. We will present only the differences with those models.

Horizontal velocity is expressed as:

(1)

With

(2)

where S is the ice-sheet surface altitude, ρ the ice density. g the acceleration from gravity, x and y the horizontal coordinates, z the vertical coordinate, positive upwards with the origin at sea level. B _AT(z) is temperature-dependent parameter appearing in the relation between deformation rale and stress. A tuning coefficient sf is included in B _AT(z).

The equation describing the time evolution of the ice thickness is a diffusion equation with a highly non-linear diffusivity coefficient:

(3)

(4)

where H is the ice thickness, M the mass balance and F the bottom melting. The ice-sheet extent and thickness are generated by the model.

The novelty of the model is that S _a and ∫ S _adz are calculated only once and the calculation is optimized (see below). The vertical velocity, which is derived from the incompressibility condition, has also been expressed as a function of ∫ S _adz.

The temperature field must be known in order to calculate S _a and ∫ S _adz. It is derived from the heat equation, with the following conditions at the boundaries: at the ice-sheet surface, the temperature depends on climate and surface elevation (see next section), while at the ice/bedrock interface, the temperature depends on the geothermal heat flux and the type of base (lower than or at the melting point).

The bedrock adjustment is described by means of a diffusion equation similar to that used by Reference Letréguilly, Reeh and Huybrechts.Letréguilly and others (1991b).

Surface Temperature

The surface temperature is used for the mass balance as well as the upper-boundary condition for the ice temperature. It has three components:

(5)

T _a0 is the present atmospheric mean annual temperature distribution at the ice-sheet surface, parameterized from a map published by Reference OhmuraOhmura (1987). The second component takes into account the changes due to altitude variations of the ice-sheet surface, using a lapse rate. The third component is the time evolution depending on the climate. For this, we use a temperature record derived from the δ ¹⁸O measurements of the GRIP ice core (Reference DansgaardDansgaard and others, 1993) as climate forcing during 250 000 years. The record spans the last glacial interglacial cycle entirely, and most of the one before that. Doubts exist as to the validity of the record prior to 85 000BP (Reference Grootes, Stuiver, White, Johnsen and Jouzel.Grootes and others, 1993). However, even if the details of the record are uncertain prior to 85 000 BP, it is useful for the experiment to simulate the ice-sheet evolution during two full climate cycles. For our model study, we filtered some of the high-frequency variations with a spline smoothing. The original record and the smoothed one are shown in Figure 4a.

Fig. 4. Ice-sheet-evolution experiment, a. Temperature record derived from the δ¹⁸O measurements performed on the GRIP deep ice core, at Summit, central Greenland (thin line), with the spline smoothing (thick line) that was used as climatic-variations input for the model. b, c, d, e and f show the evolution of selected key variables: b. Mass balance; c. Ice volume; d. Minimum ice thickness; e. Number of points with ice; f. Basal temperature. b and f are averaged over the number of points of the ice sheet.

Mass Balance

The mass-balance distribution is composed of accumulation, ablation, calving and bottom melting. Since each of these respond differently to changes in temperature, they are parameterized separately.

Accumulation

We used the present accumulation distribution (Reference Ohmura and Reeh.Ohmura and Reeh, 1991), corrected for climatic temperature variations and changes in the ice-sheet surface slope according to the equation:

(6)

where acc(t) and acc(0) are accumulation rates at time t and the present time, respectively, and α and α ₀ are the surface slopes of the ice sheet at time t and at the present time, respectively. The exponential factor in this equation comes from a study of past accumulation rates of the GRIP ice core, central Greenland (Reference Dahl-Jensen, Johnsen, Hammer, Clausen, Jouzel. and PeltierDahl-Jensen and others, 1993). The second factor has been introduced to take into account the changes in the orographic component of the precipitation when the ice sheet is very different from the present one. Because the Greenland bedrock topography is in some places mountainous and very irregular, the 0.001 and 1.5 coefficients ensure that the factor varies between reasonable limits. During glacial conditions, for an ice sheet only slightly larger than the present one, the orographic correction term does not deviate much from 1. It becomes more important for climate warming when a large part or all of the ice sheet may disappear. The orographic correction then shifts the areas of high accumulation towards more sloping grounds.

For a climatic temperature decrease of 10°C. representing the Last Glacial Maximum, and no change in the ice-sheet surface, accumulation rates calculated from Equation (6) are 46% of the present value. For a climatic increase of 5°C, corresponding to the maximum of the previous Interglacial, accumulation rates are 148% of the present value, all other things kept constant. The amplitude of accumulation changes as a function of temperature changes is higher than that used with the Huybrechts model for the Greenland experiments (Reference Huybrechts, Letréguilly and Reeh.Huybrechts and others, 1991; Reference Letréguilly, Huybrechts and Reeh.Letréguilly and others, 1991a,Reference Letréguilly, Reeh and Huybrechts.b).

Numerical Method and Special Features of the Model

Ice temperature and velocity must be calculated through the depth of the ice. which means solving the model equations on a three-dimensional network of grid points. We did so by means of finite differences, on a network of 41 × 71 grid points in the horizontal with a spacing of 40 km. and 11 evenly spared grid points in the vertical. The bedrock and ice-sheet-surface topography are derived from the 20 km spaced data set that was used for the experiments with the Huybrechts model (Reference Huybrechts, Letréguilly and Reeh.Huybrechts and others, 1991; Reference Letréguilly, Huybrechts and Reeh.Letréguilly and others. 1991a, Reference Letréguilly, Reeh and Huybrechts.b). In order to avoid boundary problems at the bedrock and the ice-sheet surface with the three-dimensional variables, such as temperature and velocity, the vertical coordinate was scaled to the ice thickness. This also has the advantage of making it possible to greatly simplify the writing of the velocity and temperature (paper in preparation by C. Ritz).

Another particular feature of the model concerns the integrals S _a and ∫ S _adz in Equations (1) and (4). Since these are used in the calculation of ice-thickness evolution as well as the velocity field, which in turn is used for the temperature field, special attention was given to its calculation. In a finite-difference model, all variables are assumed to behave linearly between the grid points, so that the temperature can be expressed as a linear function of depth between the grid points. Instead of using an ordinary numerical integration, B _AT(z′) (S – z′)³ was expressed, by means of a Taylor’s expansion, in such a way that the integral has an analytical solution for the interval between two grid points (paper in preparation by C. Ritz). This method is more exact, and ensures that no deformation is lost in the lowermost layer of grid points because of the discontinuity at the ice/bedrock interface. Healing caused by deformation in the lowermost layer is also not completely taken into account in the basal grid points, and must be calculated separately and added to the basal heat flux. Reference HuybrechtsHuybrechts (1990, etc.) solves these problems in a manner different from ours, by using a network with unevenly spaced vertical grid points, concentrating near the base.

The time step for the ice-thickness evolution and bedrock isostasy is 5 years, and for the mass balance, temperature and velocity 50 years. A model run takes about 30 min of CPU time for a 200 000 year experiment on a Cray-C90 computer.

Steady-State Experiments

Steady-state experiments are performed by imposing a set of initial values of the climate and model parameters under study, and letting the model run 50 000 years with unchanged parameters. The first experiment was to ensure that the model could reproduce the present ice sheet when the present climatic conditions were prescribed, with the present state as initial ice-sheet configuration. Values of the model parameters are given in Table 1. Because the enhancement factor (sf) and geothermal heat flux (G _hf) are poorly known, tests were made to evaluate how uncertainties in these parameters affect the ice sheet for the present climatic conditions. Results of these tests are summarised in Table 2. None of the changes influences the extent of the ice sheet: only the shape and ice thickness are affected. The ice sheet is more sensitive in the enhancement factor, sf, than to any other parameter. Experiment 0 gives the best match with the present ice-sheet topography (see Fig. 1. and is our steady-state reference for the present climate. In all further experiments, the values sf = 3 and G _hf = 0.042 W m⁻² were used (Table 1). About 50% of the base is at the melting point, in areas away from ice divides.

Table 1. Values of the parameters used in the model

Table 2. Summary of the steady-stale experiments performed to test the sensitivity of the least-known parameters of the model, for the present climatic conditions. No. is the experiment number. Model parameters: sf is the softening-enhancement factor of the flow parameter: G_hf is the geothermal heat flux. Model outputs: np is the number of grid points covered with ice; h_max and h_mean, are the maximum and the average thickness of the ice sheet, respectively: T_b*, the mean basal temperature below the pressure melting; Vol the ice sheet volume

Fig. 1. Maps of the present Greenland ice-sheet thickness digitized an a 40 km grid. Contour lines (thin lines) represent ice-surface elevation with a contour interval of 0.5 km, starting at 1 km. The thick line indicates she ice-sheet margin; the dark gin hatching indicates ice-free land above sea level. Left, actual ice sheet: right, steady-state modelled ice sheet.

Experiments were also done for two glacial climates and two warmer climates (Table 3). Experiment 0 was used as initial configuration of the ice sheet, and the model was run for 50 000 years under the new climatic conditions. For ΔT = –10°C (Fig. 2a), representing full glacial conditions, the ice sheet reaches the coast nearly everywhere. The ice sheet is thinner than that in experiment 0, due to the lower accumulation. Correspondingly, the basal temperatures are much colder and are below the melting point almost everywhere. The initial temperature of the series used for the evolution experiment is ΔT = –4.5°C (Fig. 2b), so a steady-state run was performed for that temperature. The resulting ice sheet is not much different from the full glacial one.

Table 3. Summary of the steady-state experiments performed for different temperature perturbations ΔT. sf=3, G_hf = 0.042. No., np, h_max, h _mean, T _b*, Vol are the same as in Table 2

Fig. 2. Sensitivity of Greenland steady-state model to four different climatic parameterizations. Contours represent surface elevation.

The ΔT of 5°C (Fig. 2c) and 7°C (Fig. 2d) represent drastic climate warmings. They were chosen in preference to a more realistic warming because they are necessary for initiating a sizeable retreat of the ice sheet. The ice-margin retreat is more pronounced in the north and northeast where the accumulation is the lowest, and in the southwest where the altitude of the bedrock is relatively low. The highest mountain range in Greenland runs along the southern half of the east coast, and another, much lower range along the northern half of the west coast. Because of their higher altitude, these areas (presently under ice) are less sensitive to a warming. The ice sheet appears to be globally more stable than has been found in previous works (Reference Huybrechts, Letréguilly and Reeh.Huybrechts and others, 1991; Reference Letréguilly, Huybrechts and Reeh.Letréguilly and others, 1991a). This is because of the feedback between mass balance and ice-sheet elevation: for a climate warming; the ablation and accumulation both increase. If the ablation increase dominates over the accumulation increase, ice will melt and the surface elevation will decrease. Since temperature at the ice-sheet surface increases for a lowering altitude, the ablation will increase further, creating a strong feed-back. A collapse of the Greenland ice sheet can then he initiated for climates only moderately warmer than the present. The accumulation parameterization used in similar experiments with the Huybrechts model (Reference Letréguilly, Huybrechts and Reeh.Letréguilly and others, 1991a) led to such a situation. However, our accumulation parameterization, as given in Equation (6), produces a much higher accumulation increase for warmer climates, and for a temperature change of up to 5°C it almost compensates for the ablation increase, so that the ice sheet decreases only slightly.

Two further steady-state experiments were conducted (experiments 7 and 8), to show the sensitivity of the ice sheet to the accumulation: Equation (6) was modified so that the accumulation could not increase with the surface temperature. The resulting ice sheets for climate warmings of 3° and 5°C are shown in Figure 3. They are obviously completely different from those in experiments 5 and 6. For a temperature increase ΔT = 3°C (experiment 7), the ice sheet breaks into two smaller parts. Because of the increased ablation and margin retreat, the ice advection is no longer large enough to maintain the connection between the two domes of the ice sheet. For a temperature increase ΔT = 5°C (experiment 8), the ice sheet disappears almost completely, and only a few isolated ice caps on the highest mountains remain. Table 4 summarises the results obtained for these two experiments.

Table 4. Summary of the steady-state experiments performed for different temperature perturbations ΔT, with the modified accumulation. sf = 3, G_hf = 0.042. No., np, h _max, h _mean, T _b*, Vol are the same as in Table 2

Fig. 3. Maps of the steady-stole modelled ice sheet for two warmer climates (ΔT = 3°C and ΔT = 5°C) with the accumulation limited to the present distribution.

The accumulation parameterization for experiments 7 and 8, although not identical, is close to one used in Reference Letréguilly, Huybrechts and Reeh.Letréguilly and others (1991a). A comforting result is that the ice sheet is very close to that obtained with the Huybrechts model, in which the numerical schemes and flow parameters are different. Most models reproduce the present ice sheet under the present climatic conditions well: however, they usually differ when it comes to simulating the ice sheet under unknown climates.

Evolution Experiment

The ice-sheet configuration of experiment 4 was used as initial state, and the time-dependent temperature forcing derived from the GRIP ice core was used as climate input. The results of the evolution experiment are shown in Figure 4. The accumulation parameterization given in Equation (6) was used.

The changes in the ice-sheet volume and area are not very dramatic. Changes in area (Fig. 4e) depend mostly on ablation, while changes in the maximum ice thickness (Fig. 4d) depend mostly on accumulation. The ice volume depends on both accumulation and ablation: large changes such as those of 130 000 years and the present result from the margin retreat due to ablation increase, while smaller changes such as those that occurred during the glacial periods result from the ice-thickness changes. Margin moves are then limited because the ice sheet reaches the coast nearly everywhere. The basal temperature (Fig. 4f) seems to follow the surface temperature forcing, but with a lag of a few thousand years, and much damped: the surface temperature change from the glacial to the Interglacial is 10°C, while it is only 3°C at the base. The colder basal temperatures of the Glacial Period may be expected to create harder basal ice, leading to a thickening of the ice sheet. In our experiment, what happens may seem at first glance to be the opposite, because the maximum ice thickness depends mostly on the accumulation-rate changes. However, the fact that the general trend of the maximum ice thickness (Fig. 4d) during the last Glacial Period is relatively constant, while the general trend in surface temperature decreases during the same period, can probably be explained by the increasing stiffening of the basal ice.

The resulting ice-sheet topography for the present is 100 m thicker in the centre than that of experiment 0, and basal temperatures are 4°C lower. They are still influenced by the Ice Age conditions.

As with the steady-state experiments, the results of the evolution run differ from those of a previous study with the Huybrechts model (Reference Letréguilly, Reeh and Huybrechts.Letréguilly and others, 1991b). The temperature time series used as climate forcing differs, but the main difference comes from the accumulation parameterization. In the previous study, accumulation was not allowed to rise above the present accumulation rates for warmer climates. The two accumulation parameterizations represent extreme situations, reality probably being somewhere in between. For glacial climates, accumulation rates can be inferred from the ice cores, but for a warmer climate there is no history to guide us as to what the accumulation increase should be. It is then difficult to choose a parameterization. General circulation models may provide an answer when the particularities of the polar climates are better accounted for and their spatial resolution is sufficiently refined.

Conclusion

The ice-sheet model presented in this work is able to reproduce the present topography of the Greenland ice sheet well. Future improvements to the model will include sea-level change and heat conduction in the bedrock. In this study, limited available computing time made it difficult to take heat conduction into account, because the model would have taken much more than 50 000 a to approach steady state. Considering the large areas of the base reaching the melting temperature, a sliding law is needed. Since the ice viscosity depends on the age and formation conditions of the ice, we additionally plan to use a flow law taking ice transitions into account. To do that, the ageing of the ice will have to be computed. Finally, we plan to do some experiments with the finer. 20 km grid, in order to have a more detailed simulation of the ice sheet.

Experiments were conducted to test the model’s sensitivity to some poorly known parameters. They show that the modelled ice sheet is only slightly sensitive to a 20% variation in the geothermal heat flux, but that the ice thickness is very dependent on the tuning flow parameter, sf.

The main result of the study concerns the sensitivity of the modelled ice sheet in a warmer climate to the type of accumulation parameterization used: for a temperature increase of 5°C, one parameterization results in a margin retreat of 0–120 km (three grid points), the rest of the ice sheet being little changed, while the other results in a complete collapse of the ice sheet. Unfortunately, little evidence is available predicting what the accumulation in warmer climates should be.