Repeat Subglacial Lake Drainage and Filling Beneath Thwaites Glacier

Active subglacial lakes have been identified throughout Antarctica, offering a window into subglacial environments and their impact on ice sheet mass balance. Here we use high‐resolution altimetry measurements from 2010 to 2019 to show that a lake system under the Thwaites Glacier undertook a large episode of activity in 2017, only four years after the system underwent a substantial drainage event. Our observations suggest significant modifications of the drainage system between the two events, with 2017 experiencing greater upstream discharge, faster lake‐to‐lake connectivity, and the transfer of water within a closed system. Measured rates of lake recharge during the inter‐drainage period are 137% larger than modeled estimates, suggesting processes that drive subglacial meltwater production, such as geothermal heat flux or basal friction, are currently underestimated.


Introduction
The vast majority of ice in the Antarctic ice sheet drains from the continent to the ocean through fast-flowing ice streams and glaciers (Rignot et al., 2011). The presence of meltwater at the bed reduces basal stress, allowing the ice masses to sustain high velocities in some regions (Alley et al., 1986;Kamb, 2001). The movement of water has also been linked to transient glacier flow acceleration (Stearns et al., 2008) and to enhanced melt at the grounding line (Le Brocq et al., 2013;Wei et al., 2020). Therefore, the presence, location, and movement of water at the ice-bed interface are likely significant controls on the mass balance of Antarctica (Bell, 2008). The transport of water from upstream regions to downstream zones was once thought to be a steady-state process (Parizek et al., 2002); however, satellite observations indicate that the movement of subglacial water might be episodic (Gray, 2005;Wingham et al., 2006). Observations of localized height anomalies have been interpreted as subglacial water moving in and out of subglacial lakes causing a response at the surface of the glacier. Subglacial lakes located within the interior of the ice sheet are thought to be in a steady state with only localized impact on ice flow (Siegert et al., 2005), while lakes located in fast-flowing regions could temporarily alter Antarctic mass balance by modulating the amount and location of subglacial water through episodic drainage events (Siegfried & Fricker, 2018).
Active subglacial lakes have been identified throughout Antarctica with satellite altimetry and ice-penetrating radar (Smith et al., 2009;Wright & Siegert, 2012). Observations of surface elevation changes indicate that subglacial lakes are hydraulically connected (Fricker & Scambos, 2009;Wingham et al., 2006) and often exist in groups beneath Antarctic ice streams. During the ICESat-1 mission, operating from 2003 to 2010, no subglacial lakes were observed under the Amundsen Sea Sector of the Antarctic Ice Sheet, which was attributed to inadequate measurements due to cloudy conditions . Analysis of ice-penetrating radar identified a region of high specularity under the Thwaites Glacier, interpreted as evidence of the presence of water in distributed channels at the base of the ice sheet (Schroeder et al., 2014). In 2017, four large connected active subglacial lakes were discovered under the Thwaites Glacier from analysis of swath processed CryoSat-2 data, which indicated that the lakes drained simultaneously between June 2013 and January 2014 . Examination of the modeled subglacial melt production in the region suggested that the lakes should have a refill and drainage recurrence interval between 5 and 83 years, depending on whether the recharge scenario involved local melt production only or melt generated across the larger upstream catchments.
Here we use CryoSat-2 altimetry to produce elevation time series which extends the record of lake activity to mid-2019 and describes a second drainage event in 2017. The occurrence of two drainage events within a short timeframe allows us to explore the impact of drainage activity on the evolution of the subglacial system and to quantify subglacial melt supply providing rare insights into subglacial processes and basal melt generation.

