Cascading water underneath Wilkes Land , East Antarctic ice sheet , observed using altimetry and digital elevation models

We describe a major subglacial lake drainage close to the ice divide in Wilkes Land, East Antarctica, and the subsequent cascading of water underneath the ice sheet toward the coast. To analyse the event, we combined altimetry data from several sources and subglacial topography. We estimated the total volume of water that drained from Lake CookE2 by differencing digital elevation models (DEM) derived from ASTER and SPOT5 stereo imagery acquired in January 2006 and February 2012. At 5.2 ± 1.5 km3, this is the largest single subglacial drainage event reported so far in Antarctica. Elevation differences between ICESat laser altimetry spanning 2003–2009 and the SPOT5 DEM indicate that the discharge started in November 2006 and lasted approximately 2 years. A 13 m uplift of the surface, corresponding to a refilling of about 0.6 ± 0.3 km3, was observed between the end of the discharge in October 2008 and February 2012. Using the 35-day temporal resolution of Envisat radar altimetry, we monitored the subsequent filling and drainage of connected subglacial lakes located downstream of Cook E2. The total volume of water traveling within the theoretical 500-km-long flow paths computed with the BEDMAP2 data set is similar to the volume that drained from Lake CookE2, and our observations suggest that most of the water released from Lake Cook E2 did not reach the coast but remained trapped underneath the ice sheet. Our study illustrates how combining multiple remote sensing techniques allows monitoring of the timing and magnitude of subglacial water flow beneath the East Antarctic ice sheet.


Introduction
It has been known for decades that liquid water can exist and accumulate in lakes at the base of the Antarctic ice sheet (e.g. Siegert and Dowdeswell, 1996). Geothermal heat flux, ice deformation and ice bed friction provide heat that is trapped by the thick insulating ice sheet and can melt a few millimetres of basal ice each year (e.g. Clarke, 2005;Pattyn, 2010). The water at the base of the ice sheet lubricates the ice-bedrock interface and the diminution of friction creates distinctive topographic features on the ice-sheet surface such as low slopes areas above the lakes (Rémy et al., 1999;Rémy and Legrésy, 2004).
More recent observations show that this water can migrate rapidly under the ice sheet and induce movements of the surface as ice is lifted or lowered by floatation when lakes fill in or drain (Clarke, 2006). The advent of new satellite observation techniques has allowed the measurements of these surface elevation changes (e.g. Gray et al., 2005;Wingham et al., 2006;Fricker et al., 2007), which led to a better general understanding of the dynamic nature of subglacial hydrological systems (Fricker and Scambos, 2009). The impact of subglacial water on ice-sheet dynamics is still being investigated (Bell et al., 2007;Stearns et al., 2008;Schoof, 2010) and the modelling of the phenomenon is the subject of active researches (e.g. Johnson and Fastook, 2002;Goeller et al., 2013). Smith et al. (2009) provided the first continent-wide inventory of active subglacial lakes in Antarctica from repeat ICE-Sat laser altimetry. Their inventory includes the drainage of a lake, Cook E2  2008 along several kilometres on two crossing ICESat tracks. They estimated that 2.7 km 3 of water was discharged during this event, the largest drainage observed so far in Antarctica. The several-tens-of-metres drawdown of the surface during this drainage is similar in magnitude to the surface drawdown due to Jökulhlaups (e.g. Björnsson, 2003) and more than 10 times the second largest drawdown reported by Smith et al. (2009).
Here, we revisit this event in more detail using laser altimetry and by differencing multi-temporal digital elevation models (DEMs) from stereo imagery. Those additional data sets allowed us to recalculate the total volume of water that drained from the Lake Cook E2 and document the changes that this lake has undergone since the end of its drainage. Our values are compared to those reported in a similar recent study using CryoSat-2 radar altimetry (McMillan et al., 2013). Smith et al. (2009) also pointed out a possible connection between Lake Cook E2 and another lake, Cook E1 , located 100 km downstream. To further investigate the effect of the large amount of water that drained out of Lake Cook E2 , we analysed time series of ice-sheet elevations derived every 35 days from the Envisat radar altimeter.

Data
In this study, we used a combination of altimetry measurements derived from a radar altimeter (Envisat RA-2), a laser altimeter (ICESAT GLAS) and two optical stereo imagers (SPOT5 HRS and Terra ASTER). An overview of those satellite data is given in Table 1.

