Surface meltwater drainage and ponding on Amery Ice Shelf, East Antarctica, 1973–2019

Abstract

Surface melting on Amery Ice Shelf (AIS), East Antarctica, produces an extensive supraglacial drainage system consisting of hundreds of lakes connected by surface channels. This drainage system forms most summers on the southern portion of AIS, transporting meltwater large distances northward, toward the ice front and terminating in lakes. Here we use satellite imagery, Landsat (1, 4 and 8), MODIS multispectral and Sentinel-1 synthetic aperture radar to examine the seasonal and interannual evolution of the drainage system over nearly five decades (1972–2019). We estimate seasonal meltwater input to one lake by integrating output from the regional climate model [Regional Atmospheric Climate Model (RACMO 2.3p2)] over its catchment defined using the Reference Elevation Model of Antarctica. We find only weak positive relationships between modeled seasonal meltwater input and lake area and between meltwater input and lake volume. Consecutive years of extensive melting lead to year-on-year expansion of the drainage system, potentially through a link between melt production, refreezing in firn and the maximum extent of the lakes at the downstream termini of drainage. These mechanisms are important when evaluating the potential of drainage systems to grow in response to increased melting, delivering meltwater to areas of ice shelves vulnerable to hydrofracture.

1. Introduction

Surface meltwater drainage systems have been regularly observed in summer on the Antarctic ice shelves since the exploration of these regions began (Mellor and McKinnon, 1960; Bell and others, 2017). Meltwater flows laterally, creates channels and fills surface depressions (Winther and others, 1996; Lenaerts and others, 2017; Kingslake and others, 2017; Dell and others, 2020) and in some cases drains englacially into or through the ice shelf (e.g. Bell and others, 2017; Dow and others, 2018; Stokes and others, 2019; Dunmire and others, 2020; Schaap and others, 2020) or into the firn (e.g. Hubbard and others, 2016; Montgomery and others, 2020). The change in local load from water movement and ponding on the shelf can cause flexural stresses large enough to generate fractures (MacAyeal and others, 2015). Meltwater can flow into and enlarge these fractures, a process referred to as hydrofracturing (Weertman, 1973; Rack and Rott, 2004; van der Veen, 2007; MacAyeal and others, 2015; Banwell and others, 2019; Robel and Banwell, 2019; Lai and others, 2020). These processes have been linked to increased meltwater production that caused the collapse of several Antarctic Peninsula ice shelves (Scambos and others, 2004; Banwell and others, 2013). Kuipers Munneke and others (2014) proposed the refreezing of meltwater in firn pore space as a precursor to ice-shelf collapse through its impact on near-surface permeability and the ability of water to form supraglacial lakes. Based on a relationship between surface melting and summer air temperature derived from satellite data and a regional climate model, Trusel and others (2015) predict a doubling in Antarctic surface melting by the end of this century. As surface drainage systems relocate water long distances across ice shelves, it is important to identify the controls on their multi-year growth under warming conditions.

Multi-year changes in meltwater lakes have been studied in several locations around Antarctica over the satellite period (1972–present). Surface depressions refill with meltwater in austral summer and then drain englacially, into or through the ice shelf, or refreeze in austral winter (Leppäranta and others, 2013; Langley and others, 2016). Satellite observations highlight large interannual variability in the maximum size of drainage systems (e.g. Liang and others, 2012; Langley and others, 2016; Lenaerts and others, 2017; Dell and others, 2020; Banwell and others, 2021)). It is unclear if this variability can be explained by or is linked to seasonal meltwater production or whether more complex glacio-hydrological coupling is at play. For example, large refrozen ice masses have been found advected downstream of lakes on Larsen C Ice Shelf (Hubbard and others, 2016) and Amery Ice Shelf (AIS) (Phillips, 1998), reducing firn permeability over large areas. If this encourages subsequent surface meltwater transport, as has been observed in Greenland (MacFerrin and others, 2019), complex, multi-year coupling between meltwater production, drainage and refreezing is possible.

To examine the potential for such coupling, we present a detailed examination of the seasonal and interannual evolution of meltwater drainage on AIS, East Antarctica, an ice shelf with a large surface drainage system that forms regularly. We document the growth and evolution of the drainage system using multispectral and synthetic aperture radar (SAR) satellite imagery covering nearly five decades from 1972 to 2019. We focus on seasonal and interannual evolution – including the production, drainage and wintertime refreezing of meltwater – as well as the comparison between observed drainage system extent and melt production estimated from a regional climate model from January 1979 to December 2018. We discuss possible explanations for the observed low correlation between modeled melt production and drainage system extent and their potential implications for the future evolution of meltwater drainage systems around Antarctica.

2. Study area

AIS (Fig. 1) is the largest ice shelf in East Antarctica, measuring 62,620 km² in area (Foley and others, 2013) and draining 16% of the East Antarctic ice sheet (Fricker and others, 2002). Located from 73°S to 69°S in the Prydz Bay embayment between the Prince Charles Mountains and Princess Elizabeth Land, AIS is further north than most Antarctic ice shelves. AIS experienced little to no net thickness change and little velocity change between 1968 and 2007, but experienced one period of abrupt changes in surface elevation (1996–2003) potentially linked to surface melting and firn densification (King and others, 2009). Between 1999 and 2017, the Australian Antarctic Division maintained an automatic weather station on AIS (Amery G3, Allison and Heil (2001); Fig. 1; 70° 53′31″S, 69° 52′21″E) that made hourly observations. Over this period of observation, the mean summer (December–February) 4 m (above the surface) temperature was − 8.6°C, and taking the maximum temperature of each summer and finding the average of these annual maxima, we estimate a mean annual-maximum 4 m temperature of 4.7°C.

Fig. 1. Landsat 8 true-color (Bands 2, 3 and 4) image of AIS, from 17 January 2019. The CTL is indicated by the red box. The two main nunataks referred to in the text (Clemence Massif and Fisher Massif), the location of the Amery G3 Automatic Weather Station, and the location of two observed englacial drainage events are labeled. Grounding lines are shown in magenta (Depoorter and others, 2013). The inset indicates the location of AIS in East Antarctica, image from NASA's Blue Marble: Antarctica.

