Introduction

Subglacial lakes have been identified at 77 locations within Antarctica from airborne radio-echo-sounding (RES) data covering over 50% of the ice sheet (Reference Oswald, and Robin.Oswald and Robin, 1973; Reference Robin,, Drewry, and Meldrum,Robin and others, 1977; Reference McIntyre,McIntyre, 1983; Reference Siegert,, Dowdeswell,, Gorman, and McIntyre,Siegert and others, in press). Seismic data indicate that, close to Vostok station, the water depth of a large 230 km by 50 km subglacial lake is about 500 m (Reference Kapitsa,, Ridley,, Robin,, Siegert, and Zotikov,Kapitsa and others, 1996). The existence of large volumes of water beneath the East Antarctic ice sheet has significant consequences for analyzing the thermal regime of the ice sheet. The presence of subglacial lakes indicates that, in the region of the ice sheet over and adjacent to the lake, the temperature of the ice-sheet base is at the pressure-melting point.

Basal ice-sheet temperatures are controlled by a number of parameters including ice thickness, ice-sheet surface temperature and accumulation rate, heat transported through horizontal advection of ice. the basal heat gradient (the sum of geothermal heal flux and heat produced from basal sliding) and heat derived from internal ice deformation. The basal heat gradient is dependent on the flow at (or effectively close to) the ice-sheet base, and its magnitude will therefore increase with distance from the ice divide as ice-sheet velocity increases. This paper aims to determine numerical information on the spatial distribution and magnitude of heat derived at the ice-sheet base due to (i) geothermal sources, (ii) internal ice deformation and basal sliding (or, collectively, horizontal ice-motion terms) and (iii) heat transported by the horizontal advection of ice.

An inspection of the spatial distribution of lakes over the Antarctic continent (Fig. 1) indicates that approximately 43% of those found are located within 100 km of an ice divide in the direction of ice flow (Reference Siegert,, Dowdeswell,, Gorman, and McIntyre,Siegert and others, in press). Furthermore, a large concentration of Antarctic subglacial lakes (∼78%) lies beneath major ice domes at Dome C, Ridge B and Titan Dome (Fig. 1). However, there are seven subglacial lakes (representing 9% of the known total), at Terre Adélie and west of the Transantarctic Mountains, that are situated over 400 km along the line of ice flow from the ice divide (Fig. 1). The direction of ice motion was determined using information on ice-sheet surface morphology (Reference Drewry,Drewry, 1983), assuming that the direction of ice flow is perpendicular to ice-surface contours.

Fig. 1.

(a) Locations of subglacial lakes determined from the SPRI–NSF–TUD airborne RES data (after Siegerl and others, in press). Ice divides are illustrated as dashed lines. Abbreviations to place-names are as follows: DC, Dome C; GVL, George V Land; HD, Hercules Dome; OL, Oates Land, RB, Ridge B; S. Sovetskaya station; SP, South Pole: TA, Terre Adélie; TAD, Talos Dome; TID, Titan Dome; V, Vostok station; WM, Whilmore Mountains. (b) Calculated basal temperatures from a three-dimensional thermomechanical ice-sheet model of the Antarctic ice sheet (from Reference Huybrechts,Huybrechts, 1992). The geothermal heat flux used to determine the basal temperatures was set at 54.6 m W m⁻².

Recent three-dimensional thermomechanical glaciological modelling of the Antarctic ice sheet indicates that, under contemporary environmental conditions, and using a geothermal heat flux of 54.6 m W m⁻², several regions of the ice-sheet base, corresponding to the proximity of the calculated ice divides, are at the pressure-melting point (Reference Huybrechts,Huybrechts, 1992). However, although the glaciological model used by Huybrechts is one of the most sophisticated of its kind, the temperature in several regions of the modelled ice-sheet base, where subglacial lakes are known to exist, was calculated to be below the pressure-melting value (Fig. 1b). For example, around Ridge B, where up to six subglacial lakes have been identified, Huybrechts’ glaciological model calculated the temperature of the ice-sheet base at around −5°C (Fig. 1b).

