Drainage networks, lakes and water fluxes beneath the Antarctic ice sheet

Abstract

Antarctica Bedmap2 datasets are used to calculate subglacial hydraulic potential and the area, depth and volume of hydraulic potential sinks. There are over 32 000 contiguous sinks, which can be thought of as predicted lakes. Patterns of subglacial melt are modelled with a balanced ice flux flow model, and water fluxes are cumulated along predicted flow pathways to quantify steady-state fluxes from the main basin outlets and from known subglacial lakes. The total flux from the continent is ~21 km³ a⁻¹. Byrd Glacier has the greatest basin flux of ~2.7 km³ a⁻¹. Fluxes from subglacial lakes range from ~1 × 10⁻⁴ to ~1.5 km³ a⁻¹. Lake turnover times are calculated from their volumes and fluxes, and have median values of ~100 a for known ‘active’ lakes and ~500 a for other lakes. Recurrence intervals of a 0.25 km³ flood range from ~2 months to ~2000 a (median ≈130 a) for known ‘active’ lakes and from ~2 to ~2400 a (median ≈ 360 a) for other lakes. Thus, several lakes that have recently been observed to fill and drain may not do so again for many centuries; and several lakes that have not, so far, been observed to fill and drain have the potential to do so, even at annual to decadal timescales.

INTRODUCTION

The magnitude and variability of subglacial water fluxes beneath glaciers, ice caps and the margins of the Greenland ice sheet affect their dynamics, by controlling water storage, subglacial water pressures and the morphology of subglacial drainage pathways (Kamb and others, 1985; Iken and Bindschadler, 1986; Björnsson, 1998; Mair and others, 2003; Bartholomew and others, 2010; Schoof, 2010; Banwell and others, 2013). Beneath the Antarctic ice sheet, however, the characteristics of drainage pathways, their spatial and temporal patterns of water flux and the effects of water movement on ice-sheet dynamics are far less well known.

The low subglacial hydraulic gradients beneath the Antarctic ice sheet mean that water ponding is prevalent. Subglacial lakes and their interconnecting drainage pathways form a dynamic hydrological system beneath the ice sheet (Gray and others, 2005; Wingham and others, 2006; Fricker and others, 2007; Pattyn, 2008; Carter and others, 2009; Peters and others, 2009; Smith and others, 2009; Siegert and others, 2014). This system may affect the ice-sheet dynamics at a large scale, with several major subglacial lakes marking the onset zones of fast ice flow (Siegert and Bamber, 2000; Bell and others, 2007; Wright and Siegert, 2012). Centennial scale changes in water routing may cause long-term variations in ice stream velocities (Alley and others, 1994; Joughin and others, 2002; Vaughan and others, 2008). Decadal to annual scale fluctuations in water movement between lakes causes medium to short term variations in vertical deformation (Gray and others, 2005; Wingham and others, 2006; Fricker and others, 2007; Fricker and Scambos, 2009; Smith and others, 2009; McMillan and others, 2013; Siegfried and others, 2014), and there is currently one study showing that water movement at these timescales also affects horizontal ice velocities (Stearns and others, 2008).

The drainage system beneath the Antarctic ice sheet and the associated mechanisms of water flow are difficult to observe due to the large ice thicknesses involved. Here we report on the characteristics of Antarctica's subglacial drainage system that can be inferred from steady-state hydrological and ice dynamics modelling. We build on methodologies used by, amongst others, Wright and others (2008), Le Brocq and others (2009), Pattyn (2010) and Livingstone and others (2013). We use the latest 1 km gridded surface and bed data to map the subglacial hydraulic potential (φ) field, the location of major drainage catchments, and the locations of the main drainage pathways within each catchment. We compare the locations of ‘sinks’ in the φ field with the positions of known subglacial lakes to investigate the extent to which features of the φ surface may be used to predict subglacial lake locations. We then couple basal melt rate calculations, derived from a spatially varying estimate of the geothermal heat flux and a balanced ice flux flow model to the drainage catchments and pathways inferred from the φ field, to provide first order calculations of subglacial water fluxes flowing beneath the ice sheet. These are used to investigate the implications for sub-ice shelf processes, the residence times of water in the lakes, and the recurrence intervals (RIs) of floods between subglacial lakes.

METHODS

Subglacial hydraulic potential and sinks

DEMs of the surface and bed of Antarctica at 1 km horizontal resolution were acquired from the Bedmap2 dataset (Fretwell and others, 2013). The datasets use the GL04C geoid as an absolute reference for elevation and the Polar Stereographic projection based on the WGS84 ellipsoid, with true scale at 71°S (Fretwell and others, 2013). The two DEMs were clipped using the accompanying grounded ice and exposed bedrock extent masks respectively. The clipped DEMs were used to calculate the φ across the continent from: φ = ρ _i gz _s + (ρ _w – ρ _i) gz _b, where ρ _i = 910 kg m⁻³ is the density of ice, ρ _w = 1000 kg m⁻³ is the density of water, g = 9.8 m s⁻² is the acceleration due to gravity, z _s is the ice surface elevation, and z _b is the bed elevation (Shreve, 1972). Sinks in the φ surface were identified and filled, and the number, area, depth and volume of all the contiguous filled sinks were calculated using the algorithm of Arnold (2010). These contiguous filled sinks can be thought of as predicted lakes.