Many studies have reported surface meltwater streams and lakes on AIS (Mellor and McKinnon, 1960; Swithinbank, 1988; Hambrey and Dowdeswell, 1994; Phillips, 1998; Fricker and others, 2002; Glasser and Scambos, 2008; Fricker and others, 2009; Kingslake and others, 2017; Stokes and others, 2019; Dirscherl and others, 2020; Moussavi and others, 2020; Dell and others, 2020; Lai and others, 2020), with the earliest of these recounting observations of meltwater from early 20th century explorers (Mellor and McKinnon, 1960). Phillips (1998) describes the radar altimeter waveform and radar backscatter data from the austral summers of 1992/93 and 1993/94, showing evidence of refrozen meltwater in a surface trough aligned with ice flow. Multi-spectral satellite imagery later showed that this trough hosts the largest ephemeral lake on AIS (e.g. Kingslake and others, 2017). Due to its location at the downstream terminus of an extensive drainage network that forms in AIS's central drainage catchment, we refer to this lake as the ‘Central Terminal Lake (CTL)’ (Fig. 1). Multiple such individual drainage systems form regularly, many of which terminate in their own large terminal lakes. Despite the long history of observations, there has been no in-depth assessment of the interannual variability of AIS's surface meltwater system.

3. Data and methods

3.1. Datasets

We used visible (optical) and SAR satellite imagery, satellite-derived products and model output. We list these datasets below.

3.1.1. Satellite imagery

We used a variety of satellite data for various time periods: (i) optical images from MODIS Terra/Aqua 2002–19, courtesy of NASA's Level-1 and Atmosphere Archive and Distribution System Distributed Active Archive Center (LAADS DAAC). (ii) Landsat 1, 4, 5, 7 and 8 optical images from 1972 to 2019, courtesy of the U.S. Geological Survey (Table S1). There were no usable Landsat 1, 4 and 5 images from 1975 to 1987, 1991, 1993 to 2000; in some years only one usable, low cloud-cover Landsat image was available. (iii) WorldView optical images from 2009 to 2019, courtesy of the Polar Geospatial Center (images © 2019 DigitalGlobe, Inc.). Figure 2 depicts the temporal availability of these data.

Fig. 2. The temporal coverage of usable, low-cloud cover (images with enough visible ice shelf surface areas to determine the presence of meltwater ponding) Landsat 1, 4–8 (red), MODIS Terra/Aqua (orange), and Sentinel-1 SAR (yellow) imagery, as well as the years with observations of melt lakes on AIS (light blue), and the observations of CTL (dark blue, see Fig. 1).

3.1.2. Satellite-derived products

We used the Reference Elevation Model of Antarctica (REMA), a mosaicked DEM with an 8 m spatial resolution created using stereographic pairs of WorldView imagery (Howat and others, 2019).

We analyzed a time series of HH-polarization SAR imagery over AIS from 2015 to 2019 using the Google Earth Engine Sentinel-1 Ground Range Detected SAR product, which is the European Space Agency's (ESA's) Copernicus Programme's C-band Sentinel-1 imagery, courtesy of the ESA, pre-processed to remove noise and converted from backscatter intensity to decibels (Gorelick and others, 2017). These images were downloaded at 10 and 30 m resolution as GeoTIFFs (code courtesy of Stef Lhermitte) (Table S2).

3.1.3. Atmospheric model output

We complemented these image-based datasets with output from the Regional Atmospheric Climate Model (RACMO 2.3p2; van Wessem and others, 2018). RACMO 2.3p2 is forced by the climate reanalysis product ERA-Interim to simulate atmospheric circulation and the surface energy balance across Antarctica from January 1979 to December 2018.

3.2. Mapping surface melt extent and drainage

The locations of the downstream extents of the terminal lakes vary interannually by tens of kilometers in satellite imagery between 1973 and 2019. We refer to the location of the downstream margin of the CTL, expressed as a latitude, as the ‘maximum downstream extent’ of the drainage system. We defined time periods as follows: spring/early summer (November–December), midsummer (January–February) and end-of-summer/early fall (February–March); since austral summers begin in one calendar year and end in the next, we refer to summers by the year they began followed by the year they terminated (i.e. 1991/92 or 2014/15).

3.2.1 Surface catchments

We used REMA to estimate surface drainage catchments on AIS. After resampling REMA to 40 m (ignoring voids), we used the MATLAB package TopoToolbox (Schwanghart and Scherler, 2014) to calculate drainage catchments in the area of meltwater ponding and to isolate the central catchment, which contains the CTL (Figs 1, 3, and 4a, shown in green).

Fig. 3. AIS's surface drainage basins. Each color represents a different drainage basin, computed using the 40 m resampled REMA (Howat and others, 2019) and drainage basin delineation from TopoToolbox (Schwanghart and Scherler, 2014). The drainage basin in which CTL forms is shown in green. The background image is a Landsat 8 true-color image from 17 January 2019. Grounding lines are shown in black (Depoorter and others, 2013).

Fig. 4. (a) Stream networks mapped from Landsat imagery over 1987/88–1991/92 (green) and 2012/13–2014/15 (pink). The drainage basin of the CTL is shown in green. Grounding lines are shown in gray (Depoorter and others, 2013). (b) Mapped maximum area CTL margins from summer images, vertically offset to show multi-year variation. The margins are digitized from Landsat 7 and 8 true-color imagery from (top) 11 March 2004, 27 January 2005 and 11 February 2006, and (bottom) 13 February 2013, 11 February 2014 and 29 January 2015. (c) The number of summertime images of maximum measured lake area that have surface water in each Landsat pixel, indicating the relative frequency of water coverage (see Table S1 for images).

3.2.2. Drainage system formation

We took two approaches to approximate the seasonal formation of the drainage system on AIS: mapping streams from optical Landsat imagery composites from two multi-year periods, and an automated lake detection method to map meltwater ponds.

1. (i) We selected Landsat images of consecutive years of summertime meltwater drainage during two multi-year periods (3–4 years) ~ 20 years apart. We take this multi-year approach because stream networks are seasonally transient phenomena with water flowing into and out of different sections of the network, each image captures only a part of the whole network. Spatially variable cloud cover further limits our ability to use a single image (or even multiple images from a single year) to fully characterize the distribution of the drainage network. We manually digitized the visible stream networks in multiple true-color images from consecutive summers to obtain the most complete map of channels from each multi-year-long time period. We used 4 years of imagery from the austral summers of 1987/88, 1988/89, 1989/1990 and 1991/92 (Landsat 4, Bands 1, 2 and 3) and 3 years from the austral summers of 2012/13, 2013/14 and 2014/15 (Landsat 8, Bands 2, 3 and 4). We combined the mapped drainage networks into two composite drainage maps corresponding to the two multi-year periods.

