Brief communication: Heterogenous thinning and subglacial lake activity on Thwaites Glacier, West Antarctica

A system of subglacial lakes drained on Thwaites Glacier from 2012–2014. To improve coverage for subsequent drainage events, we extended the elevation and ice-velocity time series on Thwaites Glacier through austral winter 2019. These new observations document a second drainage cycle in 2017/18 and identified two new lake systems located in the western tributaries of Thwaites and Haynes glaciers. In situ and satellite velocity observations show temporary < 3 % speed fluctuations associated with lake drainages. In agreement with previous studies, these observations suggest that active subglacial hydrology has little influence on thinning and retreat of Thwaites Glacier on decadal to centennial timescales.

1 Introduction

Although subglacial lakes beneath the Antarctic Ice Sheet were first discovered more than 50 years ago (Robin et al., 1969; Oswald and Robin, 1973), they remain one of the most enigmatic components of the subglacial hydrology system. Initially identified in ice-penetrating radar data as flat, bright specular reflectors (Oswald and Robin, 1973; Carter et al., 2007), subglacial lakes were thought to be relatively steady-state features of the basal hydrology system with little impact on the dynamics of the overlying ice on multi-year timescales. The advent of high-precision repeat satellite observations in the 1990s, however, revealed an entirely new class of active subglacial lakes that fill and drain on annual to decadal timescales and possibly affect the flow of the overlying ice (e.g., Gray et al., 2005; Wingham et al., 2006; Fricker et al., 2007; Smith et al., 2009). Under the central trunk of Thwaites Glacier, in particular, satellite radar altimetry revealed a large (∼ 4 km³ volume discharge), connected subglacial lake drainage event from 2012–2014 (Smith et al., 2017). Initial subglacial lake recharge estimates suggested lake drainages of this magnitude should occur every 20–80 years. To better constrain refill and discharge time, we extended the Thwaites Glacier velocity and altimetry record to include the most recent drainage events. We also expanded the spatial coverage to include Haynes Glacier and the western tributary of Thwaites Glacier (Fig. 1). Here we describe the recent subglacial lake behavior in these regions and discuss the impact of these subglacial lake systems on slip velocity at the ice–bed interface.

Figure 1 Location map of Thwaites Glacier and subglacial Thwaites lakes. (a) Average ice speed between 2015–2019 omitting period when lakes were active (color) plotted over Moderate Resolution Imaging Spectroradiometer (MODIS) image mosaic (Haran et al., 2014). Thwaites Glacier, Thwaites Lake 124 (Thw₁₂₄), Thwaites Lake 142 (Thw₁₄₂), Thwaites Lake 170 (Thw₁₇₀), Haynes Glacier (HG) lake, western Thwaites (WT) lake and GNSS sites (LTHW and UTHW) are labeled. Thwaites lakes are named by their approximate distance from the grounding line. (b) LTHW and (c) UTHW GNSS position plotted over time (color) with contoured mean velocity between 2015–2019.

2 Methods

We used in situ Global Navigation Satellite System (GNSS), satellite synthetic aperture radar (SAR) and satellite radar altimetry data to derive velocity and elevation-change time series.

2.1 Ice velocity and vertical displacement

We used the speed anomalies recorded by two long-term on-ice GNSS receivers, LTHW and UTHW, deployed on Thwaites Glacier from 2009 to present to augment Eulerian velocity products derived from Sentinel-1A and Sentinel-1B synthetic aperture radar (SAR) imagery of the Thwaites Glacier catchment collected between 2015–2019 (Fig. 1). Processing of GNSS data follows the workflow of Christianson et al. (2016). We first created a position time series with sub-5 cm uncertainties in all dimensions by calculating geodetic positions every 30 s relative to two fixed rock sites located 200 km (Backer Island) and 300 km (Howard Nunatak) away from the LTHW site, using differential carrier phase positioning as implemented in the Track Software (Chen, 1998). We then constructed velocity time series from these geodetic solutions using 3 d Savitzky–Golay filtered moving averages (Press et al., 2007). Finally, we subtracted the SAR-derived Eulerian speed at each GNSS position to solve for the Lagrangian velocity anomaly relative to the mean 2015–2019 velocity, omitting periods of lake activity (reference velocity field in Fig. 1).

