1. Introduction
The vast majority of ice in the Antarctic Ice Sheet drains from the continent to the ocean through fast-flowing glaciers and ice streams (Rignot and others, Reference Rignot, Mouginot and Scheuchl2011). The high velocities of these features are theorised to be maintained by the presence of meltwater at the ice-sheet bed, which reduces basal friction (Tulaczyk and others, Reference Tulaczyk, Kamb and Engelhardt2000). The movement of subglacial water has been linked to transient glacier flow acceleration (Stearns and others, Reference Stearns, Smith and Hamilton2008), nutrient mixing (Vick-Majors and others, Reference Vick-Majors2020) and enhanced melt at the grounding line (Le Brocq and others, Reference Le Brocq2013; Wei and others, Reference Wei2020). This implies that the location and movement of subglacial water are first-order constraints controlling Antarctica's mass balance (Bell, Reference Bell2008). Therefore, constraining subglacial water is vital to understanding the ice-sheet's response to climate change (Karlsson and others, Reference Karlsson2021).
Water transport from the interior of Antarctica was once thought to be a steady-state process (Parizek and others, Reference Parizek, Alley, Anandakrishnan and Conway2002). It is now known that subglacial water collects in hydrological sinks, which store and release water in episodic events (Gray, Reference Gray2005; Wingham and others, Reference Wingham2006; Siegfried and Fricker, Reference Siegfried and Fricker2018). Satellite altimetry can detect these features by searching for a localised elevation change signal at the ice-sheet's surface. This behaviour is interpreted as the movement of basal water in and out of ‘active’ subglacial lakes (Smith and others, Reference Smith, Fricker, Joughin and Tulaczyk2009; Wright and Siegert, Reference Wright and Siegert2012). These active subglacial lakes are often located beneath Antarctica's fast-flowing regions and could temporarily alter ice-sheet mass balance by modulating the amount and location of water during episodic drainage events (Siegfried and Fricker, Reference Siegfried and Fricker2018).
Following a drainage event, subglacial lakes often recharge by collecting a proportion of subglacial water flux and thus regain volume. Various factors influence the change in subglacial lake volume, including inflow from upstream subglacial melt production, inflow from the discharge of upstream subglacial lakes and outflow into the downstream subglacial hydrology system. Precise quantification of both positive and negative contributions to subglacial lakes is infeasible, given the unknown and unobservable nature of both inflow and downstream release. Nevertheless, the net flux into a subglacial lake is responsible for the detected change in its volume as estimated by satellite altimetry. As such, observing the behaviour of subglacial lakes during their recharge period can provide insight into subglacial melt production and hydrology at steady-state conditions, while the drainage phase is sensitive to transient dynamics. One example of this link can be illustrated in four subglacial lakes beneath the Thwaites glacier. These lakes underwent a drainage event in 2013 (Smith and others, Reference Smith, Gourmelen, Huth and Joughin2017) and again in 2017 – with a clear recharge period between these two events (Malczyk and others, Reference Malczyk, Gourmelen, Goldberg, Wuite and Nagler2020). Malczyk and others (Reference Malczyk, Gourmelen, Goldberg, Wuite and Nagler2020) extracted recharge rates during the networks recharge phase and compared these observations against modelled values generated by Smith and others (Reference Smith, Gourmelen, Huth and Joughin2017). The modelled values relied on a single melt map and were routed through an older bed topography realisation. Observation-driven recharge rates were significantly greater than those produced from the modelled output, which implied that either subglacial melt production under the ice sheet was underestimated or the presence of inaccuracies in the modelled subglacial network (Malczyk and others, Reference Malczyk, Gourmelen, Goldberg, Wuite and Nagler2020). However, such an analysis was isolated to Thwaites and has not been extended to Antarctica.
Direct observations of the subglacial system are difficult due to the ice-sheet's thickness – with in situ observations involving drilling experiments (Lukin and Vasiliev, Reference Lukin and Vasiliev2014; Tulaczyk and others, Reference Tulaczyk2014; Priscu and others, Reference Priscu2021). As such, estimates of subglacial melting rates across Antarctica are primarily constrained by models considering the impact of geothermal heat flux, vertical dissipation and frictional heating (Joughin and others, Reference Joughin2009; Pattyn, Reference Pattyn2010). With few in situ observations of the subglacial system across Antarctica, it is currently challenging to validate the results produced by such models. Here we use CryoSat-2 altimetry to produce time-dependent volume changes for subglacial lakes across the Antarctic ice sheet. From these, we extract recharge rates, which act as a proxy for the lower bound on subglacial melt production. These values are compared against the theoretical recharge rates derived by running a routing algorithm forced by estimates of subglacial melting. This approach provides us with the unique opportunity to compare direct observations of the subglacial network against the predicted behaviour.
2. Method
2.1. Deriving observation-driven subglacial lake recharge rates
Time-dependent ice-surface elevations were generated using swath processing of CryoSat-2 level L1b SARIn data acquired between 2010 and 2021 over all known active subglacial lakes (Fricker and others, Reference Fricker, Scambos, Bindschadler and Padman2007; Smith and others, Reference Smith, Fricker, Joughin and Tulaczyk2009; Siegfried and Fricker, Reference Siegfried and Fricker2018; Malczyk and others, Reference Malczyk, Gourmelen, Goldberg, Wuite and Nagler2020) situated within the CryoSat SARIn mask. In contrast to the commonly used point of closest approach (POCA), swath-processed SARIn exploits the entire radar waveform and can result in a one to two orders of magnitude increase in data density compared to the POCA alone (Gray and others, Reference Gray2013; Gourmelen and others, Reference Gourmelen2018).
To determine the spatial extent of subglacial lake activity during our study period, average surface elevation change rates were computed using a plane-fitting algorithm (McMillan and others, Reference McMillan2014) applied to swath-processed SARIn data. 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 and others, Reference Foresta2016). The regression was fitted iteratively to the data, omitting elevations differing more than three std dev. 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 localised elevation change (>0.5 m a−1) after a Gaussian low-pass filter was applied to the map (Fricker and others, Reference Fricker, Scambos, Bindschadler and Padman2007; Flament and others, Reference Flament, Berthier and Rémy2014; Smith and others, Reference Smith, Gourmelen, Huth and Joughin2017; Malczyk and others, Reference Malczyk, Gourmelen, Goldberg, Wuite and Nagler2020).
The temporal behaviour of the lake systems was determined by creating a surface elevation time series from 2010 until 2021. We used an adapted version of the point-to-point method (Malczyk and others, Reference Malczyk, Gourmelen, Goldberg, Wuite and Nagler2020) (see Text S1) first outlined in Gray and others (Reference Gray2015, Reference Gray2019) over our lake outlines to determine time-dependent elevations at a 30 d resolution, with elevations averaged over a 90 d search radius. This approach provides N − 1 time series realisations, with N being the number of time steps, which allows for calculating a mean time series and the associated error. The standard error of the mean of the time series realisations gives the statistical error at each time step. To isolate the behaviour of each lake relative to the catchment, we removed the background thinning or thickening 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.
Volume change through time was derived by integrating elevation change against the area of each lake mask. We approximate the volume budget of subglacial water flux by the volume corresponding to the surface elevation change. Subglacial lake 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 respect to the standard error of the regression slope.
We assume that the observed volume changes which correspond to our recharge rates are strictly due to the accumulation of subglacial meltwater, excluding other influencing factors. Nevertheless, we acknowledge certain challenges associated with this assumption. For instance, ice flow divergence, the accumulation of blowing snow in the upwind direction and changes in basal traction can all lead to a surface response that may be misinterpreted as recharge (Sergienko and others, Reference Sergienko, MacAyeal and Bindschadler2007). Moreover, the identification of the water source contributing to a lake's recharge presents inherent difficulties. It remains plausible that water leaking from an unobservable upstream lake could contribute to the recharge of a downstream lake.
2.2. Deriving modelled subglacial lake recharge rates
2.2.1. Subglacial melt rates
Basal melt rates, m r, (m a−1) for Antarctica were calculated using Eqn (1):
where G represents the geothermal heat flux (J a−1 m−2), F h basal frictional heating (J a−1 m−2), V d vertical conduction (J a−1 m−2), ρ i the density of ice (kg m−3) and L i the latent heat of fusion (J kg−1) (Cuffey and Paterson, Reference Cuff and Paterson2010).
We consider four different geothermal heat flux estimates for Antarctica, magnetically derived (Maule and others, Reference Maule, Purucker, Olsen and Mosegaard2005; Martos and others, Reference Martos2017) and seismically derived (Shapiro and Ritzwoller, Reference Shapiro and Ritzwoller2004; Shen and others, Reference Shen, Wiens, Lloyd and Nyblade2020).
The flow of the Antarctic Ice Sheet was assimilated using the higher-order ice-flow model STREAMICE (Goldberg, Reference Goldberg2011), with bed topography and ice thickness from BedMachine Antarctica V2 (Morlighem and others, Reference Morlighem2020) and observations of ice surface velocity (MEaSUREs InSAR-Based Antarctica Ice Velocity Map, Version 2 (Rignot and others, Reference Rignot, Mouginot and Scheuchl2017)). The basal friction and ice stiffness parameters were inferred at 5 km resolution via an adjoint method (Utke and others, Reference Utke2008; Goldberg and Sergienko, Reference Goldberg and Sergienko2011), obtaining values for basal shear stress and velocity. The product of these terms gives the basal frictional heating, F h.
Vertical heat conduction was estimated using the gradient of englacial temperature from Van Liefferinge and Pattyn (Reference Van Liefferinge and Pattyn2013) and the thermal conductivity of ice (Cuffey and Paterson, Reference Cuff and Paterson2010).
We examine the influence of two additional basal melting products (Joughin and others, Reference Joughin2009; Van Liefferinge and Pattyn, Reference Van Liefferinge and Pattyn2013) in conjunction with the melting rates generated using the aforementioned method within the Thwaites glacier drainage basin – a region specifically selected due to better observational constraints on the total water volume mobilised. It is worth noting that the Van Liefferinge and Pattyn (Reference Van Liefferinge and Pattyn2013) estimate incorporates an average of the Maule and others (Reference Maule, Purucker, Olsen and Mosegaard2005) and Shapiro and Ritzwoller (Reference Shapiro and Ritzwoller2004) Geothermal Heat Flux (GHF). Furthermore, the basal melting estimates by Joughin and others (Reference Joughin2009) do not encompass the entire Thwaites drainage basin. We address this limitation by utilising a composite melting map derived from the average of our four estimates. The subglacial flow paths across the Thwaites glacier and the corresponding recharge rates for Thw124 are derived by employing the FD8 routing approach (Le Brocq and others, Reference Le Brocq, Payne, Siegert and Alley2009), as outlined in section 2.2.2.2.
2.2.2. Subglacial routing
To identify likely subglacial drainage pathways and fluxes, we calculate hydraulic potential (ϕ) as,
where ρ i represents the density of ice, ρ w the density of water, Z s the surface height of the ice sheet and Z b the bed topography, assuming basal water pressure is at overburden everywhere (Shreve, Reference Shreve1972). We used surface height and bed elevation from Bed Machine Antarctica V2 (Morlighem and others, Reference Morlighem2020).
If the basal water pressure at the bed is equal to the overburden pressure, then water will tend to flow along pathways perpendicular to the hydropotential, from high to low. Therefore, our potential maps allow us to determine the pathway that subglacial water will likely follow. However, localised closed basins within the hydropotential at short spatial scales prevent water flow and thereby block water transport over distance. Subglacial water can escape these basins by seeping through subglacial till, flowing through low-pressure channels, or being transported by valleys that our current surveys cannot resolve (Smith and others, Reference Smith, Gourmelen, Huth and Joughin2017). To address this short coming, our hydropotential map is conditioned by artificially filling closed depressions to yield a map where long-distance water transport is possible. Such artificial filling assumes no volume is loss to any hydrological sinks and that flow is perfect from source to grounding line.
We estimate the flux of water moving through the subglacial network by combining our routing maps with our subglacial melting rates. From this flux map, we can estimate the yearly supply of water entering a subglacial lake. In the case of a lake network, the predicted supply of meltwater to downstream lakes does not consider water released from upstream lakes or other hydrological sinks.
We outline the three subglacial routing methods and their corresponding depression-filling algorithms used in our study below. For each method, we also describe how recharge rates are calculated.
2.2.2.1. D8 routing algorithm
In the first step, a depression-filling algorithm identifies localised basins and contiguous flat areas in the hydropotential map and then fills them using a grey-weighted distance transformation algorithm (Schwanghart and others, Reference Schwanghart, Groom, Kuhn and Heckrath2013). A deterministic routing scheme, the so-called D8-scheme, is then applied to the conditioned hydropotential map using the TopoToolBox flow model (Schwanghart and Scherler, Reference Schwanghart and Scherler2014). This scheme transports a cell's water to the neighbouring cell with the lowest hydropotential. As this method is non-diffusive, each cell belongs to a single catchment, and thus it is possible to derive a map of regions contributing to the inflow of water towards a subglacial lake. Modelled recharge rates for each lake are calculated by summing the subglacial melting rates over its corresponding catchment area.
2.2.2.2. Thin-film-based subglacial water flow (FD8)
This method, proposed by Le Brocq and others (Reference Le Brocq, Payne, Siegert and Alley2009), differs from that of the D8 algorithm in its routing approach. The method first identifies localised basins within the hydropotential map and then fills them using a ‘flood-filling’ algorithm. This algorithm assigns the average of the four neighbouring cells to the sink site and repeats until localised sinks are no longer detected (Le Brocq and others, Reference Le Brocq, Payne, Siegert and Alley2009). As this method is computationally expensive, particularly over large catchments, after 2000 sweeps, the remaining sinks are filled by raising all sink cells by a small but random amount until long-distance flow can be achieved. This introduces a minor stochastic element into the method, but in practice, this randomised sink filling has minimal impact on the output flux and routing.
The approach routes water by incorporating a deterministic scheme, where water is distributed among the neighbouring cells based on the magnitude of the surrounding hydropotential, thus similar to the D-infinity approach often used in hydrological routing (Tarboton, Reference Tarboton1997). As a result, the Le Brocq and others (Reference Le Brocq, Payne, Siegert and Alley2009) routing model has the potential to produce multiple flow paths, which is impossible under the D8 scheme. When calculating recharge rates for a subglacial lake, it is essential to consider that the network might be branched, and therefore there may be multiple inflows and outflows to the lake. Each cell's contribution to a lake might vary depending on the magnitude of the underlying FD hydropotential map. Therefore, recharge rates are calculated by identifying cells on the edge of a subglacial lake and the corresponding proportion of flux directed into the lake. Once all cells have been identified, all the inflowing fluxes are summed to produce our modelled recharge rate.
2.2.2.3. Stochastic D8
The stochastic D8 approach combines a classic D8 routing scheme with an uncertainty quantification via Monte Carlo ensemble runs. The employed D8 scheme does not require a depression-filling algorithm to pre-process the hydropotential map. Instead, it employs a breach algorithm which routes the water from the minimum of a depression uphill to the lowest overflow point, from where it is then routed according to the standard D8 procedure again.
Uncertainty quantification is achieved by performing Monte Carlo simulations, wherein the model input fields are systematically varied within plausible ranges. It is important to note only the mean values of the ensemble of runs are presented. The model inputs subject to uncertainty include the bed elevation, surface elevation and the flotation fraction, which represents the ratio between subglacial water pressure and ice overburden pressure used in the hydraulic potential. Since only mean values are reported, uncertainty is not considered for basal melt. To introduce variations in each run, a Gaussian random field is added to each input parameter. The Gaussian random field is characterised by its amplitude (σ), expressed as a std dev., and correlation length. For most variables, the correlation length is set to 10 km, while for the surface elevation, it is set to 1 km to account for the high spatial resolution of the underlying REMA dataset. The amplitude (σ) is determined as follows: for the bed elevation, it is derived from the spatially varying field errbed of the BedMachine dataset, which represents an estimate of uncertainty. However, as this estimate is closer to a maximum uncertainty than a std dev., errbed is divided by 4.5 to obtain an appropriate value for σ. For the surface elevation, σ is set to 0.6 m (Howat and others, Reference Howat, Porter, Smith, Noh and Morin2019). For the flotation fraction, which has a mean value of 0.97, σ is assigned a value of 0.01. For each run, the inflow into the studied lakes is recorded using the standard D8 scheme.
2.3. Comparing modelled and observation-driven recharge rates
During the recharge phase of a subglacial lake, the influx of water into the lake is assumed to be solely driven by melt production in the upstream catchment. Therefore, the recharge rates obtained from this phase can be used as a proxy to estimate the volume of subglacial meltwater that needs to be directed towards the lake in order to account for the observed changes in volume (Fig. 1a). However, it is important to acknowledge that certain external factors can influence the inflow of water into a lake. For instance, water released from upstream subglacial lakes can contribute additional flux, supplementing the annual melt production. In such cases, it is beneficial to consider the subglacial lakes as a collective entity (Fig. 1b). By summing up the recharges of these lakes, a more stringent constraint on the inward flux can be obtained compared to considering each lake individually. This approach remains valid only if the subglacial lakes are hydrologically connected and receive recharge through the same drainage pathway, with no other significant pathways supplying the lakes or subsections of the hydrological system. Under these conditions, the subglacial melt between the lakes is expected to be negligible compared to the volume mobilised throughout the rest of the catchment. A more stringent constraint can be achieved if it can be established that the subglacial system is ‘blocked’ downstream from a subglacial lake (Fig. 1c), thereby restricting or completely halting the downstream flux. The cause of this blockage could be attributed to low permeability or the absence of well-developed efficient channels. In this scenario, there would be less ambiguity regarding the downstream flux. Consequently, the summation of the lakes' recharges would approximately match the predicted inflow into the network.
Figure 1. An illustration of the mechanics influencing subglacial lake recharge in three different settings. (a) The recharge components of a singular lake. An upper bound of flux into the lake can be estimated with routing models, and the total change of the system can be summarised with estimates of volume change derived from altimetry. However, the flux leaving the lake is unknown. (b) The recharge components of a connected lake network. The outwards flux of an upstream lake feeds a downstream feature. Therefore, the sum of volume changes must be less than the predicted flux into the network. (c) The recharge components of a connected lake network under the assumption that downstream discharge is blocked. Under such a scenario, the sum of volume changes must be comparable to the predicted flux into the network.
3. Results
3.1. Subglacial lake observations
Out of the 140 currently identified active lakes in Antarctica (Smith and others, Reference Smith, Fricker, Joughin and Tulaczyk2009, Reference Smith, Gourmelen, Huth and Joughin2017; Siegfried and Fricker, Reference Siegfried and Fricker2018; Livingstone and others, Reference Livingstone2022), we focus on the 40 lakes that fall within CryoSat's SARIn mask mode (Fig. 2). The remaining lakes are excluded from our analysis due to the inapplicability of our method to those cases.
Figure 2. Distribution of active subglacial lakes across the Antarctic. Yellow circles represent subglacial lakes discovered using IceSat-1 or CryoSat-2 altimetry, while red circles represent the subset of these lakes that display some recharge activity during our study period. The background map represents subglacial flux and flow paths across Antarctica derived using the FD8 routing approach. The highlighted orange region represents the SARIn data mode coverage for CryoSat-2 and acts as a boundary for our method.
Among the 40 previously documented subglacial lakes within our study area, we have identified 11 lakes distributed across six different drainage systems exhibiting recharge activity between 2010 and 2021. For a detailed rationale for each lake's inclusion, refer to Text S2. Additionally, we have discovered two previously unknown subglacial lakes, expanding our understanding of the region. The following sections present our findings concerning the 13 lakes demonstrating drainage and recharge activity throughout the study period. Detailed information regarding the volume gains and corresponding recharge rates during each lake's recharge period can be found in Table 1. Furthermore, Table S1 provides comprehensive documentation of volume changes during the drainage activity of each lake under investigation. The time series of surface elevation change for the lakes not discussed below can be found in the Supplementary material (see Figs S1–S6), along with individual flux maps for each routing method (Figs S7–S12).
Table 1. Volume gains and recharge rates over each lake during the recharge period
3.1.1. Thwaites glacier subglacial lake system
The Thwaites Lake network encompasses a set of four hydrologically connected subglacial lakes, which Smith and others (Reference Smith, Gourmelen, Huth and Joughin2017) discovered and are located on the central part of the Thwaites glacier in West Antarctica. Ice velocities over the lake network average 170 m a−1.
3.1.1.1. Subglacial lake activity and recharge
Our time series of surface elevation change over Thwaites glacier capture the 2013 drainage event, as discussed by Smith and others (Reference Smith, Gourmelen, Huth and Joughin2017), and the 2017 event, as discussed by Malczyk and others (Reference Malczyk, Gourmelen, Goldberg, Wuite and Nagler2020). Subglacial lakes Thw70, Thw124, Thw142 and Thw170 drained between April 2013 and May 2014, with a second period of activity between March 2017 and March 2018 (see Table S1 for volume changes). Between these two drainage events, the lakes steadily regain volume via subglacial melt production (Malczyk and others, Reference Malczyk, Gourmelen, Goldberg, Wuite and Nagler2020), with corresponding volume gains and recharge rates documented in Table 1.
3.1.1.2. Connectivity and modelled recharge rates
The simultaneous drainage and water transfer between the lakes indicates hydrological connectivity. Our model of subglacial flow paths (Fig. 3a) suggests varying degrees of hydrological connection between the lakes, depending on the routing approach employed. Regardless of the routing scheme used, the primary pathway passes through Thw142 and Thw124 before connecting with a second branch located north of Thw124. All three routing schemes link Thw170 to the downstream lake network. However, under the D8 scheme the main drainage pathway does not flow through Thw170, resulting in significantly lower modelled recharge rates (Table 2). Thw70 is connected to the main pathway under the FD8 and stochastic D8 approaches, whereas, under the D8 scheme, water is routed east of Thw70, resulting in its disconnection from the upstream network.
Figure 3. Rates of surface elevation change, bed elevation, watering routing and time-dependent volume changes for subglacial lakes in the Thwaites Lake Region. (a) Bed elevation (from BedMachine (Morlighem and others, Reference Morlighem2020)) and surface elevation changes. Thw70, Thw142 and Thw170 display rates of surface elevation change from 2010 to 2016, while Thw124 displays rates of surface elevation change from 2014 to 2020. The rates of elevation changes represent the surface response on the ice-sheet surface to changes in water volume at the bed. White outlines represent the 2013 drainage event lake masks from Malczyk and others (Reference Malczyk, Gourmelen, Goldberg, Wuite and Nagler2020) for Thw70, Thw124, Thw142 and Thw170. The map insert illustrates the location of the lake region. The blue lines represent flow paths derived using the D8 routing approach, the green from FD8 and orange from stochastic D8. These represent the probable hydrological flow paths under the three different schemes. (b) Mean volume changes for each lake. The shaded region represents a 95% confidence interval. The dotted lines represent periods of recharge activity. See Supplementary Figure S1 for original elevation change time series and background thinning component.
Table 2. Averaged modelled recharge rates for each routing approach (D8, FD8 and Stochastic D8) of all melting maps, and range of modelled recharge rates obtained across all melting maps
Table 2 provides the modelled recharge rates for each lake, illustrating that the Stochastic D8 approach yields the highest flux into the network. The recharge rates differ based on the choice of GHF model, with the Shapiro and Ritzwoller (Reference Shapiro and Ritzwoller2004) GHF predicting the largest recharge rate, while the Martos and others (Reference Martos2017) GHF predicts the smallest. For Thw124, the GHF and frictional heating components contribute 72 and 28%, respectively, towards the total meltwater mobilised towards the lake, while vertical conduction removes 22% of this value.
3.1.1.3. Hydrologically linked network
Thw124, Thw142 and Thw170 are all part of the same drainage pathway, with minimal major branches supplying water to the system. Therefore, it is likely that these three lakes receive recharge from the same upstream meltwater source, and that while recharging, upstream lakes are acting as sink for the lakes further downstream, which needs to be accounted for when comparing observation-driven and modelled recharge. Considering these three subglacial lakes as a set when evaluating subglacial melt production is therefore needed, with the cumulative observation-driven recharge rates serving as a better-constrained lower bound. Under this assumption, a modelled recharge rate of at least 0.81 km3 a−1 is required to account for volume gains within the lake network.
3.1.1.4. Comparison of melt products over the Thwaites Lake network
We use the FD8 routing approach to estimate the flux entering Thw124. The modelled flux can vary significantly depending on the subglacial melt rate estimate used (Fig. 4). Our subglacial melt rate estimate provides the largest recharge rate, of 1.32 km3 a−1, which is 0.51 km3 a−1 greater than the sum of the observed rate of volume change at Thw124, Thw142 and Thw170. Recharge rates calculated using the melt estimate from Van Liefferinge and Pattyn (Reference Van Liefferinge and Pattyn2013) are 0.81 km3 a−1. While the subglacial melt rate from Joughin and others (Reference Joughin2009) predicts a recharge rate of 0.77 km3 a−1. Here 0.46 km3 a−1 of the recharge rate is sourced from the portion of the catchment covered by the melt estimate from Joughin and others (Reference Joughin2009), with the remaining 0.31 km3 a−1 sourced from our new subglacial melt estimate over the remaining portion of the catchment.
Figure 4. Modelled recharge rates at Thw124, compared against observation-driven recharge rates, were derived using the FD8 routing scheme and forced with an average composite of our melting maps (blue), the Van Liefferinge and Pattyn (Reference Van Liefferinge and Pattyn2013) melt and the Joughin and others (Reference Joughin2009) melt. The black dashed line represents the minimum subglacial water required to account for observation-driven recharge rates. The dark green bar represents the proportion of water forced into Thw124 using the Joughin and others (Reference Joughin2009) melt alone, while the light green represents the remaining proportion caused by the average composite of our melt maps.
3.1.2. Mercer and Whillans subglacial lake system
The Mercer and Whillans lake network comprises 11 subglacial lakes distributed along the Mercer and Whillans ice streams, situated west of the Transantarctic Mountains. It is worth noting that two subglacial lakes fall outside the SARIn collection zone and, therefore, are omitted from this study. Among the nine subglacial lakes within the study area, seven are within 100 km of the grounding line of the Ross Ice Shelf. Ice velocities over the network average 90 m a−1.
3.1.2.1. Subglacial lake activity and recharge
Our time series of surface elevation change over the Mercer and Whillans subglacial system (Fig. 5b) captures the 2011–14 drainage event discussed by Siegfried and Fricker (Reference Siegfried and Fricker2018). Upper Subglacial Lake Conway (USLC) drained between July 2010 and September 2013. Second to drain was Subglacial Lake Mercer (SLM), which rapidly discharged between May 2012 and April 2015. Finally, Subglacial Lake Conway (SLC) drained steadily from November 2012 until January 2014.
Figure 5. Rate of surface elevation change, bed elevation, water routing, melting rate and time-dependent elevations/volume change for lakes in the Mercer and Whillans region. (a) Bed elevation (from BedMachine (Morlighem and others, Reference Morlighem2020)) and surface elevation change rates. SLE displays rate of surface elevation change from 2010 to 2021, SLM from 2015 to 2018, SLC from 2014 to 2019 and USLC from 2013 to 2020. The remaining lakes show rate of surface elevation change from 2010 to 2021. The map insert illustrates the location of the lake region. The cyan line represents routing derived from the D8 approach, green from FD8 and orange from stochastic D8. (b) Mean volume changes for each lake. The shaded region represents a 95% confidence interval. The dotted lines represent periods of recharge activity. See Supplementary Figure S2 for original elevation change time series and background thinning component.
Our observations indicate a second, previously undocumented, episode of lake activity starting in 2018. SLM activates from May 2018 until January 2019 before regaining volume and draining between December 2019 and August 2020. SLC is active from June 2019 until July 2020. Table S1 documents volume changes for both these periods of lake activity.
In between these two periods of lake activity, SLM and SLC experience recharge. USLC recharges over the rest of the study period following the termination of its drainage activity in 2013. Subglacial Lake Engelhardt (SLE) displays a sharp increase in volume between October 2013 and March 2015 before regaining volume at a constant rate. Volume gains and recharge rates for each lake can be found in Table 1.
3.1.2.2. Connectivity and modelled recharge rates
Our modelled subglacial flow paths (Fig. 5a) illustrate a strong hydrological connection between USLC, SLC and SLM, irrespective of the routing approach. This pathway meets with another flow path east of SLM before flowing towards the grounding line. However, it is worth noting that the FD8 scheme predicts a significantly weaker subglacial connection than the other two methods (see Fig. S8). SLE has a drainage path distinct from the abovementioned pathway, flowing from the northeast towards the west before meeting the grounding line.
Table 2 provides the modelled recharge rates for each lake, illustrating that the D8 approach yields the highest flux into the network. The recharge rates differ based on the choice of GHF model, with the Maule and others (Reference Maule, Purucker, Olsen and Mosegaard2005) GHF predicting the largest recharge rate, while the Martos and others (Reference Martos2017) GHF predicts the smallest. For SLM, the GHF and frictional heating components contribute 80 and 20%, respectively, towards the total meltwater mobilised towards the lake, while vertical conduction removes 18% of this value.
3.1.2.3. Hydrologically linked network
USLC, SLC and SLM all exist on the same drainage pathway, with minimal major branches feeding the system and as such, assuming all three lakes recharge from the same meltwater supply is a safe assumption. This pathway must mobilise at least 0.74 km3 a−1 to account for the total volume gain within the three lakes. The D8 and stochastic D8 schemes can account for the observed volume changes, while the FD8 scheme cannot.
3.1.3. Slessor subglacial lake system
The Slessor network comprises a group of seven subglacial lakes situated in the central and upstream regions of the Slessor glacier, which is a constituent of Coats Land and feeds the Filchner Ice Shelf. Four of the seven lakes lie beyond the SARIn collection zone and, thus, are excluded from the analysis. The remaining trio of lakes is positioned ~150 km from the grounding line where ice flows at an average of 146 m a−1.
3.1.3.1. Subglacial lake activity and recharge
We have identified a localised region of surface elevation change (Fig. 6a) covering the extent of subglacial lake Slessor 2, albeit with a different shape reported by Smith and others (Reference Smith, Fricker, Joughin and Tulaczyk2009) and Siegfried and Fricker (Reference Siegfried and Fricker2018). Consequently, we updated the subglacial lake mask accordingly. Our time series of surface elevation change capture the 2014 drainage event observed by Siegfried and Fricker (Reference Siegfried and Fricker2018). Slessor 2 drains between April 2014 and October 2014 before recharging over the remaining proportion of the study period.
Figure 6. Rates of surface elevation change, bed elevation, water routing, melting rate and time-dependent elevations and volume changes for Slessor 2. (a) Bed elevation (from BedMachine (Morlighem and others, Reference Morlighem2020)) and rates of surface elevation change. Slessor 1 and Slessor 2 both display rates of surface elevation change from 2010 to 2021. The map insert illustrates the location of the lake region. The cyan line represents routing derived from the D8 approach, green from FD8 and orange from stochastic D8. (b) Mean volume changes for each lake. The shaded region represents a 95% confidence interval. The dotted lines represent periods of recharge activity. See Supplementary Figure S3 for original elevation change time series and background thinning component.
3.1.3.2. Connectivity and modelled recharge rates
Our three routing approaches indicate that the primary drainage pathway passes through Slessor 2 before connecting with Slessor 1 and flowing towards the grounding line (Fig. 6a). Table 2 provides the modelled recharge rates for Slessor 2, illustrating that the stochastic D8 approach provides the highest flux into the lake and FD8 the lowest. The Maule and others (Reference Maule, Purucker, Olsen and Mosegaard2005) GHF model predicts the largest recharge rate, while the Martos and others (Reference Martos2017) predict the smallest. GHF and the frictional heating components contribute 59 and 41% towards the total meltwater mobilised towards the lake, while vertical conduction removes 33% from the subglacial system.
3.1.4. Lambert subglacial lake system
The Lambert network comprises one subglacial lake located in the central region of the Lambert glacier, which is positioned in East Antarctica and feeds the Amery Ice Shelf. The lake is positioned ~100 km from the grounding line where ice velocities average 121 m a−1.
3.1.4.1. Two new lakes
We have identified two regions of localised elevation change to the north and south of previously documented subglacial lake Lambert 1 (Fig. 7a). We interpret these regions as subglacial lakes based on several lines of evidence. Firstly, the regions are positioned prominently along the modelled subglacial hydrology network (Fig 7a). Secondly, they exhibit a high degree of localisation with rates of elevation change far exceeding the surrounding region. Finally, the off-lake elevation change signal significantly differs from the on-lake signal (Fig. S4). In light of these observations, we have named these new lakes Lam80 and Lam110, reflecting their respective distances from the grounding line and following the naming convention established after the discovery of lakes at Thwaites (Smith and others, Reference Smith, Gourmelen, Huth and Joughin2017). These lakes form a hydrologically connected lake system, of which Lambert 1 is a member. It is worth noting that no activity was observed at Lambert 1 over the study period (Fig. S4).
Figure 7. Rates of surface elevation change, bed elevation, water routing, melting rate and time-dependent elevations/volume change for lakes in the Lambert region. (a) Bed elevation (from BedMachine (Morlighem and others, Reference Morlighem2020)) and rates of surface elevation change. Lam110 and Lambert 1 display rates of surface elevation change from 2015 to 2018. Lam80 displays rates of surface elevation change from 2010 to 2021. The map insert illustrates the location of the lake region. The cyan line represents routing derived from the D8 approach, green from FD8 and orange from stochastic D8. (b) Mean volume changes for each lake. The shaded region represents a 95% confidence interval. The dotted lines represent periods of recharge activity. See Supplementary Figure S4 for original elevation change time series and background thinning component.
3.1.4.2. Lake activity and recharge rates
Lam110 demonstrates two distinct periods of activity within our study period. The first period spans from September 2010 to December 2014, during which the lake gradually discharged its volume. The second period extends from March 2017 until the end of the study period and similarly exhibits a gradual discharge of volume. Between these two events, Lam110 enters a phase of recharge. Lam80 exhibits a consistent increase in volume throughout the entire study period, suggesting a continuous period of recharge for this subglacial lake.
The SARIn collection mode of CryoSat-2 employed in our study limits our ability to investigate potential upstream sources that may be supplying water to Lam80 and Lam110. Consequently, it remains plausible that the observed volume gains in these lakes could be attributed to water transfer originating from an unidentified upstream subglacial lake discharge. However, it is equally plausible that these volume gains might solely result from subglacial melting, signifying subglacial lake recharge. Due to the absence of a method to ascertain the origin of this water, we assume that the observed signals in both lakes are exclusively driven by meltwater production. We acknowledge that the presence of unknown hydrological processes further upstream introduces additional uncertainty into the system.
3.1.4.3. Connectivity and modelled recharge rates
There is a clear subglacial connection between Lam110 and Lambert 1, with water routed from the south towards the north and the grounding line. Under the FD8 and stochastic D8 schemes, Lam80 is disconnected from the other two subglacial lakes, with water routed towards the east of the lake. Under the D8 scheme, Lam 80, Lam110 and Lambert 1 are all hydrologically connected.
The D8 routing approach provides the highest flux into the subglacial system and the FD8 approach the lowest. The Shapiro and Ritzwoller GHF predicts the largest recharge rate, while the Martos and others (Reference Martos2017) predict the smallest. GHF and the frictional heating components contribute 71 and 29% towards the total meltwater mobilised towards the system, while vertical conduction removes 44% from the subglacial system.
3.1.5. David subglacial lake system
The David network comprises six subglacial lakes beneath the central and upstream regions of David glacier in Victoria Land. Four subglacial lakes lie beyond the SARIn collection zone and are omitted from consideration. The remaining two subglacial lakes are located ~140 km from the Drygalski Ice Tongue and both these lakes experienced drainage activity in 2004 (Smith and others, Reference Smith, Fricker, Joughin and Tulaczyk2009). Ice velocities over the network average 26 m a−1.
3.1.5.1. Subglacial lake activity and recharge rates
The time series of surface elevation change for David s1 (Fig. 8b) reveals a consistent linear increase, with oscillations that appear to have a seasonal component. The nature of these oscillations is unknown. One possible source is from surface processes such as surface mass balance, firn compaction or changes in the scattering depth of the radar signal. However, it is worth noting that these oscillations are strictly limited to the lake's boundaries (Fig. S5). One would expect the surface processes listed above to occur inside and outside the lake area equally unless the lake depression causes specific conditions affecting surface conditions seasonally. As such, the exact cause of these oscillations remains uncertain. Superimposed on these oscillations is a consistent trend of increased surface elevation, suggesting recharge following a drainage event in 2004 (Smith and others, Reference Smith, Fricker, Joughin and Tulaczyk2009). We have documented the corresponding volume gains and recharge rates associated with this recharge event in Table 1.
Figure 8. Rate of surface elevation change, bed elevation, water routing, melting rate and time-dependent elevations/volume change for Cook E2. (a) Bed elevation (from BedMachine (Morlighem and others, Reference Morlighem2020)) and rates of surface elevation change. David 1 and David s1 display rates of surface elevation change from 2010 to 2021. The map insert illustrates the location of the lake region. The cyan line represents routing derived from the D8 approach, green from FD8 and orange from stochastic D8. (b) Mean volume changes for each lake. The shaded region represents a 95% confidence interval. The dotted lines represent periods of recharge activity. See Supplementary Figure S5 for original elevation change time series and background thinning component.
3.1.5.2. Connectivity and modelled recharge rates
A distinct hydrological connection is evident between David 1 and David s1 (Fig. 8a) wherein water flows from the southern region towards the northwest. The FD8 and stochastic D8 routing approaches both direct water through David s1, while the D8 approach diverts water northward, bypassing the lake. Notably, the stochastic D8 routing approach yields the highest flux into the subglacial system, while the D8 approach yields the lowest. The Maule and others (Reference Maule, Purucker, Olsen and Mosegaard2005) GHF predicts the largest recharge rate, while the Martos and others (Reference Martos2017) predict the smallest. GHF and the frictional heating components contribute 71 and 29% towards the total meltwater mobilised towards the system, while vertical conduction removes 29% from the subglacial system.
3.1.6. Cook subglacial lake system
The Cook network comprises two subglacial lakes beneath the Cook Glacier's upstream region in Victoria Land. The most upstream subglacial lake, Cook E2, is positioned below very slow-flowing ice and is <30 km from the ice divide with the David Glacier. Cook E2 experienced a drainage event in 2008, with the lake discharging between 2.7 and 6.4 km3 (McMillan and others, Reference McMillan2013; Li and others, Reference Li, Lu and Siegert2020). Ice velocities over the network average at 4 m a−1.
3.1.6.1. Subglacial lake activity and recharge rates
The time series of surface elevation change for Cook E2 (Fig. 9b) demonstrates a consistent linear increase, with oscillations that appear to have a seasonal component. As with David s1, these oscillations are strictly limited to within the lake's boundaries (Fig. S6). The nature of these oscillations is unknown. Superimposed on these oscillations is a consistent trend of increased surface elevation, suggesting recharge following a drainage event in 2008 (McMillan and others, Reference McMillan2013). We have documented the corresponding volume gains and recharge rates associated with this recharge event in Table 1.
Figure 9. Rate of surface elevation change, bed elevation, water routing, melting rate and time-dependent elevations/volume change for lakes in the David region. (a) Bed elevation (from BedMachine (Morlighem and others, Reference Morlighem2020)) and rates of surface elevation change. Cook E1 and Cook E2 display rates of surface elevation change from 2010 to 2021. The map insert illustrates the location of the lake region. The cyan line represents routing derived from the D8 approach, green from FD8 and orange from stochastic D8. (b) Mean volume changes for each lake. The shaded region represents a 95% confidence interval. The dotted lines represent periods of recharge activity. See Supplementary Figure S6 for original elevation change time series and background thinning component.
3.1.6.2. Connectivity and modelled recharge rates
In light of Cook E2's proximity to an ice divide, its upstream subglacial network is considerably smaller than other catchments considered in this study. Consequently, only a negligible proportion of water originating from upstream passes through the lake (see Fig. S12). Modifying the routing methodology has little impact on the predicted magnitude of water that feeds Cook E2 or the subglacial network, except for minor localised variations.
All three routing approaches yield similar fluxes into the system (Table 2). Half of this flux can be attributed to subglacial melt within the lake basin alone. There is little variation between the four subglacial melt estimates. GHF accounts for most of the observed melt feeding Cook E2, with the frictional heating and vertical dissipation components contributing <0.1%.
3.2. Subglacial lake recharge rates
Except for Cook E2, the inflow of subglacial meltwater into all studied subglacial lakes and interconnected networks across Antarctica is accounted for by at least one of the routing approaches, regardless of the subglacial melt rate used (Fig. 10). The stochastic D8 scheme is the most successful approach, as it can explain observation-driven recharge under all heat flux settings for all but one of the observed lakes; Cook E2. The FD8 scheme can account for recharge in four of the six studied regions but fails to account for recharge in lakes within the SLC network and Cook E2. D8 can explain our observations in three of the six regions but falls short in accounting for recharge at Thw70, Thw170 (Fig. S14) and David s1, primarily because the modelled main pathway bypasses these lakes. Cook E2 is the only lake where the observation-driven recharge rates significantly exceed the modelled rates, regardless of which routing approach and subglacial melt rate estimate are used.
Figure 10. Modelled recharge rates, forced with four different melting approaches, compared against observation-driven rates of recharge derived from CryoSat-2 altimetry. The black dashed line represents the minimum subglacial water required to close the water budget.
4. Discussion
4.1. Comparison of observation-driven and modelled recharge rates
Except for Cook E2, the modelled subglacial lake recharge rates consistently equalled or surpassed the observation-driven rates under at least one routing approach and all basal melt models for all regions (Fig. 10). This outcome was anticipated, considering that an ample supply of subglacial water is essential for facilitating lake recharge, which highlights the validity of modelling recharge. In regions with a single recharging lake, such as Slessor 2, or with non-hydrologically linked lakes, such as SLE, modelled recharge rates were generally several times greater than the observation-driven values. The proportion of water not retained for recharge is anticipated to be discharged downstream, eventually reaching the grounding line, or being absorbed into hydrological sinks.
In contrast to singular lakes, regions that formed a network where observation-driven recharge rates were closer to the model-predicted inflow, such as the Thwaites and SLC networks. This suggests that these singular lakes provide only a loose lower bound on subglacial melt production, as observation-driven rates are typically a fraction of the magnitude of water mobilised within the catchment. Furthermore, outward fluxes are unknown. However, this ambiguity on outward fluxes can be partially alleviated when considering a set of connected lakes, which recharge along the same network and have no other additional water sources. In such a scenario, only the unknown outward flux of the most downstream lake impacts the total budget. Therefore, under the assumption of comparable lake leakiness, a stricter bound on subglacial melt production can be provided by taking the sum of recharge rates across a set of connected subglacial lakes, as opposed to considering each lake in the network individually.
The Thwaites Lake network presents a unique case that offers an exceptional opportunity to investigate subglacial melt production due to distinctive characteristics observed at Thw124, suggesting a ‘restricted system’ where the outward flux from the most downstream lake is less ambiguous (i.e. Fig. 1c). Three key observations provide support for this interpretation. Firstly, despite recharging from the same hydrological pathways, the recharge rates at Thw124 are six times higher than those at Thw142 and Thw170. Consequently, Thw124 must discharge less flux than its upstream counterparts. Secondly, during 2017 activity, Thw124 fills while the upstream lakes drain, contrary to all three lakes draining in 2013. This implies a modification to the hydrological system post-2013 activity. Finally, Malczyk and others (Reference Malczyk, Gourmelen, Goldberg, Wuite and Nagler2020) calculated the water budget for the network, demonstrating that the volume gain at Thw124 is the product of combined discharge from Thw142 and Thw170, along with annual melt production. This implies that Thw124 had to restrict a significant portion of downstream discharge during its recharge phase; otherwise, the water budget would not have closed. It is important to note that Malczyk and others (Reference Malczyk, Gourmelen, Goldberg, Wuite and Nagler2020) determined that the downstream flux at Thw124 must be at least 0.17 km3 a−1, indicating that the lake does not act as a hydrological roadblock but instead restricts the downstream flux. Nevertheless, the reduced ambiguity in the outward flux at Thw124 makes the Thwaites network ideal for comparing melt production against observations.
4.2. Cook E2
Cook E2 is the only subglacial lake where none of the modelled recharge rates can account for the observation-driven recharge rate, regardless of the routing approach and subglacial melt estimate. Our observation-driven recharge rates are 770% greater than the average of our modelled rates. The behaviour of lake activity at Cook E2 following the 2007 drainage event has been well documented. Li and others (Reference Li, Lu and Siegert2020) suggest a recharge rate of 0.052 km3 a−1 – which is within our estimate's uncertainty range (Table S6). In conjunction with Li's findings, our observations indicate that the lake is filling much more quickly than current estimates of subglacial melt production suggest. We have three hypotheses which could explain this disparity. Firstly, the modelled subglacial melt rates within the Cook catchment are underestimated. Secondly, the subglacial network or the catchment basin is incorrect, and more water must be mobilised towards the lake. Finally, another water source could be mobilised.
We can calculate the expected subglacial melt rate needed to account for our observation-driven recharge rates by dividing the recharge rate by the Cook E2 drainage basin size. This implies an average melt rate of at least 0.054 m a−1, 580% greater than current estimates. The required geothermal heat flux is estimated by rearranging Eqn (1) and forcing it with the average values for frictional heating and vertical conduction over the Cook E2 drainage basin. Doing so leads to an estimated heat flux of 526 mW m−2 (see Text S3 for a calculation breakdown), centralised over the lake. In contrast, Antarctica's largest documented heat flux is 285 ± 80 mW m−2 (Fisher and others, Reference Fisher, Mankoff, Tulaczyk, Tyler and Foley2015), implying that the elevated subglacial melt needed to account for our observations is highly implausible unless locally elevated geothermal heat fluxes are present within the deep trough that contains Cook E2 (Li and others, Reference Li, Lu and Siegert2020).
Increasing the spatial extent of the subglacial network feeding Cook E2 also comes with its issues. Cook E2 sits near an ice divide, so the network feeding the lake is topographically restricted. Under current subglacial melt rates, the network would have to increase in area by a factor of at least eight to account for the observation-driven recharge rates, which is unrealistic given the lake's position. A relative error of ~10% in the basal topography is unlikely to cause a significant change to where water is routed, as this would have been picked up in the stochastic D8 approach.
Li and others (Reference Li, Lu and Siegert2020) note that subglacial water could be mobilised from basal till layers contributing to the subglacial network (Christoffersen and others, Reference Christoffersen, Bougamont, Carter, Fricker and Tulaczyk2014; Siegert and others, Reference Siegert2017). This impact is difficult to assess and validate, but the location of the Cook E2 within a deep trough would facilitate the presence of large sediment deposits.
Finally, it is conceivable that a potential bias exists between the volumetric changes occurring at the bed and the corresponding surface response. Li and others (Reference Li, Lu and Siegert2020) utilised radio-echo sounding to ascertain that the area of Cook E2 spans ~46 km2, a considerably smaller measurement compared to our estimate of 103 km2. As such, the bias between lake area and surface response could lead to an overestimation of recharge.
An additional noteworthy observation concerning Cook E2 pertains a deviation from the recharge rate power-law association with lake volume (Fig. S15) as established by Livingstone and others (Reference Livingstone2022). This association posits that larger lakes generally experience greater recharge than smaller lakes.
4.3. Impact of routing approach
The choice of routing scheme is the primary source of variation within our modelled recharge rates. For example, at Lam80, the D8 approach predicts an average inflow of 4.28 km3 a−1, 3300% greater than the FD8 estimate of 0.13 km3 a−1. For comparison, the maximum range between varying the GHF within the region is 1.07 km3 a−1. The significant variations in predicted recharge rates between the approaches reflect where subglacial water is routed. Altering the routing approach can cause minor to medium-scale changes regarding where water is directed, significantly impacting the derived recharge rate. This is reflected in drainage pathways at Lam80, where under D8 the primary pathway flows directly through the lake, allowing recharge from water mobilised from the whole catchment, while under other schemes the flow path is diverted away from the lake (Fig. 7a). This results in the potential of extreme disparity in estimating recharge rates between methods.
Certain routing approaches exhibit limitations in accurately capturing observation-driven recharge rates within specific catchments. For instance, the D8 routing approach proves inadequate in predicting recharge rates for David S1, Thw70 and Thw170, as it diverts water around these features. Although the FD8 approach establishes hydrological connections among the three lakes within the SLC network, it fails to mobilise sufficient water from upstream sources to account for the observed volume gains.
In contrast, the stochastic D8 method emerges as the sole routing approach capable of accommodating the observation-driven recharge rates for all lakes across Antarctica, except for Cook E2, under all subglacial melting rate estimates. This effectively validates the lower bounds on all four GHF maps. The deterministic nature of the D8 and FD8 approaches renders them more susceptible to errors arising from inaccuracies in bed topography. This can result in minor to medium-scale variations regarding where water is routed. The stochastic D8 approach provides a means to address these uncertainties by introducing artificial variability and averaging over many simulations. This approach effectively negates the impact of these uncertainties in the routing. Incorporating a stochastic element into deterministic routing methods might be helpful for the future study of subglacial flow paths.
Our three approaches' disparity is unexpected, as the mechanical differences between D8 and FD8 schemes are not profound. In particular, the differences between the D8 and stochastic D8 methods are unexpected, as they both incorporate a D8 system. The primary difference between the two methods is that the stochastic D8 scheme includes a stochastic element by adjusting the hydropotential map over a predefined error field. This could imply that the underlying mechanics of the method does not link the differences between our three approaches – but is instead due to hydropotential map variation. Stochastic modelling in Greenland has illustrated that small changes within the hydropotential significantly impact subglacial flow paths (Mankoff and others, Reference Mankoff2020; MacKie and others, Reference MacKie, Schroeder, Zuo, Yin and Caers2021). Our results from running the stochastic D8 approach confirm this behaviour across Antarctica. We believe that small but significant changes to the hydropotential are introduced during the hole-filling process in each routing methods algorithm. As a result, each method produces a unique hydropotential map, which can result in differing large-scale water flow. An example of this behaviour can be observed within the Thwaites network, as shown by Figure S13. Under each scheme's respective hole-filling approach, the large-scale differences between the methods are profound: the FD8 scheme suggests a hydrological connection between all four lakes, while D8 implies a degree of disconnection. However, when both methods are routed over identical filled depression hydropotential maps, the flow paths and corresponding fluxes become nearly identical. This behaviour highlights the importance of small-scale changes in modelling subglacial flow paths. As such, we recommend a transition from deterministic methods that are suspectable to the propagation of errors to stochastic or physically based routing schemes.
4.4. Impact of geothermal heat flux
The differences in derived recharge rates based on the subglacial melt rate estimate using different GHF products can be substantial, with the larger estimate often being twice the lower estimate (Fig. 10). The water feeding the lakes is mobilised by collecting water over the catchment. For subglacial lakes in upstream regions, ice velocity beyond these lakes is low, so geothermal heat flux becomes the primary driver of subglacial melting. As such, minor changes to the heat flux lead to a significant increase in water being mobilised towards a lake when considering the vast area occupied by subglacial catchments.
Modelled recharge rates for the Thwaites network can vary significantly depending upon the subglacial melt rates forcing the FD8 routing scheme (Fig. 4). The subglacial melt rate estimate from Joughin and others (Reference Joughin2009) cannot account for the observation-driven recharge rates within the Thwaites network with a deficit of 0.04 km3 a−1. This subglacial melt product did not include estimates of elevated geothermal heat flux that were hypothesised to exist from ice-penetrating radar observations (Schroeder and others, Reference Schroeder, Blankenship, Young and Quartini2014). As such, the Joughin and others (Reference Joughin2009) melt rate was theorised to underestimate melt volumes (Smith and others, Reference Smith, Gourmelen, Huth and Joughin2017), a result confirmed by our analysis. This underestimation of subglacial melting rates could be a factor as to why modelled recharge rates derived by Smith and others (Reference Smith, Gourmelen, Huth and Joughin2017) could not account for rates derived from altimetry (Malczyk and others, Reference Malczyk, Gourmelen, Goldberg, Wuite and Nagler2020).
The subglacial melting product proposed by Van Liefferinge and Pattyn (Reference Van Liefferinge and Pattyn2013) presents a modelled recharge rate that precisely aligns with the observed rates within the Thwaites network, as depicted in Figure 4. This agreement substantiates the theoretical framework outlined in Figure 1c. Such a finding assumes that all annual subglacial water production within the ice stream is effectively mobilised towards the lakes with nothing lost to other hydrological sinks. Moreover, in this scenario, Thw124 must function as a perfect hydrological sink, impeding downstream discharge. Consequently, the absence of downstream discharge hinders channel growth during the recharge phase of the lake. This phenomenon could potentially explain the feature's resistance to activation during the volume transfer in the 2017 drainage event.
The discrepancy observed between the Van Liefferinge and Pattyn (Reference Van Liefferinge and Pattyn2013) melt estimates and the Maule and others (Reference Maule, Purucker, Olsen and Mosegaard2005) and Shapiro and Ritzwoller (Reference Shapiro and Ritzwoller2004) estimates is rather surprising, considering that the former incorporates an average of the Maule and others (Reference Maule, Purucker, Olsen and Mosegaard2005) and Shapiro and Ritzwoller (Reference Shapiro and Ritzwoller2004) GHFs (Pattyn, Reference Pattyn2010). As a result, we would anticipate a similar pattern of recharge in comparison to our own estimates. While the precise explanation for this disparity remains speculative, it is possible that it reflects the influence of adjusting geothermal heat flux values based on local topography, as highlighted by Colgan and others (Reference Colgan2021) or constraining stricter basal friction dissipation through model inversion.
4.5. Limitations of observing active subglacial lakes
We have collected high temporally constrained observations of subglacial lake activity across Antarctica using swath-processed surface elevation changes. However, our method prevents us from studying subglacial lake activity further inland due to CryoSat-2 SARIn data collection mode constraints. Swath processing is impossible under CryoSat-2 Low-Resolution Mode (LRM), and our time series approach diverges without regular temporal data. Therefore, we cannot constrain the behaviour of the other 96 known active lakes outside of the SARIn collection mode. It is noteworthy that some of these lakes exhibit hydrological linkages to the catchments we have examined (see Text S4). Under such circumstances, it becomes plausible to consider the possibility that unobservable upstream lakes may be undergoing recharge processes, thereby regulating the downstream flow of water. In this context, our current routing approach may inadvertently lead to an overestimation of the water supply directed towards an observable lake.
Moreover, knowledge of recharge activity at these upstream lakes would be valuable, as basal frictional heating has a diminishing impact further inland as ice flow reduces, so melt would be primarily constrained by geothermal heat flux. Comparing modelled against observation-driven recharge rates for such lakes might allow for discriminating different geothermal heat flux products.
The upcoming radar altimeter CRISTAL mission (Kern and others, Reference Kern2020) will provide SARIn mode coverage over the entire Antarctic Ice Sheet, thus providing complete coverage of subglacial lake activity across the ice sheet at high spatial resolution. In addition, as one of Copernicus' Sentinel missions, CRISTAL will provide the community with long-term monitoring of the subglacial environment. Other missions, such as IceSat2 (Markus and others, Reference Markus2017), high-resolution stereo imagery such as that used in REMA (Howat and others, Reference Howat2022) or InSAR missions such as Sentinel 1 or TanDEM-X, can be used to sense surface expression of active subglacial lakes, and as such can be part of an subglacial lake observing system strategy (e.g. Sandberg Sørensen and others, Reference Sandberg Sørensen2023).
Subglacial melting has primarily been inferred by models and rare point measurements at drilling locations (i.e. Lukin and Vasiliev, Reference Lukin and Vasiliev2014; Tulaczyk and others, Reference Tulaczyk2014; Priscu and others, Reference Priscu2021). Our method presents a unique and novel way to explore processes in the subglacial environments and relate remote-sensing observations to subglacial melting. However, our approach has been unable to distinguish between the four different GHF estimates. The potential exists, but some assumptions and limitations concerning our modelled and observation-driven recharge rates must first be addressed. Uncertainty within the subglacial network is the largest source of variation for our modelled recharge rates. Different routing approaches can result in large-scale flow path change in regions with high relative bed uncertainty (i.e. Thwaites), which results in significant uncertainty in recharge rates. If the bed topography can be better constrained to prevent these large-scale flow path divergences, our modelled recharge rates will be significantly more accurate.
Our primary source of uncertainty arises from our assumption of a one-to-one relationship between surface elevation changes and volumetric changes at the bed – an invalid assumption as shown by Sergienko and others (Reference Sergienko, MacAyeal and Bindschadler2007). The impact of this assumption is that our recharge rates underestimate the actual filling rate, providing a weaker constraint on melt production than if this relationship could be established. We recommend that research be performed to determine the connection between changes at the bed and the corresponding response on the ice-sheet surface. Such a relationship is vital in reducing the uncertainty attached to altimetry-derived time-dependent volume changes and recharge rates. By gaining a deeper understanding of the interactions between subglacial processes and their surface expressions, we can enhance the accuracy and reliability of future assessments and further constrain subglacial melt production.
5. Conclusion
Here we have used observations of surface elevation change from CryoSat-2 satellite altimetry to study the drainage cycles of active subglacial lakes in Antarctica and quantify rates of recharge. We also derive expected rates of recharge using a range of modelled subglacial melt and routing algorithms. We have detected 13 lakes that exhibit evidence of recharge following drainage events, from which we have extracted recharge rates and compared these with modelled values.
Under a stochastic routing approach, the modelled estimates successfully account for the observed values for each studied lake, validating the lower bound of four geothermal heat flux estimates. There is an exception for Subglacial Lake Cook E2, where the modelled recharge rates are notably lower than those derived from altimetry data. This suggests the possibility of either enhanced melt or the release of groundwater from basal till layers within the region.
We demonstrate that considering subglacial lakes as a network rather than individually provides a better constraint on melt production. The Thwaites Lake network is a valuable system for validating subglacial melting rates, as Thw124 limits downstream discharge during its recharge phase. Modelled recharge rates cannot account for observed recharge under the Joughin and others (Reference Joughin2009) estimate. This suggests that subglacial melt production in the region may be greater than previously thought.
It is important to note that constraining melt production using active subglacial lakes has its limitations. Non-hydrological factors have the potential to contribute to the filling of surface depressions, adding a degree of uncertainty to our altimetry-derived recharge rates. Moreover, we cannot document inland active lakes under CryoSat-2. Therefore, we recommend further research to develop methods for establishing accurate volume changes at the bed during lake activity and the ongoing documentation of drainage activity of inland active lakes.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2023.70.
Acknowledgments
This work was performed at the University of Edinburgh under grants from the European Space Agency's project 4DAntarctica (ESA: Grant 4000128611/19/I-DT), and from the PROPHET project, a component of the International Thwaites Glacier Collaboration (ITGC). This work received support from National Science Foundation (NSF: Grant 1739031) and Natural Environment Research Council (NERC: Grants NE/S006745/1, NE/S006796/1, and NE/T001607/1) (ITGC Contribution No. 024). G. R. M. acknowledges a NERC PhD Studentship (NERC: Grants NE/L002558/1). The authors wish to thank ESA for providing open access to CryoSat-2 data and M. Morlighem for open access to Bedmachine. The authors are grateful to the editor, Frank Pattyn, and to three anonymous reviewers, whose comments have significantly improved the manuscript.
Data and code availability
Our lake masks, elevation and volume change time series are freely available from 4D Antarctica (https://4d-antarctica.org/products/). CryoSat-2 swath data are available through the CryoTEMPO EOLIS product (https://cryotempo-eolis.org/product-overview/) and via the cs2eo data platform (https://cs2eo.org/). 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/3). The D8 routing algorithm is freely available on GitHub (https://github.com/wschwanghart/topotoolbox). The stochastic D8 routing algorithm is freely available on GitHub (https://github.com/mauro3/WhereTheWaterFlows.jl).
Author contributions
G. M. and N. G. designed the work. G. M. analysed the results with input from N. G. and M. Werder. M. Werder provided code for the stochastic routing model. M. Wearing and D. G. provided data for heating from basal friction and vertical heat conduction. G. M. wrote the manuscript, with input from all co-authors.
Open access