Subglacial melt rates

Basal melt rates across the continent are derived using a three-dimensional (3-D) balanced ice flux flow model, which is described more fully elsewhere (Leysinger Vieli and others, 2007, 2011; Hindmarsh and others, 2009). The model calculates spatial patterns of basal melt across the continent by combining the geothermal heat flux with the frictional heat produced by basal shear strain, and accounting for horizontal heat advection. The approach assumes that all viscous heating occurs at the ice-sheet base, rather than englacially. Spatial patterns in the basal shear strain rate are calculated by multiplying the balance velocity by the driving stress. The model uses the shallow ice approximation to derive the internal deformation and velocity field. Here, the model is applied at a 10 km resolution using ice-sheet surface and ice thickness data from Bedmap1 (Lythe and others, 2001).

The other key boundary conditions are the spatially varying geothermal heat flux values derived from satellite magnetic data (Fox Maule and others, 2005) and the surface accumulation rate pattern based on field measurements (Arthern and others, 2006). Sensitivity experiments examining how well the model reproduces time-varying accumulation rate patterns inferred from radar layers, show that model outputs are relatively insensitive to grid size and the geothermal heat flux, but more sensitive to surface accumulation patterns (Leysinger Vieli and others, 2011). The ‘Results and Discussion’ section below provides a qualitative assessment of the errors in our basal melt flux calculations by comparison with other similar studies. The calculated basal melt rates are re-gridded to 1 km to match the φ grid.

Drainage catchments, pathways and water flux

We calculate the location and direction of the flow pathways across the sink-filled φ surface, and accumulate the steady-state water flux along the pathways, where each cell is assigned the sum of its own basal melt flux and that of all upstream cells with a higher φ (Tarboton and others, 1991). The accumulated fluxes can be expressed either as depth per time or, if multiplied by the area of the grid cell, as volume per time. The resulting grid is also used to define the individual subglacial drainage catchments beneath the ice sheet. The calculations are used to quantify the water discharge at primary drainage outlets along the coast and at known subglacial lake locations.

Lake turnover times and flood RIs

An examination of calculated water fluxes near known subglacial lakes provides some insight into water residence timescales within the lakes, or the relationship of potential RI versus outburst flux if the lakes drain episodically. To focus our analysis on the most potentially dynamic lakes, we concentrate on those that lie within 10 km of a predicted lake that accumulates at least 0.1 m a⁻¹ of melt. There are 232 such lakes, 60% of the total. For each of these lakes, we calculate the fluxes at the predicted lake outflow grid cell.

The volumes of the 232 lakes are estimated by cumulating the volumes of all the filled sinks in the φ surface that lie within a 20 km radius of the known lake locations. These are likely to be underestimates of the true lake volumes since the radar data used to derive much of the Bedmap2 bed DEM detect the water surface and not the underlying bed surface. An estimate of the steady-state lake turnover times for these 232 lakes is derived by dividing the estimated lake volumes by the calculated steady-state water fluxes.

The mean flood volume from the ‘active’ lakes observed by ICESat was ~0.25 km³ (Smith and others, 2009). Assuming this represents a typical flood volume, we calculate the flood RIs for each of the 232 lakes defined above by dividing this flood volume by the calculated steady-state water fluxes.

RESULTS AND DISCUSSION

Subglacial hydraulic potential and sinks

The pattern of φ across Antarctica (including sinks) is relatively smooth and ranges between 39 and –0.4 MPa (Fig. 1). Given Bedmap2 surface and bed DEM errors (Fretwell and others, 2013), uncertainty in the φ values are calculated to range between ±0.27 MPa for most of the ice sheet (where ice thickness uncertainty is ±30 m) and up to ±1.34 MPa for the mountainous regions (where ice thickness uncertainty is ±150 m). These errors have implications for the calculation of subglacial drainage pathways and catchments as discussed further below. Values of φ peak within the interior of East Antarctica where ice is thickest and decline towards the coasts. There are many (N = 407 115) individual sink cells forming 32 380 contiguous sinks (predicted lakes) ranging in area from 1 to 14 700 km², ranging in maximum depth from 0.004 to 870 m and ranging in volume from 3.67 × 10⁻⁶ to 4684 km³ (Fig. 2). The lower values for area and for volume are fixed by the 1 km² DEM resolution. Many predicted lakes are big enough to appear visibly at the continent-wide scale shown in Fig. 1. Predicted lakes occur predominantly in regions where bed elevations are relatively low (filled sink cell elevations: mean, $\bar x = - 223.5\,{\rm m}$ , SD, σ = 701.3 m; continent-wide bed elevation $\bar x = 92.8\,{\rm m}$ , σ = 677.3 m). The total area of all predicted lakes is 520 007 km² (representing 4.3% of the grounded ice sheet) and their total volume amounts to 7323 km³, a sea level equivalent of 0.021 m assuming a conversion factor of 2.8 × 10⁻⁶ m km⁻³ (Vaughan, 2006). The total area is greater than that found by Wright and others (2008) (1.7% of the grounded ice sheet), who used older versions of the surface and bed DEMs to those used here, which were gridded to 5 km.