Figure 2 Surface elevation-change time series over the Thwaites Glacier lakes showing the 2017 drainage cascade from (a) vertical displacement computed from integrated vertical displacement rates (V_z) from Sentinel-1 SAR data and (b) swath-processed radar altimetry in a polar stereographic projection (EPSG:3031). Water volume (km³) associated with observed vertical displacement is labeled for each lake. (c) Time series of uplift rates (V_z) from SAR LOS results (colored dots, left abscissa); locations marked in (a) and (b) and horizontal speed from GNSS observations (right abscissa). Solid lines represent period over which SAR vertical displacements (V_z) were integrated to produce the vertical displacements shown in (a). Dotted lines represent the quarters of gridded CryoSat-2 data differenced to create (b).

Following the methods of Joughin (2018, updated 2019) and Joughin et al. (2010), we constructed a speckle-tracked velocity time series of Thwaites Glacier from 2015 through the austral winter of 2019 using SAR data collected by the European Union's Copernicus Sentinel-1A and Sentinel-1B satellites and processed by the European Space Agency (ESA). We also computed the component of motion in the satellite line-of-sight (LOS) direction (Gray et al., 2005; Friedl et al., 2020). The bulk of this signal is due to relatively steady horizontal displacements, but it is also influenced by potentially more temporally variable vertical displacements. Thus, we computed the mean LOS component for the full time series and subtracted it from each individual estimate. Since the horizontal and surface-parallel flow components are relatively steady, the residual line-of-sight estimate should largely be due to vertical motion, which we corrected for incidence angle effects to produce an approximate vertical displacement rate (V_z). Because we only subtracted the mean, a small component of the horizontal velocity may map into the vertical velocity due to the glacier's acceleration; however, this contribution is generally in the noise during periods with no lake activity and small relative to the vertical signal during times of active lake change (Fig. 2). To more tightly constrain the timing of the drainage events, we spatially interpolated the time series of Sentinel-1-derived V_z to fill gaps in coverage and integrated the result during a period of filling/draining (see solid vertical bars in Fig. 2c) to produce estimates of net uplift and subsidence shown in Fig. 2a.

2.2 Elevation and lake volume change

We also extended the previous time series of ESA CryoSat-2 radar altimetry data (Smith et al., 2017) through austral winter 2019, as shown in Figs. 2b and 3. Elevation models were derived by fitting surfaces of elevation change to CryoSat-2 swath-processed elevation retrievals and points-of-closest-approach relative to a reference elevation model from the first quarter of 2011 (Smith et al., 2017). The fitting procedure minimized an objective functional that considered data misfit, spatial gradients in the constructed reference elevation model, elevation-change rate fields, temporal gradients in elevation-change rate and the magnitude of model bias parameters. In this scheme, three expected elevation statistics are used to choose weight parameters that regularize the least-squares fit. The elevation statistics, $E\left(\frac{{\partial }^{\mathrm{2}}{z}_{\mathrm{0}}}{\partial x^{\mathrm{2}}}\right)$, $E\left(\frac{{\partial }^{\mathrm{3}}z}{\partial x^{\mathrm{2}}\partial t}\right)$ and $E\left(\frac{{\partial }^{\mathrm{2}}z}{\partial t^{\mathrm{2}}}\right)$, represent expected values for spatial and temporal derivatives of the reference elevation model, z₀ and the time-dependent height-change field, z. The values chosen for this study are $E\left(\frac{{\partial }^{\mathrm{2}}{z}_{\mathrm{0}}}{\partial x^{\mathrm{2}}}\right)=\mathrm{6.7}\times {\mathrm{10}}^{-\mathrm{8}}$, $E\left(\frac{{\partial }^{\mathrm{3}}z}{\partial x^{\mathrm{2}}\partial t}\right)=\mathrm{6}\times {\mathrm{10}}^{-\mathrm{9}}\mathrm{m}{\mathrm{yr}}^{-\mathrm{1}}$ and $E\left(\frac{{\partial }^{\mathrm{2}}z}{\partial t^{\mathrm{2}}}\right)=\mathrm{1.0}{\mathrm{m}}^{\mathrm{2}}{\mathrm{yr}}^{-\mathrm{2}}$ and tighten the spatial variations in the least-squares elevation time series, $E\left(\frac{{\partial }^{\mathrm{2}}{z}_{\mathrm{0}}}{\partial x^{\mathrm{2}}}\right)$ and $E\left(\frac{{\partial }^{\mathrm{3}}z}{\partial x^{\mathrm{2}}\partial t}\right)$, compared to the original Smith et al. (2017) implementation by factors of 5 and 10, respectively. These radar altimetry measurements complement SAR observations of integrated vertical displacement, which we use together to understand the character of new lake drainage activity.

