Radar Sounding Confirms a Hydrologically Active Deep-Water Subglacial Lake in East Antarctica

Lake CookE2, upstream of Cook Glacier in East Antarctica, is an “active” subglacial lake that experiences episodic discharge and recharge of basal water. Although around 130 active lakes are known to exist, the majority are not able to be identified by ice-sounding radar techniques, suggesting they are ephemeral and/or distributed stores of small amounts of water rather than permanent significant singular features. However, airborne radar data from Lake CookE2 reveal a bright and flat ice-bed interface, providing clear evidence of deep (>10 m) water surrounded by elevated topography. The data show the lake area is ∼46 km2; three times less than a previous estimate (145 km2) from Ice, Cloud and land Elevation Satellite (ICESat) satellite altimetry, suggesting a bias in identifying subglacial lake area from surface depressions. Using time-series altimetry from ICESat, Cryosat-2, and the Reference Elevation Model of Antarctica, we re-estimate the lake discharged ∼2.73 km3 of water (or ∼59.6 m in lake level) between February 2006 and October 2008. Subsequently, the ice surface over the lake rose steadily and experienced a mean uplift of ∼9 m between January 2011 and November 2016, indicating continuous recharge with total volume increase of ∼0.42 km3. The lake is recharging at a rate of ∼1.1 m/year, which means it could take another ∼39 years to reach the lake level that triggered the previous discharge.

