Introduction

Beneath the 4 km thick East Antarctic ice sheet there is a complex system of over 70 lakes (Reference Siegert, Dowdeswell, Gorman and McIntyreSiegert and others, 1996) that may host subglacial ecosystems. Little is known of the flow of water into these lakes or their longevity. The strongest evidence that has been presented for mass exchange between a subglacial lake and the overlying ice sheet is the freezing of Vostok lake water to the base of the ice sheet (Reference JouzelJouzel and others, 1999; Reference Siegert, Kwok, Mayer and HubbardSiegert and others, 2000; Reference Bell, Studinger, Tikku, Clarke, Gutner and MeertensBell and others, 2002; Reference Tikku, Bell, Studinger and ClarkeTikku and others, 2004). Yet to support isolated ecosystems that might exist in subglacial lakes (e.g. Reference Karl, Bird, Bjorkman, Houlihan, Shackelford and TupasKarl and others, 1999; Reference PriscuPriscu and others, 1999) an influx of water is necessary to provide nutrients. Here we present evidence for melting at the base of the ice sheet introducing water into subglacial Lake Concordia (Fig. 1).

Fig. 1. Location of Lake Concordia. (a) East Antarctica with location of map for (b) denoted by gray shaded area. Survey map box denotes location of 1999-2000 aerogeophysical survey (Reference Studinger, Bell, Buck, Karner and BlankenshipStudinger and others, 2004). (b) Ice surface topography and identification of Lake Concordia in vicinity of Dome C. The ice surface topography, in meters above sea level, is from the European Remote-sensing Satellite-1 (ERS-1) grid (Reference Bamber and BindschadlerBamber and Bindschadler, 1997) with a 5 m contour interval. Lake Concordia is outlined with a thick black line. The thin dashed black lines are flight-lines from the 1999-2000 aerogeophysical survey; the thin white lines denote identification of the lake in the radar data from both the Italian and US surveys. Diamonds denote the extension of the Lake Concordia subglacial valley to the south. Circles are lakes identified previously (Reference Siegert, Dowdeswell, Gorman and McIntyreSiegert and others, 1996). The star denotes the location of the European Project for Ice Coring in Antarctica (EPICA) core at the ?, west of Dome C, indicates the location of the 25 cm a ¹ InSAR summit of Dome C (Reference Jouzel, Orombelli and LoriusJouzel and others, 1996). The solid black li interferogram velocity estimates (Reference Legrésy, Rignot and TabaccoLegresy and others, 2000).

We use aerogeophysical data to define the setting and describe the basal processes of Lake Concordia (74.0° S, 125.15° E), an 800 km² subglacial lake near Dome Concordia (Dome C), East Antarctica (Fig. 1). To measure ice thickness and lake extent and to identify melting, freezing and internal deformation (vertical strain-rate) processes, we use 60 MHz ice-penetrating radar data collected by the Italian Antarctic Programme and the US Support Office for Aerogeophysical Research. We also use airborne gravity data (Reference Studinger, Bell, Buck, Karner and BlankenshipStudinger and others, 2004b) to estimate the lake volume.

Prior to this analysis Vostok lake was the only large subglacial lake for which basal processes had been analyzed, (e.g. Reference JouzelJouzel and others, 1999; Reference Siegert, Kwok, Mayer and HubbardSiegert and others, 2000, 2001;Reference Souchez, Petit, Tison, Jouzel and VerbekeSouchez and others, 2000; Reference Jean-Baptiste, Petit, Lipenkov, Raynaud and BarkovJean-Baptiste and others, 2001; Reference Bell, Studinger, Tikku, Clarke, Gutner and MeertensBell and others, 2002). Lake Concordia provides valuable additional insight into subglacial lake dynamics. In this paper, we demonstrate that Lake Concordia is characterized by both melting and freezing, resulting in active mass exchange between the ice sheet and the lake. We therefore consider it an excellent candidate site for future subglacial lake exploration.

Lake Setting