Fig. 1. Subglacial hydraulic potential (Pa) of the Antarctic ice sheet derived from 1 km² surface and bed elevation grids contained in the Bedmap2 dataset. Also shown (in blue) are the locations of the filled sinks in the subglacial hydraulic potential field (the predicted lakes) and the positions of the 379 known lakes (yellow triangles) (Wright and Siegert, 2012). Background image is the MODIS Mosaic of Antarctica (Haran and others, 2014).

Fig. 2. Frequency histograms of (a) area, (b) maximum depth and (c) volume for the 32 380 contiguous filled sinks (predicted lakes) shown in Figure 1. Note the logarithmic scale on the x-axes.

Subglacial melt rates

The modelled basal melt rate distribution is shown in Fig. 3. Melt rates reach a maximum of 430 mm a⁻¹, with the greatest values occurring beneath the thick, fast-moving ice streams. Zones of negative melt (i.e. freezing) occur primarily near the coasts where ice is thin but also occasionally inland where ice is slow moving. Melt rates reach a minimum (i.e. freezing rates reach a maximum) of −185 mm a⁻¹ where ice from East Antarctica abuts against the Transantarctic mountain range. The mean melt rate under the ice sheet is 2.0 mm a⁻¹. This is less than the 3.5 mm a⁻¹ computed in a similar study by Llubes and others (2006) and less than the range of 3.2–6.5 mm a⁻¹ (with a mean of 5.3 mm a⁻¹) calculated by Pattyn (2010). Llubes and others (2006) used a simpler ice flow model to that used here, and a simple geothermal heat flux pattern based on orogen distributions. Pattyn (2010) used a more complex ice flow model, incorporating an ice stream/ice shelf component in ice streams and over large subglacial lakes, and used a range of geothermal heat flux patterns (including that derived by Fox Maule, 2005) and surface accumulation rate distributions (but not that produced by Arthern and others, 2006) in a suite of model runs. Comparing our results to those of others, shows that melt rate magnitudes and patterns are sensitive to the precise flow model and boundary conditions used. The implication of this comparison is that we are unlikely to be overestimating basal melt rates but that we may be underestimating them, although more sensitivity testing would be required to test this further.

Fig. 3. Basal melt rates (m a⁻¹) derived from the 3-D balance flux ice flow model. Also shown are the locations of the 379 known lakes (yellow triangles) (Wright and Siegert, 2012) and the 12 595 hydrological catchments (thin grey lines). Background image is the MODIS Mosaic of Antarctica (Haran and others, 2014).

Drainage catchments and networks

There are 12 595 individual basins ranging in size from 1 to 1.18 × 10⁶ km² draining water to the ocean. These are depicted in Fig. 3. The 54 largest basins by area are shown with an ID in Fig. 4 and their areas are given in Table 1. The five largest basins, in size order, are those feeding Byrd Glacier, Lambert Glacier, Recovery Ice Stream, Foundation Ice Stream and Budd Glacier. Networks of drainage pathways within all basins trend towards the coast in the general direction of the ice surface gradient (Fig. 5). There is a close association between the predicted lakes and the locations of the main drainage pathways (defined as where accumulated melt flux >0.1 m a⁻¹) (Fig. 5). For example, Fig. 6 shows the relationship between the major drainage pathways and predicted lakes for the Recovery Lakes A–D (Fig. 6a) and Lake Vostok in East Antarctica (Fig. 6b). Pathways often begin in or reach a confluence at a predicted lake because the reduction in φ associated with sinks is large enough to concentrate flow from upstream. In cases where several drainage pathways enter a predicted lake, only one outflow pathway occurs, and this creates a confluence.

Fig. 4. The 54 largest drainage basins for the main discharge points along the Antarctic Coast. Catchment IDs are listed in Table 1, which provides the hydrological characteristics of each basin, and are based on those in Rignot and others (2008). The black areas contain the remaining 12 526 smaller basins shown in Fig. 3. Background image is the MODIS Mosaic of Antarctica (Haran and others, 2014).

Fig. 5. Map of water flux (m a⁻¹) in the main drainage pathways generated from accumulating the basal melt rate field (Fig. 3) along the subglacial hydraulic potential field (Fig. 1). Only those pathways with an accumulated melt flux >0.1 m a⁻¹ are shown. Also shown (in magenta) are the filled sinks (predicted lakes) in the subglacial hydraulic potential field. Background image is the MODIS Mosaic of Antarctica (Haran and others, 2014).

Fig. 6. Map of the main drainage pathways (as in Fig. 5) (light blue) and predicted lakes (as in Figs 1 and 5) (mid blue), showing details of the areas around: (a) Recovery Lakes A–D; and (b) Lake Vostok. Only those pathways with an accumulated melt flux >0.1 m a⁻¹ are shown. The positions of known lakes (Wright and Siegert, 2012) are also shown (yellow triangles) as are the outlines for the Recovery Lakes A–D and Lake Vostok (green). Background image is the MODIS Mosaic of Antarctica (Haran and others, 2014). Inset maps show the location of the main maps within Antarctica.

