The European Project for Ice Coring in Antarctica (EPICA) is a European initiative to drill two deep ice cores in East Antarctica. In the first 5 year phase (to December 2000) an ice core is being obtained from Dome Concordia (Dome C) (central East Antarctica, see Fig. 1), an area characterized by low-shear ice flow and a snow-accumulation rate of about 35 mm w.e. yr−1. It is expected to provide a climate archive covering five entire glacial-interglacial cycles (500 kyr). A further objective is to conduct field studies in the Atlantic sector of East Antarctica as a prerequisite for drilling in Dronning Maud Land (DML) in the second phase of EPICA. In contrast to Dome C, the planned DML core will be focused on rapid climate changes within the last glacial-interglacial cycle (100 kyr) in a region under the influence of South Atlantic precipitation sources and with a much larger accumulation rate.
The position of the DML drill site has not yet been established. In accordance with the aims of the coring, criteria which can be assessed by dynamic/thermodynamic ice-sheet modelling are:
(i) no melt in the entire ice column in order to avoid disturbances to the stratigraphy by meltwater;
(ii) a minimum ice age of 130 kyr in order to cover the entire recent climate cycle;
(iii) a large accumulation rate in order to resolve rapid climate changes on time-scales of decades; and
(iv) from a technical point of view, a small shear strain rate in the borehole (cf. Reference Jouzel, Hammer, Miller, Orombelli, Peel and StaufferJouzel and others, 1994).
In an earlier study (Reference Calov, Savvin, Greve, Hansen and HutterCalov and others, 1998), an initial attempt was made to determine an appropriate site for the DML core by modelling the present state of DML using a paleoclimatic simulation of the entire Antarctic ice sheet through two climate cycles, and investigating the model output with respect to coring criteria. This led to a preliminary proposal for a drill site at 73°57’ S, 03°35’W. However, due to severe limitations of the input data available at that time and the coarse grid resolution of the model (109 km) this was only seen as the first step towards an appropriate analysis for the EPICA-DML core. Here we present and discuss a strongly improved simulation based on high-resolution bedrock data in DML, a new accumulation dataset for Antarctica and, on the modelling side, a high-resolution subgrid for DML coupled to the previous coarse grid for the entire ice sheet. The drill-site proposal is adjusted according to these results, and the computed distribution of temperature, age, horizontal velocity and shear deformation is shown for the corresponding ice column.
2. Modelling of Dronning Maud Land
2.1. Large-scale ice-sheet model
The large-scale thermomechanical ice-sheet model SICOPOLIS (Simulation COde for Polythermal Ice Sheets), here used for simulations of the entire Antarctic ice sheet, is based on the continuum-mechanical and thermodynamical theory of polythermal ice sheets in the shallow-ice approximation (Reference Fowler and LarsonFowler and Larson, 1978; Reference HutterHutter, 1982, Reference Hutter1993; Reference GreveGreve, 1997b). It computes three-dimensionally the temporal evolution of ice extent, thickness, velocity, temperature, water content and age as a response to external forcing, specified by: (i) the mean annual air temperature above the ice; (ii) the surface mass balance, which is ice accumulation (snowfall) minus ablation (melting, evaporation); (iii) the mean global sea level; and (iv) the geothermal heat flux from below, imposed 5 km below the ice-bedrock interface in order to account for thermal inertia effects in the lithosphere. Isostasy is treated by a simple local-lithosphere-relaxing-asthenosphere model (Reference Le Meur and HuybrechtsLeMeur and Huybrechts, 1996), which balances the ice load locally with the buoyancy of the solid lithosphere in the liquid asthenosphere, and accounts for the asthenosphere viscosity by an isostatic time lag TV. Ice shelves are ignored, and the ice thickness is set to zero where the ice sheet reaches the coastline. Further general model features are described in greater detail by Reference GreveGreve (1997a) and Reference Greve, Weis and HutterGreve and others (1998).
The model domain for the entire Antarctic ice sheet consists of a 5341 km × 4796 km rectangle in the stereographic plane with standard parallel 71° S, spanned by Cartesian coordinates x, y (the vertical coordinate is z). Horizontal resolution is 109 km, so that the domain is discretized by 50 × 45 gridpoints. Vertical resolution is 51 gridpoints in the cold-ice column, 11 gridpoints in the temperate-ice column (if existing) and 11 gridpoints in the lithosphere column. Time-steps are 2.5 years for the dynamic evolution (ice topography, velocity) and 50 years for the thermodynamic evolution (temperature, water content, age).
2.2. Nesting for Dronning Maud Land
Reference Calov, Savvin, Greve, Hansen and HutterCalov and others (1998) used simulation results computed on the 109 km grid for assessing the present conditions of DML to determine a suitable drill site for the EPICA-DML core. As this coarse resolution only led to 10 × 9 gridpoints in the investigated part of DML, it is desirable to carry out simulations at a distinctly higher resolution. This may be achieved by increasing the resolution for the entire ice sheet at the cost of extremely long CPU times. For that reason, we pursued the alternative idea of nesting, that is the Antarctic ice sheet as a whole is still modelled on the 109 km grid, but for DML a refined subgrid is introduced (Reference SavvinSavvin, 1999). On this subgrid, the same model equations as on the coarse grid are solved for the entire ice sheet. Dirichlet-type boundary conditions for the subgrid domain are defined by interpolating the computed coarse-grid fields for surface and bedrock topography, temperature, water content and age onto the margin points of the subgrid. In turn, the subgrid returns the computed topography to the coarse grid at positions where coarse-grid and refined-subgrid points fall together (dynamic feedback). The inclusion of thermodynamic feedbacks from the refined subgrid to the coarse grid was unfavourable for numerical stability and was therefore not employed (Reference SavvinSavvin, 1999).
The nested domain chosen is a 872 km × 436 km rectangle completely covered by the high-resolution ice-thickness data discussed in section 2.3 (Fig. 1). Horizontal resolution is 10.9 km (refinement factor 10), so the region is now discretized by 81 x41 gridpoints. The time-steps for the refined DML computations are 0.1 year for the dynamic evolution and 1 year for the thermodynamic evolution.
It is clear that a prerequisite for accurate subgrid computations is the accuracy of the boundary conditions provided by the coarse grid. In this context, the 109 km coarse-grid resolution was sufficient because the simulated overall ice topography (not shown) agrees well with the Reference DrewryDrewry (1983) data. Therefore, and so as not to increase the required CPU times drastically, no increased coarse-grid resolution was applied.
The horizontal subgrid resolution of 10.9 km is less than three times the maximum ice thickness in DML. This is on the edge of applicability of the shallow-ice approximation because normal stress deviators become increasingly important on horizontal scales comparable to the ice thickness. Therefore, in a similar study for the Summit region of the Greenland ice sheet, higher-order algorithms were employed for the nested region which included the normal-stress effects (Reference Greve, Mugge, Baral, Albrecht, Savvin, Hutter, Wang and BeerGreve and others, 1999; Reference BaralBaral, 1999). The inclusion of such algorithms for the nested DML domain will be done in the future.
2.3. Bedrock topography
As a part of the DML reconnaissance within the first phase of EPICA, airborne radio-echo-sounding (RES) measurements of the ice thickness at an average resolution of 20 km were made and evaluated by the Alfred Wegener Institute (AWI) (Reference Steinhage, Nixdorf, Meyer and MillerSteinhage and others, 1999; Figs 1 and 2b). This is of particular importance because the classical maps of the Antarctic ice sheet by Reference DrewryDrewry (1983) are based on rather sparse data in DML, so their accuracy does not meet the requirements for the EPICA-DML project. With the new high-resolution ice-thickness data, a much more precise bedrock map was constructed first by interpolating them on the numerical 10.9 km grid for DML and then subtracting them from Reference DrewryDrewry’s (1983) already gridded surface topography (Fig. 2a). The isostatically relaxed bedrock topography for vanishing ice load, z = b 0(x,y), which is required for the simulations, was finally computed under the assumption that the present bedrock, z = b(x,y), is close to local isostatic equilibrium with its ice load, so that b 0(x,y) = b(x,y) + (p/pa) × H(x,y) (where H is ice thickness, p is ice density and pa is asthenosphere density). The result is shown in Figure 2c.
2.4. Snow accumulation
Giovinetto, Zwally and Bentley (personal communication, 1997; publication in preparation) have compiled a new accumulation dataset which includes approximately 1700 data points for the entire Antarctic ice sheet. This is a distinct improvement over the 550 data points on which the simulation of Reference Calov, Savvin, Greve, Hansen and HutterCalov and others (1998) was based. For the simulations here, the dataset of Giovinetto, Zwally and Bentley is interpolated onto the coarse 109 km numerical grid for the entire Antarctic ice sheet as well as onto the refined 10.9 km DML grid. Figure 3 depicts the resulting accumulation pattern in DML, which shows a clear increase from inland towards the coast.
2.5. Simulation set-up
Similar to Reference Calov, Savvin, Greve, Hansen and HutterCalov and others (1998), the present state of the Antarctic ice sheet and DML in particular was modelled as the result of a paleoclimatic simulation from 242 200 yr BP to today, initialized by a previous 100 000 year steady-state run for the conditions at 242200 yrBP. The simulation is driven by surface-temperature and sea-level histories derived from the Vostok ice core Reference JouzelJouzel and others, 1993, Reference Jouzel1996) and the SPECMAP sea-level record (Reference Imbrie, Berger, Imbrie, Hays, Kukla and SaltzmanImbrie and others, 1984), respectively, and for the present reference state of the surface temperature the parameterizations given by Reference HuybrechtsHuybrechts (1993) are used. The accumulation input is based on the modern dataset discussed in section 2.4, modified by a time-dependent factor coupled linearly to the surface-temperature deviation from today (Reference Calov, Savvin, Greve, Hansen and HutterCalov and others, 1998). Surface melting is modelled by a degree-day approach with the parameters β snow = 3 mmw.e. d–1K−1 (snowmelt), β ice = 8 mm w.e.d–1 k−1 (ice melt), P max = 60% (firn –saturation rate) and σ stat = 5°C (standard deviation for statistical air-temperature fluctuations) (Reference ReehReeh, 1991). Further model parameters are listed in Table 1.
3. Results and Implications for the DML Drill Site
The simulated present state of DML is depicted in Figures 4–8. First, comparison of the simulated ice surface and thickness (Fig. 4) with their data counterparts (Fig. 2) shows good agreement for the surface and excellent agreement for the thickness, for which many of the fine structures are even very well reproduced as a consequence of using the high-resolution input data with the refined subgrid in DML. The main test for the model performance, however, is the surface topography, the good agreement achieved being an essential prerequisite for the suitability of this study in terms of providing model support for the EPICA-DML project. A conspicuous discrepancy between the modelled and the measured ice surface occurs in the northwestern part of the region, where the 2 km surface contour shows a large bulge towards the coast which is much less pronounced in the data. Further, the divide area in the southwest is not reproduced completely. Nevertheless, these discrepancies are not so large that they might invalidate future results.
The gross direction of the ice flux (Fig. 5a), which points exactly down the surface slope due to the applied shallow-ice approximation (Reference HutterHutter, 1983), is from southeast to northwest (coastward), and the mass flux increases in the same direction. The very pronounced northward drainage basin around 0°E is sustained by the presence of the temperate ice base and, further north, even a basal temperate layer (Fig. 6). Apart from this, cold ice prevails in the region, the basal temperatures being typically a few degrees below pressure melting. The occurrence of bands of temperature contours in the plot, which run along the gridlines, indicates a numerical problem. On the other hand, the age solution (see below), which is based on the same type of advection-diffusion-production equation, does not show such features, so the bands may be a combination of a physically real large gradient and a numerical locking of this gradient to the gridlines.
At glacial-maximum conditions (17.7 kyr BP, Fig. 5b), the simulated ice flux shows smaller absolute values due to the lower ice temperatures; however the general shape of the flow field is the same as for the present. This result is important for the interpretability of any ice core in DML, because strongly varying flow directions over time may induce stratigraphic disturbances such as the strongly inclined layers and overturned folds observed in the Greenland Ice-core Project (GRIP) and Greenland Ice Sheet Project Two (GISP2) cores (Reference Alley, Gow, Johnsen, Kipfstuhl, Meese and ThorsteinssonAlley and others, 1995; Reference Greve, Mugge, Baral, Albrecht, Savvin, Hutter, Wang and BeerGreve and others, 1999). However, the simulated preservation of the flow-field shape is strongly linked to the assumed spatial constancy of the accumulation pattern through time (section 2.5) which is not guaranteed in reality so this result is rather uncertain.
The distribution of the age at 85% depth, given in Figure 7, shows a mean gradient toward the northwest, and older, deep ice appearing naturally inland. Whereas the simulated ages in the southeast are generally > 300 kyr, and thus at least three climatic cycles are covered, close to the coast they are ≤ 100 kyr, so that ice columns in this part of the region do not necessarily reach back to the late Eemian.
Together with the mass flux, the horizontal shear deformation Sh at the base, defined by
(where vx and vy are velocities in the x and y direction, respectively), increases mainly from the southeast to the northwest, with variations of two orders of magnitude (Fig. 8). However, the fine structure of the ice velocity leads to distinct local variations of the basal-shear deformation, so it can vary by an order of magnitude within some 10 km.
In agreement with the aims of the EPICA-DML core, the accumulation field (Fig. 3) makes it desirable to drill as close to the coast in the northwest as possible (criterion (iii) of the Introduction). However, this is limited by criterion (h) which requires coverage of an entire glacial-interglacial cycle. According to Figure 7, this cannot be guaranteed in the northwesternmost part of the considered region.
A reasonable compromise for this tradeoff may be the location at 73°59′ S, 00°00′ E. This position is marked in all the DML plots, and its main properties are shown in Figure 9 and Table 2. The accumulation rate is still rather large for Antarctic conditions (12.6 cm ice equiv. yr−1), and the age at 85% depth (215.3 kyr) is large enough by far to contain the last glacial-interglacial cycle. Furthermore, the basal temperature (–10.9°C) is more than 9°C below pressure melting, which leaves sufficient tolerance to fulfill criterion (i) in the whole ice column even if the local geothermal heat flux is somewhat larger than the 54.6 mW m−2 we apply here (Table 1). The shear deformation is also relatively small (0.047 yr 1 at the base, less above) compared to the maximum values within the DML region, so that criterion (iv) is equally met.
Figure 4a also shows the geographic origin of the ice in the column at 73°59’ S, OOW E. It was computed by tracing the ice particles of the 51 vertical gridpoints back along the three-dimensional velocity field until they crossed the ice surface (Mügge, 1998). To this end, the particle-tracing algorithm was run diagnostically with the SICOPOLIS output for the present time-slice only, implying that the evolution of the velocity and topography was neglected. Full coupling between SICOPOLIS and the particle-tracing algorithm with consideration of the time-dependent velocity field has not yet been implemented, this will be done in the near future. Evidently, the ice-origin line calculated in this way extends approximately 320 km upstream, so the entire 2060 m thick column was accumulated within the DML domain.
Compared to the previous drill site proposal (Reference Calov, Savvin, Greve, Hansen and HutterCalov and others, 1998; 73°57′S, 03°35′W), the site considered here (73°59′S, 00°00′E) is 109 km further east. The old position still satisfies the four criteria; however, its basal temperature is actually –5.7°C (4.8°C below pressure melting), whereas Reference Calov, Savvin, Greve, Hansen and HutterCalov and others (1998) reported –24.7°C. As the changes of the accumulation input are small, this huge shift is clearly a consequence of drastically improved knowledge of the basal topography in DML in terms of the gross shape as well as the fine structure, and, directly linked to this, the high numerical resolution which makes efficient use of the fine-structure information. For that reason, the reliability of our new results is distinctly improved.
4. Conclusion and Outlook
The present dynamic and thermodynamic state of a 872 km × 436 km region in western Dronning Maud Land (DML) was computed as the result of a paleoclimatic simulation through two glacial-interglacial cycles using the three-dimensional polythermal ice-sheet model SICOPOLIS. For DML, recent high-resolution radio-echo-sounding ice-thickness data were used and a nested model with a grid spacing of 10.9 km was applied. The computed maps of basal temperature, age and shear deformation, together with the accumulation distribution from interpolated field data, contain important information for determining a suitable drill site for the EPICA-DML ice core, the major objective of which is to obtain a high-resolution climate archive for the last glacial-interglacial cycle. Main results are:
(1) The accumulation rate increases from southeast to northwest (toward the coast), which supports a near-coastal location for the EPICA-DML core.
(2) The near-basal age decreases in the same direction, so that locations too close to the coast do not cover the entire last glacial-interglacial cycle.
(3) Except for a temperate tongue in the north and some minor scattered temperate patches, the basal ice in DML is below pressure melting.
(4) The shear deformation varies by two orders of magnitude in DML with a general coastward increase.
Based on these results, a drill-site at 73°59′ S, 00˚00′ E seems favorable. However, the optimum position for the EPICA-DML core cannot be determined by ice-sheet modelling alone; considerations of atmospheric circulation (Atlantic origin) and chemistry (interpretability of ice-core parameters) must also be included (Reference Jouzel, Hammer, Miller, Orombelli, Peel and StaufferJouzel and others, 1994). Ongoing studies of the accumulation distribution in the area and the inclusion of an improved surface-topography dataset (Reference BamberBamber, 1994; Reference Bamber and HuybrechtsBamber and Huybrechts, 1996), showing more fine structure than the classical Reference DrewryDrewry (1983) data used here, will also result in more refined boundary values. Nevertheless, the results of this study represent a sound basis for further discussion.
We thank M. B. Giovinetto, University of Calgary, Alberta, Canada, for providing the accumulation dataset of the Antarctic ice sheet, and H. Miller and D. Steinhage, Alfred Wegener Institute for Polar and Marine Research, Bremerhaven, Germany, for providing the radio-echo-sounding, ice-thickness data of DML. The comments of A.J. Payne and an anonymous referee helped improve the quality of this paper.
This work is a contribution to the "European Project for Ice Coring in Antarctica" (EPICA), a joint European Science Foundation/European Commission (EC) scientific programme, funded by the EC under the Environment and Climate Programme (1994–98) contract ENV4-CT95-0074 and by national contributions from Belgium, Denmark, France, Germany, Italy, The Netherlands, Norway, Sweden, Switzerland and the United Kingdom. Support for this work from the Deutsche Forschungsgemeinschaft under project No. Hu 412/19–3 is gratefully acknowledged.