ICESat measurements
ICESat measurements were taken from release 531 available from the NSIDC (http://nsidc.org/data/icesat/index.html) and span campaigns 2A to 2F, from October 2003to October 2009. ICESat measurements were first converted from the TOPEX/Poseidon to the WGS-84 ellipsoid. ICESat altitudes, originally defined above EGM2008, were corrected for the elevation difference between EGM2008 and EGM96, EGM96 being the reference for the DEMs (see below).

Digital elevation models from stereo imagery
DEMs of the area surrounding Lake Cook E2 were constructed from stereo pairs of satellite optical images acquired in 2006 by the Advanced Spaceborne Thermal Emission and Reflection Radiometer, ASTER (e.g. Kääb, 2002), and in 2012 by the "Haute Resolution Stéréoscopique" (HRS) sensor onboard SPOT5 (e.g. Bouillon et al., 2006). The date of the image acquisitions and their resolution are provided in Table 1.

BEDMAP2
We use the BEDMAP2 subglacial topography and ice-sheet surface DEMs (Fretwell et al., 2013) to compute the hydrologic potential at the base of the ice sheet and to infer the flowpath of subglacial water draining from Lake Cook E2 .

Envisat measurements
We use elevation measurements from the Envisat radar altimeter (RA-2) during 85 repeat cycles from September 2002 to October 2010.
The main frequency of the altimeter is within Ku-band at 13.75 GHz. Microwaves at this frequency penetrate through the snow pack. However, here we do not mix Envisat data with data acquired at different wavelengths of the electromagnetic spectrum. Thus, we can overlook the issue of calibrating elevation between different missions. The variation of penetration depth with time is corrected as in Flament and Rémy (2012). See Sect. 3.3 for further details on the processing of Envisat data.

DEM differencing
Comparing multi-temporal DEMs is now commonly used for measuring elevation changes of glaciers (e.g. Gardelle et al., 2013) or at the steeper ice-sheet margins (e.g. Howat et al., 2007), where other altimetry techniques do not perform as well or lack sufficient spatial coverage. To our knowledge, DEM differencing has only been used once to characterize the drainage of a subglacial lake, in response to the rapid thinning of Crane Glacier in the Antarctic Peninsula (Scambos et al., 2011). This is probably because DEMs from optical imagery are still regarded as being of limited use for the study of the flat and uniformly white central part of ice sheets. The lack of identifiable features on their surface (e.g. Bindschadler et al., 2011) often precludes the accurate computation of correlation between the images of the pair, a necessary step in the computation of DEMs from stereo images.
Here, we benefit from some recent improvements achieved by the French mapping agency (IGN, National Institute of Geographic and Forest Information) in the processing of stereo pairs on low contrast areas to derive a nearly complete DEM of the Lake Cook E2 area from SPOT5-HRS imagery acquired on 9 February 2012, about 3 years after the drainage ended. Those improvements consisted in (i) reducing the size of the final correlation window to match small-scale variations in radiometry present in the two images of the stereo pair; and (ii) accepting DEM pixels with a lower correlation score (CS, ranging from 0 to 100). All pixels where the CS is greater than 30 are now retained, whereas the standard IGN algorithm retained only pixels with a CS greater than 50. Those modifications improved the DEM coverage. Gaps in the SPOT5 DEM only covered 1 % of the study area vs. 51 % using the standard processing of the SPIRIT project (Korona et al., 2009). Importantly, it was verified using ICE-Sat data that the newly retained DEM pixels had the same accuracy as the pixels retained in the standard algorithm. The ice-sheet surface topography prior to the start of the lake drainage was derived from an ASTER stereo pair acquired on 26 January 2006. The ASTER DEM and the corresponding correlation mask were calculated using PCI Geomatica 10.1 (Toutin, 2008). Unreliable pixels, covering 60 % of the area of interest, were masked in the ASTER DEM. Both DEMs are computed from stereo imagery without ground control points and, thus, can have significant biases. Horizontal biases (e.g. Nuth and Kääb, 2011) cannot be estimated here due to the lack of stable terrain but can be neglected because of the very low slopes of the area of interest. Given the mean slope for our study area (0.1 • ), a typical planimetric shift of 30 m would lead to an elevation error of only 5 cm. Vertical biases in the pre-and post-drainage DEMs were estimated and corrected using ICESat. We extracted the ASTER and SPOT5 DEM elevations at the location of each ICESat footprint using bilinear interpolation. Compared with all the data from laser period 3E (mean date: 12 March 2006; before the drainage of Lake Cook E2 ), a vertical bias of −7.2 m (standard deviation: SD = 11.7 m; N = 1039) was found and corrected for the ASTER DEM. The SPOT5 DEM was also compared to laser period 3E but only outside the lake area (a first guess of the lake area was obtained by differencing the non-adjusted ASTER and SPOT5 DEM). A mean bias of 11.6 m (SD = 2.1 m; N = 8246) was found. There was a tilt in the SPOT5 DEM in the northsouth direction which we corrected by fitting a plane to the elevation differences with ICESat. The tilt is probably due to errors in the measurements of the position and attitude (roll/pitch/yaw) of the satellite. After correction of this tilt, the standard deviation of the elevation difference is reduced to 1.6 m (laser period 3E, N = 8246). This vertically adjusted SPOT5 DEM was then compared to ICESat data from other laser periods outside of the lake area. The absolute mean difference was always lower than 0.5 m, indicating the stability of the region surrounding the lake. The standard deviations of the elevation differences were similar for all laser periods (range: 1.50 to 1.67 m).
Reliable values in the ASTER DEM covered only 40 % of our area of interest. To obtain a complete pre-event topography, we combined those reliable values from the ASTER DEM with elevation from the SPOT5 DEM (away from the lake) and from the pre-drainage ICESat laser campaigns. First, altitude contour lines were generated with an altitude step of 5 m from the vertically adjusted SPOT5 DEM. They were used only away from the lake. Then, in the Lake Cook E2 area, we also generated 5 m altitude contour lines from the vertically adjusted ASTER DEM. However, those contour lines are not complete due to data voids in the ASTER DEM. Thus, to obtain continuous contour lines, we interpolated manually between those discontinuous contour lines using the altitude from the pre-drainage ICESat laser campaigns as a constraint. To do so, we extracted from pre-drainage ICE-Sat data the locations where ICESat tracks had the same elevation as the DEM contours. Each ASTER contour line was then extended manually to pass through ICESat points with the same elevation. The complete altitude contour map was then gridded to obtain a pre-event DEM without data voids. The pre-and post-event DEMs are then differentiated to obtain the map of elevation differences and to infer the total volume loss during the drainage.

Subglacial water flowpath inferred from BEDMAP2
Following Shreve (1972), we assumed that the hydrologic potential at the ice-sheet base is a function of the bedrock elevation and the thickness of the overlying ice column. From the BEDMAP2 elevation and ice thickness maps, we compute the hydrologic potential (H p ) as where ρ w and ρ i are respectively the densities of water and ice, and h b and h s the elevations of the bed and the ice-sheet surface. We thus obtain a map of hydrologic potential with a 1 km grid spacing from which we can derive flowpaths. This grid spacing used here (1 km) is the one provided in the BEDMAP2 products and is denser than the actual resolution of the bed elevation data in our study area. The coarser resolution is taken into account as we use the 1 km grid with a routing algorithm that computes the hydropotential gradient from the hydrologic potential at all points within a 5 km radius disc, which is closer to the actual resolution of the BEDMAP2 data set in most of our study area. This algorithm is described below.
With the assumption that water will flow following the steepest gradient of the hydrologic potential, we can compute possible flowpaths for the water draining from Lake Cook E2 . We use a gradient descent algorithm, starting from Lake Cook E2 . The slope of the hydrologic potential is computed by fitting a plane on a disc of 5 km radius and taking steps of 5 km in the direction of steepest descent. In order to account for the uncertainty in the BEDMAP2 data set (up to several hundreds of metres, see Fig. 2d), we use a Metropolis-Hastings algorithm (Hastings, 1970). We allow the step direction to vary around the direction of steepest descent, proportionally to the accuracy provided with the map of the bedrock, that is, the greater the uncertainty in the bedrock elevation, the greater the angle by which the descent direction is allowed to vary. Our area of interest is not perfectly covered by BEDMAP2, with large data gaps so that over large areas, the bed elevation is interpolated from radar measurements performed at a few tenths of km (see Fig. 3 in Fretwell et al., 2013). While the upper part of our region of interest is relatively well surveyed, in places, uncertainty in the bed elevation can reach more than 300 m (Fig. 2). With such uncertainty, the derived flowpath is expected to spread in the lower reaches of our study area, close to the coast.
We then draw 10 000 flowpaths to obtain a large set of possible flowpaths. The envelope of all those flowpaths is used to delineate the area where we will search for subglacial lakes (Sect. 3.4).

Processing of Envisat data
The along-track processing applied to the Envisat radar altimetry data has been presented in detail elsewhere (Flament and Rémy, 2012) and is only briefly summarized here: elevation measurements are collected in bins of 1 km along the satellite track and processed by a least square fit. We fit 10 parameters that account for local topography (mean altitude, slope and curvature, i.e. a surface controlled by 6 parameters), echo shape (backscatter, leading edge width and trailing edge slope) and a linear temporal trend. This set of parameters forms the least square model. The processing provides an estimate of the best fitting coefficient for each parameter in this model and the residuals. The fitted coefficients are used to correct elevations taken at different dates and locations and make possible a comparison of relative elevations along time. The residuals, that is, the elevation differences between the best fitting model and the actual data, contain all elevation changes that cannot be explained by the simple least square model. In this case, the residuals include non-linear temporal variations caused by a subglacial flood and noise. Elevation time series are computed from the residuals and the elevation trend; at time t, the elevation is given by where ε(t) is the residual at time t, dh/dt is the linear elevation trend and t0 is an arbitrary reference time. Here, we set t0 = 2006.85 because it is the date of the last measurement before the drainage of Lake Cook E2 (see below).
Noise in the processed series is computed as the standard deviation of the difference between the time series and the same time series smoothed with an 11-point moving average. In other words, the higher frequencies are considered as noise and removed. This noise level is very dependent on the kilometre-scale topography and can vary from under one decimetre in a flat and horizontal area to several metres over rugged terrain (crevasses, faster flowing areas, etc.). In what follows, time series are considered of acceptable quality if this noise is below 1 m.
Envisat cannot resolve the drainage of Lake Cook E2 because the walls of the trough are too steep (see discussion, in Sect. 5.2). If we use only measurements from the predrainage period, before the dramatic change of the topography, the simple quadratic model for the surface topography performs very well and the noise level is lowered to 0.1 m (Fig. 3). Thanks to this pre-drainage elevation time series, we are able to detect accurately the first large perturbation of surface elevation, around 2006.85 (November 2006), our best estimate for the start of the drainage.

Automatic identification of active lakes
We fit each Envisat elevation time series with a Gaussian curve with varying start dates, widths and amplitudes. This allows us to automatically identify time series affected by "anomalous elevation changes" due to the flood following the drainage of Lake Cook E2 (i.e. locations where subglacial lakes filled and then drained). The algorithm is applied to a zone encompassed within the envelope of all water flowpaths derived from BEDMAP2 (Sect. 3.2). The lake-detection algorithm assigns a start date, an amplitude and a flood duration to each time series along with the correlation to the fitted Gaussian curve. If this correlation is higher than a given correlation threshold (ranging between 0.6 and 0.9), the time series is considered as a potential lake. The screening of the potential lakes is then refined using the following criteria: -The start date of the event, defined as the peak date of the Gaussian curve minus one standard deviation, must not be too early. The earliest allowed start date is 2006.85 (see Sect. 3.3).
-The noise level in the Envisat time series must be below 1 m.
-The series must have enough samples, a minimum of 70 measurements out of 85.
After detecting those active lakes, the cumulative volume of water is then obtained by using the full elevation time series, not the fitted Gaussian curve. The latter was only used for the automatic detection of active lakes. The elevation series are smoothed with an 11-point moving average and set to 0 on 2006.85 to provide a common reference date. To convert elevation time series to volume, we approximate the area covered by each time series at 11.85 km 2 because the area of interest (Fig. 1b) covers 50 000 km 2 and is sampled by 4220 Envisat elevation time series. Doing so, we assume that Envisat provided a representative sampling of elevation changes and we thus account for the fact that Envisat did not sample all active lakes due to the separation between the satellite ground tracks. Elevation time series of all detected lakes downstream of Lake Cook E2 are converted to volume and summed to test whether the volume of water that transited through those lakes is similar to the volume that drained out of Lake Cook E2 .
Other surface processes such as snowfalls and firn compaction also contribute to surface elevation changes (e.g. Helsen et al., 2008). Those processes lead to a regional volume change that must be separated from the lake drainage. In order to compensate for this and separate the lake signal from other causes of elevation change, we compute the mean elevation time series for "non-lake points" within the flowpath area and consider that this average "non-lake" elevation series accounts for surface processes. We then remove the "non-lake" series from the "lake" series to obtain corrected time series. The results of this correction are presented in Sect. 4.2.

Drainage of Lake Cook E2
We use a combination of ASTER and SPOT5 DEM differencing to derive the geometry of the surface depression, together with differences between ICESat and the SPOT5 DEM to constrain the timing of the lake drainage and reestimate the total volume of water discharged.
By chance, ICESat tracks cross almost at the locus of maximum surface subsidence (Fig. 4) and thus provide a temporally well-resolved view of the drainage (Fig. 5). In the centre of the main lake, the total surface lowering reached 70 m between laser campaign 3G (November 2006) and 3K (October 2008). The map of elevation differences also shows that, in addition to this main lake, a secondary lake drained simultaneously leading to a maximum surface lowering of 30 m. Since October 2008, a surface uplift of 13 m is observed in the central part of the main depression.
Following others (Smith et al., 2009;Fricker and Scambos, 2009), we assumed that the volume changes measured at the ice-sheet surface equal the volume of water that drained from the lake. We are aware that this hypothesis is questionable (Sergienko et al., 2007), but do not have more information to refine the estimate. The ice sheet is flowing slowly in our study area (Fig. 2a) such that dynamic effects (such as stick-slip) are unlikely to affect surface elevation. The total drainage is thus computed as the product of the "lake area" by the mean elevation difference within this area. The lake area is defined as the closed region where the surface lowering is greater than 3 m. The sensitivity of our result to this arbitrary choice is tested by estimating the lake discharge when the threshold for minimum surface lowering is set to 1.5 and 4.5 m. With a threshold at 3 m, the lake area is 219 km 2 and the mean surface lowering 20.6 m, leading to a volume loss of 4.52 km 3 of water. This volume varies by only 0.10 km 3 (referred to as the "lake-area" error) when the upper (

Fig. 5. Time series of elevation difference between the various
ICESat laser periods and the 9 February 2012 SPOT5 DEM. In black, the mean elevation difference outside of Lake Cook E2 . In blue, the maximum elevation difference between ICESat and the SPOT5 DEM. Note that, due to non-perfect repetitiveness of ICE-Sat tracks, the blue dots do not correspond to the same location for each laser campaign (Fig. 4b) which partly explains why the noise seems larger within Lake Cook E2 than outside. and lower (4.5 m) limits for the surface lowering are used to determine the extent of the lake area (Fig. 4c). Another source of uncertainty affecting our estimate of the volume loss is the error in the differential DEM. A first way to compute this DEM error is by considering the standard deviation of the elevation differences between each DEM and ICESat data (1.6 m for SPOT5, 11.7 m for ASTER). The error on the elevation change over the lake area is then obtained by dividing 11.8 m (the square root of the sum of the squares of 1.6 and 11.7 m) by the square root of the number of independent pixels in the Lake Cook E2 area using an auto-correlation length of 1 km (Bretherton et al., 1999;Nuth and Kääb, 2011). The correlation length cannot be estimated specifically for this study because of the lack of stable terrain but 1 km is a conservative value because shorter correlation lengths have previously been used with similar data set (0.5 km in Berthier et al., 2010; 0.2 km in Howat et al., 2008). The resulting elevation change error is ±0.22 m. However, this formal error does not take into account the spatial variations of the biases in the DEMs. A larger error (±5 m) was obtained for the northern Antarctic Peninsula by using a null test, i.e. examining elevation differences between ASTER and SPOT5 DEMs acquired the same day (Berthier et al., 2012). This more conservative error estimate is preferred here because it includes spatially varying biases between the DEMs and better reflects the uncertainties due to the manual interpolation of data voids in the ASTER DEM. This ±5 m error on the elevation change, multiplied by the area of the lake, leads to a total DEM error of ±1.10 km 3 of water. These "lake-area" and DEM errors cannot be considered independent and, hence, the total uncertainty is the arithmetic sum of the two, leading to a total uncertainty of 1.20 km 3 . This volume of 4.52 ± 1.20 km 3 represents the net volume loss between early 2006 and early 2012 and thus, underestimates the total volume of water discharged during the event itself given that the ice-sheet surface rose between March 2009 and February 2012 (Fig. 5). Using the elevation differences between ICESat laser period 2E and the SPOT5 DEM (Fig. 4c), we estimate that the volume change was +0. latter value is difficult to estimate, so a conservative 50 % error is assumed. Thus, in total, the discharge is estimated to amount to 4.52 ± 1.20 + 0.64 ± 0.32 = 5.16 ± 1.52 km 3 .

Subsequent cascade of lakes
The map of elevation changes for the Wilkes Land derived from Envisat radar altimetry features a complex pattern downstream of Lake Cook E2 , with spots of very localized uplift in a region of overall negative elevation trend (Fig. 6).
Most active lakes can be identified on this map due to their positive elevation trend because the surrounding area is either stable or slightly losing elevation. The positive elevation trend does not necessarily imply that the lakes permanently stored water, as can be seen on the series of Lake Cook E1 (Fig. 6b). Lake Cook E2 drainage and the subsequent passage of water through downstream lakes happened during the second half of the Envisat observation period so that a trend fitted through the full time series is positive, even though the surface elevation returned to the pre-drainage level after the passage of water.
Starting from Lake Cook E2 , active lakes are aligned along a curved line first oriented northward before bending to the north-west approximately at the middle of the path. This apparent flowpath, linking regions of anomalous surface elevation changes, is in good agreement with the likely flowpath for water draining from Lake Cook E2 inferred from BEDMAP2 (Fig. 7). The start date of the filling of those active lakes, as detected by our lake identification algorithm, is also in good agreement with the natural flow direction of subglacial water. The start date increases as hydrologic potential decreases, that is, as the distance to the coast decreases.
As described in Sect. 3.4, we computed a correction for surface elevation change due to regional processes (mostly firn volume change) based on "non-lake" points. Uncorrected times series of lake volume evolution are shown in Fig. 8a while the corrected ones are in Fig. 8b. The uncorrected volume stored in active lakes downstream of Lake Cook E2 (Fig. 8a) starts to vary before the drainage with a small annual cycle, some inter-annual variations and a minimum volume during 2004. The corrected time series (Fig. 8b) are almost flat before the start of the drainage of Lake Cook E2 which gives good confidence in the validity of the correction.
The downstream flood clearly starts in early 2007 and the flood volume grows till the end of 2008. We chose to focus on the correlation threshold of 0.8, using the 0.6 and 0.9 thresholds to estimate the sensitivity of our cumulative water volume to this variable. The peak volume, up to 3.7 km 3 , is reached towards the end of 2008. This volume varies between 4.3 and 3.0 km 3 when the alternative correlation thresholds (0.6 and 0.9) are used instead (Fig. 8). After the peak, between late 2008 and late 2011, the volume of water stored in the downstream lakes seems to remain stable. This suggests that water is permanently stored in the downstream lake system.

Usefulness of DEM derived from stereo imagery in central Antarctica
A nearly complete surface topography of the Lake Cook E2 area was derived from a SPOT5 stereo pair with a precision greater than ±2 m at the 1-sigma level after adjustment against ICESat data. This is a promising result for the central plateau of the East Antarctic ice sheet where DEMs derived from stereo imagery have been of little use up to now and thus have often been replaced by photoclinometry (e.g. Bindschadler et al., 2011). Because of their limited coverage and the associated cost, DEMs from SPOT5 or other sensors with similar stereo capabilities cannot currently replace frequent passages of the laser and radar altimeters to continuously monitor elevation changes of the whole ice sheet but can be helpful to better describe some discrete events, such as large subglacial lake drainages. The DEM and the associated imagery revealed that the lake shape is not elliptical, which explains most of the difference with the estimate based on ICESat measurements alone (Smith et al., 2009). Our estimate of the volume of water that drained from Lake Cook E2 using DEM differencing can be compared to a recent alternative estimate based on ICESat and CryoSat-2 data (McMillan et al., 2013). In their study, the volume difference computed between preand post-drainage DEMs is scaled with ICESat measurements to estimate the lake volume during each ICESat campaign. The pre-event topography is obtained by gridding the sparse ICESat data (assuming a smooth surface), whereas the post-event DEM is measured by CryoSat-2. They estimated that 6.4 km 3 of water drained from Lake Cook E2 . This alternative study does not provide an error bar. Their volume is higher than ours but falls within the bounds of our estimate of 5.16 ± 1.52 km 3 . Their lake area (260 km 2 , no error bar) is also similar to ours (between 192 and 267 km 2 ).
The DEM differencing technique presented here makes it possible to observe elevation changes with a vertical precision of about 1 to 2 m when two SPOT5 DEM are compared (e.g. Berthier et al., 2012). Subglacial lake drainages resulting in surface elevation changes of smaller amplitudes than Lake Cook E2 may thus be studied using this technique. The main limitation is the need to find suitable satellite stereo pairs acquired before the event. Large-scale acquisition programs dedicated to glaciers and ice sheets, such as SPIRIT (Korona et al., 2009) or GLIMS (Raup et al., 2007), are instrumental to build such historical archives of images.

Inability of classic radar altimetry to capture the drainage of Lake Cook E2
The small footprint of ICESat laser altimeter (around 70 m) is better suited than the ∼ 10 km footprint of Envisat radar altimeter for observing the drainage of Lake Cook E2 . The The  Lakes reported by Wright and Siegert (2012) are plotted as red circles and labelled. The black dashed line is the search area, i.e. the envelope of all possible water flowpaths draining out from Lake Cook E2 (Sect. 3.2). Hydrologic potential contours computed from BEDMAP2 data are superimposed every 500 000 kg m −1 s −2 (thin black lines). The drainage basin of Lake Cook E2 (5 % probability to reach the lake; see Sect. 5.5) is plotted in dashed blue line, and the lake contour is in red. All active lakes (correlation > 0.9, see Sect. 3.4) are reported along with a colour code indicating the start date of the flood at each lake point. Light grey lines show the ground tracks of Envisat 35-day repeat cycle.
problem had already been noticed by Fricker et al. (2010) who made similar findings in lakes on the Siple Coast. The size of Lake Cook E2 being commensurable with the diameter of the Envisat altimeter footprint, the radar altimeter is not able to retrieve any echo from the bottom of the trough left by the lake drainage. The difference in elevation shown in Fig. 9 is in fact a misinterpreted difference in radar range: as the lake trough deepens, the point sending back the first echo to the satellite is not at nadir any more but has shifted towards the rim of the depression, which yields a longer range. The minimum distance from a nadir point at the centre of the depression to the lake margin is about 5 km (e.g. for the time series in Fig. 9). For a 5 km shift and a satellite altitude of 800 km, the range difference is 15.6 m. This is in remarkable agreement with the apparent drawdown of about 15 m observed in Fig. 9, especially if we take into account that the echo is expected to be strongly distorted over such nonflat topography. Because of this limitation, inherent to classic radar altimetry, no accurate estimate of the lake drainage can be done from Envisat measurements. Our along-track processing only fits a quadratic model to all measurements to account for the surface topography and cannot correct the large topographic changes due to the lake drainage. It explains the large noise level (1.4 m) in the time series (Fig. 9). Downstream of Lake Cook E2 , elevation changes are less abrupt and can be monitored using Envisat. Envisat coverage is temporally denser than ICESat (regular 35-day cycle for Envisat vs. 2-3 measurements a year with ICESat). This compensates for the slightly higher noise in the time series of the radar altimeter. ICESat has a single point accuracy of 0.1-0.15 m (Shuman et al., 2006;Brenner et al., 2007), whereas noise in the Envisat series after along-track processing and in 682 T. Flament et al.: Cascading water underneath Wilkes Land, East Antarctic ice sheet Fig. 8. Chronology of the drainage of Lake Cook E2 and the subsequent storage/release of water in downstream lakes. (a) uncorrected volume of water downstream of Lake Cook E2 from elevation change at lake points. (b) Same volume but corrected for the regional elevation change measured at non-lake points (see Sect. 3.4). In both panels, curve "from Cook E2 " is the volume of water that drained from Lake Cook E2 , estimated from DEM differencing and ICESat measurements; 0.6-0.9 curves: volume of water stored in active lakes downstream of Lake Cook E2 , These curves are labelled according to the different correlation thresholds used for detecting the lakes (see Sect. 3.4).

Fig. 9.
Relative elevation time series from Envisat radar altimetry over Lake Cook E2 at two adjacent locations along track. The two time series are taken at the same location as series in Fig. 4; only the processing differs. Here, all data are used in the processing, whereas in Fig. 3 only the data acquired before the start of the drainage were used. Note the larger noise level here. the area of interest here is about 0.2 m. Thanks to the higher repetition rate of Envisat, events large enough to stand out of the noise can be more precisely dated.

Ability to close the subglacial water budget
The water budget, the difference between the volume of water that drained out of Lake Cook E2 and the volume of water that transited through the downstream lakes is closed within errors bars. Yet the error on the volume stored downstream is high and there is a high sensitivity of our estimate to the correlation threshold of the fitted Gaussian curve. Several sources contribute to the error: i. The "lake signal" is mixed up with elevation changes from other phenomena acting on the ice-sheet surface. Snowfall variability produces elevation change signals on interannual timescales (Magand et al., 2007;Lenaerts et al., 2012), comparable to the lake drainage timescales. We proposed a correction of the firn volume change based on volume change at "non-lake" points. This correction has the advantage of being selfconsistent with the rest of the Envisat data set. We preferred this correction to another one based on firn modelling because both data sets (radar altimetry and modelled firn densification) could not be cross calibrated yet (Ligtenberg et al., 2012).
ii. Our estimate is dependent on the validity of the area, defined as the envelope of all water flowpaths out of Lake Cook E2 , which is scanned to detect all active lakes. In a sensitivity test and to account for possible error in the flowpaths inferred from BEDMAP2, we searched for lakes within a 40 km buffer zone outside of the flowpaths boundary. Very few points were identified as lakes and the overall water volume remained virtually unchanged.
iii. It is possible that part of the water flowed rapidly through well-carved channels and did not leave a persistent imprint on the surface.
iv. Faint signals could go undetected. This could be a large contribution if water was stored as a thin layer over a large area.
v. Finally, some water may have reached the coast towards mid-2008 and exited our investigation area before Lake Cook E2 completely drained. That could explain why the budget is not closed after mid-2008.
However, at the peak date more than 70 % of the volume drained from Lake Cook E2 is found in the downstream active lakes which implies that most of the water flowing underneath the ice sheet had a direct impact on the surface elevation, and that this impact was large enough to be detected by Envisat radar altimetry. We also infer that most water released by Lake Cook E2 remained trapped in downstream lakes. However, this interpretation is subject to caution due to the high level of uncertainty of lake detection in the lower reaches of the area, where ice flow is faster and surface slope steeper.

Comparison to other known drainage events
The event we studied here is not unique but it is the largest directly observed drainage by volume and distance travelled. The drainage of Lake Cook E2 can be compared to the Adventure Trench subglacial drainage reported by Wingham et al. (2006) and further investigated by Carter et al. (2009). The total water volume drained in the Adventure Trench event was estimated between 1.8 and 2.2 km 3 , but surface elevation change was only in the range of a few metres. The distance between the two most remote connected lakes was about 290 km, slightly smaller than the distance of 380 km we found between Lake Cook E2 and its most remote connected subglacial lake. With the increasingly dense spatial coverage provided by space missions (altimetry, SAR and optical imagery) we can expect that more of such events will be discovered.
In other places of the Wilkes Subglacial Basin, the bedrock was marked by a palaeo-flood of much larger amplitude (Jordan et al., 2010). Analysing subglacial landforms, Jordan et al. (2010) estimated that a lake containing 850 km 3 of water drained. This supposed lake was located slightly south-west of Lake Cook E2 and the outburst flood probably flowed through the deep bedrock channel where active lakes Cook W1 and Cook W2 lie (Fig. 2c). This channel is parallel to and separated from the flowpaths showed in Fig. 7. Recently, Livingstone et al. (2013) computed the hydrologic potential from BEDMAP2 and concluded that many lakes are still to be discovered, mostly because of the sparse coverage of the bedrock by radio echo sounding. They also proposed that major lakes should drain when an ice sheet retreats because the surface slope locally increases (see also Scambos et al., 2011). Here, we were not able to measure any change in surface elevation or any increase in slope over Lake Cook E2 using the Envisat pre-drainage time series (Fig. 3) which could explain the triggering of lake drainage. The precision and duration of our data set is probably not sufficient for this task.

Possible causes for the uplift of the ice surface since 2009
As shown earlier (Sect. 4.1 and Fig. 5), the centre of the Lake Cook E2 depression has risen by about 13 ± 1.6 m between the last ICESat measurement in late 2009 and the acquisition of the DEM in February 2012. Several processes can potentially have contributed to the surface rising: water refilling the subglacial cavity, accumulation of windblown snow or local ice-flow convergence. By the same method used to estimate the downstream subglacial flowpaths from Lake Cook E2 , we can determine the flowpaths converging toward Lake Cook E2 . Starting from a pixel at the rim of the lake depression, we systematically explore the lake surroundings by drawing a hundred flowpaths going down the hydrologic potential. If more than a given percent of these flow paths reach the lake, the pixel is considered part of the drainage basin of Lake Cook E2 and all of its eight neighbours are added to the list of pixels to be tested. The exploration stops when no pixel is left in the list. We consider the most likely drainage basin as the zone with a more than 50 % chance to reach the lake (530 km 2 ) plus the lake area (220 km 2 ). The areas with 95 and 5 % chance to reach the lake give our uncertainty margin, respectively, 450 and 1030 km 2 .
The resulting water collection area is thus about 750 ± 300 km 2 . The small size of this collecting basin can be explained by the proximity of the ridge leading to Talos Dome. The observed refilling of 0.64 ± 0.32 km 3 in 3 years would imply a basal melt rate between 0.3 m a −1 (refilling of 0.32 km 3 and a basin of 1050 km 2 ) and 2 m a −1 (refilling of 0.96 km 3 and a basin of 450 km 2 ). These values are at least two orders of magnitude higher than expected basal melt rates in this region (Pattyn, 2010).
The hypothesis of the depression filling by windblown snow is difficult to test. One could expect that snow transport close to a ridge is small because the katabatic winds are unlikely to be strong there. However, the large topographic anomaly left by the drainage could perturb nearsurface winds and contribute to the accumulation of snow. We cannot formally exclude this hypothesis.
The last possibility, filling of the trough by local ice flow convergence (e.g. Aðalgeirsdóttir et al., 2000), could be the dominant phenomenon. The depression "refills" at a rate of  Scambos, personal communication, 2013). This subset corresponds to the dashed line box in Fig. 1b. The white box is the zone covered by Fig. 4. The white thin contour is the −3 m contour extracted from the map of elevation difference from DEMs, our preferred estimate of the lake area (Fig. 4). 0.2 km 3 a −1 . Given its 70-km-long rim and a 2700-m-thick ice column, this ice flux corresponds to a 1 m a −1 change in average ice flow velocity. This is only 20 % of the velocity, 5 ± 5 m a −1 at this location (Fig. 2a, Rignot et al., 2011).
Our observation period is too short to infer anything about a cyclic behaviour of Lake Cook E2 . However, the MODIS Mosaic Of Antarctica, acquired over austral summer 2003-2004, i.e. before the start of the drainage, reveals a faint feature on the ice-sheet surface. We interpret this feature as a modulation of shading induced by an undulation of the surface and made visible by the low Sun. This pre-drainage pattern of shading is similar to the lake contours (Fig. 10a). This surface feature could have been created by the accumulation of water preceding the drainage (locally relaxing basal shear stress), or by a previous drainage.

Conclusion
A major drainage event was detected by radar and laser altimetry and analysed in detail with DEMs derived from stereo imagery in the Wilkes Land, East Antarctica. This dramatic drainage lasted 2 years, occurred under 2700 m of ice, led to a surface lowering of up to 70 m, affected a total area of about 220 km 2 and released over 5 km 3 of water in the subglacial hydrological network. The total volume drained is larger than reported in earlier similar events underneath the Antarctic ice sheet. However, the water discharge flux remained 1 to 2 orders of magnitude smaller than during jökulhlaups in Iceland (Björnsson, 2003). Surface features were present on optical imagery prior to this event. This is another hint that this region is hydrologically active and suggests that Lake Cook E2 could be subject to recurring drainage. It will be of great interest to monitor the evolution of this deep depression in the future using ongoing (CryoSat-2) or future altimetry missions (SARAL/AltiKa on the same orbit as Envisat or ICESat-2, Sentinel-3 on different orbits).
The study of such local and rapid events is only possible with both sufficient temporal and spatial resolution and continuous coverage.
ICESat laser altimetry yields better observations of the drainage at Lake Cook E2 than Envisat radar altimeter. The latter, because of its large footprint, is not able to observe the bottom of the trough left after the water has flowed away. Still, the absence of cross-track coverage from ICESat precludes an accurate computation of the volume drained from the lake. The dense temporal sampling of Envisat allows a more precise constraint for the onset of lake drainage. Envisat radar altimetry also revealed a complex pattern of elevation changes downstream of Lake Cook E2 and we interpret them as the surface expression of transient storage of water in a succession of subglacial lakes. We suggest that most of the water released during the drainage of Lake Cook E2 remained trapped under the ice sheet and did not reach the ocean.
Outlet glaciers feeding the Cook Ice Shelf are flowing over the deep Wilkes subglacial basin, where much of the bedrock lies below sea level (Fretwell et al., 2013). Marine ice sheets are suspected to be the most prone to collapse in a warming climate (e.g. Weertman, 1974;Durand et al., 2009). Future studies should aim to determine whether this large amount of water influenced the flow of ice, in particular for the eastern tributary of Cook Ice Shelf, and modify the basal melt rates of this ice shelf (Jenkins, 2011).