Table 1. Hydrological characteristics of the 54 largest drainage basins

Accumulated melt flux represents the total melt added up for the entire basin. Discharge represents the accumulated melt flux multiplied by the cell size. Discharge per unit area represents the discharge divided by the total basin area, and gives a measurement of the average discharge for each basin cell. ID corresponds to numbered basins shown in Fig. 4 and to the glacier basin codes used by Rignot and others (2008). 1 JUT = Jutulstraumen; 2 AST = Princess Astrid coast glacier; 3, 4, 5 & 10 RAG = Princess Ragnhild coast glaciers; 6 STA = Stancomb-Wills Glacier; 7 ROB = Robert Glacier; 8 RAY = Raymond Glacier; 9 SHI = Shirase Glacier; 13 BAI = Bailey Ice Stream; 15 SLE = Slessor Ice Stream; 16 SCY = Scylla & other glaciers; 17 SUP = Support Force Glacier; 19 FOU = Foundation Ice Stream; 20 REC = Recovery Ice Stream; 21, 28 & 29 AME = American Highland glaciers; 22, 26 & 30 PHI = Philipi & other glaciers; 24 RUT = Rutford Ice Stream; 31 LAM = Lambert & other glaciers; 33 INS = Institute Ice Stream; 34 FOU/MOL = Foundation/Moller Ice Streams; 36 PIG = Pine Island Glacier; 37 SCO = Scott Glacier; 38 AMU = Amundsen Glacier; 39 MER = Mercer Ice Stream; 40 DEN = Denman Glacier; 41WHI = Whillans Ice Stream; 42 BEA = Beardmore Glacier; 43 KAM = Kamb Ice Stream; 44 NIM = Nimrod Glacier; 45 TWG = Thwaites Glacier; 46 BUD = Budd Glacier; 47 BIN = Bindschadler Ice Stream; 48 & 51 TOT = Totten Glacier; 52 BYR = Byrd Glacier; 54 MAC = MacAyeal Ice Stream; 55 & 58 MOS = Glaciers flowing to Moscow University ice shelf; 59 MUL = ; 60 MOS = ; 62 & 63 FRO = Frost Glacier; 64 DAV = David Glacier; 65 DIB = Dibble Gacier; 66 MER = Mertz Glacier; 67 NIN = Ninnis Glacier; 68 & 69 COO = Glaciers flowing into Cook ice shelf.

As mentioned above, the Bedmap2 surface DEM, and particularly the bed DEM, contain errors which produce associated errors in the φ field. Previous work has shown that small changes in the surface topography or moderate changes in the bed topography can alter the distribution of φ and therefore drainage network patterns and the size and shape of drainage catchments. For example, Wright and others (2008) found that raising three grid cells in the surface DEM by 10 m at the southern end of the Adventure Subglacial Trench switched water flow from the Byrd Glacier catchment to the Mulock Glacier catchment. Subsequently lowering another four surface DEM cells by 10 m at the northern end of the Trench caused another switch which resulted in water entering the Cook Glacier system (Wright and others, 2008). Changes in the elevation of bed DEM cells of approximately ten times those given above would have produced similar results. A full continent wide sensitivity test investigating how our calculations of subglacial drainage catchments and pathways are dependent on surface and bed DEM uncertainties is beyond the scope of this paper. Given this sensitivity found by Wright and others (2008), it is likely that future improvements in surface and bed data would lead to some changes in the location of drainage catchments and pathways, and their calculated water fluxes reported here.

Known subglacial lakes, sinks and drainage pathways

Of the 379 known lakes (Wright and Siegert, 2012), 361 occur in zones with melting bed conditions, but 18 occur in places that are predicted to be frozen (Fig. 3). This is likely due to errors in the mapping of the geothermal heat flux, and the calculations of strain heating and heat advection in the ice flux flow model. The number of known subglacial lakes (N = 379) is two orders of magnitude less than the number of predicted lakes (N = 32 380) (Fig. 1). Many of the predicted lakes from the sink filling algorithm may not be real, but artefacts of the interpolation routine associated with producing the surface and bed DEMs. Conversely, some may represent real lakes, but lakes that have not yet been mapped using airborne radar or satellite InSAR/altimetry methods and which do not, therefore, appear in the inventory. As with the predicted lakes, the 379 known lakes are biased towards areas of locally low bedrock elevation (lake elevations: $\bar x = - 290.1\,{\rm m}$ , σ = 615.2 m; continent-wide bed elevations: $\bar x = 92.8\,{\rm m}$ , σ = 677.3 m). Known lakes frequently occur in or near predicted lakes (distance between known lakes and the centre of the nearest predicted lake: $\bar x = 6.30\,{\rm km}$ , σ = 6.32 km) (Fig. 1). These associations are most apparent in East Antarctica where many large lakes have been identified. For example, the largest predicted lake which occurs in a topographic low and acts as a confluence for several drainage pathways, overlaps with the known location of Lake Vostok (Fig. 6b). Known lakes also tend to occur on or close to the main drainage pathways (defined as where accumulated melt flux >0.1 m a⁻¹); the distance between known lakes and the nearest drainage pathway is: $\bar x = 1.52\,{\rm km}$ , σ = 2.12 km.