Figure 3 Time series of GNSS velocity anomalies at UTHW and LTHW corrected for advection using the Eulerian velocity products and CryoSat-2 lake elevation change averaged over each lake area. See Fig. 1 for site locations and abbreviations. Also plotted are LTHW GNSS clockwise direction change relative to 2010 flow direction (purple). The dark grey shaded periods indicate intervals when the LTHW GNSS accelerated significantly (99 % confidence), while the light grey periods indicate when the Thw${}_{\mathrm{124},\mathrm{142},\mathrm{170}}$ lakes are active. When the largest lake fills in 2017, the LTHW GNSS closest to the lake accelerates and flows towards the lake.

3 Results: new observations of lake activity

A complete chronology of progressive thinning and lake activity across Thwaites Glacier from the extended CryoSat-2 time series is shown in the video in the Supplement (Movie SV1). These new observations reveal that the upper Thwaites Lakes, Thw₁₇₀ and Thw₁₄₂, drained in 2017, filling Thw₁₂₄ (Figs. 1 and 2). The SAR-derived elevation-change data show that the largest lake, Thw₁₂₄, filled by 1.9 km³ during the 2017 drainage, roughly balancing the volume that drained from Thw₁₄₂ (0.6 km³) and Thw₁₇₀ (1.4 km³). The quarterly CryoSat-2 results show less clear evidence of water budget balance (Figs. 2 and S2), which may be due to the degree of smoothing used in producing the time series. From CryoSat-2 elevation-change data, between 2015 and before the 2017 drainage event, the areas inside the Thw₁₂₄, Thw₁₄₂ and Thw₁₇₀ lake outlines increase in elevation, which is strong evidence of filling (Figs. 3 and S2).

The extended elevation time series (Figs. 3 and S2) also reveals the fill-drain cycle of two new lake systems: one in the western shear margin of Haynes Glacier and another in the western tributary of Thwaites Glacier (Fig. 1). From these combined observations, the western Thwaites tributary lake (WT) drained by 1.1 km³ in 2013 and the Haynes Glacier lake system (HG) drained by 0.2 km³ in 2017 (Figs. S1 and S2 in the Supplement). Complete fill-drain cycles of the Haynes Glacier lakes and the western Thwaites Glacier lake are not observed in the existing altimetry record, and the Haynes Glacier lakes do not discernibly refill after draining in 2017 (Fig. S2). The western Thwaites tributary lake, however, fills significantly at a rate of 0.1 km³ yr⁻¹ after draining in 2014.

4 Discussion

Cascading lake drainages have been observed under many Antarctic ice–stream systems (e.g., Wingham et al., 2006; Fricker et al., 2007). The positions of all identified lakes beneath Thwaites Glacier, including the new lakes in the Haynes Glacier shear margin and western tributary of Thwaites Glacier, appear to be controlled primarily by the bed and associated surface geometry (Smith et al., 2017). There are large topographic ridges at the bed with corresponding expressions at the surface that are oriented orthogonal to flow and likely act as hydraulic baffles trapping water and sediments (Holschuh et al., 2020), causing hundred-kilopascal-scale deviations in basal traction (Joughin et al., 2014; Fig. S5a in the Supplement). The weak basal shear stress in these till-draped basins combines with large-scale catchment topography to promote variations in ice thickness and surface slope that form large hydropotential lows (Fig. S4 in the Supplement; Smith et al., 2017; Holschuh et al., 2020). In these hydropotential lows, the lakes remain disconnected from their neighbors as they fill until cascading drainages driven by the upstream lakes interrupt the background fill rate in the cycle. Densely sampled SAR vertical displacement rates from 2017 (V_z in Fig. 2c) demonstrate this process, capturing the Thw₁₇₀ drainage that initiated a combined drainage with Thw₁₄₂ into Thw₁₂₄ (Fig. 2).