Surface Elevation and Volume Change Estimates
Time-dependent elevations were generated using swath processing of CryoSat-2 level L1b SARin data acquired between 2010 and 2019. In contrast to the commonly used point of closest approach (POCA), swath processed SARin exploits the full radar waveform to resolve substantially more elevation than that of the POCA (Gourmelen et al., 2018;Gray et al., 2013;Hawley et al., 2009).
To determine the spatial behavior of subglacial lake activity, average rates of surface elevation change were computed using a plane-fitting algorithm (McMillan et al., 2014) applied to swath processed SARin data from late-2014 to mid-2019. Due to the dense elevation field provided by swath processing our region was gridded at a 500-m posting, with each cell incorporating a search radius of 1.5 km to lower map noise. Within each pixel time-dependent elevations were obtained by fitting a weighted hyperplane against easting, northing, and time, with a time-dependent coefficient retrieved from the regression representing the linear rate of surface change (Foresta et al., 2016). The model was fitted iteratively to the data, omitting elevations differing more than three standard deviations away from the model fit until no further outliers were detected. These maps were used to create masks encapsulating lake activity, which we define as a region with significant localized elevation change (>0.5 m yr −1 ) relative to the background signal (Flament et al., 2014;Fricker et al., 2007;Smith et al., 2017).
The temporal behavior of the lake system was determined through the creation of a surface elevation change timeseries between 2010 and 2019. We used an adapted version of the point-to-point method (see Text S1 for a detailed description) outlined in Gray et al. (2015Gray et al. ( , 2019 over our lake outlines to determine time-dependent elevations at a 45-day resolution, with elevations averaged over a 45-day search radius. To isolate the behavior of each lake relative to the catchment we removed the background thinning signal. This was achieved by deducing time-dependent elevations, as per the method above, using a 5 km exclusion area around each lake and subtracting it from the lake's signal. Note that, although the input dataset is similar, both the spatial and temporal approach followed here should lead to slightly different spatio-temporal smoothing compared against previous estimates of 2013 activity. ; see Text S2). Volume change through time was derived by integrating elevation change against the area of each lake mask. For any time-dependent volume change, an approximate statistical error is given by the average of the standard deviations divided by the square root of the number of timeseries realizations. Our volume estimates typically have a standard error of ±0.02 km 3 . This is a lower bound on our uncertainty, as the method of spatial and temporal sampling are likely to introduce additional uncertainty. In the absence of an alternative approach, we approximate the volume budget of subglacial water flux by the volume corresponding to the surface deflation, although we acknowledge that this assumption leads to additional uncertainty (Sergienko et al., 2007;Smith et al., 2017). Recharge rates were calculated by applying linear regression against volume change and time during the inter-drainage period, with the resulting rate representing the annual water supply to each lake. Our recharge rates have an uncertainty range derived by calculating a 95% confidence interval with the standard error of regression slope. These rates were compared against modeled local and total melt supplies ( Table 1 in Smith et al., 2017).

Hydraulic Potential Mapping and Estimating Subglacial Water Flow
To identify likely subglacial flow routes, and therefore determine the possible location of subglacial channels, we mapped hydraulic potential (see Text S3), forced using BedMachine ice thickness and bed elevation data, assuming water pressure everywhere at overburden (Morlighem et al., 2020). Closed depressions within our hydropotential map were filled to represent the large-scale basal flow pattern, as discussed by Smith et al. (2017). We applied a D8 routing scheme to our edited hydropotential grid to calculate the predicted motion of water throughout the glacier bed (Schwanghart & Scherler, 2014).
To gauge whether changing lake height might have an impact on the hydraulic gradients which drive water transport, background hydraulic gradients between the lakes were calculated by averaging hydraulic pressure change within each lake mask normalized by the distance between lakes. This was calculated from a hydraulic potential map without depressions filled. The background gradients represent potential gradients between each lake calculated from bed topography and ice thickness alone. During the drainage events, basal water pressure was assumed equal to ice overburden pressure, and we estimated the change in potential gradients between the lakes by normalizing the change in lake height against distance along predicted flow routes. It is worth noting that as the lake's water level rises or drops, basal water pressure is likely to exceed or be less than overburden, which could introduce additional uncertainty to our change in gradients.
We allowed for the possibility that subglacial melt generated by dissipation in the subglacial network could contribute to the basal water budget. An assumption of Röthlisberger (R-) channels was made, allowing for calculation of channel characteristics and melt production (Schoof, 2010), forced with average discharge rates (see Text S4).

Divergence Mapping
A divergence map was derived to determine the impact of ice flow divergence on surface elevation change during the inter-drainage period. Monthly mean ice surface velocity maps of Thwaites Glacier, gridded at 200 m, were derived from 6 and 12 day repeat pass Sentinel-1A and Sentinel-1B Synthetic Aperture Radar (SAR) data acquired in Interferometric Wide (IW) swath mode using offset tracking (Nagler et al., 2015). We created a velocity composite from January 2014 to March 2017 by averaging the monthly velocity maps over the same period. Maps with less than 50% coverage of the lake region were omitted from this calculation. A divergence map was produced according to Alley et al. (2018), forced using our velocity composite and BedMachine ice thickness data from October 2018 (Morlighem et al., 2020). The length scale used to produce the maps is adaptive and based upon ice thickness multiplied by a factor of eight.

2013 and 2017 Lake Drainage Activities
Our timeseries of surface elevation change over the Thwaites subglacial lake system (Figure 1) captures the 2013 activity discussed by Smith et al. (2017) and indicates that lake activity commenced in succession. It appears that the most upstream lake, Thw 170 , was first to activate in early April 2013, draining until early January 2014 with a total volume loss of 0.45 ± 0.03 km 3 . Second in the procession was Thw 124 , draining from mid-May 2013 until May 2014 with a total water loss of 3.83 ± 0.11 km 3 . Thw 142 activated from early July 2013 draining until October 2013 with an average volume loss of 0.55 ± 0.03 km 3 . Last in succession was Thw 70 which drained from mid-August 2013 until May 2014 with a total water loss of 0.90 ± 0.06 km 3 . Note that this succession is nearly identical to the one proposed by Smith et al. (2017) except for Thw 124 , which we find to drain later in the sequence. This disagreement is attributed to the temporal smoothing effect of the solution proposed by Smith et al. (2017). We observe a second episode of lake activity upstream of Thw 70 from early-2017, indicating a previously unobserved episode of lake drainage. Thw 142 and Thw 170 deflated from mid-March 2017 until mid-March 2018, with a total volume loss of 0.89 ± 0.05 km 3 and 1.91 ± 0.06 km 3 , respectively, draining significantly more water than during 2013 activity. Over the same period Thw 124 inflated by 3.20 ± 0.06 km 3 and settled 5.2 m higher than prior to 2013 drainage. Thw 70 shows no evidence of either drainage or recharge during 2017, remaining at a near-constant elevation.

Increase in Lake Area
During the 2017 drainage event the upstream lakes change sized compared against the extent of activity in 2013 ( Figure 2). Thw 170 and Thw 142 lake area expanded by 55 km 2 and 64 km 2 , respectively, both to the east of the 2013 boundaries. Thw 124 area decreased by 120 km 2 during 2017 activity, while maintaining the general shape of the 2013 boundary.

Recharge Period
Following the termination of lake activity in 2014, the lakes steadily regained volume over the inter-drainage period. Thw 70 , Thw 124 , Thw 142 , and Thw 170 gained 0.90 ± 0.06, 1.44 ± 0.11, 0.29 ± 0.03, and 0.46 ± 0.03 km 3 over 5.1, 2.6, 3.3, and 3.1 years, respectively. Notably, by mid-2017, Thw 170 regained the elevation that was lost during the 2013 drainage event, while elevation for the other lakes increased but remained below pre-2013 levels. We believe that this volume gain was predominantly caused by recharge through subglacial water transport, rather than ice flow divergence or blowing snow. Divergence has a negligible role, as our observations suggest an impact of no more than 0.1 m yr −1 (see Figure S1). If blowing snow had an impact, we would expect a bias in elevation change towards the prevailing wind direction-a result that we do not observe (see Figures S2 and S3). These volume changes correspond to a yearly recharge rate of 0.18 ± 0.02, 0.57 ± 0.05, 0.10 ± 0.01, and 0.14 ± 0.03 km 3 yr −1 at Thw 70 , Thw 124 , Thw 142 , and Thw 170 , respectively. It is worth noting these rates potentially reflect the balance between positive contributions from subglacial melt production and leakage from upstream lakes and negative contribution from leakage of the lakes into the downstream system.

Water Budget
Our observations of volume change for 2017 (Figure 3) can be used to determine the behavior of the subglacial drainage system during the drainage activity. Both the 2013 and 2017 drainage events, as well as the predicted drainage pathway (Figure 2), demonstrate the connectivity of the lake system. We assume that all discharged water from the two upstream lakes directly contributes to the rapid recharge of Thw 124 . A total of 2.80 km 3 of water from the two upstream lakes contributed to fill Thw 124 , accounting for 87.5% of the observed volume gain. The inclusion of Thw 124 's recharge rate in the budget leads to a predicted volume gain of 3.37 km 3 : 0.17 km 3 larger than the observed increase of 3.20 km 3 . This excess water falls within the uncertainty range attached to our volume change estimates, which suggests that the filling at Thw 124 is a product of the drainage of the two upstream lakes and background melt production.
As the lakes are connected subglacially, any change to water levels impacts the hydraulic gradients between the lakes. Background hydraulic gradients from Thw 170 to Thw 142 , Thw 170 to Thw 124 , and

Subglacial Water Flow
During 2013 activity Thw 124 , Thw 142 , and Thw 170 displayed dynamic volume change over 240, 180, and 150 days, respectively, while in 2017 each lake was active for approximately 300 days (Figure 1). Water fluxes also show contrasting behavior between the 2013 and 2017 events. In 2013 the rate of volume change was roughly symmetrical on either side of the peak discharge, while 2017 displays clear asymmetry-with post-peak discharge spanning three times the duration of pre-peak discharge. Using a simple R-channel assumption for modeling drainage pathways between the lakes, we found that average discharge rates from Thw 170 during 2017 activity lead to a channel with cross-sectional area of 14.6 m 2 and a radius of 3.0 m. Water within this channel would flow at a mean velocity of 3.5 m s −1 , causing melt at the channels side walls at a rate of 0.25 m 3 s −1 . Assuming this melting rate was sustained over the period of lake activity, an additional 0.008 km 3 of water would be injected into the subglacial system-negligible compared against the water mobilized from the lakes. Should channel behavior be dictated by the combined discharge from Thw 170 and Thw 142 , we would expect a channel with a cross-sectional area of 22.1 m 2 and radius of 3.8 m. Water would flow at a mean velocity of 3.8 m s −1 , which corresponds to a melting rate of 0.42 m 3 s −1 . A channel of this size would contribute 0.013 km 3 of water to the subglacial system during 2017 activity.

Water Budget
The net volume change of the ice sheet must be conserved from principles of mass conservation. Therefore, observing the flux of each lake allows us to infer the movement of water throughout the subglacial system. Thw 124 appears to be the downstream limit of the drainage event in 2017. Hence, we assume that water discharged from upstream lakes combined with Thw 124 recharge rate and drainage-related melting of the channels side was responsible for the observed volume gain. Under this assumption we expect to see a volume gain of 3.37 km 3 , which is 0.17 km 3 greater than the actual volume gain observed at Thw 124 -but falls within the uncertainty range. Background recharge rates are required to close the water budget, as without this component the water supply into Thw 124 would be too low to explain the observed volume gain. Unlike in 2013, where the four lakes drained out of the system, there is no evidence of subglacial lake activity downstream of Thw 124 in 2017 (Figure 1).

Notable Differences in Behavior Between Drainage Events
Our observations highlight a marked difference in the rate and evolution of water movement between lakes during the 2013 and 2017 events. First, during the 2013 drainage event lake activity followed a cascading pattern spanning 6 months, while in 2017 all lakes activated nearly simultaneously (Figure 1). Second, volume Despite the apparent simultaneous activity in 2017, we suspect the lakes activity followed a cascading pattern, similar to what occurred in 2013, but with a significantly faster transfer of water that our timeseries temporal resolution could not resolve. In such a scenario, Thw 170 would activate first and trigger Thw 142 , with discharged water filling Thw 124 . This scenario is further evident considering that Thw 170 drained at similar volumes in 2013 and 2017, which suggests the lake overcame its hydro-potential barrier during both drainage events. Therefore, Thw 170 could be considered the trigger to lake activity within the system and might be responsible for controlling future drainage events.
The rapid transfer of water in 2017 might have been made possible due to the development of a more efficient drainage system, likely following 2013 activity. Several possible mechanisms could be responsible for this change. For instance, discharge rates from the 2013 drainage event might have caused the formation of a channelized system between the lakes. While such channels would likely shrink due to creep closure, they may not have fully closed due to the transfer of water between the lakes during the inter-drainage period. This could precondition the system, allowing for rapid channel expansion during 2017 activity. Alternatively, discharge from the 2013 event might have caused sediment mobilization and potential channel erosion, leading the development of a more efficient drainage system (Brisbourne et al., 2017;Kirkham et al., 2019). Our reasoning is speculative as there is insufficient evidence available to either validate or negate these hypotheses. Nevertheless, this change in efficiency indicates complex behavior of the subglacial system, which is deserving of further study.
Discharge rates displayed a clear symmetrical ramp up and descent pattern in 2013 (Figure 1c). Conversely, discharge rates in 2017 spiked rapidly before tailing off over 6 months. All lakes drained in 2013 indicating an open system, while in 2017 the subglacial system could be considered closed, with a limit at Thw 124 which collected water and prevented significant discharge downstream. The influx of water into Thw 124 would have increased the pressure head within the lake, while the pressure head of upstream features was decreased due to the lower water levels. Potential gradients between the features decreased by 20%, which might have been sufficient to prolong the discharge of upstream water. Hence, the prolonged tail of discharge rates in 2017.

Observations Regarding Thw 124
The conditions and triggers of lake drainage are still poorly understood. While Thw 170 likely acted as a trigger for both 2013 and 2017 events, Thw 124 displayed contrasting behavior. Thw 124 drained in 2013 but not in 2017, despite larger upstream discharge and the fact that the lake water level, following termination of 2017 activity, exceeded that of pre-2013 drainage (Figure 1). This suggests that Thw 124 could be acting like a roadblock, collecting water while preventing significant discharge downstream. The shutdown of downstream drainage could have taken place as early as mid-2014 during the onset of lake refill (Figure 1), in particular because Thw 124 refill took place at a rate threefold higher than at any of the other lakes ( Figure 3) and that modeling does not predict such a significant difference in refill rates between lakes .
The mechanism behind the change in behavior of Thw 124 is uncertain. As hypothesized earlier, the significant discharge rates from 2013 activity might have modified the hydraulic properties of the system, which might have formed a barrier to flow downstream. Alternatively, it could be related to the mode of channel formation. Recent modeling suggests that the accumulation of water within lake basins steepens the hydraulic gradient and allows greater flux downstream, which melts channels that can trigger drainage (Dow et al., 2016(Dow et al., , 2018. Prior to 2013 activity Thw 124 appeared to be at hydrostatic equilibrium, whereby flow into the lake is equal to flow out, evident by the sustained overall lake volume (Figure 1). This outwards flow might have melted small channels immediately downstream of the lake. During 2013 activity, when water from upstream reached Thw 124 , hydraulic gradients would steepen. This would force water over the downwards slope, melting larger channels and likely triggering drainage. As the drainage event tails off, discharge rates decrease causing creep closure within the channels. During the inter-drainage period, the behavior of the system changes which limits the amount of water discharged from Thw 124 , which would hamper the creation of channels. The influx of water from the 2017 event would increase hydropotential gradients between the lake and downstream, driving water over the reverse slope. However, inefficient downstream channels might have prevented drainage.

Recharge Rates
Our annual recharge rates are significantly above estimations based on simulations of melt generation and water routing, whether considering melt production over local or regional catchments . There is a degree of variability between lakes, with rates at Thw 70 , Thw 124 , and Thw 170 significantly above predicated recharge, while Thw 142 is close to the high-end estimate of Smith et al. (2017) (Figure 3b). The variability between predicted and observed recharge rates can be explained by lake-to-lake water transfer, along with uncertainty in the subglacial network. However, the overall discrepancy between our observed rates and modeled values suggests estimates of melting rates at the bed are likely substantially underestimated. It is likely modeled subglacial melt production was underestimated, at least in part, because modeled melt did not incorporate the elevated geothermal heat flux that has been suggested based on radar observations located within the Thwaites lakes' catchments (Schroeder et al., 2014). In particular, Schroeder et al. (2014) suggests there may localized hot spots where heat flux is greater than the background. Alternatively, the catchment and water routing used in Smith et al. (2017) may not be representative of the true conditions at the bed, impacting the accuracy of their derived recharge rates. Our altimetrically derived recharge rates seemingly imply the second scenario of Smith et al. (2017): whereby lakes regain discharged water using within catchment-scale melt production and water supplied from upstream. Assuming Thw 70 , Thw 142 , and Thw 170 recharge at the same rate as per the inter-drainage period we expect each lake to regain its pre-2013 levels in 3.2, 11.5, and 13.6 years respectively. Estimates derived from modeled total melt production instead imply a recharge time of 8.1, 9.6, and 43.4 years for Thw 70 , Thw 142 and Thw 170 respectively. Given Thw 170 likely triggered the 2013 and 2017 drainage events, and will likely trigger future events, the significantly shorter recharge time implies the drainage cycle of the system is shorter than previously thought.

Conclusions
In mid-2013 a system of interconnected subglacial lakes under the central part of the Thwaites Glacier drained . Our altimetry measurements reveal a second period of lake activity in 2017, with discharged water tracked throughout the system. Both events are compatible with a cascading transfer of water, initiated by the most upstream lake. Observations reveal significant differences between the 2013 and 2017 drainage events. Unlike 2013 activity, in 2017 a downstream lake acts as a limit for the movement of subglacial water. This lake displayed rapid recharge, forced with discharge from the upstream lakes, increased in volume by 3.20 km 3 and settled 5.2 m higher than pre-2013 levels. Across the lake system, discharge is 29% greater in 2017 than in 2013 with lake sizes expanding by 119 km 2 . During 2013 activity each lake initiated within a 6-month period, while in 2017 the lakes activated within 45 days of each other. These characteristics point towards an increase in efficiency of the active subglacial system in 2017. Observations during the inter-drainage period indicate that lake recharge rates are 137% higher than modeled estimates. This implies that subglacial lakes recharge using melt supplied from local and upstream sources and that geothermal heat flux and basal friction produce more melt water than currently predicted.

Data Availability Statement
Our rate of change maps, lake masks, and lake timeseries are freely available from 4D Antarctica (https://4dantarctica.org/products/). The CryoSat-2 satellite altimetry data are freely available from the European Space Agency (https://earth.esa.int/web/guest/data-access). The BedMachine ice thickness and bed elevation data are freely available from the National Snow and Ice Data Centre (https://nsidc.org/data/ NSIDC-0756/versions/1). The ice velocity products are based on Copernicus Sentinel-1 data made available through the European Space Agency; the products are available upon request (http://cryoportal. enveo.at). The wind velocity and direction data are freely available from the Physical Sciences Laboratory (https://psl.noaa.gov/data/gridded/data.ncep.reanalysis.derived).