A general source of error in ice-sheet models is derived from the smoothed bedrock topography that is used as model input. As a consequence, the ice thickness that is calculated by glaciological models, when the present ice-surface elevation is determined, will be prone to similar errors. In the case of Reference Huybrechts,Huybrechts’ (1992) model, the Antarctic bedrock-topography input was in the form of a finite-difference grid composed of 1600 km² sized (40 km by 40 km) cells, adapted from maps compiled by Reference Drewry,Drewry (1983). However, due to the smoothing that occurs when a continuous topography is represented as a series of individual cells, there are discrepancies between the ice-thickness array, representing the finite-difference grid, and direct point measurements from both Drewry’s maps and information gathered from the raw RES data.

The ice-sheet thickness is an important parameter in the calculation of temperature at the ice-sheet base. It is therefore useful to calculate the thermal regime of the ice sheet above each known subglacial lake, using the RES-measured ice thickness as input to the thermal equations. Subsequently, the accuracy of our numerical results will not suffer as a consequence of bedrock smoothing.

Basal Ice-sheet Temperature Calculation

We assume that, since a large proportion of subglacial lakes are known to occur at or near an ice divide (Reference Siegert,, Dowdeswell,, Gorman, and McIntyre,Siegert and others, in press), a simple Robin-type steady-state thermal model is applicable (Reference Robin,Robin, 1955). Consequently, we do not account for temperature variations at the base of the Antarctic ice sheet caused by climate change during, for example, the early Holocene (Reference Whillans,Whillans, 1978; Reference Jenssen,, Campbell, and Robin,Jenssen and Campbell, 1983). The equation used to calculate the basal temperature is:

where

and