2. (ii) We used a simple automatic lake detection algorithm based on the automatic lake detection method of Moussavi and others (2020), supplemented by manual mapping to map the maximum lake area for the 24 years in which evidence for melt production is observed in optical imagery. We selected Landsat 4, 7 and 8 images with visible meltwater, but were limited to dates with cloud-free imagery across all four Landsat tiles that cover AIS. We processed each image to yield top-of-atmosphere reflectance (Pope and others, 2016; Moussavi and others, 2020). In each image, we computed the ratio, r _RB, of the blue band B (Landsat 4: 450–520 nm; Landsat 7: 441–514 nm; Landsat 8: 452–512 nm) to the red band R (Landsat 4: 630–690 nm; Landsat 7: 631–692 nm; Landsat 8: 636–673 nm). We used the threshold for Landsat 8 published in Moussavi and others (2020) and used their unmodified, published code to remove nunataks, clouds, satellite artifacts and lakes <5000 m² (five Landsat 4–8 multispectral pixels). The imagery from Landsat 4 (1988–92) and 7 (2000–13) could not be analyzed with Moussavi and others’ (2020) method directly because the wavelengths of the bands in the earlier imagery differ from Landsat 8's bands, and their method was tuned using Landsat 8 data. We chose to use a similar band ratio threshold method to automatically detect lakes, but this method was less reliable in removing clouds, shadows and satellite artifacts. The thresholds for Landsat 4 and 7 were determined by finding the value that yielded the closest match to lake perimeters mapped manually (by delineating the lake margin pixels in ArcMap) in red-green–blue imagery. We used r _RB = 1.02 for Landsat 7. For Landsat 4 images, we found the normalized difference of the blue and red bands ((R−B)/(R+B)) with a threshold of 0.017 more effective.

Following Moussavi and others (2020), to minimize erroneously selected pixels in Landsat 4 and 7 images, we removed lakes <5000 m² in area and any lakes that overlapped mapped nunataks. Additionally, in Landsat 4 and 7 images, we inspected the lake classification and manually removed pixels observed to be clouds or satellite artifacts. Lakes that were observed through clouds were mapped manually, as the simple algorithms failed in these cases. We computed the areas of each detected lake, except for lakes detected between 2003 and 2013, when Landsat 7's broken scan-line corrector (SLC) prevented this. Lastly, selecting from our processed images for each year one image with peak summertime lake area and minimal cloud cover, we created a cumulative frequency map by summing annual lake maps to show how often each depression fills with water (Fig. 4c).

3.2.3. CTL area and downstream extent

We expanded our temporal sampling by manually mapping the margins of the CTL (Figs 1 and 3) in Landsat, MODIS Terra/Aqua and WorldView-2 optical imagery for all austral summers when water was observed in this basin (1973/74, 1987/88, 1988/89, 1989/90, 1991/92 and 2003–19). We used Landsat panchromatic band images when available, and red band images otherwise. We used the daily MODIS imagery (when available post-1999) and Worldview-2 optical imagery to ensure that the lake margin we mapped in Landsat imagery represented the maximum lake extent of that austral summer. In one case (2008), when no Landsat imagery was available to map the CTL at its maximum size, we mapped the lake margin directly from MODIS imagery. This imagery allows us to capture the lake margins which we then used as an indicator of the basin's evolution (Fig. 4b) during two periods with consecutive years of imagery showing lake formation, 2003/04–2005/06 and 2012/13–2014/15. From the mapped lake margins, we measured area change and maximum downstream extent, which corresponds to the northernmost point of the lake, for each available image. We estimated the uncertainty in margin position in this method of manually mapping lake margins to be three Landsat 4–8 pixels (90 m).

3.3. Estimating meltwater depths and volumes

3.3.1. Lake depths

We used Moussavi and others’ (2020) light-attenuation method to estimate water depth using cloud-free Landsat 8 images from the 2012/13–2018/19 austral summers. As recommended by Moussavi and others (2020), we used the average of the depths obtained using the red band (which Pope and others (2016) estimate has an RMSE = 0.43 m) and the panchromatic band (RMSE = 0.45 m). We used Moussavi and others’ (2020) light attenuation coefficients and reflectance of infinitely deep water. To estimate total meltwater accumulation, we integrated the depths across the ice shelf. To estimate the volume uncertainty, we took a conservative approach and used the integrated product of all estimated water depth pixels’ area and their uncertainty value (which is the average of the RMSE values reported for each Landsat 8 band by Pope and others (2016), i.e. 0.44 m).

3.3.2. RACMO2.3p2 meltwater volume estimates

We used output from the RACMO 2.3p2 to provide an estimate of meltwater production on AIS from January 1979 to December 2018 (van Wessem and others, 2018). The gridded output has a spatial resolution of ~27 km, and we used daily outputs. Previous work comparing in situ surface mass balance observations with RACMO output showed a good agreement in Greenland (Noël and others, 2018) and Antarctica (Lenaerts and others, 2017; van Wessem and others, 2018; Zhang and others, 2018; Agosta and others, 2019).

To estimate seasonal total melt input to the CTL, we integrated the daily surface melt volumes at 40 m resolution over the REMA-derived catchment of the lake (Section 3.2.1, Fig. 3, catchment shown in green) using MATLAB's triangulation-based linear interpolation function, ‘griddata’. Note that because RACMO models water as refreezing instantaneously in our area of interest, which is unrealistic, particularly in bare ice areas such as AIS's ablation zone, we used the RACMO modeled surface melt volume without subtracting modeled refreezing. The fraction of melt that remains in snow and firn pore space and refreezes is neglected, the implications of which we discuss later. We then computed cumulative and total summer melt volumes by integrating them through each summer (December–February). While the most informative comparison is between melt volume and lake volume, we can only calculate lake volumes from Landsat products for which light attenuation coefficients have been empirically calculated and tested, i.e. Landsat 8. As volume-to-volume comparisons are limited to the Landsat 8 period, we also compared the CTL's downstream extent and lake area with modeled melt production to maximize the number of seasons included in the comparison. When correlating modeled melt volumes with Landsat-derived volumes and extent, we performed weighted linear regressions and reported the goodness-of-fit r ² and p-values.

3.4. Lake freezing

3.4.1. A simple model of lake freezing