Lake Concordia, the size of the Dead Sea (~800km²), is clearly defined by ten radar profiles in which the basal reflectors are very flat and bright (Fig. 2), a distinguishing characteristic of subglacial lakes (e.g. Reference Oswald and RobinOswald and Robin, 1973). The basal reflector has a slope of <0.5% in the center of the lake, with no coherent reflectors detected below the base of the ice sheet (Fig. 2c-e). The observed hydrostatic equilibrium over the center of the lake indicates that the ice is afloat (a 10:1 ratio of ice surface slope to lake surface slope), with no basal shear stress. Near the periphery of the lake there are two areas where the lake reflector changes character and we infer the water depth to be very shallow and not in hydrostatic equilibrium with potential weak basal shear stress or grounding. Over the southeastern third of the lake the basal reflector is very flat, but diffuse shallow parabolas are evident below the ice sheet (Fig. 2a and b). Along the northeastern end of the lake the basal reflector is more rugged, and often discontinuous, though no diffuse parabolas are imaged beneath the ice sheet (Fig. 2e). The lake resides in a northeast-southwest-trending subglacial valley, which extends 120km towards Dome C (Reference Tabacco, Passerini, Corbelli and GormanTabacco and others, 1998; Reference Rémy and TabaccoRemy and Tabacco, 2000). Since the northern shoreline of the lake is only loosely defined by the existing data, the observed surface area of 800 km² may be an underestimate.

Fig. 2. Ice-penetrating radar profiles over Lake Concordia. (a) X08, (b) X09, (c) X10, (d) X11 and (e) Y55. The locations of the profiles are shown in Figure 1. The extent of the lake is identified by bounding upward-pointing arrows. The ice surface topography is denoted by a thin blue line coincident with the 0 km demarcation on the ice-thickness scale; three internal layers are identified by red lines with white dots; zones of deformation imaged in the upper part of the ice sheet are indicated by green arrows. TWT indicates two-way travel time (μs), with 10 μs = 842.5m in ice.

The ice over Lake Concordia is, in general, thicker than the mean regional ice thickness of 3700 m and varies from 4125 m in the southwest to 3850 m in the northeast (Fig. 3), similar to the mapped ice thickness over Vostok lake (41503740m) (Reference StudingerStudinger and others, 2003). The ceiling of Lake Concordia is steep, sloping 275 m upward from the southwest to the northeast over a distance of 40 km (Fig. 2e). A sloped ice ceiling is an important contributor to the circulation of subglacial lakes because ice thickness controls the pressure-dependent freezing temperature at the lake surface. The resulting lateral temperature gradients lead to lateral density gradients that can drive a weak large-scale circulation (Reference Mayer and SiegertMayer and Siegert, 2000; Reference Wüest and CarmackWuest and Carmack, 2000; Reference WilliamsWilliams, 2001; Reference Mayer, Grosfeld and SiegertMayer and others, 2003).

Fig. 3. Total ice thickness over Lake Concordia. The outline of the lake is given by the heavy dashed black-and-white line. The gridded maps are interpolated from data on the flightlines, with the same identifications as in Figure 1.

Constraining the water depth of Lake Concordia is necessary to estimate the lake volume. In the southern end of the lake, diffuse reflectors beneath the bright lake reflector indicate water depths on the order of several tens of meters, based on the two-way radar travel time between the lake surface and the top of individual reflectors (Fig. 2a and b). The absence of diffuse reflectors beneath the lake reflector in the lake center (Fig. 2c and d) indicates deeper water. Gravity data can be used to estimate the water depth of subglacial lakes (e.g. Reference Studinger, Bell, Buck, Karner and BlankenshipStudinger and others, 2004a). The northeastern quadrant of Lake Concordia is characterized by a 20±3mGal free-air gravity low. A residual 11 ±3mGal low remains over the lake when the regional gravity anomaly is calculated assuming 2670 kgm^-3 for bedrock density, and 900 kg m^-3 for ice density. For this bedrock density a 200-300 m deep water cavity reproduces the gravity low over the lake, which yields a maximum lake volume of 200 ± 40 km³.

The flux of material between the East Antarctic ice sheet and Lake Concordia is controlled by the accumulation rate at the ice-sheet surface, and the melt/freeze rate at its lower boundary. There are no direct or interferometric synthetic aperture radar (InSAR) measurements of ice-sheet velocities, or measurements of accumulation rate for Lake Concordia. The closest ice-sheet velocity measurements are estimated using InSAR data and range from 0 at the summit of Dome C to 25 ± 1cm a^-1 along the northeast-southwest-trending edge of the interferogram from ~74.6° S, 123.5° E to 75.0° S, 122.0° E (Reference Legrésy, Rignot and TabaccoLegresy and others, 2000) (see Fig. 1 for location). Lake Concordia lies about 80 km northeast of the northeastern end of this edge. The regional ice-sheet surface gradient suggests that the dominant component of ice flow is east-northeast.

Basal Melting and Freezing