INTRODUCTION Lake CookE2, located at the head of Cook Glacier in Wilkes Land (Figure 1), is an "active" Antarctic subglacial lake. Evidence for this is drawn from the significant surface depression that occurred during 2006(Smith et al., 2009, interpreted as a sudden loss of basal water. The ice surface lowered by ∼70 m in just 2 years (Mcmillan et al., 2013), due to the largest discharge of basal water measured in Antarctica thus far. Surface depressions with such an amplitude are easily captured by satellite altimetry and/or InSAR techniques (Flament et al., 2014). So far, 129 active subglacial lakes have been identified in Antarctica according to this interpretation, and they featured in the latest inventory of Antarctic subglacial lakes (Wright and Siegert, 2012). Lake CookE2 was discovered from ice-surface altimetry more than a decade ago by Smith et al. (2009), but there has been no measurement of the ice-water interface to confirm it as a deepwater lake as opposed to an ephemeral store of distributed water. Most "active" subglacial lakes are not associated with specular radar reflections that characterize deep-water lakes; many exist as a consequence of a complex interaction between ice flow dynamics over subglacial hills, and the temporary build-up of water on their lee sides (Siegert et al., 2014;Wright et al., 2014), resulting in a spatially variable region of radar reflectivity at the ice sheet bed. According to Ice, Cloud, and Elevation Satellite (ICESat) measurements, the area above Lake CookE2 that experiences ice-surface elevation change is ∼145 km 2 determined by Smith et al. (2009). Within this zone, the ice-surface slope can reach ∼2 • , suggesting significant hydrological forcing by the ice overburden. Indeed such surface slope is enough to drive water up a hill with a slope of over 20 • . Subglacial discharge in active lakes can reduce lubrication between ice base and bed materials, causing faster ice flow toward the margin (Stearns et al., 2008;Paolo et al., 2015). Ice surface elevation change on the order of a few meters have been shown to modify basal water pathways, leading to the spatial redistribution of water downstream (Wright et al., 2008). Hence, substantial discharges on the order estimated for Lake CookE2 could influence both the pathway of basal water and the flow of ice above.
Airborne radar sounding is a valuable means of mapping topography and conditions at the base of the ice sheet. For example, radar echoes from subglacial water bodies >10 m in depth are commonly characterized by bright, flat and specular features (Siegert and Ridley, 1998;Carter et al., 2007;Christianson et al., 2012). Radar echoes also allow for the detection of subglacial sediments (Siegert, 2000), basal water channels (Jeofry et al., 2018) and distinct geomorphological features (Rose et al., 2015). Here, we use Multichannel Coherent Radar Depth Sounder 5 (MCoRDS5) radar measurements collected during a NASA Operation IceBridge mission in November 2017 to better constrain the extent of Lake CookeE2. By comparing the SwSARIn data from Cryosat-2 and the highresolution elevation model, we derive long-term time-series of ice-surface height change for the lake between 2011 and 2016.

MATERIALS AND METHODS
In brief, we use data from the MCoRDS5 radar to examine the ice-bed interface over the Lake CookE2 region, to first determine whether the lake exists and then to measure its topographic setting. We also combine a new ice-surface elevation model, REMA (see below), with ICESat and Cryosat-2 SwSARIn data to generate a long-term time-series of ice-surface height change from 2003 to 2016, to understand how the lake has behaved in recent years.
Multichannel Coherent Radar Depth Sounder 5 airborne radar data were acquired over the Lake CookE2 on November 29, 2017 during an IceBridge mission. The radar frequency over this region is 300 MHz. Three transects, 009, 010, and 011 (black and red lines in Figure 1), were obtained over Lake CookE2. We use the MCoRDS5 Level 1B data from Centre for Remote Sensing of Ice Sheets (CReSIS), University of Kansas 1 , which include necessary geophysical parameters such as geodetic coordinates, range travel-time, aircraft elevation, surface interface travel-time and bottom interface travel-time. Radar echoes from a clear icewater interface typically show bright and flat returns that are distinct compared with echoes from an ice-rock interface that undulate and have variable but generally weaker strengths. This distinction has been used to identify pooled subglacial water of at least 10 m in depth across the ice sheet (Robin, 1975;Gorman and Siegert, 1999;Ilisei et al., 2018).
Antarctic surface elevation models have been generated from recent satellite altimetry missions (DiMarzio, 2007;Helm et al., 2014;Zwally et al., 2014), with a resolution of 500 m for most regions. Recently, a new elevation model, the Reference Elevation Model of Antarctic (REMA), was generated by stereo-photogrammetry based on sub-meter resolution optical satellite imagery (Howat et al., 2019). The spatial resolution of this elevation model used here is 8 m and the elevation error is within 1 m.
Ice-surface height measurements from two satellite altimetry missions, ICESat (between February 2003 and October 2009) and Cryosat-2 (between November 2010 and December 2016), are combined and then differenced with REMA to construct a 13-year record of topography-free ice-surface height change. ICESat was the first satellite laser altimetry mission to measure surface elevation over land and seas, with a precision of ∼15 cm (Siegfried et al., 2011). Cryosat-2 is a radar altimetry satellite, launched in April 2010. The main aim of Cryosat-2 is to measure changes in ice-surface height and sea ice thickness over the Polar regions. An onboard altimeter operates in three modes: low resolution, SAR and SARIn. The SARIn mode covers icesheet margins, providing interferometric measurements over steeply sloping regions. A new SARIn elevation product, named SwSARIn data, was generated recently (Gourmelen et al., 2017a) to provide denser elevation measurements by exploiting full radar waveforms (Foresta et al., 2016;Gourmelen et al., 2017b).
By combining an accurate reference elevation model and dense satellite altimetry measurements, we can construct a longterm time-series of ice-surface height change over Lake CookE2. We use an adapted method based on previous studies (Moholdt et al., 2010;Siegfried and Fricker, 2018). For ICESat data, we collect all available elevation data from 17 campaigns within the lake outline. The REMA is interpolated at each ICESat point to obtain a height difference. Then we filter out potential outliers in height differences based on the three-sigma rule, as follows: (1) we calculate the mean and standard deviation of remaining height differences at the beginning (according to the rule, height differences that exceed the triple standard deviation from the mean of height differences are assumed as outliers and then removed); (2) we calculate the mean and standard deviation again using filtered height differences (this step is iteratively implemented until no outlier is detected); and (3) the mean of the final height differences is used to represent the mean height difference with respect to REMA. Measurements in each campaign are processed following the same workflow.
For Cryosat-2 measurements, we construct monthly icesurface height differences in four steps as follows. First, we confine our study area within a new outline for Lake CookE2 (based on ice-penetrating radar -see below; red dashed line in Figure 1), to generate an ice-surface height change time-series for the lake. SwSARIn measurements within the lake outline are collected for each month. Second, we collect all elevation measurements for each month. If the number of measurements is less than an empirical value of ∼8000 in a particular month, the measurements are discarded and further processing is not be conducted. This is because limited numbers of measurements are likely to introduce a bias in the ice-surface change derivation, especially if data filtering is applied. After initial examination, the REMA is interpolated at each SwSARIn point collected during a month period to obtain the "topography-free" icesurface height change with respect to the reference elevation model. Three-sigma rule is iteratively used to filter out unreliable height differences in each group. The mean of the filtered height differences is used to represent height changes relative to the elevation model over the lake area. The process is iteratively implemented for each month. Third, we use a −1.5 m (Cryosat-2 minus ICESat) (Mcmillan et al., 2013) estimate to account for elevation bias between the two satellite altimetry missions. Fourth, we combine height differences with REMA between 2003 and 2016 to construct a 13-year record of ice-surface height change over Lake CookE2. Uncertainties in ice-surface change from ICESat and Cryosat-2 are also estimated. The elevation error of REMA is typically less than 1 m (Howat et al., 2019). The accuracy of ICESat measurements over ice sheet is within ±14 cm (Shuman et al., 2006). The accuracy of SwSARIn measurements is estimated to be ∼2 m. Therefore, the accuracy of ice surface change from ICESat and Cryosat-2 is ∼1 m and ∼2.2 m, respectively.
We estimate the long-term change in ice-surface height using Cryosat-2 data by a quadratic model as follows: h = a 1 x + a 2 y + a 3 x 2 + a 4 y 2 + a 5 xy + a 6 t + a 7 where h is the ice-surface height measurement, t is the measurement time, a 1 ∼ a 5 are coefficients of the terrain model, and a 6 and a 7 denote linear model coefficients of ice-surface change during study period. x and y are coordinates in the Polar Stereographic South projection.

Airborne Ice Penetrating Radar
We present three radar transects across Lake CookE2 (Figure 2). In each, a bright and smooth bed feature can be observed, typical of an ice-water interface, surrounded by undulating returns and overlain by ∼1700 m of ice. The unbroken nature of the basal reflector indicates that the water must be at least 10 m deep (Gorman and Siegert, 1999), confirming Lake CookE2 as a deep-water subglacial lake. Indeed the lake depth is likely much deeper than 10 m when one considers the steep side-wall slopes bordering the lake, if we extrapolate them toward the lake center. Based on this simple approach, it is possible that the lake is several tens, if not hundreds, of meters deep. Several active subglacial lakes identified by satellite altimetry show no clear evidence of being classic deep-water lakes with bright and smooth radar returns at ice-bed interface (Siegert et al., 2014;Wright et al., 2014;Humbert et al., 2018). These should not be considered as pooled subglacial water bounded by steep-sided topography; rather they occur through the complex interplay between ice flow over topography and the hydrological potential that is created. Pooled water can accumulate on the lee sides of major obstacles, due to a hydropotential minimum established there but, as the water builds and the surface elevation increases, this minimum is reduced. Through this effect, once the hydropotential minimum is reduced to a critical level a discharge can occur. In the case of Lake CookE2, however, we observe a rare example of an active subglacial lake that appears similar to those at the ice-sheet center; i.e., a deep subglacial lake bounded by topography.
In the radar transects, only the middle of the uplifted area contains a bright and smooth base, revealing that the lake outline as identified formerly by satellite altimetry alone is significantly wider than the actual lake beneath. During the discharge, ice directly below the lake lowered significantly due to the loss of basal water on the order of tens of meters. Forced by gravity, the surrounding ice would then flow toward the lake interior, resulting in a depressed surface surrounding the actual lake. Therefore, ice surface lowering is not only caused by collapse of ice immediately above the lake, but also by increased ice flow into the former lake area.
To identify the outline of Lake CookeE2, Smith et al. (2009) used elevation-change anomalies greater than 0.1 m to delineate the lake outline, allowing an estimate of ∼145 km 2 . By differencing pre-drainage and post-drainage surface elevation change, Mcmillan et al. (2013) estimated the lake area to be ∼260 km 2 , including a hook-shaped area adjacent to the main body of the lake (southwest blue area in Figure 1). We see no radar echoes characteristic of deep-water in this hook-shaped zone, hence we think this area should not be assumed as part of the main lake. Using elevation change thresholds of 1.5, 3, and 4.5 m, Flament et al. (2014) derived areas of 267 km 2 , 219 km 2 , and 192 km 2 , respectively. Here, we constrain the lake outline through the radar echoes. From the radar records, the lake's length and width is ∼12.2 km and ∼4.1 km (based on transects 010 and 011, respectively). We estimate the lake's area to be ∼46 km 2 . Previous studies estimate the lake area to be 149∼267 km 2 ; around 3-6 times the value we derive.

Temporal and Spatial Characterization of Ice-Surface Height Change
Using ICESat and Cryosat-2 data, we construct a time-series of ice-surface height change over Lake CookE2 between September 2003 and November 2016 (Figure 3a). From February 2006 to October 2008, the ice surface lowered ∼59.6 m, indicating (assuming conservation of volume) ∼2.73 km 3 of water was discharged from the lake. Between January 2011 and November 2016, the ice surface rose steadily at a rate of ∼1.1 m/year. The ice-surface began to rise in early 2011, ∼2 years after a drainage event had finished, and continued to do so in the subsequent 6 years. If we assume this rate of recharge is maintained in future, it would take another ∼39 years for the lake to reach the level associated with its previous discharge. Using the lake area derived by ice-penetrating radar, we estimate that the volume of water in CookE2 increased by ∼0.42 km 3 in this time. At the same time, we find significant annual variability in ice-surface heights; lowest around January and greatest around May-July in each year.
Previous research has confirmed an association between the annual change in ice-surface heights, observed by satellite altimetry, with surface accumulation of snow in Antarctica based on firn-densification modeling (Ligtenberg et al., 2012). Although the temporal variations in surface accumulation described by Ligtenberg et al. (2012) and our annual oscillations in surface elevation change match well, the respective amplitudes are much larger in our study; on the orders of meters over Lake CookE2. As described in section "Materials and Methods, " uncertainty in Cryosat-derived ice surface change is ∼2.2 m whereas uncertainty in the annual variability determined by Ligtenberg et al. (2012) is <3 cm. Although it is possible that our measurements relate to annual surface accumulation changes over the surface depression above the lake, we acknowledge that further work is needed to determine whether this is the case.

Basal Water Flow Sensitivity to Ice-Surface Height Change
Using the latest ice surface and bed data from BedMachine Antarctica, we calculate the hydraulic potential (Figures 4A,B) for Lake CookE2 and its locale (Morlighem et al., 2019). Lake CookE2 lost a great deal of water during the drainage event between 2006 and 2008, leading to a widespread ice surface depression ( Figure 4C). As surface elevation dictates the flow of water by 10x the basal slope, such change is likely to lead to modification of basal water pathways. To reveal the response and sensitivity of basal hydrological networks to the observed ice-surface changes, we calculate basal flow paths by a subglacial hydraulic potential method, which assumes that the water pressure is equal to the ice overburden pressure (Shreve, 1972;Livingstone et al., 2013). A previous study found that basal flow paths are sensitive to surface height changes of just a few FIGURE 2 | Three airborne radar transects across Lake CookE2. They are sequentially numbered by 20171129_04_009 (a), 20171129_04_010 (b), and 20171129_04_011 (c). In each image, red lines denote the tracked ice-bed interface out of lake outline determined by ICESat; blue lines denote ice-bed interface within lake outline determined by ICESat; and white lines denote lake length identified by radar and the ICESat repeat-pass method (Smith et al., 2009). For clarity, the bed interpretation is offset 50 m vertically from the basal reflection.  meters at certain sites in the ice sheet (Wright et al., 2008). We compare the basal routing before and after the discharge event ( Figure 4B); calculated by BedMachine bed data, surface data and ice-surface DEM generated by ICESat (DiMarzio, 2007). Widespread ice-surface uplift (Figure 3b) led to modification to the basal flow paths, confirming that they are very sensitive to ice-surface height change (Figures 4B,C).
Before the discharge, Lake CookE2 contributed water to downstream hydrological networks. When the lake began to refill after the drainage, basal flow paths on the east side of the lake combined with a tributary downstream of Lake CookE2 (Figure 4B). A comparison between these two states shows that the water flow changed as a consequence of the discharge event ( Figure 4D). Hence, subglacial lake dynamics appear to have an influence on the basal hydrology outside of the immediate lake boundary.

Source of Water in Lake CookE2
Given that Lake CookE2 is receiving substantial levels of water over long periods, we can ask where the water originates. To frame the problem we calculate the water storage balance in Lake CookE2. From 2011 to 2016, water in the lake increased at a rate of ∼0.052 ± 0.008 km 3 /year. Conceptually, water from a variety of sources can contribute to water increase in Lake CookE2, such as basal melt, groundwater in till layers and hydrologic reservoirs. The mean basal melt rate over Lake CookE2 is ∼6-7 mm/year (from ice-sheet modeling) potentially contributing ∼0.005 km 3 /year of water to the lake (Pattyn, 2010). If we assume the water is all produced by basal melting at this rate, the upstream catchment would require an area of ∼186 km 2 . This area is only ∼3 times larger than the lake's, indicating Lake CookE2 is close to the upper regions of a hydraulic catchment. We note recent work that suggests groundwater in basal till layers can contribute to subglacial hydrology networks (Christoffersen et al., 2014;Siegert et al., 2018), meaning that the lake could draw water from its immediate surroundings. If groundwater made a contribution to the lake, this would reduce the area of melting upstream necessary to fill the lake.

CONCLUSION
We have confirmed, using ice-penetrating radar data, that Subglacial Lake CookE2 in East Antarctica (an active subglacial lake) is similar in its physiography to many stable deep-water lakes that are contained by steep-sided topography. Combining ICESat and Cryosat-2 SwSARIn data from 2003 to 2016, and the high resolution elevation model REMA, we constructed a longterm time-series of topography-free ice-surface height change. We found that the surface elevation over the lake kept rising at a rate of ∼1.1 m/year following a drainage event. We used this rate of recharge to calculate that it would take another ∼39 years for the lake to reach the level associated with its previous discharge. We also investigated the spatial distribution of ice-surface change over 2011-2016, confirming consistency with the new lake area outline defined by ice-penetrating radar, and that previous assessments of lake area from ICESat are overestimated by as much as six times. By comparing basal water flux before and after a drainage event, we showed that water flow is sensitive to ice-surface height change in Lake CookE2, and that water discharged by active subglacial lakes can perturb subglacial hydrology at a regional scale. In addition, we revealed that periodical filling in Lake CookE2 is associated with water input from an upstream hydrologic catchment that is up to three times the lake surface area.

DATA AVAILABILITY STATEMENT
ICESat DEM, elevation data and BedMachine data are available at National Snow and Ice Data Center; Cryosat-2 SwSARIn data are available at British Antarctic Survey; REMA data are available at Polar Geospatial Center in University of Minnesota; MCoRDS5 data are available at Center for Remote Sensing of Ice Sheets. The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://github.com/ yanli2020stone/Lake-CookE2-outline.

AUTHOR CONTRIBUTIONS
YLi and YLu designed the experiment and wrote the manuscript. MS improved the writing and provided advice on estimate of new lake area, ice-surface height change, and sensitivity of subglacial hydrology to ice-surface height changes. All authors contributed to the article and approved the submitted version.