We observed a spatio-temporal pattern of changes in SAR backscatter in 10 m resolution imagery from late summer/fall 2017 (9 February 2017 to 22 April 2017) over the CTL, and interpreted this pattern as an indication of an ice lid forming and thickening until the lake was frozen through completely (Sellmann and others, 1975; Jeffries and others, 1993; Miles and others, 2017; Engram and others, 2018). To test this interpretation, we developed a simple thermal model of ice lid growth, which balances the heat conducted through the lid with the latent heat released during freezing, to estimate the depth of water H frozen in time t _H. We neglected latitudinal gradients in surface temperature over the lake as output from RACMO shows no evidence of this. Equating steady heat loss through conduction to the atmosphere per unit area and the rate of energy released by the growth of the ice lid thickness yields

(1)$${ dh\over dt} = {k\over \rho_{\rm i} L} {T_{\rm w} - T_{\rm a}( t) \over h},\; $$

where t is the time, h is the ice lid thickness, ρ _i is the density of ice (918 kg m⁻³), L is the latent heat of fusion of ice (334,000 J kg⁻¹), k is the thermal conductivity of ice at − 15°C (2.34 Wm⁻¹ k⁻¹), T _w is the water temperature (assumed 0°C), and T _a(t) is the air temperature (Cuffey and Paterson, 2010). Integrating yields

(2)$$H^2 = {2k\over \rho_{\rm i} L}\int_0^{t_H} ( T_{\rm w} - T_{\rm a}( t) ) dt$$

We took two approaches to approximate T _a. First, we integrated Eqn (2) numerically using the daily-mean surface temperature simulated by RACMO in the grid cell with the most coverage of the terminal lake. We integrated between the onset of visible freezing (9 February 2017) and the last visible SAR backscatter change (22 April 2017). Second, we used the average value of this temperature time series, $\bar {T_{\rm a}}, -25.8 ^\circ$C, and computed H using Eqn (2) with t _H again defined by the SAR observations.

3.4.2. Observed lake freezing rates from SAR imagery

Based on our interpretation of the SAR backscatter changes over the CTL as evidence of complete freeze-through and the thermal model used to support this interpretation, we examined AIS-wide SAR data to determine freeze-through across the ice shelf. We analyzed SAR imagery over AIS between 6 February 2017 and 15 February 2017, between 6 May 2017 and 15 May 2017, and between 6 September 2017 and 15 September 2017. We averaged the images in each of these three periods to calculate the average backscatter intensity in each pixel and to produce a single image across different ground tracks during the austral summer, fall and spring. We generated these composite images in Google Earth Engine at 30 m resolution.

4. Results

4.1. Overview of the AIS supraglacial drainage system

Of the 45 summers between 1972 and 2018, 26 had at least one satellite image of the southern AIS that was sufficiently cloud-free to assess the state of the supraglacial drainage system. Of these, there were 24 summers with visible, surface meltwater ponding (Fig. 2). In general, visible meltwater lakes and streams started to appear between September and November in optical imagery. Water accumulated first in the linear depressions or troughs, in the blue ice zone near the southern grounding zone (Fig. 1). Next, lakes formed in troughs in areas downslope: on AIS's eastern margin, northeast of Clemence Massif, and northeast of Fisher Massif (Fig. 1). Lakes varied in shape from thin and elongate, typically in surface troughs, to oval and ring-shaped in flatter areas near the eastern and western grounding lines. While there was some temporal variation in the shapes of individual basins, surface depressions filled with meltwater to form similar lakes to previous summers (Fig. 4b) while advecting with ice flow.

Most lakes were connected via drainage channels, forming separated drainage networks, several of which terminated in large, elongated lakes. In most cases, in situ melt and water delivered by drainage channels filled these terminal lakes over several weeks. The largest of these terminal lakes (CTL) formed in a trough aligned with ice flow in the center of AIS (Kingslake and others, 2017; Stokes and others, 2019) and reached a median length of ~50 km between 1973 and 2019. Our maximum observed length of this lake was 137 km in February 1992, and SAR imagery from 15 August 1993 indicates this feature extended even further when it over-topped into the adjacent trough (Phillips, 1998, their Figure 3).

In at least two cases, we observed a transition in the appearance of the water in the CTL from dark blue to light blue/gray at the downstream end of the basin where it protruded into the northern snow-covered region. We interpreted this color change as a transition from surface flow to near-surface flow through permeable snow or firn and the downstream edge of the light blue/gray region as a wetting front. One example was captured in daily MODIS imagery (26 January 2006 to 15 January 2006; Supplementary Video 1).

4.2. Temporal variability in drainage networks

We found a broad agreement in the areas and the lengths of the mapped channel networks from two time periods (1987/88–91/92 and 2003/04–05/06; Section 3.2.2) and the REMA-derived catchments (Section 3.2.1). Figure 3 displays the AIS's catchment structure derived from REMA. The largest lakes and the streams connecting them were mostly aligned parallel to and rarely crossed the REMA-defined catchment divides. Highlighted in green in Fig. 3 and in Fig. 4a is a thin, central drainage catchment that contains the CTL. To the west and east are two other catchments, which in most years remain disconnected from the central catchment. Of the 1268 mapped stream segments (defined as portions of streams that connect two confluences or connect the upper and lower limits of a drainage system to each other or a confluence) only 44 cross the REMA-derived catchment divides (Figs 3 and 4a). Although channels rarely cross catchment boundaries, we observed high interannual variability in the downstream maximum extent of the drainage networks.

We characterized melt seasons into four categories based on the areal coverage of the lakes on AIS, calculated using the automated lake detection method: low (areal coverage of $< 0.01\percnt$), moderate (0.01–0.1%), high (0.1–1.1%) and very high ($\gt1.1\percnt$) melt. Across these categories, we documented variability in the state of the drainage system in terms of the extent to which the lakes are connected by streams and whether the CTL forms:

Low melt (1972/73 and 1999/2000): We observed no lakes in these 2 years. Instead, we observed darkening in surface channels relative to the white snow and ice, but no indication of surface flow through channels or large-scale ponding.

Moderate melt (1988/89, 2000/01, 2001/02, 2002/03, 2006/07, 2008/09, 2009/10, 2010/11, 2011/12 and 2015/16): A small number of lakes formed in linear troughs near the southern grounding line, around Clemence Massif, near Fisher Massif, and in isolated basins along the eastern grounding line (Fig. 1). Few melt streams connected the lakes during these years, and the CTL did not form. The average areal coverage of meltwater ponds on AIS in the austral summers of 2010/11 and 2015/16 (moderate melt years with uninterrupted coverage of AIS, without SLC-error artifacts) was $0.07\pm 0.05\percnt$ of the total area of AIS.