Over a subglacial lake where the ice sheet is afloat, the basal shear stress vanishes. Therefore the horizontal velocity is expected to be uniform through the ice sheet above the lake (e.g. Reference Mayer and SiegertMayer and Siegert, 2000), and the vertical strain rate is expected to be constant with depth. In contrast, for a grounded ice sheet frozen to the bedrock, the horizontal velocity decreases with depth and vanishes at the base, while the variation of vertical strain rate with depth is controlled by the basal shear and temperature structure. Where the ice sheet is afloat above a subglacial lake, the vertical velocity, v_z , will vary linearly with depth controlled by the basal melt/freeze rate and the surface accumulation rate. In this case, vertical velocity is given by v _z (z) = — m — {[(b — m)z]/H}, where m is the melt/freeze rate, b is the surface balance rate, H is the total ice thickness, z is the elevation with respect to the base of the ice sheet, and —(b — m)/H is the vertical strain rate (Fig. 4).

Fig. 4. Schematic vertical ice-sheet velocity profiles for grounded and floating ice. Positive vertical velocities represent upward flow. The vertical velocity of the floating ice is assumed to vary linearly with elevation such that v _z(z) = –b _z/H, where b is the ice equivalent surface balance and H is the ice thickness.

At the lake surface, the vertical component of ice velocity may be downward reflecting melting, upward reflecting freezing, or neither (Fig. 4). Melting or freezing at this interface is controlled by the balance between the influx of heat from the lake, q _L, and the outflux of heat, q _B, to the base of the ice column. In terms of water volume, the melting or freezing rate in dimensions of velocity is given by m = (q _L – q _B)/p _w L where p _w is the density of water and L is the latent heat of melting. The melting temperature is pressure- dependent (e.g. Reference PatersonPaterson, 1994, p.211-213;Reference Wagner, Saul and PrussWagner and others, 1994) so that T _m(p) = T _m(p ₀) – C _T(p – p ₀), where T _m is the melting temperature, p is pressure, p ₀ is a reference pressure and C _T is the pressure-melting coefficient. For pure water and ice phases C _T = 7.42 × 10^-8K Pa^-1, and for air- saturated water C _T = 9.8 × 10^-8 K Pa^-1. The pressure influence on melting temperature leads to freezing and melting processes associated with ascending or descending water circulation in the lake, which we shall refer to as ‘glacio- hydraulic’. The heat flux q _L results from the combined effects of geothermal flux q _G and a glaciohydraulic heat flux q _V associated with water circulation. Reference RothlisbergerRothlisberger (1972) described the physics of glaciohydraulic processes as they apply to water flow in ice-walled conduits. Subsequently, other workers have described the operation of glaciohy- draulic freezing processes beneath glaciers (e.g. Reference Alley, Lawson, Evenson, Strasserand and LarsonAlley and others, 1998) and ice shelves (e.g. Reference Jenkins and BomboschJenkins and Bombosch, 1995). Both glaciohydraulic melting and freezing are likely to be important in subglacial lakes.

For a quantitative understanding of glaciohydraulic melting or freezing, the specific internal energy of water, relative to its freezing temperature, is u = c _w(T _w – T _m), where c _w is the specific heat capacity of water and T _w is the water temperature. For this discussion, we assume that the water is isothermal so that, with reference to the melting temperature, the rate of change of specific internal energy can be written du/dt = –c _wdT _m/dt = c _w C _Tdp/dt. Because the pressure change is associated with the upward or downward flow of water (rather than a temporal change in the pressure of the water body), dp/dt = –p _wgv_wz where g is the acceleration due to gravity and v _wz is the vertical component of water velocity (positive upward). Combining the expressions for du/dt and dp/dt leads to the conclusion that ascending or descending water flow is associated with a volumetric heat source of strength Φ = —p _W gc _w C _T v _wz. Thus, relative to the melting temperature, ascending isothermal water experiences cooling, and descending isothermal water experiences warming.

Reference Wüest and CarmackWuest and Carmack (2000) recognized that water circulation along a sloping ice-water contact can produce a glaciohydraulic heat flux q _V, of positive or negative sign, that can activate melting and freezing processes at the lake surface. The spatial patterns of q _V and of the regions of melting and freezing are related to the topography of the lake surface. Regions of extensive melting and extensive freezing can coexist provided that the area-averaged flux (q _v) ≈ 0.