The predicted lakes are often close in size to the calculated areas of known lakes taken from the literature (Supplementary Materials, Table S1). In the Supplementary Materials Table S1, all predicted lakes within a 20 km radius of known lakes are compared with the calculated areas of known lakes. As examples, the predicted lake around Lake Vostok (Fig. 6b) has an area of 13 811 km², about the same as the calculated known lake area of 14 000 km² (Wright and Siegert, 2012). The predicted lake corresponding to the 90°E Lake is 1830 km² compared to the known lake area of 2000 km² (Wright and Siegert, 2012). The predicted lake for Academy12 Lake is 297 km² compared to the calculated known lake size of 207 km² (Smith and others, 2009) and that for Whillans1 (Subglacial Lake Engelhardt) is 339 km² compared to the known lake size of 362 km² (Fricker and Scambos, 2009).

In other instances, the predicted lakes within 20 km of known lakes underestimate the calculated areas of the lakes (Supplementary Materials, Table S1). For example, the calculated area of CookE2 Lake is 260 km² but the predicted lakes within 20 km amount to just 9 km². However, this figure increases to a more comparable 213 km² if the search radius is increased to 40 km. Similarly, Lakes A–D at the onset of the Recovery ice stream (Fig. 6a) have a combined area of 13 330 km² (Bell and others, 2007) but the predicted lakes within 20 km sum to just 563 km². This increases to 7128 km² with a larger search radius of 50 km. The discrepancies between the predicted lake areas and the calculated areas of known lakes are likely to be due largely to inaccuracies in the ice surface, and particularly the bed topography data, from which the φ field and predicted lakes were derived. This is a bigger problem for the Recovery ice stream onset area, where Bedmap2 ice thickness data were derived from a satellite-based gravity model, than for other regions, where ice thickness was mapped by airborne radar survey (Fretwell and others, 2013).

The obvious associations between known subglacial lakes, predicted lakes and major drainage pathways supports previous assessments that lakes are not transient features but exist for long periods of time (Vaughan and others, 2008; Wright and others, 2008). As the surface topography changes, so too will the gradients of φ, the locations of lakes, and the positions of drainage pathways. Such changes will typically occur on centennial to millennial timescales, although they may occur more rapidly, for example due to sudden steepening of outlet glaciers in response to ice shelf breakup (Scambos and others, 2011).

The size, shape or exact location of lakes cannot be inferred precisely from the sink characteristics, but large sinks, in areas of good bedrock source data, will greatly increase the probability that one or more lakes will lie within a few km of the large sink location. This finding could help guide the selection of target regions for future lake surveys. For example, numerous sinks occur in the western portion of Dronning Maud Land (i.e. ~10°E, 75°S; Fig. 1), which has reasonable coverage in terms of ice thickness in the Bedmap2 data but does not feature in the lake inventory dataset.

Subglacial water fluxes to the coast

The 12 595 individual basins discharge a total of 21.4 km³ a⁻¹ of meltwater, equivalent to 6.0 × 10⁻⁵ m a⁻¹ of global sea level assuming a conversion factor of 2.8 × 10⁻⁶ m km⁻³ (Vaughan, 2006). This is approximately a third of the average value, 65 km³ a⁻¹, calculated in a suite of model runs in a similar study by Pattyn (2010), who used a different ice flow model, and a variety of different input fields compared to those used here. These differences highlight the difficulties in predicting basal melt fluxes beneath the ice sheet given the assumptions made in different flow models, and the uncertainties in topographic, basal geothermal and surface accumulation data. The water fluxes of the 54 largest basins by area are given in Table 1. These 54 basins discharge 18.8 km³ a⁻¹, which is 87.9% of the total flux, equivalent to 5.3 × 10⁻⁵ m a⁻¹ of global sea level. This basal melt flux is an important component of the Antarctic ice-sheet mass balance. It is usually overlooked in the input/output method of ice-sheet mass-balance calculation, although it is accounted for in the altimetry and gravimetry approaches. This is relevant because comparisons of Antarctic ice-sheet mass balance calculated by the different methods show that the input/output method produces estimates that are more negative than the other two methods, although within the errors of the different techniques (Rignot and others, 2011; Shepherd and others, 2012). Including the basal melt flux in the input/output method produces mass-balance estimates that are even more negative, increasing the discrepancy between estimates derived by this method and those derived by the other two methods. The basal melt water fluxes from the entire continent and from the primary basins represent 31 and 25% respectively of the average annual mass balance of the ice sheet for the period 1992–2011, which was –71 ± 53 km³ a⁻¹ (Shepherd and others, 2012).