High melt (1989/90, 2003/04, 2004/05, 2007/08, 2012/13, 2013/14, 2014/15, 2016/17, 2017/18 and 2018/19): In these years, melt lakes formed in the same locations as in moderate melt years. In addition, lakes formed in depressions along the center line of AIS, west of Clemence Massif. We observed a continuous network in the central catchment, and in most of these years, the CTL forms. The areal coverage of meltwater ponds in 1989/90, 2013/14, 2014/15 and 2017/18 (years with imagery without SLC gaps) was $0.62\pm 0.36\percnt$. In 2016/17 and 2018/19, areal coverage was higher than other high melt years, $1.04\pm 0.44\percnt$.

Very high melt (1973/74, 1987/88, 1991/92 and 2005/06): In these years, lakes filled depressions in locations not inundated even in high melt years. The terminal lakes were significantly extended, particularly in the central catchment, where the visible drainage network grew to ~ 140 km in length, in 2005/06, and possibly ~160 km in 1991/92 (though heavy cloud cover makes this uncertain). Additionally, the central drainage network expanded into regions of unsaturated firn, and a lake near the eastern grounding line appeared to drain rapidly into the ice shelf (see Section 4.4). In all very high melt years, quantification of areal coverage was prevented by a lack of cloud-free coverage across AIS (1973/74, 1987/88 and 1991/92) or by the SLC-error (2005/06). The available 2005/06 imagery indicated areal coverage by lakes $> 1.18\percnt$ and, although clouds prevented quantification of coverage in the other very high melt years, inspection indicated that they had a similar coverage.

4.3. Lake depths and volume

Focusing on summer 2018/19, a high melt year during Landsat 8 coverage, we estimated lake depth based on Landsat 8 imagery from 24 January 2019 (Fig. 5). This yielded a mean water depth of 0.97 m, with a SD of 0.56 m (Fig. 5, inset). The 99.95 percentile of all AIS water depths in 2018/19 was 3.84 m (Fig. 5). These average water depths were similar to water depths on other ice shelves, for example Banwell and others (2014) estimated a mean water depth of 0.82 m on Larsen B Ice Shelf from Landsat 7 imagery from a pre-collapse date, 21 February 2000. The CTL was relatively shallow with a 99.95th percentile depth of 1.42 m in 2016/1, and 2.33 m in 2018/19. This lake was shallower at the sides of the trough than in the center and shallowed slightly to the north (downstream) end (Fig. 6).

Fig. 5. Meltwater depths from Landsat 8 (using method of Moussavi and others (2020)) from 24 January 2019. Grounding lines are plotted in black (Depoorter and others, 2013). (Inset) Histogram of water depths.

Fig. 6. (a) The observed time to freeze-through (days) based on the timing of the transition from the Sentinel-1 SAR backscatter associated with initial freezing (Fig. 7a) to that of freezing through (Fig. 7c) of AIS's CTL. (b) Landsat 8 light attenuation water depth estimates (Section 3.3.1) from 27 January 2017. (c) REMA DEM showing the elevation above sea level of the large central trench, using a hillshade effect for better visualization. The CTL margin from the Landsat 8 27 January 2017 image is plotted in black in each panel.

Integrating the depths estimated from each 30 × 30 m pixel in the Landsat 8 imagery from the end of the melt season for each year between 2013 and 2018 provided estimates of maximum meltwater volume over AIS for each year (Table 1). We used this method to quantify variability in surface meltwater volumes in the CTL (see Section 4.5) and across AIS over the relatively short period covered by Landsat 8. In the table caption, we indicate in which scenes we expect a significant loss of meltwater volume to freeze-over. We computed (i) the average water depth across AIS (by taking the mean of all depth estimates), (ii) the average maximum depth (by computing the maximum depth of each lake on AIS and taking the mean of these maxima) and (iii) and the total meltwater volume in AIS (by summing the depth estimates and multiplying by the pixel area, 900 m²). Note that estimated depths on AIS were small compared to the estimated uncertainty from Pope and others (2016). Nonetheless, we assumed that the Landsat 8 light-attenuation method provided a useful measure of relative volumes, even though absolute uncertainty was high.

Table 1. Peak summer meltwater depths and volumes calculated using the Landsat 8 light-attenuation method across AIS (Moussavi and others, 2020)

*For 2015/16 we replaced an absent Landsat 8 image tile from 1 February 2016 (20160201128112) with an image with the same spatial coverage from 26 February 2016 (20160226127112), which may underestimate water volume in lakes that were frozen over in the later scene.

The highest volume estimate was from summer 2016/17 (0.83 ± 0.40 × 10⁹ m³), which was consistent with this year being a high melt year. The lowest volume estimate was the summer of 2015/16 (0.12 ± 0.08 × 10⁹ m³), though the missing Landsat column 112 tile on this date was substituted with a melt estimate from the next available cloud-free image with the same spatial coverage, on 26 February 2016. This substitution may have led to an underestimation due to the loss of visible, surface water pixels as lakes froze over. There was less interannual variability in average depths than in total meltwater volume: the highest average depth was 0.93 m (2016/17) and the lowest average depth was 0.61 m (2015/16).

4.4. End of melt season lake drainage and freezing

We interpreted Landsat imagery of cracks appearing on the surface of lakes from the end of the summer (mid-February to mid-March) as floating ice lid formation. We also saw evidence of drainage into and possible through the ice shelf near the grounding line in daily MODIS imagery from 10 January 2006 to 25 January 2006 (2006 drainage event location shown in Fig. 1). A large lake (~12 km) or several connected lakes shrank by $93\pm 5\percnt$ over a period of 7 days until a water body only 2.4 ± 0.1 km wide was visible. This drainage event was also observed by Ice, Cloud and Land Elevation Satellite laser altimetry (Fricker and others, 2009).

