1. Introduction

This paper presents studies of the influence of varied geothermal heat flux beneath Antarctica on the overlying ice sheet. The work is based on experiments that have been performed with Greve’s (1995) polythermal ice-sheet model which ran successful simulations for the Greenland ice sheet; its numerical code SICOPOLIS has now been adopted for the southern polar continent.

A particular property of this model is the different treatment of cold and temperate zones in the ice sheet; most of the ice volume may consist of the first type of ice which is regarded as an incompressible single-component material, while the second can be found in bottom regions as a binary mixture of ice at its local pressure-melting point and water, e.g. when the pressure is high enough to melt the ice. The energetic behaviour of cold ice is determined by its temperature field; since the temperature in the temperate zones is locally constant, here the production and advection of water content take over this role. A non-material surface called CTS (cold-temperate transition surface) separates both regions from each other; it migrates depending on the conditions supported by the ice on both sides. More detailed information about this topic and the theory of polythermal ice sheets can be read in the papers by Hutter (1982, 1983, 1993) and Greve (1995) and Greve and Hutter (1995). The numerics have been developed by Greve (1995) and a short description is given in section 3.

Interactions with the atmosphere and the lithosphere are stated in the following form: the surface temperature is parameterized as a function of geographical latitude and height; the accumulation function is represented by a temporally constant snowfield which was supplied from precipitation model output of the ECHAM3 T42-GCM belonging to the Deutsches Klimarechenzentrum Hamburg. Melting at the ice surface is not taken into account, because temperatures are so low that this phenomenon virtually does not occur. The lithosphere adjusts to the ice load following a simple model of vertical moveable rock columns sinking in the viscous asthenosphere, and the primary topography has been calculated from the present-day rock-bottom topography and ice thickness after Drewry (1983). The geothermal heat flux is kept constant in space and time; it has been varied between 42 and 105 m W m⁻² in the experiments described below. As we are only interested in steady states, thermal inertia of the lithosphere is not taken into account.

Our results demonstrate a strong dependence of the characteristic ice-sheet properties on the basal thermal regime which in turn seems to be very sensitive to variations in the geothermal heat flux, as has already been demonstrated for the Greenland ice sheet by Greve and Hutter (1995).

2. Theory and Numerics

The theoretical background of the polythermal ice model used in our study can be found in Greve (1995); we will not write down here the entire set of equations, which build the physical structure of the problem, for reasons of brevity. Briefly, it consists of balance equations for mass, momentum and energy for the region of cold ice, combined with constitutive assumptions for example the material law for viscous fluids containing Glen’s flow law and a rate factor depending on the homologous temperature. The energy balance can be evolved in an equation for the temperature field

(T temperature, Tʹ homologous temperature, v_x,y,z velocities in the x,y and z directions, ρ density of ice, c specific heat of ice, k beat conductivity of ice, E enhancement factor, A(Tʹ) rate factor for cold ice, f(σ) creep function using Glen’s flow law with n = 3 and σ shear stress). Ice-temperature changes by horizontal and vertical advection, by vertical heat conduction and by heat production based on internal friction. The shallow-ice approximation, which assumes typical height scales to be very much smaller than typical length scales, eliminates lateral heat conduction.

Temperate ice is described by similar equations to those for cold ice, yet valid for the mixture of ice and water. A balance equation for the water content ω is added. As temperature equals pressure melting

(T _m pressure-melting point, T ₀ melting point at zero pressure, β Clausius–Clapeyron gradient and h z coordinate of the free surface), the water content takes over the role of the energetic variable determining the material behaviour. Its temporal evolution is given by

(L latent beat of ice, A _t(ω) rate factor for temperate ice, f _t(ω) = f(σ) creep function for temperate ice and D(ω) water-drainage function). As in Equation (1) water content varies as a result of horizontal and vertical advection and heat production due to internal friction; furthermore, there is a small Clausius–Clapeyron correction (first two terms in the second line) and the contribution of water drainage. Note that for temperate ice the rate factor in the flow law strongly depends on the water content.