The controls on lake filling are less clear. From the altimetry observations of the Thw₁₇₀ fill cycle, the average fill rate is ∼ 0.16 km³ yr⁻¹ (Fig. S4). This agrees with the fill rate (∼ 0.14 km³ yr⁻¹) we calculate by routing modeled basal meltwater production (Joughin et al., 2009) down the glaciostatic hydropotential gradient (Shreve, 1972) into Thw₁₇₀ but requires inflow of all melt water produced upstream into the Thw₁₇₀ lake basin (Fig. S4). The glaciostatic hydropotential also routes water around Thw₁₇₀ into downstream lakes Thw₁₄₂ and Thw₁₂₄, but the fill rates associated with these flow paths (∼ 0.44 and ∼ 0.27 km³ yr⁻¹, respectively) are much larger than the fill rates derived from surface height change (Fig. S4). These discrepant observations may reflect limitations of the static hydropotential assumption and the modeled catchment meltwater budget, but they also suggest current bed-elevation models do not resolve small-scale (< 1 km) bed topography important for routing subglacial water into the upper most lake, Thw₁₇₀.

4.1 Lake impact on ice flow and coupled drainage morphology

The inland SAR and GNSS observations show a general pattern of acceleration at the LTHW and UTHW sites, consistent with an increase in driving stress due to inland propagation of thinning caused by ungrounding and loss of ice-shelf buttressing (Rignot et al., 2014; Joughin et al., 2014). Figure 3 shows the velocity anomaly at LTHW after subtracting the 2010 glacier velocity (340 m yr⁻¹). At LTHW, this secular trend is punctuated by two signals associated with the Thwaites Lakes drainage events, first in September 2012 and again in May 2017 (Fig. 3). During the 2012 drainage documented by Smith et al. (2017), surface velocities initially spike by 2 % over a several-day period but then decline by 3 % over the following 6 months. Loss of receiver power interrupted this record in March 2013. When the receiver began telemetering data again in 2015, the relative change in position suggests the speed anomaly in 2014 remained below what would be expected from the 2010–2012 trend. From January 2016 to May 2017, the LTHW receiver continued to accelerate at a rate of 4 m yr⁻². As Thw₁₄₂ and Thw₁₇₀ drained in the austral winter 2017, filling Thw₁₂₄, there was a nearly stepwise 1 % increase in glacier speed at LTHW (Fig. 3). Speed remained elevated at the LTHW site after Thw₁₂₄ stopped filling, coinciding with a 2^∘ shift in ice-flow direction to the grid-north (clockwise), toward Thw₁₂₄ (Fig. 3); however, this speed change is imperceptible in distributed velocity maps before and after Thw₁₂₄ filled in 2017 (Fig. S3 in the Supplement). The UTHW site also exhibits subtle velocity fluctuations with a magnitude less than 1 % of the mean velocity (105 m yr⁻¹) after correcting for the spatial ice-velocity gradient. These fluctuations are small relative to the background acceleration we observe at the UTHW site, 0.75 m yr⁻².

The transient velocity anomalies that depart from the trend observed at LTHW are likely related to subglacial lake dynamics. During lake filling in 2017, the areal extent of the Thw₁₂₄ lake increases. This increase should reduce traction at the lake boundary, both at the margins (near the LTHW site) of the lake and inside the lake as unresolved topographic pinning points are submerged. In 2017, changes in basal traction occur almost immediately, causing a step increase in velocity (Fig. 3). As the lakes drain, traction is restored as ice regrounds, reducing basal slip, which we observe in the months following the 2012 lake drainage near Thw₁₂₄.

The positive GNSS acceleration observed in September 2012 (Fig. 3) suggests that Thw₁₂₄ began to drain in 2012 before the quarterly resolved subsidence associated with the lake drainage becomes distinguishable in the CryoSat-2 surface elevation time series (Fig. 3). The cause of the acceleration we observe at Thw₁₂₄ before the bulk drainage cannot be unambiguously attributed to a discrete set of processes with these data, but the finite duration of the velocity increase (∼ 10 d) suggests that a distributed drainage system may have been established at the downstream edge of the lake preceding bulk drainage. Glaciostatic hydraulic water routing indicates that Thw₁₂₄ would likely drain at the grid-south edge of the lake near the LTHW GNSS receiver (Fig. S3). Enhanced lubrication outside the low-drag Thw₁₂₄ basin as the lake begins to drain likely increases local slip and drives the subtle change in ice-flow direction that we observe in the austral winter of 2012 before the peak drainage in 2013, when flow direction shifts back to the mean flow direction between 2010–2012. Non-steady effective pressure likely also affects the basal shear stress as the drainage system initially forms, then empties and closes. Similar to lake drainages under the Siple Coast ice streams (Siegfried et al., 2014, 2016), changes in basal slip speed due to lake activity are small (∼ 1 %–2 %) relative to the average sliding speed.