We interpreted a spatio-temporal pattern of change in SAR backscatter from the CTL (Fig. 7, Supplementary Video 2) as evidence of the lake freezing through to its bed. As SAR may penetrate several meters into dry snow and ice surfaces (Page and Ramseier, 1975), radar backscatter changes due to subsurface effects have been commonly observed. Previous work on Canadian and Alaskan terrestrial lakes has described the distinct appearance in radar imagery of lakes completely frozen through to their beds and lakes not frozen through to their beds (Leconte and Klassen, 1991). The cause of this difference has been confirmed experimentally to be the presence or absence of a radar reflective ice–water interface at the bottom of an ice layer, which returns a high backscatter due to the sharp dielectric contrast between freshwater ice and liquid water (Atwood and others, 2015). Thus, observed transitions from high to low backscatter in supraglacial and terrestrial lakes have previously been interpreted as evidence of lakes freezing through to their beds (Sellmann and others, 1975; Jeffries and others, 1993, 1994; Miles and others, 2017; Engram and others, 2018; Dunmire and others, 2020). SAR backscatter over the CTL was relatively high (0.5–2.0 dB) immediately after ice-lid formation (Fig. 7a; the timing of which is based on interpretation of Landsat 8 and MODIS imagery; 2 February 2017 to 9 February 2017). In each SAR image, obtained 6 days apart, the size of the high backscatter area decreased as the surrounding area of low backscatter migrated inwards (Fig. 7b; Supplementary Video 2). We observed this pattern of inward migration simultaneously across the surface of the lake. By midwinter, the lake had nearly uniform, relatively low backscatter values of − 6 to − 2 dB (Figs 7c and d). If the base of the ice lid was approximately flat, the inward migration of the high-to-low transition was consistent with the bathymetry of the lake (Figs 5 and 6), which is deepest in the center. We therefore interpreted this inward pattern in SAR as the CTL freezing through at the sides first, and as progressively deeper areas froze through completely, this high-to-low backscatter transition migrated inward. Engram and others (2018) estimated the water depth below which this high-to-low transition becomes indiscernible at 4 m depth, much deeper than the ~2.33 m maximum depth in 2018/19 estimated above (see Section 4.3).

Fig. 7. (a–c) Sentinel-1 SAR backscatter images from (a) 9 February 2017, (b) 29 March 2017 and (c) 28 April 2017 showing the inward pattern of migration of the high-to-low backscatter transition. (d) Backscatter histograms for the CTL, 9 February 2017 (blue) and 28 April 2017 (red).

Based on this interpretation, in 2016/17, 66±6 days elapsed between the initial freezing of the lake's surface and the lake freezing through completely (Fig. 6). Using our simple thermal model (Eqn (2)) with a constant temperature equal to the average RACMO surface temperature from that period, $\bar {T_{\rm a}} = -25.8^\circ {\rm C}$, a freeze-through time of t _H = 66 ±6  days corresponds to a maximum depth, H = 1.48 ± 0.09 m. Alternatively, integrating Eqn (2) using the daily surface temperatures yields H = 1.49 ± 0.09 m (Fig. 8). Note that the depths estimated using a temperature-dependent thermal conductivity in the freeze-through model deviate from those presented in Fig. 8 by < 0.05 m. These depth estimates agree within uncertainty with the 99.95th percentile of the 27 January 2017 Landsat 8-derived depths (1.42±0.44 m), the maximum measured CTL water depth of 2016/17 (Fig. 8).

Fig. 8. Histogram: we show the histogram of depths estimated from light attenuation from a 27 January 2017 Landsat 8 image, and the corresponding 99.95th percentile depth, 1.42 m (blue dashed line). We plot the range of depths calculated from a simple freezing model (Eqn (2)) using RACMO-derived skin temperatures, T _air(t) and t _H = 66 ± 6 days from the Sentinel-1 SAR image backscatter time series (9 February 2017 to 22 April 2017). This method yields a depth estimate of 1.39–1.57 m when using an average of T _a over the period of freezing (black) and 1.40–1.57 when using daily T _a(t) (red).

In order to investigate freeze-through across AIS, we inspected SAR imagery from February, May and September 2017 (Fig. 9, Section 3.4.2). In many locations across AIS, we saw evidence of the same progressive freeze-through described above for CTL. In some cases, this continued until the lake was fully frozen, as it did at CTL (Fig. 9, yellow and blue insets), but in others, central areas of high backscatter remained in the May and September imagery, which we interpret as these lakes retaining liquid water beneath their ice lids throughout the winter (Fig. 9, cyan and red insets). Tracking these lakes in SAR imagery from September 2017 to September 2019, we found that they appear completely frozen (i.e. their backscatter values are relatively low and uniform) after any melt season in which liquid water does not refill a particular basin.

Fig. 9. Ice-shelf wide composite SAR backscatter intensity images from summer (February), fall (May), and spring (September) are created by averaging the backscatter values of Sentinel-1 SAR 30 m resolution images between 6–15 February 2017, 6–15 May 2017, and 6–15 September 2017. Here, we show the February composite image with insets displaying the backscatter intensity of a few lakes across the three seasonal images, with some of the lake margins mapped from the 3 February 2017 Landsat 8 image. All SAR images are shown with the same gray scale. We interpret lakes that appear as solid, relatively lower (darker) backscatter as frozen through. The upper lake in the cyan inset retains some exposed liquid water in the February image, visible as a very dark backscatter in the upper lake. We interpret lakes with an area of high (bright) backscatter within their margins as having a floating lid of ice above liquid water.

4.5. Lake area, maximum downstream extent and modeled meltwater production

The maximum downstream extents of the main drainage networks varied interannually, primarily due to changes in the northward growth of the large terminal lakes. These lakes formed in troughs aligned with ice flow, which have smaller along-shelf slopes than across-shelf slopes. Focusing on the CTL's catchment (Fig. 3), we compared seasonal melt production (December–February) in its REMA-derived catchment to its maximum mapped area, its maximum downstream extent (i.e. northernmost point), and, when Landsat 8 imagery was available in 2013–19, its volume. We found a weak positive correlation between modeled seasonal melt in the catchment and maximum mapped lake area (Fig. 10a). Including years when no lakes were observed in the regression yields r ² = 0.43, $p-value= 9 \times 10^{-4}$ (not shown). Excluding these years yields r ² = 0.35, $p-value = 0.04$ (Fig. 10a). In several cases, we observed very large lakes when RACMO predicted moderate (e.g. 1991/92) or low melt production (e.g. 1987/88). Conversely, in many years with modeled melt production higher than in 1987/88, the CTL was either absent or small. Similarly, we found a weak positive correlation between modeled melt production and the latitude of the central lake's maximum downstream extent (r ² = 0.30, p − value = 0.007 with nonlake years included and r ² = 0.24, p − value = 0.1 when they are not).