The interface separating the two ice regions (CTS) shows no jump in temperature, i.e. the temperature on both sides must be at the melting point. In the case of melting conditions, there is a volume flux of ice directed from the cold region into the temperate region, leading to a continuous water content and temperature gradient: in the case of freezing conditions (opposite ice-flow direction), however, these two quantities can be discontinuous.

The most obvious quantity representing the dynamic behaviour of the ice sheet is the ice thickness H which is simply the difference between the z-coordinates of ice surface h and ice base b; its temporal evolution is given by

with the mass flux q defined as

(

accumulation–ablation function for the ice surface and basal melting rate derived from the mass jump condition for water w), integrated from the ice base b to the ice surface h. The ice thickness grows or shrinks depending on the divergence of the mass flux, on snowfall and on basal melting, which in turn is a function of the heat-conduction gradient at the ice–lithosphere interface and of basal friction.

As for the numerics, the three-dimensional model is integrated with finite differences. The horizontal grid-point distance is 100 km, while in the vertical σ coordinates for cold ice (51 points with densification towards the bottom), temperate ice (11 equidistant points) and lithosphere (11 equidistant points) are used. The time step is 10a for the calculation of the velocities and the topography, and 100a for integration of temperature and water content. The latter equations are discretized implicitly for the z derivatives and explicitly (upwind scheme) for the x and y derivatives, i.e. the horizontal advection terms. For the ice-thickness-evolution equation, an ADI scheme (alternating implicit and explicit discretization in the horizontal directions) is applied. The external forcing has already been described in the introduction. At the ice base, the sliding law is of Weertman-type (Calov, 1994), if the base is temperate and no-sliding conditions are assumed in the case of a cold base. The associated discontinuity that appears at the transition between sliding and adhesion does not cause any problems in the numerical solution, because it is smeared out by the finite grid resolution.* The most important physical quantities are listed in Table 1.

3. Analysis of Heat-Flux Data

The “new global heat-flow compilation” derived from on-line computer devices, a description of which has been given by Pollack and others (1993), contains some heat-flux data gained by different research workers. The positions of these measurements lie quite close together near Ross Island. The data vary over a wide range (Q _geoth = 63 ... 164 m W m⁻²); an explanation may be the rough topography of this region, which extends from the Ross Sea floor to the Transantarctic Mountains. With respect to our model grid, the maximum distances between the measurements are about 150 km in x and 300 km in y directions, which gives us almost single-point information. We decided to ignore the two peak values of 164 and 142 m W m⁻² and obtain an average from the remaining seven values of Q _geoth = 70.2 m W m⁻² with a standard deviation of 6.4 m W m⁻², which is obviously much higher than the value Q _geoth = 42 m W m⁻² commonly used in ice-sheet modelling, e.g. for Greenland. As we have no more detailed information about the thermal conditions at the ice–lithosphere interface, we assume the heat flux is spatially constant beneath the Antarctic continent and take our average of 70 m W m⁻² as standard. Therefore, experiments have been performed with geothermal heat fluxes of 42 (Greenland), 50, 60, 70 (standard), 80, 90 and 105 m W m⁻²(standard + 50%)

4. Model Runs with Various Geothermal Heat Fluxes

Because of the above-mentioned hint that heat flux beneath Antarctica may be higher than the commonly assumed Q _geoth = 42 m W m⁻², we carried out a series of experiments, starting with an initial ice-free rock topography over10⁵ years, until a steady state is approximately reached. The atmospheric forcing fields remain unchanged; only variations in geothermal heat flux will be considered.

Figure 1 shows the free ice-surface topography based on measurements, taken from Drewry (1983); the topographies resulting from the model runs are depicted in Figure 2. The associated homologous temperature at the ice base and the distribution of temperate ice above the base can be seen in Figure 3. Other characteristic quantities such as total ice volume V _tot, maximum height of the free surface h _max, maximum ice thickness H _max and maximum thickness of the temperate-ice layer H _t,max are given in figure 4a–c, supplemented with observed data where possible.

Fig. 1. Observed free-surface topography after Drewry (1983). Equidistunce of the contours is 500 m.