Regardless of the exact drainage mechanism, the locations of the lakes are governed by ice-flow response to the underlying bed topography, which promotes hydropotential basins that form as ice flows over ridges. Once a connected drainage begins, differences in water pressure in the conduits between lakes promote efficient drainage down the hydraulic gradient. During drainage, each lake is likely in local equilibrium with the ends of the conduits that directly connect to it; however, the large hydropotential differences between lakes (∼ 1 MPa) cannot equilibrate along the entire length of the drainage path, which implies that the pressure difference over the length of the conduit is likely more important in determining whether water flows into or out of the lakes than the small variations in the hydraulic potential in the lakes as they fill and drain. Conversely, drainage of upstream lakes may disrupt the steady-state drainage morphology by temporarily increasing the hydraulic conductivity of the conduits, bringing additional water into the lower lakes and allowing drainage from adjacent lakes that share the same downstream conduit. The processes governing changes in subglacial hydraulic connectivity are poorly understood but likely depend on the local dynamic hydropotential and the evolution of conduit morphology, which we do not try to infer from our observations of the fill-drain cycle. We only note that these processes likely contribute to the variability in lake fill-drain levels observed elsewhere over multiple fill-drain cycles (Fig. S2.; Siegfried et al., 2014, 2016), which cannot be explained by the evolving glaciostatic hydropotential alone.

4.2 Implications for basin-wide change

The velocity changes in Fig. 3 that we attribute to lake drainage events represent the dominant, albeit small (≲ 3 %), inland sub-annual velocity variability. The lakes sequester water and thus likely play some not well-understood role in maintaining distributed low effective pressures that control long-term flow rates, but the fill-drain cycles have little transient effect on the flow behavior of the overlying ice. On decadal timescales important for understanding ice-sheet behavior and contribution to sea level, the lakes do not appear to control inland ice-flow variation. Static inversions for bed resistance before and after the 2017 lake drainage event (see Sect. 4 in the Supplement for details) are not sensitive to the subtle surface velocity changes we measure with GNSS (Fig. S5), and the lakes identified by Smith et al. (2017) have no discernible effect on ice velocity at the UTHW site. These new observations suggest that the observed speedup at the grounding zone of the main trunk of Thwaites Glacier following the 2013 drainage (Smith et al., 2017) was associated with warming ocean conditions following anomalous Amundsen Sea wide ocean cooling from 2012/13 (Christianson et al., 2016). These warm ocean conditions likely enhanced sub-ice-shelf melt and led to increased ungrounding and acceleration. Our observations and model experiment (Sect. 4 in the Supplement; Shapero et al., 2020) invalidate proposed geoengineering solutions that seek to drain large volumes of water from beneath Amundsen Sea Embayment glaciers to increase basal resistance (Moore et al., 2018). These results further demonstrate that capturing the details of lake fill-drain cycles, and at least some elements of the associated basal hydrology system, may not be that important for modeling Thwaites Glacier's contribution to sea level on decadal to centennial timescales.

5 Conclusions

We document the temporal change in velocity and elevation far from the grounding zone in response to the steepening of Thwaites Glacier and three distinct systems of active lakes: one on the main Thwaites Glacier trunk, another in the western shear margin of Haynes Glacier and a third in the westernmost tributary of Thwaites Glacier. At the LTHW GNSS site, over 100 km from the grounding line, ice velocity has accelerated at a nearly constant rate over the last decade. This background acceleration was interrupted in 2012 by the connected drainage of lakes Thw₁₂₄, Thw₁₄₂ and Thw₁₇₀ and, in 2017, by the partial filling of Thw₁₂₄ via drainage of Thw₁₄₂ and Thw₁₇₀. Our observations suggest that the transport of ∼ 2 km³ of water beneath Thwaites Glacier, which represents approximately half the annual basal meltwater production for the entire Thwaites catchment (Joughin et al., 2009), has only a small and transient effect on glacier speed relative to ongoing thinning driven by ocean melt.