Fig. 10. (a) Comparison of summertime surface melt volume spatially integrated over the REMA-derived drainage basin of the CTL and the mapped maximal area of the CTL (r ² = 0.35, $p-value=0.04$). The slope of the weighted best-fit line (shown in red) is 0.287. (b) Comparison between RACMO 2.3p2 surface melt volume spatially integrated over the central basin and water volume estimates from light-attenuation in Landsat 8 imagery. Only water depths collected before freeze-over (inferred from optical imagery, after 31 January) are included. The slope of the weighted best-fit line (shown in red) is 0.47, r ² = 0.58, $p-value= 9.8 \times 10^{-5}$.

We compared melt volumes estimated from summertime Landsat 8 images to the accumulation of modeled melt over the period between the start of the melt season, taken as 2 September, and the date of each image (Fig. 10b). We found a weak positive correlation when using all available Landsat-derived volumes (not shown; r ² = 0.41, p − value = 1.5 × 10⁻⁵). When we excluded volumes derived using imagery from after the onset of ice-lid formation (approximated to be 31 January in each year), we saw a stronger positive correlation between Landsat-derived and RACMO-calculated melt volumes (Fig. 10b; r ² = 0.58, p-value = 9.8 × 10^-5). The slope of the linear fit was 0.47, implying that the Landsat 8-derived volumes were on average less than half of RACMO-derived volumes.

5. Discussion

Streams and lakes have formed on AIS for decades (Phillips, 1998; Kingslake and others, 2017; Stokes and others, 2019). Streams and elongated lakes form along flow stripes and surface troughs associated with the convergence of upstream flow units (e.g. Mellor and McKinnon, 1960; Swithinbank, 1988; Phillips, 1998; Fricker and others, 2002, 2009). We have extended this previous work and added an analysis of how this water drainage system varies interannually. We mapped the stream networks and analyzed how hydrologic routing and drainage catchments conform to topography as described by the REMA DEM. Agreement between the spatial distribution of REMA-derived basins and mapped channel networks indicates that ice-shelf surface topography exerts a first-order control on ice-shelf surface drainage (Figs 3 and 4a). Moreover, between 1972 and 2019, we observe meltwater-filled depressions repeatedly forming in similar locations on the ice shelf. The locations of these depressions are likely to be controlled on the large-scale by ice flow (e.g. MacAyeal and Sergienko, 2013). Three main drainage catchments co-exist on AIS and do not appear to exchange significant meltwater.

We saw evidence of drainage through lake-bottom fractures, but this appears to be rare (see Section 4.4), so we assumed insignificant volume is lost to englacial drainage into the underlying firn or into the ocean. Likewise, we selected end-of-summer imagery when the lake area reached its peak and before significant freeze-over, but we expect a fraction of meltwater was still neglected from the estimate due to freeze-over. Water volumes from Landsat 8 light attenuation estimated over the whole AIS from 2012/13 to 2018/19 are between 0.11 and $0.83\times10^9\, {\rm m}^3$ of water, the same order of magnitude as water volumes calculated by Moussavi and others (2020). The ranking from smallest-to-largest annual peak melt volumes is the same between our estimations and those of Moussavi and others (2020), except for 2016/17, which we estimate to be the largest melt volume and which Moussavi and others (2020) estimate to be the second largest after 2018/19 (see Moussavi and others (2020), their Fig. 12). This may be due to the difference in volumes estimated from a subset of AIS, such as the one Landsat tile of AIS used in Moussavi and others (2020) versus the entire AIS area used here.

At the end of each melt season, surface meltwater on AIS freezes over where it ponds. From inspection of fall and spring SAR imagery, we hypothesize that lake freeze-through is widespread, as evidenced by spatio-temporal variation in SAR backscatter (Figs 7 and 9). Several pieces of evidence support this hypothesis: the timing of these changes at the end of the melt season, the agreement between estimates of the maximum depth of the CTL from Landsat 8 and a simple thermal model of ice lid freezing, and field observations of a layer of bubble-free ice under ~ 60 cm of snow on AIS (Phillips, 1998), which indicates in situ refreezing of meltwater in this area. This interpretation is also consistent with the established SAR backscatter proxy for lake freeze-through described in Alaskan lakes (e.g. Engram and others, 2018)) and in Greenlandic supraglacial lakes by Miles and others (2017), and is supported by thermal modeling of Antarctic supraglacial lakes by Dunmire and others (2020). However, several surface properties, such as surface roughness, grain size, and internal structures lead to changes in the SAR backscatter strength (Fahnestock and others, 1993). Therefore, while we favor our interpretation, it is not definitive. To test the hypothesis will require radiative transfer modeling and field observations.

In addition to water freezing in place, we observed evidence of rapid drainage into the ice. Fricker and others (2009) also saw evidence of multi-day drainage of meltwater into and possibly through the ice shelf, near the eastern grounding line (Fig. 1, location of 2006 event). Airborne and satellite observations of rapid lake disappearance and doline formation in 1989 (not shown; Fricker and others, 2002) and 2007 (Fig. 1, location of 2007 event) suggest this occurs in other regions of AIS (Dan Zwartz, personal communication, 3 October 2017).

The drainage networks on AIS typically terminate in large lakes that extend into the snow-covered regions downstream of the region of meltwater production during years with sufficient melt production. In contrast to the long-term, steady structure of drainage networks, the downstream growth of these terminal lakes exhibits large interannual variability in maximum downstream, northward extent. One potential driver of this variability is the amount of meltwater available each year. However, we found only a weak, but significant (p-value=  0.04), positive correlation between the maximum area of the CTL and modeled melt production in the lake's catchment (Fig. 10a). Possible explanations for this weak correlation include:

1. (1) Meltwater production modeled by RACMO may not accurately represent melt in the lake's catchment. While daily mean 4 m air temperatures recorded at a weather station (Fig. 1) correlate well (r ² = 0.70) with RACMO daily mean 2 m temperatures simulated at the same location, we lack direct observations of melt in the region needed to test this possibility.

2. (2) Even if RACMO accurately simulates melt production in the catchment, input to the CTL may not be accurately approximated by simulated melt in the lake's catchment because a significant proportion of the water may be stored in porous snow or firn. RACMO simulates zero runoff in this region, instead all meltwater refreezes locally in the model. The assumption of no runoff is unrealistic; thus, we chose to compare lake volumes to RACMO-modeled melt, which implicitly assumes no refreezing. This simplification neglects water refrozen in firn and snow pore space and neglects interannual variability in refreezing. Indeed, despite previous studies (Trusel and others, 2015; Lenaerts and others, 2017) showing that RACMO often under-predicts meltwater production in blue-ice areas, such as our study area on AIS, we find that RACMO-modeled melt volumes are systematically larger than Landsat-derived volumes in the central basin (Fig. 10b), which is consistent with significant refreezing. If the unaccounted for refreezing is a large proportion of the total meltwater produced and if it varies significantly interannually, this could partly account for the low correlation between lake area and modeled melt volumes.