Of the 54 largest basins, 17 have discharges in excess of 0.3 km³ a⁻¹ and five have discharges over 1 km³ a⁻¹ (Fig. 4; Table 1). In flux order these five are Byrd Glacier, Foundation/Moller ice streams, Lambert Glacier, Recovery ice stream and Slessor Glacier. Although the largest basin by area, Byrd Glacier, is also the largest in terms of outlet discharge, the rankings of area and discharge are not always the same (Table 1). For example, the second largest basin in terms of area (Lambert Glacier) is the third largest in terms of outlet flux and Budd Glacier is the fifth largest basin by area but is ranked 33rd by flux. Thus, outlet flux is a function not only of catchment size but also melt rates (patterns of geothermal heat flux and heat produced from ice strain). The Totten Ice Stream basin (48) is ranked 36th in terms of area and 12th in terms of flux, but it has the highest discharge per unit area of all the catchments due to the relatively high geothermal heat fluxes and strain heating, producing 8.78 mm a⁻¹ of melt averaged across the basin (Table 1). By contrast, Byrd Glacier ranked first by both area and flux is ranked just 20th in terms of flux per unit area (2.23 mm a⁻¹).

In terms of fluxes to beneath the large ice shelves, ~6.5 km³ a⁻¹ of freshwater is discharged to the Ronne–Filchner ice shelf. The Ross and Amery ice shelves receive ~4.9 and ~2.9 km³ a⁻¹ respectively. Thwaites Glacier's ice shelf receives 1.6 km³ a⁻¹ and the ice shelf of Pine Island Glacier gets 0.6 km³ a⁻¹. These discharges are comparable to the calculated subglacial water flux of Jakobshavn Isbrae in winter of ~2.0–3.3 km³ a⁻¹ (Echelmeyer and Harrison, 1990), when water will be derived mostly from geothermal and strain energy (Mernild and others, 2010).

