Evolution of Supraglacial Lakes on the Larsen B Ice Shelf in the Decades Before it Collapsed

The Larsen B ice shelf collapsed in 2002 losing an area twice the size of Greater London to the sea (3,000 km2), in an event associated with widespread supraglacial lake drainage. Here we use optical and radar satellite imagery to investigate the evolution of the ice shelf's lakes in the decades preceding collapse. We find (1) that lakes spread southward in the preceding decades at a rate commensurate with meltwater saturation of the shelf surface; (2) no trend in lake size, suggesting an active supraglacial drainage network which evacuated excess water off the shelf; and (3) lakes mostly refreeze in winter but the few lakes that do drain are associated with ice breakup 2–4 years later. Given the relative scale of lake drainage and shelf breakup, however, it is not clear from our data whether lake drainage is more likely a cause, or an effect, of ice shelf collapse.


Introduction
The Antarctic Peninsula (AP) has experienced extreme warming in the middle to late 20th Century, with air temperatures increasing by almost 2.5°C (e.g., Vaughan & Doake, 1996;, up tõ 1990 (Turner et al., 2016). This has been linked to thinning and loss of the AP's ice shelves (e.g., Morris & Vaughan, 2003;Shepherd et al., 2003). During the 1990s, the Larsen B ice shelf (LBIS) began to shrink, culminating in its eventual collapse in 2002 (e.g., Glasser & Scambos, 2008;Scambos et al., 2003). Since then, the glaciers which formerly fed into LBIS have accelerated between twofold and eightfold, following the removal of buttressing forces formerly provided by the shelf (e.g., Rignot et al., 2004;Scambos et al., 2004). Through this mechanism,~9 Gt year −1 of grounded ice has been lost to the sea, accounting for one third of all ice loss observed from the AP since 2002 (Berthier et al., 2012).
Supraglacial lakes (SGLs) form in surface depressions from ponded melt water. They are a component of the ice surface hydrological network which also includes streams and rivers (Glasser & Scambos, 2008). This network is important because it can route surface meltwater into lakes (Stokes et al., 2019), moulins (Langley et al., 2016), or off of the ice into the ocean . By exporting meltwater off of the ice, supraglacial streams can strengthen ice shelves . SGLs, however, have been implicated in ice shelf breakup (e.g., Banwell & MacAyeal, 2015). Widespread SGL coverage was observed on the Prince Gustav and Larsen A ice shelves before their failure in the mid-1990s (e.g., Cooper, 2009), and lakes on the LBIS were observed to drain just before its collapse . SGLs can drain laterally through supraglacial streams (Kingslake et al., 2015) or vertically via hydrofracture; when water-filled crevasses propagate through the full ice thickness (Krawczynski et al., 2009;MacAyeal et al., 2003). Through repeated filling and draining, SGLs induce localized flexure of the ice shelf. It is thought that this can produce fractures that can cause neighboring lakes to drain in a chain reaction, leading ultimately to ice shelf disintegration (Banwell & Macayeal, 2015).
Since much of Antarctica is fringed by floating ice shelves that provide important buttressing to grounded ice flow (Fürst et al., 2016;Goldberg et al., 2019), it is important to understand the contribution of supraglacial hydrology to ice shelf stability. Here we investigate the potential role of SGLs in the collapse of LBIS (Figure 1) by analyzing their evolution from 1979 to 2002 using a combination of optical and synthetic aperture radar (SAR) satellite imagery (supporting information Table S1).

Satellite Data
Here we use all available, high quality satellite images acquired between November and March during 1979 and 2002 over the LBIS. These data are eight optical satellite images acquired for seven separate years by NASA/USGS's Landsat series and 13 SAR images acquired for six separate years by ESA's European Remote Sensing (ERS) platforms (Table S1). The SAR signature of lakes changes annually from dark lakes against a light background in winter to bright lakes against a dark background in summer (Text S2.2). Here we choose the summer SAR image in which lakes were clearest for each of the years we have data for quantitative analysis; these were mostly acquired between mid-January and mid-February. In order to compare lake characteristics between years, ERS data are radiometrically calibrated and georeferenced (Nagler et al., 2016). SGLs are then manually delineated in ERS and automatically and manually delineated in Landsat (Texts S2.1 and S2.2). We quantify uncertainty associated with using two different sensors by digitizing lakes in a 52-km 2 sample area in Landsat and ERS images acquired 2 days apart. We find that 7% more lakes are identified in ERS images, lakes are 16% larger on average in Landsat and cover an 8% greater area (Text S2.3).
To account for the fact that lake covered area changed during our study period, and to provide a fair comparison between Landsat and ERS, we subsampled our study area to that covered by lakes in 1988, omitting regions where lakes were unclear in ERS due to speckle noise (granular interference). We use a partially cloudy Landsat image from the 9th of February 1990 in section 3.2. We account for the cloud cover by calculating the change in lake characteristics in the cloud free area with respect to the Landsat image captured on the 19th of January 1988, then applying the relative change to the cloud covered portion. Finally, we use a radiative transfer model (Text S2.3) to compare lake depth and volume between two cloud-free Landsat images acquired on the 19th of January 1988 and 21st of February 2000.

Meteorological Data
We contextualize our findings with simulated meteorology and simulated firn conditions from the Regional Atmospheric Climate Model (RACMO2). Details of the model are found in van Wessem et al. (2016) and summarized here. RACMO2 combines the High Resolution Limited Area Model with the European Centre for Medium-range Weather Forecasts Integrated Forecast System and has been adapted for use over the large ice sheets of Greenland and Antarctica. ERA-Interim reanalysis data are used for forcing at lateral boundaries and surface topography is based on a combination of the 100 m and 1 km digital elevation models from Cook et al. (2012) and Bamber et al. (2009), respectively. In previous work we found that the model systematically underestimates air temperature by 1.82°C over LBIS (Leeson et al., 2017), and correct for this here by adding 1.82°C to all temperature estimates.

Southward Spreading of Lakes
Using images from Landsat and ERS (Table S1) acquired in nine of the years between 1979 and 2002, we quantify changes in lake covered area by manually delineating the southernmost boundary of lake covered area in ArcGIS ( Figure 2). In 1979, lakes are confined to the Hektoria/Green/Evans glaciers domain north of Foyn Point (locations identified in Figure 1). By the austral summer of 1987/1988 lakes had spread to the north of Exasperation Inlet (noted in  and in 1993 had reached its southern boundary. From~1996, lakes begin to encroach upon the shelf from the south and by 2002, covered almost the entire area. This corresponds to a southward spreading rate of 2.59 km (0.02°) per year between 1979 and 2002. This is split into two different periods with a threefold speed up in spreading rate after 1993 (0.01°per year 1979-1990 and 0.03°per year 1993-2002). According to simulations performed by RACMO2 over this region, 1993 was a very high melt year (more than twice as much melting as the 1990-1999 average, Leeson et al., 2017) and a year with particularly low minimum firn air content on the portion of the shelf lost in 2002 (Figure 2l). This suggests that high melting in this year saturated the ice shelf's firn pack, thereby preconditioning it for SGL formation in subsequent years (Kuipers Munneke et al., 2014). We also note that the retreat of the ice shelf (e.g., Kulessa et al., 2014) is clear in the ERS images after 1993, with the shelf progressively shrinking down to its precollapse extent by 2000.

Changes in Lake Characteristics and Relationship With Climate
We assess the variability in lake number, size, and total area in our study area between 1988 and 2002 and investigate potential relationships to climate. The presence of clouds in some Landsat images precluded their use in this analysis, reducing our sample to seven individual years across the 23-year study period. Using these data, we found no noticeable trend in lake characteristics ( Figure 3); lake number, size, and total area within our sample did not increase. Variability was reasonably high in each case with 1993 and 2000 acting as end members for each parameter with fewer, larger lakes in 1993 (335 lakes, mean area 0.64 km 2 ) and an abundance of small lakes in 2000 (1,170 lakes, mean area 0.07 km 2 ). We attribute this to topography; small lakes coalesce in high melt years to form fewer, larger lakes. Using regression analysis we found that summer (DJF) melting simulated by RACMO2 best explains interannual variability in mean and total lake area (Text S5). Mean and total lake area are strongly correlated with DJF melting (r = 0.95, p < 0.01 and r = 0.77, p = 0.01, respectively), in good agreement with Langley et al. (2016). We note that these patterns hold despite images being acquired on different days into the melt season. This suggests that lakes reach their maximum size by mid-January and do not grow or shrink significantly between then and mid-February. Lake number is weakly anticorrelated with DJF melting (r = −0.58, p = 0.18), suggesting that other factors, e.g., the time of year, act as a more important control (see Text S5).
We compare lake depth and volume between two cloud-free Landsat images acquired toward the beginning (the 19th of January 1988) and end (the 21st of February 2000) of our study period (Figures 3d and 3e). There are notable differences in both depth and volume between these images, despite similar predicted melt amounts in each year. Lakes hold more water in 1988 with a mean volume of 1.7 × 10 5 m 3 compared to 0.6 × 10 5 m 3 in 2000. This is not surprising since mean lake area is twice as large in 1988. This is likely a result of lakes refreezing; the 2000 image was captured a full month later than the 1988 image. More interesting, however, is that in 2000, lakes are more than 50% deeper than in 1988 with average maximum depths of 1.28 and 0.8 m, respectively. This is interesting because, for a fixed surface topography, one might expect the larger lakes to be deeper. This suggests that the surface topography of the ice shelf is not fixed, and that the depressions in which lakes form deepened between 1988 and 2000. This is despite the likely onset of refreezing prior to the date of acquisition of the 2000 image and is spatially independent ( Figure S7).

Evidence for Lake Drainage
In our combined ERS-Landsat record we find evidence for lakes refreezing, vertical drainage and lateral drainage.
Prior to 2002, lakes tend to persist through the winter without draining; e.g., ERS images acquired in March 1996 and 1997 and November 1996 show lakes on the LBIS after/before the DJF melt season. A cloud-free Landsat image acquired on the 1st of March 1986 shows no lakes, rather it shows abundant snow which suggests that the surface of the lakes has frozen over, and they have subsequently been buried. Buried lakes are visible in ERS because of the capability of radar to penetrate snow (Johansson & Brown, 2013). In ERS imagery from 1996,~100 lakes located in the northernmost part of the region are visible on the 29th of February but not visible on the 24th of March (Figure 4). Since RACMO2 shows that the air temperature was consistently below zero between these dates it is reasonable to assume that these lakes fully refroze. Similar numbers of lakes in the 21st of February 2000 Landsat image appear to be refreezing/refrozen, based on their similarity to refreezing lakes in previous work (e.g., Langley et al., 2016). This is supported by the temperature data; in the 7 days prior to image acquisition, mean temperatures across most of the ice shelf were below zero. In both cases, refrozen lakes seem concentrated toward the ice margin, likely because temperatures are colder there beyond the influence of warm föhn wind (e.g., Cape et al., 2015;Leeson et al., 2017).
We see evidence for lake drainage in years before 2000, mainly on portions of the ice shelf which broke off before the main collapse. In 1996, a cluster of~10 lakes disappear mid-melt season between ERS images captured on the 25th of January and the 29th of February ( Figure 5). Similarly, a cluster of approximately six lakes disappear between the 2nd of February 1997 (ERS) and the 4th of February 1997 (Landsat). In both cases the disappeared lakes are surrounded by visible lakes, suggesting that they did not refreeze. In the case of the 1996 images, the drained lakes appear slightly brighter than the surrounding surface (signature also noted in Miles et al., 2017). We see more evidence (~20 lakes) of this signature in an ERS image captured on the 24th of March 1996 (Figure 4r), mainly on the portion of LBIS lost before 2000. The drained lakes are not particularly noteworthy in that their area and shape appear consistent with surviving lakes. This suggests that their drainage was glaciologically determined, i.e., initiated as a result of local perturbations in ice flow. We note that this period in time has been associated with LBIS speed-up in other studies (e.g., Kulessa et al., 2014) and that the drained lakes are oriented perpendicular to the direction of flow, both of which support this inference. Since these lakes are grouped, it is possible that a single glaciologically driven event may have triggered a small scale chain reaction draining neighboring lakes (e.g., Banwell et al., 2013); however, higher resolution observations in time are needed to confirm this.
In the image captured on the 21st of February 2000 we see~19 shapes characterized by light and dark patterns consistent with steep-sided topographic depressions (e.g., Figure 4s). We interpret these to be lakes which have drained, supporting previous findings (e.g., Glasser & Scambos, 2008). Since we have a contemporaneous ERS image acquired on the 27th of February 2000 we examine the locations of these lakes in SAR and identify a possible signature (Section S9). We see evidence of this signature in ERS images acquired in 2002 around the time of shelf breakup; however, the images are too noisy and the signature insufficiently clear to draw any conclusions. In our record we also have two relatively cloud-free Landsat images from the early part of the 2001/2002 melt season captured on the 15th and 31st of December 2001. These images show lakes present on the 15th of December that are not present on the 31st of December (Figure 5c). Since this coincides with downstream lakes increasing in size between the two dates, and temperatures increasing rather than decreasing, we interpret this to be evidence of lateral lake drainage (e.g., Kingslake et al., 2015).

Conclusions
We use SAR and optical satellite imagery to investigate the evolution of SGLs on the LBIS prior to its collapse and find that lakes spread southward at a rate consistent with firn air depletion, covering almost the entire ice shelf by 2002. This is consistent with previous studies which suggest both a temperaturedependent latitudinal "limit of viability" for ice shelves (Morris & Vaughan, 2003) and a firn-density based metric of ice shelf vulnerability (Alley et al., 2018). We note in particular that 1993 and 1995 were the years of lowest firn air content in our record, directly preceding the beginning of the ice shelf breakup in 1995. Notably, RACMO2 does not simulate a temperature or surface melt trend during this period (Leeson et al., 2017). This suggests that ice shelf vulnerability is cumulative, and that accurate measurements and models of meltwater retention in firn are needed in order to assess the risk of collapse of other ice shelves.
We find that, according to our data, there was no trend in lake area over the 1979-2002 study period and that lake area is best described by seasonal melt volume. This is interesting because such behavior reflects an interconnected hydrological network where lake size is a function of throughput rate rather than total abundance of water. We see evidence for some linear meltwater features that terminate at crevasses or the shelf terminus, which could provide a mechanism for meltwater export off-shelf ( Figure S11). This supports the work of Kingslake et al. (2017) and Bell et al. (2017) who characterize Antarctic ice shelf hydrology as a dynamic system where excess meltwater is exported, e.g., off the ice shelf, as opposed to being stored locally. Our findings suggest that the LBIS exhibited a combination of local storage (in lakes) and meltwater export, which may explain why it was able to support an abundant population of lakes for several decades before it finally collapsed. We also find that lakes get deeper over time. Since this holds across a wide area, and we do not see evidence of repeated widespread draining, we attribute this to enhanced ablation at the lake bed (Buzzard et al., 2018) as opposed to a viscous response to successive fill/drain cycles (Banwell et al., 2013).
At the end of the melt season, we find that most lakes refreeze but a small proportion drain. Specifically we see evidence for~30 lakes draining prior to the 24th of March 1996. These potential drained lakes are mainly located on floating ice which broke off between 1996 and 2000. The first time we see evidence for multiple lakes draining on the portion of ice shelf lost in 2002 is in 2000, where we see evidence for~20 drained lake basins. Thus, we find that, in our data set, floating ice on which SGLs drain breaks up 2-4 years after drainage is first observed. We note, however, that the number of observed lake drainages is small in both cases, relative to the large size of the broken off area. As such, it is not clear from these data whether the lake drainages were a cause, or rather an effect, of the LBIS breakup.