Melting and freezing may be estimated from aerogeo- physical ice-penetrating radar data by measuring along-flow changes in the thickness of isochronous internal layers with respect to the lake surface (e.g. Reference Siegert, Kwok, Mayer and HubbardSiegert and others, 2000; Reference Bell, Studinger, Tikku, Clarke, Gutner and MeertensBell and others, 2002). This method interprets along-flow thinning of the basal layer as melting, and along-flow thickening as freezing. This method assumes that the vertical velocity can be expressed in terms of the accumulation rate, the ice thickness and the basal melting/refreezing rate (see Fig. 4. Melting and freezing estimated from basal layer thickness changes neglect variations of the flow field over time, accumulation rate, ice thickness, longitudinal stresses, upstream effects and three-dimensional ice flow.

To estimate the basal melting and freezing over Lake Concordia utilizing basal layer thickness changes, we mapped three distinct isochronous internal layers above the lake (Fig. 2). The cumulative depth error of our vertical mapping is ±15 m. We identified and mapped the deepest resolvable layer and two layers roughly equidistantly spaced above it, ~1100, 1550 and 2150 m above the lake surface. We assume the elevation of these internal layers over the lake reflects the processes active at the base of the ice sheet and deformation within the lower ice sheet. The elevation of the internal layers over the lake is consistently 200 m lower at the western grounding zone, where the ice sheet first encounters the lake, relative to the eastern grounding zone where the ice sheet regrounds assuming an east-northeast downslope azimuth of the ice-sheet flow field (Figs 2c and d and 5a).

Fig. 5. Maps of ice thickness between different layers over Lake Concordia. The gridded maps are interpolated from data on the flight-lines, (see Fig. 2).(a) Variations in the thickness of the basal ice represented as the ice-thickness variations between internal layer 1 and the lake/ bedrock surface with respect to average basal ice thickness at the southwestern shoreline (1075 m). North-northwest and east-southeast flowlines are shown as arrows at the southern end of the lake, with basal thinning indicated by a red solid line and basal thickening by a white dashed line. The locations of grounded ice identified in the individual profiles of Figure 2 are also identified. (b) Estimated thickness variations in the basal unit from internal deformation: ice thickness between internal layers 1 and 2 scaled to thickness between layer 1 and elevation where melting ice touches the lake surface in the southeastern corner (300m elevation above the lake at the southwestern shoreline). (c) Amount of melting (red) and freezing (blue): the elevation of internal layer 1 with respect to the 1075m baseline, corrected for the estimated vertical internal deformation. The black arrows indicate the range of possible downslope ice-flow directions (north-northwest-east-southeast).

We use the thickness of the basal ice unit, bounded by the deepest internal layer (layer 1) and the lake surface, to estimate the sign and rate of basal thermal processes, assuming minimal temporal changes in ice thickness, accumulation rates and variations of the flow field. We use the thickness changes of the ice units bounded by internal layers as a measure of the vertical strain and deformation. We have examined the thickness of the units bounded by internal layers 1 and 2 and by internal layers 2 and 3.

The basal ice unit varies in thickness over Lake Concordia (Fig. 5a). We present the variations relative to the thickness of the unit along the southwestern shoreline. The basal unit thickens to the northeast, thickening by 200 m close to the northern grounding line. A region of distinct 100 m thinning of the basal unit is identified along the southeastern shoreline (Fig. 5a). Neglecting the effects of vertical strain and variations in ice flow, the changes in the thickness of the basal layer reflect melting and freezing over Lake Concordia. In this case, the basal thickening reflects accretion of ice due to freezing of lake water onto the ice sheet, and the observed basal thickening would represent a net basal accretion of 200 m ice in the northern part of the lake. The 125 m thinning of the basal unit in the southeastern corner of the lake could reflect melting of the ice sheet.

The well-defined thinning of the basal unit in the south is indicative of a region of focused melting for all possible downslope (north-northwest to east-southeast) ice-flow directions (Fig. 5a). The region of focused basal thinning about the -100 m contour of the basal unit is coincident with the disrupted basal reflectors interpreted as grounding of the ice sheet (Fig. 2a and b and 5a). We interpret the close association of the thinning of the basal layers with the disrupted basal reflectors as indicating that the basal melting is related to grounding of the ice sheet in the shallow part of the lake. The grounding may be triggering the basal melting via extreme strain heating by driving the conductive heat flux in the bottom of the ice sheet to zero (Reference Clarke, Nitsan and PatersonClarke and others, 1977). Assuming a geothermal heat flow q _G = 57 mW m^-2 as observed at the Dome C site (personal communication from C. Ritz, 2002), an isothermal basal ice-sheet condition (q _B = 0) would allow 6 mm a^-1 of melting (m = q _G/p _w L) taking p _w = 1000 k gm^-3 and L = 3.335 × 10⁵ Jkg^-1. An alternative mechanism for introducing meltwater into the lake is through glaciohydraulic processes associated with the water flow along the tilted ice ceiling of Lake Concordia, analogous to those inferred for subglacial Vostok lake (e.g. Reference Mayer and SiegertMayer and Siegert, 2000; Reference Wüest and CarmackWuest and Carmack, 2000; Reference WilliamsWilliams, 2001; Reference Mayer, Grosfeld and SiegertMayer and others, 2003). Basal thinning is also expected for more easterly ice flow in the regions west of Y55 and south of X09 (Fig. 5a) where the ice thickness is ~4100 m (Fig. 3). The ice is afloat on the western end of X09 (Fig. 2b and 5a). Both grounding-line strain heating and glaciohydraulic melting at the tilted ice ceiling are viable mechanisms for introducing water into Lake Concordia.

Influence of Vertical Strain on Melting and Freezing Estimates

The flow pattern of the ice sheet over a subglacial lake is simpler than for grounded ice, because there is no basal shear stress. With relatively small variations in ice thickness (H) and surface balance rate (b), the vertical strain rate, —(b — m)/H, is controlled by the melting/refreezing rate (m) at the ice-sheet/lake interface. With a constant vertical strain rate, the downwarping of internal layers increases linearly with depth. In our calculations we assume steady-state, two-dimensional flow, uniform ice thickness and accumulation rate, and again neglect variations over time of the flow field, accumulation rate, ice thickness, longitudinal stresses and three-dimensional ice flow.

Accounting for the observed vertical strain within the basal ice sheet can alter the volume and rate estimates of melting and freezing over the lake derived from examining thickness changes of the basal ice unit. We have used the thickness changes of the units bounded by internal layers 1 and 2 and by internal layers 2 and 3 (see Fig. 2) to estimate the vertical strain in the lower half of the ice sheet over Lake Concordia (Fig. 5b). The distributions of thickening and thinning for these two units are very similar, implying the vertical strain does not vary as a function of depth between ~1500 and 3100 m over the lake. Given the similarity of the thickness changes for these two units, we have only shown the thickness change for the unit bounded by internal layers 1 and 2 (U1-2). The thickness change of this unit is presented relative to a point along the western shoreline where the unit is 425 m thick. Little thickening is observed along the northern shoreline where the thickening of the basal unit was observed. The dominant feature of U1-2 is a region of distinctive thickening of up to 150 m along the eastern shoreline. This 20% thickening is centered in the southeast of the lake (Fig. 5b), and is roughly coincident with the area of basal thinning (Fig. 5a) and the disrupted basal reflectors (Fig. 2a and b). The maximum basal thinning is further south, between profiles X08 and X09 (Fig. 5a).

The absence of significant thickness change to U1-2 in the north lake (Fig. 5b) indicates that the interpretation of up to 200 m accretion in that region from thickening of the basal ice unit (Fig. 5a) is not significantly altered by vertical strain. We assume the downstream effects of ice-sheet regrounding on basal ice thickening over the north central lake will be minimal, as such deformation in Vostok lake appears to be fairly localized (Reference Bell, Studinger, Tikku, Clarke, Gutner and MeertensBell and others, 2002), generally <5 km from the grounding line. The interpretation of basal accretion from basal ice thickening in northern Lake Concordia appears to be robust.

The apparent thickening of U1-2 by vertical strain in the southeast will increase the estimate of net basal melting inferred from thinning of the basal ice unit bounded by internal layer 1 and the lake surface. Internal layer 1, ~1 km above the lake surface, provides the deepest observational constraint on the basal thickness variations. In the southeast the observed vertical strain, reflecting the downwarping of the internal layers, increases with depth. For example, along radar profile X09 at the southern end of the lake the downwarping of internal layer 1 near the southeastern lake margin is significantly greater than for the shallower internal layer 2 by ~85 m (Figs 2b and 6). In this simplified twodimensional case, the vertical strain is estimated as δr _a/X _a, where X _a is the thickness between layers 1 and 2 at the western shoreline and δr _a is the observed thickening between layers 1 and 2 relative to X _a (Fig. 6). Estimation of strain rate is dependent on the time required for this thickening to occur (Δt), which we do not know. We could calculate a basal melt/freeze rate (m) from the estimate for the vertical strain rate, —(b — m)/H = δr _a/X _a/Δt, assuming we knew Δt and the accumulation rate (b). Without knowing these parameters we can still estimate the effect of vertical strain in the basal ice. Assuming a constant vertical strain rate above the lake, we can calculate the elevation of an isochronous surface below layer 1, denoted layer 0, without knowing δt. An estimate for the additional component of vertical downwarping between layers 1 and 0, δr _b, assumes it is proportional to the ice thickness at the western shoreline, i.e. δr _b = X _b(δr _a/X _a) where X _b is the thickness between layers 0 and 1 at the western shoreline. This follows from the assumption of a constant vertical strain rate, i.e. (δr _b/X _b)/Δt = (δr _a/X _a)/Δt. This deep isochronous surface will show more pronounced melting than layer 1, as illustrated schematically for profile X09 (Fig. 6).

Fig. 6. Schematic of observed internal deformation and projection of estimated internal deformation in the basal ice above Lake Concordia for profile X09 (Fig. 2b). The scaling relationships between the observed and projected internal layers are explained in the text.

Our approach to estimating the effect of vertical strain in the lowermost part of the ice sheet above the lake assumes that the internal deformation in the 1000 m above the lake is proportional to the internal deformation between 1000 and 2500 m above the lake. The predicted thickening of the basal layer, assuming a constant vertical strain rate in the basal ice sheet, reaches a maximum of 150 m near the southeastern shoreline (Fig. 5b) for an isochronous surface with a basal thickness of ~-300m above the western shoreline. If the ice sheet is grounded in the southeastern lake, as interpreted from the disrupted basal reflectors (Fig. 2a and b), the vertical strain rate may not be constant, but instead could be much smaller, as seen schematically in Figure 4.

One method of estimating the effect of vertical strain on thickness changes in the basal ice is to assume a constant vertical strain rate and calculate the basal thickness of an isochronous surface below the deepest observed internal layer, as presented above. Correcting the observed changes in basal thickness of an internal layer 1075 m above the lake surface at the southwestern shoreline (Fig. 5a) for the predicted basal thickening at just ~300m above the lake surface at the southwestern shoreline and assuming a constant vertical strain rate (Fig. 5b) produces an estimate of basal thickness changes within the lowermost 450 m of the ice sheet (from the lake surface at -300 m to the thickest ice in the northeastern corner of the lake at 150 km) (Fig. 5c). The estimated deformation-corrected basal layer displays a maximum 100 m decrease in elevation along the western shoreline from south to north (Fig. 5c). The estimated deformation-corrected maximum amount of freezing is 250 m (ranging from -100 to 150 m assuming an east- northwest ice flow) in the northeastern corner of the lake, and the maximum melting is 300 m along the southeastern margin of the lake. Below we quantify the discrepancies between the estimated melting and freezing.

The distribution, rates and annual volumes of freezing and melting are derived from the slope of our estimated deformation-corrected basal thickness variation map (Fig. 5c) along the azimuth of the ice flow over the lake. Freezing is the dominant process for all downslope azimuths (north-northwest-east-southeast), with 30-50% more freezing than melting. The dominance of freezing can be interpreted as (i) an ongoing reduction in lake volume, (ii) an indication of significant transport of water into the lake in a subglacial hydrologic system that would balance the melt and freeze or (iii) an understimate of the amount of melting.

For the estimated deformation-corrected basal thickness map, the total volume of accreted ice above the lake is 0-60 km³, dependent on the ice-flow direction. This volume is ~25-30% of the estimated 200 km³ lake volume. Most of the accretion occurs over the open lake, with an average freeze-on of 125 m over the northern 400 km² of the lake for an east-northeast flow. Only 20% of the total basal thickening (~10km³) occurs over the grounded northernmost 100km² of the lake, (e.g. Fig. 2e), less than the 30-50% discrepancy between the melt and freeze. Therefore any error associated with misinterpretation of the basal thickening in the grounded zone will not change the observation that the freezing volume exceeds the melting volume. We estimate that the deformation in the lowermost 450 m of the ice sheet adds another 12-15 km³ of melt, based on the observed internal deformation (Fig. 5b). This is still insufficient to account for the total volume of net freezing. Again, only if the internal deformation in the 1000 m above the lake is significantly greater than we have assumed, and therefore adds much larger volumes of melt, or if there is a large contribution from the subglacial hydrologic system, can the excess freezing be reconciled.

To evaluate basal melting and freezing rates we assume an east-northeast downslope ice-sheet flow field with a speed of 25 cm a^-1 as estimated from the rates calculated with InSAR at a distance of ~80km from the lake in the direction of Dome C (Reference Legrésy, Rignot and TabaccoLegresy and others, 2000). The proximity of Lake Concordia to the ice divide (Fig. 1) and the associated low surface velocities will result in a significantly slower passage of the ice sheet over this lake than over Vostok lake. The ice sheet will traverse the 20 km across the northern end of Lake Concordia from west- southwest to east-northeast in close to 80 000 years based on the 25 cm a^-1 rate. In comparison the ice sheet traverses the 60 km across the southern end of Vostok lake from northwest to southeast in 20 000 years at a rate of 3 m a^-1 (Reference Bell, Studinger, Tikku, Clarke, Gutner and MeertensBell and others, 2002). The estimated melting and freezing rates at Lake Concordia display a range from -6.0 to 6.0 ± 0.9 mm a^-1, with freezing dominating the larger northern part of the lake. The melting and freezing rates are estimated as the change in surface slope of the internal deformation-corrected basal layer ice thickness (Fig. 5c) using a 25 cm a^-1 east-northeast ice-sheet flow field. The range of rate estimates is similar for all downslope ice-flow azimuths (from north-northwest to east-southeast). The error is based on the ±15 m vertical uncertainty in the identification of internal layers. The sum of the calculated melt/freeze rates over the lake yields a net freezing rate that is about 15% of the total amplitude of the melt/freeze rates. The Lake Concordia freezing rates are slower than those observed at Vostok lake of 7-38 mm a^-1 (Reference Bell, Studinger, Tikku, Clarke, Gutner and MeertensBell and others, 2002), but are consistent with geothermal constraints and rates required for the extreme strain heating of 6 mm a^-1 discussed earlier. Over Lake Concordia the observation that the volume of basal freezing may exceed the observed input of water from melting of the ice sheet raises the question of the longevity of subglacial lakes and their response to global climate change.

Climate Change and Subglacial Lake Longevity

Long-term global climate change may affect the water budget of subglacial lakes. At Vostok station, 700 km from Lake Concordia, there is a ~10-12°C surface temperature change between glacial and interglacial episodes (Reference PetitPetit and others, 2000) (Fig. 7a) extracted from the ice-core record. These surface temperature changes will propagate through the ice sheet. We modeled the basal ice temperature regime using the paleo-surface temperatures from the Vostok ice core assuming one-dimensional heat transfer. The purpose of this analysis is to isolate the influence of surface temperature changes on the basal melting rate. The redistribution of thermal energy by lake circulation and the resulting modification of the effective geothermal flux at the water-ice contact is not considered.

Fig. 7. Assumed surface forcing and calculated subglacial response based on thermal modelling of a 4100 m ice column above Lake Concordia. (a) Surface temperature variations derived from deuterium record in Vostok ice core adjusted so that the present surface temperature coincides with that at Dome C (Reference PetitPetit and others, 2000). (b) Calculated variation in basal melting/freezing rate (solid line) associated with heat flow from Lake Concordia into the base of the ice column compared with the basal melting rates associated with a geothermal heat flux of 50 mWm^_2 (dashed line). The net rate of freezing or melting resulting from these two processes is the difference between the geothermal contribution and the conductive loss. If the geothermal flux were 50 mWm^_2 the corresponding meltwater equivalent melting rate would be 4.75 mma^_1 and a cyclic variation of net melting and net freezing would result. F: freezing; M: melting; LGM: Last Glacial Maximum. The dotted lines link the glacial minimum observed in the Vostok record with the resulting basal freezing.

The heat-transfer model is based on numerical solution of the equation

where T is temperature, z is elevation above the base of the ice column, t is time, v _z is the vertical component of ice velocity, and σ_xz is the shear stress. The terms on the right- hand side of this energy-balance equation account for the effects of vertical ice advection, thermal conduction and deformational heating. The assumed physical properties of ice are p = 900 kg m^-3 (density), c = 2177.3 J kg^-1 K^-1 (specific heat capacity), K = 2.1 W m^-1 K^-1 (thermal conductivity), B ₀ = 2.755 × 10^-12 Pa^-3 s^-1 (flow-law coefficient), E = 60.7 kJ mol^-1 (creep activation energy) and n = 3 (flow-law exponent). Other physical constants are the universal gas constant R = 8.31434 J mol^-1 K^-1 and the gravity acceleration g = 9.8 ms^-2. The vertical component of ice velocity is assumed to vary linearly with ice thickness such that v _z(z) = —bzH, where b = 0.037 m a^-1 is the ice equivalent surface balance and H = 4100 m is the thickness of the ice column. The shear stress distribution is taken as σ_xz = pg(H — z) sin θ, where θ = 0.023° is the ice surface slope. The boundary conditions are taken as T(0, t) = -2.712°C (the ice melting temperature at the base of a 4100 m ice column) and T(H, t) = T _s(t), where T _s(t) is derived from the deuterium record in the Vostok ice core which spans the time interval from 422.8 kyr bp to the present (Reference PetitPetit and others, 2000). A small temperature shift is applied to adjust for the present surface temperature at Dome C relative to that at Vostok. The initial temperature distribution is assumed to vary linearly with elevation from a maximum of -2.712°C at the base of the ice column to a minimum of -61.5°C at the surface. The effects of initial transients are suppressed by starting the simulation at 800 kyr BP whereas results in Figure 7 are only presented for the past 400 kyr. From 800 to 422.8 kyr bp, when Vostok data are lacking, the surface temperature was held constant at the initial value.

The regionally averaged rate of melting or freezing at the lower boundary is controlled by the balance between the influx of geothermal heat and the conductive flux of heat into the base of the ice column as introduced previously. The conductive flux of heat into the base of the ice sheet is defined as q _B = –K∂T(0,t)/∂z. To isolate the influences of geothermal flux and conductive loss (Fig. 7b) on the melt/freeze rate m (defined above) we express them in terms of their separate contributions m _G = q _G/p _w L and m _B = q _B/p _w L where m = m _G – m _B is the net result of these opposing processes. This separation proves to be advantageous because the geothermal flux is not well constrained by observations. If, as we assume, the geothermal flux does not vary with time then the variation in the bottom melting or freezing rate is associated with variation in the conductive flux into the base of the ice column.

We found that the surface temperature cycle over the ~400000year Vostok record produces a variation in the basal melting/freezing rate of 0.1-0.2 mm a^-1 over a glacial cycle and that this cycle is delayed by 50 000-100 000 years with respect to the surface changes (Fig. 7b). We have found that the magnitude of this variation is not greatly affected by changing the form of v _z(z) or by other plausible alterations of the glaciological assumptions. Over the past 100 000 years the melting rate was a maximum about 60 000 years ago and is currently decreasing with time as the base of the ice sheet is just now sensing the surface temperatures at the Last Glacial Maximum (Fig. 7). The slow transit of the ice sheet above Lake Concordia means that the basal ice along the eastern side of the lake has experienced both glacial and interglacial conditions. The predicted changes in melting and freezing rates are below the resolution of our method. For a geothermal heat flux of ~50 mW m^-2 the variation in melting rate could yield phases of melting alternating with phases of freezing. This heat flux is plausible given the margin for error in the 57mWm^-2 geothermal heat-flux estimate from ~1000m above the bedrock at the Dome C ice core (personal communication from C. Ritz, 2002). Subglacial lakes with small volumes may therefore be ephemeral features.

Subglacial Lake Exploration

Our analysis of the ice-penetrating radar data above Lake Concordia indicates that there is melting and freezing, possibly associated with a tilted ice ceiling, consistent with the operation of glaciohydraulic processes in the lake (e.g. Reference Mayer and SiegertMayer and Siegert, 2000; Reference Wüest and CarmackWuest and Carmack, 2000; Reference WilliamsWilliams, 2001; Reference Mayer, Grosfeld and SiegertMayer and others, 2003). Together, basal processes and water circulation could drive a flux of water, biota and sediment through Lake Concordia, making the lake a good candidate for hosting life and preserving sediments. Lake Concordia may be representative of other subglacial lakes and shares many characteristics with Vostok lake, including the important role of glaciohydraulic freezing and melting processes.

Ice-penetrating radar data (Reference Tabacco, Passerini, Corbelli and GormanTabacco and others, 1998; Reference Rémy and TabaccoRemy and Tabacco, 2000) provide an understanding of the local and regional glaciological setting. The geological setting is also being addressed using aerogeophysical surveys (Reference Tabacco, Passerini, Corbelli and GormanTabacco and others, 1998; Reference Studinger, Bell, Buck, Karner and BlankenshipStudinger and others, 2004b). Lake Concordia is accessible given its proximity to the Dome C European station (e.g. Reference Jouzel, Orombelli and LoriusJouzel and others, 1996). Proximity to the Dome C ice-core site will also enable age dating of the ice above the lake. Lake Concordia therefore meets many requirements for subglacial lake exploration. The next stages of exploration are more detailed reconnaissance and mapping, such as seismic profiles to determine lake volume and sediment thickness (these are being planned; personal communication from A. Camerlenghi, 2002), and determination of the ice flow above the lake. These measurements would lay the groundwork for an initiative of in situ sampling of the ice overlying the lake, the lake water and the sediments below the lake.

Conclusion

This is the first study to document widespread melting of glacial ice feeding a subglacial lake, and to postulate that extreme strain heating and glaciohydraulic melting may be important mechanisms for introducing meltwater into Lake Concordia. We have mapped freezing, melting and vertical internal deformation of the basal ice sheet over Lake Concordia. Vertical strain and melting are associated with a grounding zone in the southeastern corner of the lake. Melting is also potentially associated with the thick- ice end of the tilted ice ceiling. The lake may be dominated by freezing and net loss of water from the lake. Lake Concordia is a promising potential site for future subglacial lake exploration.