Introduction

The deep ice cores CRIP and GISP2, which were drilled close to the highest point of the Greenland ice sheet (“Summit”) at 72°34′N 37°38′W, contain direct and indirect information on palaeoclimatic conditions such as atmospheric composition, surface temperature, snowfall rate, reaching back at least 250 000 years into the past (Reference DansgaardDansgaard and others, 1993; Reference SowersSowers and others, 1993; Reference MeeseMeese and Others, 1994). In order to obtain time series of these climatic parameters, it is necessary to provide age-depth profiles of the boreholes. This can be achieved partly by couniing annual layers downward; however, in the lower parts of ice cores, annual layers cannot be resolved (Reference JohnsenJohnsen and others, 1992; Reference MeeseMeese and others, 1994). Therefore, flow models must be used to compute age-depth profiles from the point where annual layers cease to be resolvable down to the bottom of the ice sheet.

Here, the three-dimensional transient large-scale ice-sheet model SICOPOLIS (Reference GreveGreve, 1995, Reference Greve1997) is applied to the Greenland ice sheet and a simulation is presented which covers 250 000 years of climate hisiory. It is shown that the present state of the ice sheet in the vicinity of Summit is reproduced very well and that the simulation provides reasonable age-depth profiles except for the vcry near-bottom regions.

The Model Sicopolis

SICOPOLIS (SImulation COde for POLythermal Ice Sheets) is a three-dimensional dynamic/thermodynamic ice-sheet model which simulates the time-dependent extent, thickness, velocity, temperature, water content and age for any specified grounded ice sheet as a response to external forcing. External forcing is given by (i) surface temperature, (ii) snowfall rate, (iii) sea level surrounding the iee sheet, and (iv) geothermal heat flux from below.

In this study, SICOPOLIS is applied to the Greenland ice sheet, and a palaeoclimatic simulation which covers the period from 250 ka BP (here, BP refers to thousand calendar years before present) until today, is discussed. The simulation is driven by the mean annual air-tempcrallire (T_ma) reconstruction derived from the δ¹⁸O profile of the GRIP core (Reference DansgaardDansgaard and others, 1993) and the conversion formula (Reference Johnsen, Dahl-Jensen, Dansgaard and GundestrupJohnsen and others, 1995)