Symbols used in the above equations are defined as follows: T _B is the basal temperature of the ice sheet (°C); T _s is the mean annual surface temperature of the ice sheet (°C); z is the coordinate in the vertical direction, positive upwards and zero at the ice-sheet base: h is the ice thickness above the subglacial lake (m); b is the mean annual surface accumulation of the ice sheet above the lake (m year⁻¹; k is the thermal diffusivity of ice (36.3 m² year ⁻¹); K is the thermal conductivity of ice (2.1 W m ⁻¹ °С⁻¹; and ∧_geo is the Earth’s geothermal heat flux (54 m W m⁻² unless stated otherwise).

Equation (4) represents the basal heat gradient caused by the flux of geothermal heat. This equation is relevant only at the ice divide, where horizontal ice motion is minimal. However, a way from the ice divide, extra heat terms associated with horizontal ice advection, internal ice deformation and basal sliding will result in more heat at the base of the ice sheet than that calculated in Equation (4). By adjusting the basal heat gradient to account for the extra heat terms at the ice-sheet base, we calculated an “effective” basal heal gradient, which we can then substitute in Equation (1).

The pressure-melting temperature (T _pmp) of pure ice, in °C, is given by:

which corresponds to a Clausius–Clapeyron gradient of 8.7 × 10⁻¹ °Cm⁻¹ (Reference Paterson,Paterson, 1994).

Equation (1) holds only if the basal temperature is at or below the pressure-melting point. If the theoretical basal temperature determined from Equation (1) is higher than the pressure-melting point, the basal temperature recorded is equivalent to the pressure-melting temperature.

The simple steady-state thermal model used in this study is one of several that may be used to determine the temperature at the base of an ice sheet. However, it should be noted that these models often provide conflicting temperature values at the ice-sheet base (e.g. Reference Paterson,Paterson, 1994). We consider, therefore, that the relative variation in basal temperatures determined by Equation (1), and their spatial distribution beneath the ice sheet, are as important to this study as the absolute basal temperatures calculated for ice above each lake. We note also that, since the Antarctic ice sheet is not in complete steady slate, our assumption of steady state may cause basal temperature calculations to be too high. Subsequently, due to the simple method by which we calculate the effective basal heat gradient and basal temperatures, we acknowledge that the results presented are, at best, semi-quantitative in nature.

Model Inputs

The thermal model requires inputs of mean annual surface accumulation (b), and surface temperature (T _S), basal temperature (T _B), basal heat gradient derived from the flux of geothermal heat (∧_geo), and ice thickness (h). The ice thickness above a subglacial lake can be measured directly from the time-dependent raw RES data (time-continuous “Z-scope” or spot-measurement “A-seope”), assuming a 60 MHz radio-wave velocity in ice of 1.69 × 10⁸ ms⁻¹ (Reference Bogorodsky,, Bentley, and Gudmandsen,Bogorodskiy and others, 1985). We have calculated the ice thickness above each lake from our own independent measurements of the raw Z-scope data, held at the Scott Polar Research Institute, University of Cambridge, which are accurate to within 1.5%.

We assume that the surface air temperature approximates the surface temperature of the ice sheet since, in the Antarctic ice sheet, the surface air temperature is generally within 1°C of the firn temperature at 10 m depth (Reference Loewe,Loewe, 1970). Thus, the surface temperature of the ice sheet can be obtained from maps of mean annual air temperature, determined from field measurements (Reference Robin,Robin, 1983; Fig. 2a).

Fig. 2.

(а) Mean annual ice-sheet surface temperature (°C), based on surface measurements of Antarctica. Adapted from Reference Robin,Robin (1983). (b) Mean annual surface mass-balance distribution on Antarctica. Values are in water equivalent cm year⁻¹. Adapted from Reference GiovinettO, and Bentley,Giovinetto and Bentley (1985).

Surface accumulation on the ice sheet above each lake was obtained from a map of accumulation, interpolated from direct field measurements (Reference GiovinettO, and Bentley,Giovinetto and Bentley, 1985; Fig. 2b). The reliability of the contours, interpolated between the field data points, has been enhanced in recent years, after a detailed numerical modelling study of the Antarctic climate yielded results similar to those measured in the field (Reference Furtuin, and OerlemansFortuin and Oerlemans, 1990).

The geothermal heat flux used in this study, 54 m W m⁻² unless stated otherwise, represents the approximate mean value of the Earth’s surface heat flux, and is similar to that used by Huybrects (1992) in recent glaciological modelling studies of the ice sheet.

However, the geothermal heat flux varies globally between about 40 and 90 m W m⁻² (e.g. Reference Stacey,Stacey, 1977). It is likely, therefore, that the actual value of the geothermal heat flux may vary spatially beneath the Antarctic ice sheet. In order to account for (i) the uncertainty in the value of the geothermal input to the thermal equation, (ii) heat generated at the ice-sheet base due to horizontal ice-motion terms (internal ice deformation and basal sliding), and (iii) heat transported through horizontal advection of ice, the basal heat-gradient parameter within Equation (1) is varied (thus producing an effective basal heat gradient), with other variables kept at constant values. Subsequent results thus provide a sensitivity experiment on the relation between the ice-sheet basal temperature above subglacial lakes and the effective basal heat gradient.

The calculation procedure adopted in this study is detailed below. Essentially, Equation (1) is solved for (1) T _B, in order to derive the basal temperature for when the geothermal heat flux is 54 m W m⁻², and (2) the basal heat gradient (Equation (4)), assuming that Т _B is at the pressure-melting value. Specifically, if the temperature above a subglacial lake was initially calculated to be below the pressure-melting value (when 54 m W m⁻² is used as the geothermal parameter) then extra heat, derived from geothermal and/or other effective basal heat-gradient terms, is required in Equation (1). We account for this extra heat by altering the geothermal heat input in Equation (4). to a minimum value that yields pressure melting at ice ice-sheet base in Equation (1). In performing this procedure, we calculate values for the effective basal heat gradient in terms of the equivalent geothermal heat input (i.e., the adjusted geothermal heat can be thought of as a surrogate for the effective basal heat gradient). This new basal heat gradient, particular for each subglacial lake, is then recorded.

Heat Transfer within the Antarctic Ice Sheet

Thermal model results indicate that, using a Robin-type approach (Equation (1)), the ice-sheet base above 70% of all known subglacial lakes is at the pressure-melting temperature when a geothermal heat flux of 54 m W m⁻² is used as the sole component of the ice-sheet basal temperature gradient (Fig. 3). In addition, the temperature of the ice-sheet base above 75% of subglacial lakes was calculated within 1°С of the pressure-melting point (Fig. 3). The majority of these lakes are located around the Dome C (75° S, 125° E), Ridge B (77° S, 92° E) and Vostok station (78° S, 105° E) regions of East Antarctica (Fig. 4). At Titan Dome (80° S, 140 E) and Talos Dome (73° S, 157° E), where several subglacial lakes exist, the temperature of the ice-sheet base above the lakes is calculated at, or within 1°С of, the pressure-melting temperature. However, in contrast to these relatively warm basal ice-sheet temperatures, we calculate the basal temperature above lakes located near Hercules Dome and the Whitmore Mountains to be greater than 5°C below the pressure-melting value (Fig. 4). Indeed, basal temperatures above lakes located, for the most part, away from ice divides in Terre Adélie (68° S, 137° E), George V Land (70° S, 140° E) and Oates Land (71° S, 151° Ε), are calculated at 11.2–20.4° C below the pressure-melting temperature (Fig. 4).

Fig. 3.

Histogram indicating the theoretical basal temperature distribution above Antarctic subglacial lakes, using Equation (1) to derive T_B for each lake whilst assuming a geothermal heat flux of 54 m W m⁻².

Fig. 4.

Distribution if basal ice temperatures above Antarctic subglacial lakes, calculated by solving Equation (1) for T_B, given Λ_geo = 54 m W m⁻². The positions of three flowlines (solid lines with arrows) along which several subglacial lakes exist are shown (after Reference Drewry,Drewry, 1983). Information concerning the ice-sheet thermal regime along flowlines (i) and (ii) is referred to later. Abbreviations to place-names are as in Figure 1.

By varying Λ _geo until pressure melting is attained for each lake, the effective basal heat gradient required to ensure basal melting was calculated. A surrogate geothermal heat flux was thus determined which causes, in Equation (4), a basal heat gradient equivalent to that derived from the actual geothermal-heal and horizontal ice-motion terms. From this exercise, we can determine the percentage of subglacial lakes that experience melting above them, as a function of the imposed geothermal heat flux (Λ _geo) (Fig. 5).

Fig. 5.

Basal temperature gradient against the number of lakes experiencing pressure melting above them. Tbe graph indicates how varying the geothermal heat flux affects the basal temperature of the ice sheet above subglacial lakes.

Comparison with Glaciological Modelling Information

Glaciological modelling can be used to calculate the temperature of the ice-sheet base by a procedure different to that employed in this study. It is useful, therefore, to compare our results with those determined from previous modelling investigations of the Antarctic ice sheet. We examine the results of two independent glaciological modelling studies, namely the thermomechanic studies of Reference Budd, and JenssenBudd and Jenssen (1989) and Reference Huybrechts,Huybrechts (1992). In both cases, the basal temperature of the ice sheet was calculated by a thermal model which accounted for heat derived from ice deformation and basal sliding. What remained largely unknown in these modelling investigations was the value for the geothermal heat flux from the Antarctic Plate. However, both studies investigated the ice-sheet response to variations within the spatially uniform geothermal heat-llux input.

Reference Budd, and JenssenBudd and Jenssen (1989) found that when a geothermal heat flux of up to 51.7 m W m⁻² was used, the centre of the ice sheet at Dome C, Ridge B and the South Pole experienced basal temperatures below the pressure-melting value. However, when the geothermal heat flux was 57.5 m W m⁻², basal temperatures within the region representing Dome C were at the pressure-melting value. Similarly, although Reference Huybrechts,Huybrechts (1992) generated basal melting temperatures around Dome C when the geothermal heat flux was 54.7 m W m⁻², his model did not predict the pressure-melting temperature beneath Ridge B. In our study, we predict that basal melting over subglacial lakes in the area of Ridge B should occur even when the geothermal heat flux is reduced to around 43 m W m⁻² (Fig. 6). The discrepancy between our results and those of the two glaciological models may be a consequence of bedrock smoothing that occurs within numerical modelling of ice sheets. We predict that the geothermal heat flux varies spatially beneath the Antarctic ice sheet between 37 and 64 m W m⁻². However, to date, glaciological models have not accounted for spatial variability in geothermal heat.

Both glaciological models predicted basal melting temperatures around 80° S, 150° E for values of the geothermal heat flux between around ±10% of 54 m W m⁻². This location is over 400 km from the nearest ice divide, and corresponds with the downstream end of the two ice flowlines analyzed in our study (Figs 1 and 4). The glaciological models predict pressure-melting values because of the relatively high basal temperature gradient in this region, which is caused by heat derived from basal sliding of the ice sheet. Upstream from this location, the basal temperature decreases to below the pressure-melting value (e.g. Fig. 1b). because of the absence of basal sliding and the heat derived from this process. Our analysis suggests that, along the two flowlines, the rapid variation within the basal heat flux may correspond with the onset of basal sliding (Fig. 7). We consider there to be good agreement between our independent analysis of the basal thermal properties along these two flowlines within the East Antarctic ice sheet and glaciological model results from this same region.

Summary and Conclusions

A one-dimensional vertical heat-transfer equation (Reference Robin,Robin, 1955) was used to analyze the thermal regime of the Antarctic ice sheet above known subglacial lakes. Three inputs to this thermal equation (Equation (1)) were provided by existing Antarctic field measurements (i.e. ice thickness (h), surface accumulation (b) and surface temperature (T _s)). A fourth variable, the basal ice temperature above subglacial lakes (T _B), was assumed to be equal to the pressure-melting value. Equation (1) was subsequently solved for Λ _geo, the geothermal heat input. In addition, assuming a value for the geothermal heat-flux input, Equation (1) was solved for T _B. Results from the numerical experiments indicated the following points about the thermal regime of the Antarctic ice sheet above subglacial lakes.

• 1. The temperature at the base of the Antarctic ice sheet was modelled to be at pressure-melting over about 70% of the subglacial lakes (which are located in the regions of Dome C, Ridge B and Vostok station) when using a geothermal heat flux of 54 m W m⁻².

• 2. The geothermal heat flux was calculated to determine the minimum basal heat required to produce basal melting of the ice sheet above all lakes (Equation (4)). The minimum heat flux was 40 m W m⁻² when eight lakes out of the whole data set experienced basal melting above them (Fig. 5). Conversely, the ice-sheet base over only one lake (in George V Land) remained below the pressure-melting value when the minimum geothermal flux of heat was 90 m W m⁻² (Fig. 5).

• 3. The amount of additional basal heat required to melt the ice-sheet base over a number of lakes was found to be related strongly to the geographical location of the lakes (Fig. 6).

• 4. The additional heat can be derived either geothermally (requiring a substantial spatial variability in the geothermal heat flux within Antarctica) or from horizontal basal sliding and internal ice deformation (or through a combination of the two).

• 5. Near (<100 km) to ice divides, where a high concentration (∼43%) of lakes exists, since the horizontal advection of ice should be relatively small, the geothermal heat flux calculated to ensure basal melting of the ice sheet was 37–64 m W m⁻². However, away from the ice divide at Terre Adélie, where horizontal advection may be relatively high, the effective basal heat flux required to ensure pressure melting at the ice-sheet base was 1.5–2 times the value needed at Dome C. Subsequently, heat derived from ice flow at Terre Adélie may be of the same order as the heat flux derived from the Earth.

• 6. Plotting the effective basal heat flux (required to calculate basal melting above sub-ice lakes) against the distance of the lake to the nearest ice divide, provides information on the variation in heat derived from horizontal ice-motion terms along the flowline. Two such flowlines were investigated, between the ice margin east of the Transantarctic Mountains and (a) Dome C ice divide and (b) Ridge R ice divide. In both cases, we parameterize empirically the variation of heat due to horizontal ice motion along the flowline as an exponential-type relation. The location along the flowline at which a significant rise in the effective basal heat flux occurs may correspond to the onset of basal sliding (Fig. 7).