The main difference between the floating ice tongues of Greenland and the ice shelves of Antarctica is their area. For example, the tongue of Jakobshavn Isbrae was ~105 km² prior to 1997 (Motyka and others, 2011), giving winter water fluxes per unit area of 19–31 m a⁻¹. By contrast, the areas of the Antarctic ice shelves are between one and three orders of magnitude larger than the Greenland ice tongues: Ronne–Filchner, ~427 000 km²; Ross, ~474 000 km²; Amery, ~60 000 km²; Thwaites, ~3930 km²; and Pine Island, ~3570 km² (Fretwell and others, 2013). The combination of water flux and ice shelf area give the following water fluxes per unit area to the Antarctic ice shelves, which are two to three orders of magnitude less than those beneath Jakobshavns Isbrae: Ronne–Filchner, 0.015 m a⁻¹; Ross, 0.010 m a⁻¹; Amery, 0.048 m a⁻¹; Thwaites, 0.393 m a⁻¹; and Pine Island, 0.179 m a⁻¹. These low fluxes per unit area mean that any significant melting beneath the Antarctic ice shelves associated with freshwater plumes must be highly localised leading to upward channel incision to the underside of the ice shelves, as suggested previously by others (Jenkins, 2011; Carter and Fricker, 2012; Le Brocq and others, 2013; Langley and others, 2014; Marsh and others, 2016). Such channels may provide lines of weakness which may be exploited during calving events (O'Leary and Christoffersen, 2013). In this regard, the higher flux per unit area for the Thwaites ice shelf compared to the others suggests that it may be particularly susceptible to instability.

Subglacial water fluxes from subglacial lakes

The calculated fluxes from the 232 known lakes that lie within 10 km of a predicted lake that accumulates >0.1 m a⁻¹ are shown in Fig. 7. The fluxes from these lakes range from 1.05 × 10⁻⁴ to 1.49 km³ a⁻¹ ( $\bar x = 3.45\, \times \,{10^{ - 2}}\,{\rm k}{{\rm m}^3}\,{{\rm a}^{ - 1}}$ , σ = 0.13 km³ a⁻¹). The fluxes from the ‘active’ lakes in Smith and others's (2009) inventory, which typically lie beneath the large outlet glaciers and ice streams, are generally higher than those from the other lakes that tend to lie in the ice-sheet interior (active lake fluxes: $\bar x = 6.19 \times {10^{ - 2}}\,{\rm k}{{\rm m}^3}\,{{\rm a}^{ - 1}}$ , σ = 1.76 × 10⁻¹ km³ a⁻¹; other lake fluxes: $\bar x = 5.76\, \times {10^{ - 3}}\,{\rm k}{{\rm m}^3}\,{{\rm a}^{ - 1}}$ , σ = 1.65 × 10⁻² km³ a⁻¹; Supplementary Materials, Table S2).

Fig. 7. Frequency histogram of water fluxes from 232 (60%) of the 379 known subglacial lakes. The 232 are those that lie within 10 km of a predicted lake that accumulates at least 0.1 m a⁻¹ of melt flux.

The 18 highest ranked lakes in terms of calculated water flux are all ‘active’ lakes identified by Smith and others (2009). The top four are, in order, Byrd1, Byrd S5, Recovery 8 and Byrd S6. Three lakes known to be very active beneath the Mercer and Whillans Ice Streams are ranked 5, 6 and 7; respectively Lakes Mercer, Engelhardt and Upper Conway. In contrast, most of the large known lakes in the interior of the ice sheet are ranked very low in terms of water flux (Supplementary Materials, Table S2). For example, the fluxes from the Adventure Trench Lake, 90°E Lake and Lake Sovetskaya are ranked 119, 155 and 185 respectively; Recovery Lakes A–D at the onset region of the Recovery Ice Stream are ranked 210, 43, 64 and 62 respectively.

The analysis above shows that there is a general associated between high calculated water fluxes and ‘active’ lakes that have been observed to fill and drain, and between low water fluxes and other lakes that have not been observed to fill and drain. There are anomalies, however. For example, Lakes ‘DCS/DCSx/X01d’ and ‘DCS/DCSx/X02b-X02e’ in the Byrd catchment (Wright and Siegert, 2012) are ranked 19 and 21 respectively in terms of water flux (Supplementary Materials, Table S2). We suggest that these are potentially ‘active’ lakes that have not, so far, been observed to fill/drain. Conversely, several active lakes receive relatively low fluxes; for example Whillans 2a (Lake 12) and Institute W2 are ranked 227 and 229 in terms of flux but were observed to fill and drain from ICESat altimetry (Smith and others, 2009). This may be due to errors associated with our calculations of subglacial meltrates, or φ gradients and therefore drainage pathways. Alternatively, it may demonstrate that periodic lake filling and draining is a function of processes other than simply the steady-state water flux to the lake, for example the alteration of local hydraulic gradients through the filling and draining of nearby lakes.

Lake turnover times and flood RIs

The steady-state turnover times for the 232 lakes that lie within 10 km of a predicted lake receiving a flux >0.1 m a⁻¹ are given in Supplementary Materials, Table S3 and shown in Figs 8a and 9a. For the ‘active’ lakes, the turnover times range from about a month (Bindschadler 6) to nearly 20 000 a (Recovery 1). For the other ‘non-active’ lakes, turnover times span ~2 months (Lake DCS/DCSx/X01d, Byrd Glacier) to over 400 000 a (Lake U2, Adventure Trench). The median turnover times are ~100 a and ~500 a for the ‘active’ and ‘non-active’ lakes respectively.

Fig. 8. Frequency histograms of (a) turnover times and (b) RIs of a 0.25 km³ flood, for 232 of the 379 known subglacial lakes. Lakes are split into 116 ‘active’ lakes (from Smith and others, 2009) and 116 other lakes that have not been observed to fill/drain using ICESat altimetry. Note the logarithmic scale on the x-axes.

Fig. 9. Maps showing: (a) the turnover times, and (b) the RIs of a 0.25 km³ flood, for 232 of the 379 known subglacial lakes. Lakes are split into 116 ‘active’ lakes (from Smith and others, 2009) and 116 other lakes that have not been observed to fill/drain using ICESat altimetry.

Having estimated the total volume of all 232 lakes receiving an accumulated flux >0.1 m a⁻¹, we now focus on the 116 ‘active’ lakes in that list and compare their total volumes with the maximum flood volumes observed from ICESat altimetry (Smith and others, 2009, their Table 4). The observed flood volumes vary between 0.2% and 177% of the total lake volumes. Clearly, it is impossible for a lake flood volume to be more than 100% of its total volume, therefore flood volumes that are >100% we attribute to errors, either in the calculation of flood volumes (which Smith and others, 2009 claim may be as much as 50%) or in our calculation of total volumes (which are likely to be underestimates due to inaccuracies in the Bedmap2 bed DEM). The median flood volume to total volume ratio is 17% ( $\bar x = 37\% $ , σ = 46%).

The mean flood volume from the ‘active’ lakes observed by ICESat was approximately 0.25 km³ (Smith and others, 2009, their Table 4). Assuming this represents a typical flood volume, we calculate the flood RIs for each of the 116 ‘active’ lakes. We also do this for the remaining 116 lakes that were not observed to fill or drain during the ICESat era. Results are given in the Supplementary Materials, Table S4 and depicted in Figs 8b and 9b. The calculations indicate RIs of between ~2 months and ~2000 a (median = 131 a; $\bar x = 316\,{\rm a}$ , σ = 412 a) for floods from the ‘active’ lakes and RIs of ~2 a and ~2400 a (median = 359 a; $\bar x = 518\,{\rm a}$ , σ = 465 a) for floods from the ‘non-active’ lakes. Thus, floods comparable to the average flood observed from ICESat altimetry data typically occur more frequently from the ‘active’ lakes beneath the ice streams towards the ice-sheet margins and less frequently from lakes that have not been observed to flood beneath thicker ice towards the ice-sheet interior. We calculate that 35 of the ‘active’ lakes have a RI < 10 a (Fig. 8b). This compares with 10 lakes that underwent more than one fill/drain cycle over the 5 a (statistically equivalent to 20 lakes with a RI < 10 a) observed by Smith and others (2009). Given the uncertainties in our calculations, these statistics compare well. For example, we can match the observed 20 lakes with a RI < 10 a if we increase the assumed flood volumes from 0.25 to 1 km³ or if we reduce the steady-state water fluxes to the lakes by 70%.

Some of the most active lakes are beneath Byrd Glacier (Byrd1, ByrdS5 and ByrdS6 are ranked 1st, 2nd and 4th respectively) and Mercer and Whillans Ice Streams (Lakes Mercer, Engelhardt and Upper Conway are ranked 5th, 6th and 7th respectively). Their RIs are comparable to that measured for Lake Conway of 6.5 a (Siegfried and others, 2014). Lake Recovery8 is ranked 3rd, and Lakes Nimrod2, Rutford1 and David1 are ranked 8th, 9th and 10th respectively. However, the analysis shows that several lakes that were not observed to fill or drain from ICESat altimetry have the potential to fill/drain at periodicities as low as a few years. For example, the RI for a 0.25 km³ flood is 56 a for Lake Vostok and 43, 65 and 63 a respectively for Recovery Lakes B, C and D. Of course, these ‘non-active’ lakes may not undergo rapid transfers of water, in which case these RIs are simply the steady-state turnover times for 0.25 km³ of water.

Since meltwater flux accumulates down-glacier, the lakes with the greatest steady-state flux are those close to the grounding lines of glaciers in relatively large catchments with relatively high rates of basal melt production. Such lakes also tend to be the ones with the shortest RI between lake drainage events (Fig. 9b; Supplementary Materials, Table S4). For this reason, attempts to identify short-term fluctuations in ice velocity associated with lake drainage events (e.g. Stearns and others, 2008) should focus on lakes just above the grounding lines on Byrd, Recovery, Mercer, Whillans, Nimrod, Rutford and David Glaciers/Ice Streams (Fig. 9b; Supplementary Materials, Table S4).

SUMMARY AND CONCLUSIONS

• • Using the Antarctica Bedmap2 datasets, we have mapped the subglacial hydraulic potential (φ) field, the area, depth and volume of hydraulic potential sinks (predicted lakes) and the locations of the main drainage pathways and catchments beneath the ice sheet. We derive basal melt rates from a balanced ice flux flow model and accumulate them along the drainage pathways to provide first order calculations of subglacial water fluxes beneath the ice sheet, to the ice shelves and from subglacial lakes.

• • The number of predicted lakes exceeds the number of known lakes by two orders of magnitude, and the total volume of all predicted lakes is 7323 km³, a sea level equivalent of 0.021 m. The known lakes frequently occur in or near predicted lakes. The size, shape or exact location of lakes cannot be inferred precisely from the sink characteristics, but the presence of large sinks may be a useful indicator for the presence of large subglacial lakes by increasing the probability that one or more lakes occurs close by. This finding could help guide the selection of target regions for future lake surveys.

• • The ice-sheet drains a total of 21.4 km³ of meltwater annually, equating to 6.0 × 10⁻⁵ m a⁻¹ of global sea level. This flux represents ~31% of the total mass balance of the ice sheet and is usually overlooked in the input/output method of mass-balance calculation although it is accounted for in the altimetry and gravimetry approaches.

• • Around 6.5 km³ a⁻¹ of freshwater is discharged to beneath the Ronne–Filchner ice shelf. The Ross and Amery ice shelves receive ~4.9 and 2.9 km³ a⁻¹ respectively. Thwaites Glacier's ice shelf receives ~1.6 km³ a⁻¹ and the ice shelf of Pine Island Glacier gets ~0.6 km³ a⁻¹. Thwaites has the highest flux per ice shelf area and might therefore be expected to be the least stable.

• • Steady-state water fluxes from known subglacial lakes range over five orders of magnitude between 10⁻⁴ and 10⁰ km³ a⁻¹, with the smallest fluxes from inland lakes and the largest fluxes from lakes close to the grounding lines of the major glaciers and ice streams. The varying water fluxes to the different lakes might be expected to correlate with varying rates of sediment delivery to and deposition within the lakes (Bentley and others, 2011).

• • Lake turnover times range from about a month to nearly 20 000 a for the known ‘active’ lakes and between ~2 months and 400 000 a for lakes which have not been observed to fill/drain. The median turnover times are ~100 a and ~500 a respectively. These times have important implications for life beneath the ice sheet, with longer water residence times increasing the potential for unique life forms (Priscu and others, 1999; Christner and others, 2014).

• • RIs for a flood of 0.25 km³ (the average measured by ICESat; Smith and others, 2009) range between ~2 months and ~2000 a (median ≈ 130 a) for the observed ‘active’ lakes, and between ~2 a and ~2400 a (median ≈ 360 a) for the known lakes that have not been observed to fill/drain. Thus, several lakes that have recently been observed to fill and drain may not do so again for many centuries. Conversely, many lakes that have not, so far, been observed to fill or drain may experience floods similar in magnitude to those that have been observed from the ‘active’ lakes, but less frequently. Alternatively, they may experience floods as frequently but their flood magnitudes will be less.

• • Attempts to identify short-term fluctuations in ice velocity associated with lake drainage events should focus on lakes with the shortest RIs between floods, notably lakes just above the grounding lines on Byrd and Nimrod Glaciers and on the Recovery, Whillans, Mercer and Rutford Ice Streams.