As expected, the free-surface topography is modified significantly by the geothermal heat flux. The local height shrinks with intensified heating from below and simultaneously the total ice volume decreases by about 11.5%, but the available rock-bottom area remains ice-covered. Also, smaller features are concerned, e.g. the central Argus Dome flattens and the domes in West Antarctica also seem to vanish. The maximum ice thickness H _max is reduced by more than 600 m over the range of the varied geothermal heat fluxes; in the most probable range 42–70 m W m⁻² for Q _geoth, H _max also reacts very sensitively. Since effects of atmospheric forcing can be excluded in this case, the geothermal heat flux must induce the mechanics to build up ice sheets with different volumes and slightly varying structures. The stronger the lithosphere heats the ice, the weaker the lower ice layers become; these are responsible for greater horizontal velocities which then result in a smaller ice sheet.

The amount of ice which is directly being influenced by the geothermal heat flux, namely the volume of temperate ice V _temp, astonishingly decreases by more than 40%; at first sight, this seems to be a contradiction as the area covered with temperate ice grows. Yet, a closer look reveals that areas with an overlying temperate-ice layer do not expand that much; they only change their site and the maximum thickness H _t,max (Fig. 4c) of the temperate-ice layer (this maximum is located in the widest area in Wilkes Land) shrinks for the same reasons as the whole ice mass does. Generally, zones with a temperate-ice base can be found all along the coastline. For small values of the geothermal heat flux, most of these points possess an overlying temperate-ice layer with melting conditions at the CTS, while central Antarctica shows a cold-ice base. As heating from the lithosphere increases in strength, the ring of temperate ice becomes more uniform around the continent and the inner land also acquires a temperate base. For 105 m W m⁻², almost the entire Antarctic ice base is at the pressure-melting point.

Fig. 2. Steady-state free-surface topography for the model runs described in the text with differing genthermal heat flux (lower lefthand corner Q_geoth in m W m⁻²); all other conditions are the same. Equidistance of the contours is 500 m.

Fig. 3. Isotherms of the steady-state homologous basal temperature shown as solid lines, the homologous temperature being indicated in °C. Open diamond symbols indicate positions where the basal ice is at the pressure-inciting point but with no temperate layer above; solid diamonds (solid circles) indicate positions where there is a basal layer of temperate ice with melting (freezing) CTS. Numbers in the lower lefthand corner give Q_geoth in m W m⁻²; all other conditions are the same.

Fig. 4. Model output for (a) total ice volume v _tot, the dashed line indicating the observed value; (b) maximum ice thickness H _max (solid diamonds) and maximum height of the ice shed h _max (open diamonds) with the corresponding observed values indicated by a thin solid line for H _max and a thin dashed line for h _max; (c) maximum thickness of the temperate-ice layer H _t,max; all are displayed as a function of the geothermal heat flux Q _geoth.

In comparison with reality, which is documented in Figure 1 and as single values in figure 4a and c, the model reproduces the major features and structures of the ice-sheet topography quite well. Argus Dome and the West Antarctic maxima are located at the righthand site; however, the Antarctic Plateau becomes too steep in the middle and too flat towards the coast; West Antarctica, on the other hand, is simulated very close to reality. The best coincidence of model results and measurements with regard to properties such as total ice volume and maximum ice thickness is achieved for a geothermal heat flux of about 50–60 m W m⁻². Therefore, 42m W m⁻² seems to be too low but 70 m W m⁻² is exaggerated as a spatially constant geothermal heat flux for the entire south polar continent. Still, the maximum height of the ice sheet obtained is too high; this might be an effect of the incorrect shape, which in turn can be a consequence of the rather low resolving grid and insufficient precision of atmospheric forcing fields; such problems will be removed in the future. Another idea is to improve the lithospheric behaviour by implementing a more elaborate model, because the position of the simulated height maximum need not be the same as the actual peak.

A final encouraging remark is that most of the observed ice lakes cited by Herterich (1990) from Oswald and Robin (1973), are located in regions where our model predicts the ice temperature is at least at the melting point.