3. (3) The shape and size of the drainage system feeding the lake (its catchment area) may change interannually. Although we integrate modeled melt production over a static REMA-derived catchment, the catchment could evolve in time. For example, if the flow of the ice shelf changed, adjusting surface topography, or if years with particularly high melt production led to over-topping of and incision through drainage divides, melt could be transported into the lake from outside the modeled catchment, which would not be included in the volume calculation. However, a lack of ice-dynamic change (King and others, 2009) and the broad agreement between channel networks mapped 25 years apart (Fig. 4a) indicate that the structure of the drainage system has remained approximately unchanged. MODIS imagery (not shown) revealed one example of catchment rearrangement in 2005/06 (a year when the lake grew to its second-largest area) when more westerly lakes drained into the central catchment than in previous years. Unfortunately, no sufficiently cloud-free imagery is available from other years with particularly large terminal lakes (e.g. 1987/88 and 1991/92) to map the channel networks that formed in these years. So, while we found no evidence of dynamic catchment change, we cannot rule this out. This is a potentially important mechanism because if catchment areas grow as melt production increases, this could facilitate a complex relationship between meltwater production and drainage network growth.

4. (4) Refreezing in the lake basin could cause it to evolve interannually and influence the growth of the next year's drainage network. Terminal lakes typically freeze through completely during the fall. This mechanism modifies the topography of the basin, filling low points. In some instances, refrozen ice masses are visible in later years as obstacles that subsequent summers’ meltwater must flow around, extending the drainage system farther downstream along the trough than in the previous year, even if the total melt volume is similar. Furthermore, where melt flows into snow-covered areas, it will fill pore space and refreeze. When meltwater flows into these areas the following year, the water flows over the low permeability, refrozen meltwater. As the lakes are laterally confined in elongated troughs, this would allow meltwater drainage to extend further downstream by reducing the loss of meltwater to the underlying substrate.

These latter two mechanisms involving the evolution of basin permeability and topography potentially constitute a nonlinear link between meltwater production and maximum downstream drainage extent (see also Stokes and others, 2019). If moderate to very high melt is sustained for consecutive years, these mechanisms could lead to year-on-year growth of the drainage systems. Supporting this idea, we observe in the Landsat record that years with the largest drainage networks and largest lakes are often the final year in a series of consecutive years of high melt (1991/92, 2005/06, and 2014/15; e.g. Fig. 4b). These very long lakes extend into unsaturated snow to the north of where melt extends during low or moderate melt years. Moreover, such a link between firn permeability and drainage system growth has been observed in the percolation zone of the Greenland ice sheet. Machguth and others (2016) and MacFerrin and others (2019) observed large near-surface ‘ice slabs’, leading them to conclude that refreezing of meltwater in firn pore space limits downwards percolation by decreasing firn permeability, forcing subsequent meltwater to drain along the surface. MacFerrin and others (2019) showed that these low permeability zones have expanded since 2001 in Greenland extending the area dominated by surface runoff by up to $26\pm3\percnt$, and we posit a similar effect could occur on Antarctic ice shelves.

Our ability to test if processes associated with catchment and basin changes are operating in Antarctica is limited by cloud-free satellite data availability. While the variety and volume of available satellite observations are growing rapidly, the data available to test hypotheses regarding multi-year or decadal variation in drainage networks is restricted to moderate-resolution Landsat and MODIS imagery. While these long-duration datasets are valuable resources, we lack usable clear-sky imagery coverage of the entire AIS from 1974 to 1986 and from 1992 to 1999 (Fig. 2). Additionally, the broken SLC obscures up to 50% of lake pixels in Landsat 7 imagery images, making computation of total lake areas difficult during this satellite's operation period. Given these limitations, testing the hypothesis that drainage systems respond nonlinearly to increases in melt through the processes proposed above will require extending this work to other ice shelves and potentially using more sophisticated automatic lake mapping techniques (Dirscherl and others, 2020; Moussavi and others, 2020; Dell and others, 2020).

6. Summary

We have examined the spatial and temporal patterns of surface meltwater drainage on AIS, East Antarctica from 1973 to 2019 using Landsat (1 and 4–8), MODIS and WorldView optical imagery and Sentinel-1 SAR imagery. We have shown that AIS's surface drainage networks are topographically controlled and exhibit significant interannual variability in maximum downstream drainage extent, lake size, areal coverage and meltwater volume. We focused on a large lake at the terminus of AIS's central catchment (CTL), which exhibits wintertime freeze-through and evidence of refreezing in firn pore space, processes we link to interannual variability. In the highest melt years, this lake transports water further than 100 km and extends into firn-covered areas downstream of the main drainage network. We estimated water volumes held in the AIS surface drainage network based on depths from Landsat 8; our estimate for the summer of 2016/17, a ‘high melt’ year and the highest volume observed, was 0.83 ± 0.40 × 10⁹ m³. Based on fall- and winter-time SAR imagery, we hypothesize that lake freeze-through is widespread across AIS, creating many impermeable ice masses that remain in the next melt season. We also hypothesize that a combination of melt volume and near-surface refreezing controls the maximum downstream extent of the drainage networks on the ice shelf, which is primarily expressed by the size of large terminal lakes that form in ice-flow-aligned troughs.

We examined the relationship between the area of the CTL and surface melt volume modeled by RACMO 2.3p2 over the lake's catchment and found only a weak positive correlation. This weak correlation could be explained by inaccuracy in the modeled melt input to the lake (e.g. related to ignoring refreezing in the lake's catchment), or by refreezing in the lake's basin changing its topography and permeability. This latter mechanism has the potential to cause drainage networks that terminate in firn-covered areas to expand nonlinearly in response to increases in meltwater input. As this firn-terminating drainage configuration is common around Antarctica, it will be important to understand these processes and determine how they will contribute to accumulating water in ice-shelf areas vulnerable to hydrofracture (Lai and others, 2020). This knowledge of ice shelf vulnerability will improve our predictions of future ice mass loss.