with α = –211.4°C, β = –11.8°C × (per mil)^–1 and γ = –0.1925°C × (per mil)^–2 Figure1. This relation is established back to marine-isotope stage 5d, 113ka BP; however, in this study, it is used to compute a surface-temperature history for the entire 250ka covered by the δ¹⁸O data. Before 113ka BP, the so-obtained surface-temperature hisiory is smoothed by 2000 year averaging and used as a reasonable initialization of the simulation. The impact of the uncertainty level of Equation (1) on the ice sheet can be judged by comparing computed and measured borehole temperature profiles. The discussion below (section 3) shows that the discrepancy is less than 2°C throughout the boreholes GRIP This relation iand GISP2. The effect of these slight deviations on the ice fludity due to the thermomechanical coupling is negligible when compared to uncertainties in the ice rheology itself (dust content, anisotropys established back to marine-isotope stage 5d, 113ka BP; however, in this study, it is used to compute a surface-temperature history for the entire 250 ka covered by the δ¹⁸O data. Before 113ka BP, the so-obtained surface-temperature hisiory is smoothed by 2000 year averaging and used as a reasonable initialization of the simulation. The impact of the uncertainty level of Equation (1) on the ice sheet can be judged by comparing computed and measured borehole temperature profiles. The discussion below (section 3) shows that the discrepancy is less than 2°C throughout the boreholes GRIP This relation iand GISP2. The effect of these slight deviations on the ice fludity due to the thermomechanical coupling is negligible when compared to uncertainties in the ice rheology itself (dust content, anisotropy).

Fig. 1. Air-temperature forcing ΔT_ma, sea-level forcing z_s1, ice thickness H at GRIP (solid) and GISP2 (dash-dotted), basal temperatures T_b at GRIP (solid) and GISP32 (dash-dotted).

The snowfall rate S, the temporal evolution of which has a pronounced effect on computed age-depth profiles, is coupled linearly to the air-temperature deviation

with γ _s = 0.03°C^–1 corresponding to a 75% reduction of S for glacial conditions, specified by ΔT_ma = −25°C (cf. Reference Cuffey, Clow, Alley, Stuiver, Waddington and SaltusCuffey and others, 1995; Reference Cutler, Raymond, Waddington, Meese and AlleyCutler and others, 1995). Certainly, this is a simplified approach. In reality, theexistence of a unique relation between S and ΔT_ma , or, equivalently, between S and δ¹⁸O is not granted, and spatial variabilities inchanges of S are excluded completely in Equation (2). In future work, it should therefore be aimed at replacing simpleparameterizations like Equation (2) by more detailed information obtained by climate modelling (e.g. Reference Ohmure, Wild and BengtssonOhmura and others, 1996).

For the sea-level, z_s1 , which only affects the ice-sheet margin to some extent and is not crucial tor the temporal evolution of the ice sheet in central Greenland, a piecewise linear scenario is applied (Fig. 1), which represents a simple approximation of current reconstructions (e.g. SPECMAP record, Reference Imbrie, Berger, Imbric, Hays, Kukla and SaltzmanImbrie and others (1984); New Guinea record, Reference Chappell and ShackletonChappell and Shackleton (1986)). For the geothermal effects, a temperature-gradient boundary condition is imposed 5 km below the ice base to account for thermal inertia of the lithoshere and the geothermal heat flux is kept constant at

Horizontal grid spacing is 20 km. The vertical resolution is 51 grid points in the cold-ice region, 11 grid points inthe temperate-ice region (if existing) and 11 grid points in the lithosphere. Emphasis is put on a 160 × 160 km window around Summit where high-resolution radio-echo-sounding data of the ice-Surface and bedrock topography are available (Reference Hodge, Wright, Bradley, Jacobel, Skou and VaughnHodge and others, 1990) and where the GRIP and GISP2 drill sites are located.

Results

In addition to the external forcing of ΔT_ma and z_s1 Figure1 shows the simulated temporal evolution of the ice thicknesses H at GRIP and at GISP2, H_GRIP and H_GISP2 , respectively, and the temporal evolution of the basal temperatures T_b , at CRIP and GISP2, T_GRIP and T_GISP2 , respectively. The thicknesses follow the air-temperature forcing in phase, because central Greenland is virtually not affected by ablation but reacts most sensitively to changes of the snowfall rate which is larger during warm climates (see Equation (2)). The agreement between modelled and observed present ice thicknesses is very good. The predicted GRIP column is 61 m thicker than observed (3090 m instead of 3029 m) and the predicted GISP2 column is 32 m thicker than observed (3076 m instead of 3044 m ).

The basal temperatures do not change by more than 3°C during the simulation, even though the air temperatures vary by 27.7°C. This is due to the impact of the large ice thicknesses in the Summit region, which strongly dampen surface oscillations. Surprisingly. T_GRIP and T_GISP2 are low during the warm Eemian and their values at the Last Glacial Maximum 21 ka BP exceed today’s values. Cooler basal temperature during interglacial climates is expected because ice velocity is larger, and this enhances advection of cold surface ice towards the base. The counteracting effect of enhanced strain heating is not relevant close to ice domes, consequently the advection effect is dominant and warm climates entail low basal temperatures around Summit (and vice versa). Again, the agreement for the present values is very good: the predicted value T_GRIP is −8.35°C as compared to a measured value of –8.56°C for GRIP. The measured value for GISP2 is not yet published, so comparison with the −7.90°C predicted value is not yet possible.

In Figure2 the measured (Reference Hodge, Wright, Bradley, Jacobel, Skou and VaughnHodge and others, 1990) and simulated present surface topographies and iee thicknesses in the vicinity of Summit are presented. The shape of the surface topography is reproduced very well, the most conspicuous discrepancy being the region above 3200 m elevation, which is distinctly too large. Nevertheless, the difference between simulated and measured elevations nowhere exceeds 80 m and is approximately 30 m at CRIP and GISP2. The distance betweenthe simulated and actual summit position is 12.3 km. Better agreement is not possible because of the 20 km grid resolution. The agreement for the ice-thickness distribution is equally convincing. Even though the fine structure of the data (2 km resolution) cannot be reproduced by the 20 km grid, the main features (the almost closed 3200 m contour in the southwest, the 2800m contour in the east, the general decrease towards the northeast) are very well predicted.

Fig. 2. The present vicinity of Summit. (a) Measured surface topography (Reference Hodge, Wright, Bradley, Jacobel, Skou and VaughnHadge and others, 1990), (b) measured ice thickness (Reference Hodge, Wright, Bradley, Jacobel, Skou and VaughnHodge and others, 1990), (с) simulated surface topography. (d) simulated ice thickness. Solid show the GRIP position, solid circles the GISP2 position and open squares the position of the simulated summit (only in panels с and d). Surface elevations in km a.s.l. (spacing 10 m), ice thicknesses in km (spacing 200 m).

Figure3 depicts the computed present temperature and age profiles at GRIP and GISP2, and compares them with the measured temperature profiles (Reference Cuffey, Clow, Alley, Stuiver, Waddington and SaltusCuffey and others, 1995; Reference Johnsen, Dahl-Jensen, Dansgaard and GundestrupJohnsen and others, 1995) and previous datings (Reference DansgaardDansgaard and others, 1993; Reference SowersSowers and others, 1993; Reference MeeseMeese and others, 1994), respectively. The agreement for the two temperature profiles is very good, the distinct improvement compared to Reference Greve, Weis and HutterGreve and others (in press) (the slight temperature inversion was not reproduced in that study) being due to the different δ¹⁸O-to-T_ma conversion with lower glacial air temperatures applied here (see Equation (1)). The computed age-depth profiles are in almost perfect agreement with the observed datings in the upper parts of the two ice cores where the latter were determined from annual-layer counting. In contrast, below 2000 m depth, the simulation predicts ice up to 35% younger compared to the datings by Reference DansgaardDansgaard and others (1993), Reference SowersSowers and others (1993) and Reference MeeseMeese and others (1994), which is apparently a consequence of the fully three-dimensional transient modelling applied here.

Fig. 3. Stmulated temperature-depth and age-depth profiles for GRIP and GISP2 (solid), measured temperature-depth profiles (Reference Cuffey, Clow, Alley, Stuiver, Waddington and SaltusCuffey and others, 1995; Reference JohnsenJohnsen and others. 1995) (dashed), previous age-depth profiles (Reference DansgaardDansgaard and others, 1993; Reference SowersSowers and others, 1993; Reference MeeseMeese and others, 1994) (dashed).

Below 2500m depth, however, the solution scheme for the age equation makes the predicted age-depth results imprecise for the following reason. The age field A is calculated by

where d/dt implies the material time derivativewhich incorporates ice movement. The bracketed term in Equation (3) is a purely artificial vertical diffusion which is required for numerical stability. The value of the diffusivity D_A, D_A = 5 × 10⁻⁸ m² s⁻¹, is chosen to he as small as possible but large enough to ensure stable integration. Thu, an artificial boundary condition at the ice base, z = b(x, y, t), is required:

where m_age is the thinning factor, chosen as m_age = 200 and S_mean is the mean snowfall rate on the ice sheet (for further discussion see Reference Greve, Weis and HutterGreve and others (in press)). By choosing a typical ice thickness [H] = 3000 m and a typical vertical velocity [VH] = 0.1 ma⁻¹, the time-scale for advection is

and the time-scale for (artificial) diffusion is

Obviously, t _diff ≫ t _adv, so that the influence of the artificialdiffusion term on the computed age-depth profiles is negligible except for a basal boundary layer where the results depend on the artificial boundary condition (4). By conducting additional simulations with m_age = 0, Reference Greve, Weis and HutterGreve and others (in press) demonstrated that the thickness of this boundary layer is approximately 15% of the total ice thickness in the Summit region of the Greenland ice sheet, which corresponds to the lowest 500 m of the GRIP and GISP2 cores.

Conclusion

This study demonstrates that the large-scale ice-sheet model SICOPOLIS, when driven by a realistic palaeoclimatic scenario of air temperature, snowfall rate and sea level, reproduces the present state of the Summit region of the Greenland ice sheet very well. This good reproduction justifies using the model as a means of computing age-depth profiles for the deep parts of the GRIP and GISP2 ice cores where annual-layer counting is not possible. Hereby, the main shortcoming is the necessity of including numerical diffusion in the age calculation, which makes the results in a near-bottom boundary layer unreliable. Future work must Consequently aim at avoiding this by either adopting a more stable discretisationscheme for the purely advective age Equaiion (3) or by applying a direct particle-tracing algorithm.