Filling and drainage of a subglacial lake beneath the Flade Isblink ice cap, northeast Greenland

. The generation, transport, storage and drainage of meltwater play important roles in the Greenland Ice Sheet (GrIS) subglacial system. Active subglacial lakes, common features in Antarctica, have recently been detected beneath the GrIS and may impact ice sheet hydrology. Despite their potential importance, few repeat subglacial lake ﬁlling and drainage events have been identiﬁed in Greenland. Here we examine the surface elevation change of a collapse basin at the Flade Isblink ice cap, northeast Greenland, which formed due to sudden subglacial lake drainage in 2011. We estimate the subglacial lake volume evolution using multi-temporal ArcticDEM data and ICESat-2 altimetry data acquired between 2012 and 2021. Our long-term observations show that the subglacial lake was continuously ﬁlled by surface meltwater, with the basin surface rising by up to 55 m during 2012–2021, and we estimate 138.2 × 10 6 m 3 of meltwater was transported into the subglacial lake between 2012 and 2017. A second rapid drainage event occurred in late August


Introduction
At total of 64 subglacial lakes have been identified in Greenland from airborne radio-echo sounding Livingstone et al., 2022). Most of them are stable, showing little or no evidence of volume change or input from the surface, and are located between the equilibrium line altitude (ELA) and the relatively flat, frozen-bed ice sheet interior. Only a few hydrologically active lakes that are recharged by surface meltwater have been identified from ice surface elevation change measurements Howat et al., 2015;Livingstone et al., 2019Livingstone et al., , 2022Palmer et al., 2015;Willis et al., 2015). Compared to the widely distributed stable subglacial lakes, the active subglacial lakes are affected more directly by surface meltwater, and their drainage may significantly influence glacier flow dynamics (Davison et al., 2020;Livingstone et al., 2019). Despite this importance, our understanding of Greenland's subglacial lakes has been primarily developed from theoretical studies or inferences from geophysical exploration due to sparsity of direct observations (Davison et al., 2019). The presence and movement of meltwater at the ice-bed interface are considered to significantly affect ice dynamics (Meierbachtol et al., 2013). Given the expected increases in surface meltwater production in a warming climate (Mottram et al., 2017;Sellevold and Vizcaino, 2021), it is of critical importance to understand Greenland Ice Sheet (GrIS) hydrology, especially the routing, storage, drainage and recharging of subglacial water.
Satellite remote sensing techniques have been used to monitor the subglacial lakes and detect their activities. As an indirect observation of subglacial lake activity, long-term ice surface elevation changes are usually derived from satellite altimetry (e.g., Fricker et al., 2007;Fricker, 2018, 2021). More recently, time-stamped digital elevation models (DEMs) have been utilized to reveal the detailed patterns of surface deformation (e.g., Livingstone et al., 2019;Willis et al., 2015). A few studies also use synthetic aperture radar (SAR) speckle tracking (Joughin et al., 2016;Hoffman et al., 2020) and interferometry SAR (InSAR) (Gray et al., 2005;Neckel et al., 2021) to detect ice surface displacements. However, few studies have investigated the long-term filling and drainage of subglacial lakes in Greenland. In particular, the subglacial lake volume change, water residence times and drainage are still poorly understood.
At the Flade Isblink ice cap (81.3 • N, 15.0 • W) in northeast Greenland (Fig. 1), a collapse basin in the ice cap surface about 70 m deep, created by sudden subglacial lake drainage between 16 August and 6 September in 2011, was first revealed by Willis et al. (2015). Basin surface elevation estimates with DEMs created from stereoscopic satellite imagery suggest that rapid surface uplift occurred over the 2 years following the collapse as supraglacial meltwater was transported to the ice base, refilling the subglacial lake. Although the Flade Isblink ice cap is not directly connected to the wider GrIS, its glacial setting is similar to that of the northern GrIS. It is important to investigate its behavior and impact on ice dynamics, which may lead to improvements in our understanding of subglacial lakes beneath the GrIS. In order to better understand the repeat subglacial lake filling and drainage, here we extended the surface elevation time series records to early 2021 with ArcticDEM repeat surface models and ICESat-2 altimetry data. We describe the longterm subglacial lake behavior, analyze its volume change and compare it with the surface runoff supply. We also identify a second drainage event in 2019 and explore the impact of drainage on glacier dynamics.

Surface elevation and basin volume change calculation
Surface elevations from 2012 to 2017 were first acquired from multi-temporal ArcticDEM strip data (Porter et al., 2018). The initial absolute accuracy of ArcticDEM strip data is less than 4 m in horizontal and vertical planes. Therefore, the DEM strips should be vertically co-registered before calculating elevation changes. Only a few DEM strips extend over bedrock or have ICESat footprints as ground control points in our study area, so we cannot directly coregister each of them. Instead, we first co-registered a DEM acquired on 20 April 2015 using the 3-dimensional offset val-ues provided by the metadata text file as a reference. A square window centered over the collapse basin with sides equal to twice the length and width of the basin (∼ 7.6 km) was defined. Another 1500 m buffer was set outward along the boundary of the collapse basin (Fig. 2i). Then, all the other DEMs were vertically co-registered to the reference DEM by calculating the mean elevation differences using the pixels within this window but outside the 1500 m buffer. We applied an iterative, 3-standard-deviation filter to remove outliers when estimating the elevation differences . The DEM precision was estimated from the standard deviation of the elevation differences that remained after the iterative filter. In this way, the influence of both the systematic vertical offsets and snow accumulation or melting was removed.
Besides ArcticDEM data, Advanced Land Observing Satellite (ALOS) Global Digital Surface Model "ALOS World 3D" (AW3D30) Takaku et al., 2014Takaku et al., , 2020 was also used to analyze the elevation change. The AW3D30 DEM in our study area is derived from data spanning the period 2006-2010, just before the late summer of 2011 when the deep basin formed. As above, the AW3D30 DEM was vertically co-registered to the reference DEM. Note that, the co-registered DEMs only represent "relative" ice surface heights that have eliminated systematic changes over the larger ice cap due to surface accumulation or melting and other processes, rather than the true elevation.
The surface elevation measurements from the Advanced Topographic Laser Altimeter System (ATLAS) onboard ICESat-2 were also used to extend the time series to early 2021. As a successor to the ICESat satellite mission, ICESat-2, a polar-orbiting satellite with a 91 d repeat cycle and 92 • orbit inclination, was launched in September 2018 (Markus et al., 2017). ATLAS generates six green (532 nm) laser beams in three pairs along one reference ground track, and each pair contains one weak and one strong beam. In the across-track direction, the spacing between each beam pair is ∼ 3.3 km, and each pair of strong and weak beams is separated by ∼ 90 m. There are eight tracks (four pairs) that pass through the collapse basin, with two pairs (Track 0126 pair3 and Track 0321 pair2) passing over the main basin and another two pairs (Track 1266 pair3 and Track 1107 pair2) passing over the area between the main basin and thumb basin (Fig. 1b). We only used repeat cycles 3-9 for our study because the first two are not repeats due to pointing control issues.
We used the ICESat-2 level-3a Land Ice Height (ATL06) data product, removing poor-quality elevation measurements caused by clouds or random clustering of background photons based on the ATL06 quality summary flag . Further, we checked for height consistency by calculating adjacent elevations using the along-track slope parameter and comparing the estimated to the measured elevations. Only the data where the difference between original elevations and the estimated elevations is less than 2 m were used  (Li et al., 2020). In order to reduce errors introduced by large across-track slopes, we merged the two single-beam track data for the left beam and right beam into one beam pair. A reference track was first calculated by averaging all of the single-beam tracks from both left and right ground tracks. The elevation of the reference track for each cycle was then estimated from the left and right single-beam track measurements and the across-track slope parameter (Li et al., 2020). This procedure provides four repeat-track observations for elevation change analysis.
After all of the ICESat-2 data were co-registered to the reference DEM using the method described above, the time series of elevation change over the collapse basin were estimated along the four reference tracks using both the registered ArcticDEM and ICESat-2 data. Additionally, average ice surface elevation changes were also estimated at three reference track crossovers (Fig. 2i).
Changes in subglacial water volumes were estimated by integrating elevation change over the basin area. Previous studies show that a reduction in the depth of the depression would result from the inflow of the ice around the basin (Aðalgeirsdóttir et al., 2000;Willis et al., 2015). Therefore, we expect that the basin volume change here was mainly caused by ice inflow and subglacial lake filling. Assuming the subsidence that occurs around the basin outline in a 1500 m buffer region corresponds to ice flowing into the basin, we calculated the inflow volume by integrating the surface elevation changes over the buffer area . The volume change of the subglacial lake was then estimated by differencing the basin volume change and ice inflow volume.

Catchment delineation and surface melting analysis
The catchment boundary was extracted using ArcticDEM surface elevation as follows Yang et al., 2019). First, we filled the ArcticDEM surface to create a sink-free DEM raster. Then we identified the flow directions from the slope direction on the partially filled DEM. Finally, the Basin function in ArcGIS software was used to delineate the catchment boundary.
To assess the surface meltwater dynamics, we used estimates of meltwater runoff from the high-resolution Regional Atmospheric Climate Model (RACMO2.3p2) (Noël et al., 2018). Daily runoff produced in the catchment was generated from RACMO2.3p2, which is statistically downscaled to a 1 km horizontal resolution (Noël et al., 2019). The total runoff within the catchment was calculated by summing the 1 km grid cells within the catchment boundary. Furthermore, a series of Landsat 8 Operational Land Imager (OLI) and Sentinel-2 MultiSpectral Instrument (MSI) images acquired during the 2014-2020 melt season were used to better illustrate the supraglacial lakes and streams.

Ice velocity estimate
We obtained estimates of the ice surface velocity from the MEaSUREs Greenland Monthly Ice Sheet Velocity Mosaics from SAR and Landsat dataset, Version 3 (Joughin et al., 2018). These include monthly surface velocity estimates for the Greenland Ice Sheet and periphery and are posted at a 200 m grid resolution. Due to the limited coverage of the ice velocity product in the summer, an 800 m by 800 m region located downstream of the collapse basin was chosen to evaluate changes in ice surface velocity (Fig. 6c). We calculated the mean velocity within this region to estimate the velocity time series from 2018 to 2020.

Collapse basin surface elevation change
After the basin surface rose by up to 38 m during 2012-2014 , the elevation of the entire basin continued to increase during the ArcticDEM period (2012-2017) (Fig. 2i). The surface of the main and thumb basins uplifted by up to 65 and 50 m, respectively, while the southern part of the collapse basin only had a maximum uplift of ∼ 10 m.  Fig. 2a-h). The ice surface elevation then showed a sudden decrease in 2019, followed by a gradual increase from January 2020. Profiles CC and DD show that the elevation changed gradually while the surface maintained approximately the same shape.
Combining ArcticDEM and ICESat-2 data, we estimate changes in surface elevation at three crossovers (Fig. 3). Elevation at the south edge of the collapse basin (crossover G) continuously increased by ∼ 10 m from 2012 to 2021. At the shallow saddle between the main basin and the thumb basin (crossover F), the surface rose at a faster rate of ∼ 5 m yr −1 during 2012-2021, with a sudden subsidence of ∼ 2 m between 20 June and 18 September in 2019. The main basin (crossover E) had the most rapid surface uplift of ∼ 9 m yr −1 from May 2012 to April 2019. After continuously increasing for the 8 years after the basin first collapsed in 2011, the surface of the main basin subsided by more than 10 m between 19 April 2019 and 16 January 2020. Afterward, the elevation increased again at a rate of ∼ 5 m yr −1 . The elevation increased dramatically in the melt season during 2014-2016 (Fig. 3 inset). During the melt season in 2014 and 2015, the surface of the main basin rose ∼ 3 m at a rate of ∼ 33 and ∼ 28 m yr −1 , respectively. In 2016, the elevation gained ∼ 7 m between 8 July and 4 September. This rate of elevation increase of ∼ 49 m yr −1 is about half of the observed rapid surface uplift during the 2-week period in 2012 .

Subglacial lake volume change and surface meltwater runoff
We define the volume of the collapse basin to be the volume between the pre-collapse ice surface and the post-collapse ice surface. Time series of volume change of the collapse basin and subglacial lake during the period of 2012-2017 are shown in Fig. 4. Between 3 May 2012 and 5 May 2013, the volume of the collapse basin decreased by 47.5 × 10 6 m 3 , with ∼ 55 % (26.3 × 10 6 m 3 ) of the changes occurring as a result of surface uplift caused by increasing subglacial lake volume and the remainder due to rapid infilling by ice flow. After 2013, however, the rate of ice inflow slowed, accounting for a small portion of the basin volume loss. Basin volume showed notable changes corresponding with rapid surface uplift in the 2014-2016 melt seasons. In the 2014 melt season, the basin lost a total volume of 4.2 × 10 6 m 3 between 7 July and 7 August, with the majority of the loss (3.0 × 10 6 m 3 ) due to influx of surface meltwater to the subglacial lake. During the 2016 melt season, the volume of the surface basin decreased by 17.0 × 10 6 m 3 between 26 July and 4 September, and ∼ 97 % (16.6 × 10 6 m 3 ) of the volume change was due to subglacial lake refilling. Over the entire 5-year period, the collapse basin lost 176.0 × 10 6 m 3 of volume. About ∼ 21 % (37.8 × 10 6 m 3 ) of the loss was due to ice inflow, and the remaining 138.2 × 10 6 m 3 was the result of subglacial lake refilling by surface meltwater.

Discussion
Few active subglacial lakes have been observed in Greenland Howat et al., 2015;Livingstone et al., 2019;Palmer et al., 2015;Willis et al., 2015). This may be partly because subglacial lakes under the GrIS are nearly 8 times smaller than in Antarctica  and, therefore, may not be resolved by altimetry observations due to sparse track density. Alternatively, the surface of the GrIS margin is typically steeper than in Antarctica, making the depressions in hydraulic potential required for lake formation less likely to occur . Finally, efficient subglacial drainage systems formed in the melt season may release the stored water, preventing subglacial lake formation. Here we investigate an active subglacial lake located under the Flade Isblink ice cap, which is on the periphery of, as well as separated from, the northern GrIS. The lake is similar in size to the active lakes found beneath the ablation zone of the GrIS. Willis et al. (2015) first discovered the sudden subglacial lake drainage event under the Flade Isblink ice cap during the autumn of 2011. A collapse basin was formed on the surface of the ice cap, and the surface rose over the next 2 years due to recharging of the subglacial lake. Our estimates of the collapse basin and subglacial lake volume change between    , who reported a similar amount of volume change of 46.5 × 10 6 and 29.6 × 10 6 m 3 , respectively. Additionally, we also concur that volume change caused by ice inflow accounted for a large portion of basin volume loss over the first 2 years (2012-2014) of our investigation period. The rate of influx declined as the depression became shallower in the following years, decreasing its contribution to basin volume change.
Surface meltwater may drain into crevasses or moulins every melt season and lead to a rapid elevation increase in a short period. Contrasting with the north-flowing meltwater that mainly drained into crevasses on the southern margin of the collapse basin in 2012 (the polygon in Fig. 5a), much of the meltwater accumulated locally in a supraglacial lake at the southern part of the basin during the 2014-2016 melt season (Fig. 5c). Changes in supraglacial hydrology may have been due to the burial of the crevasses and the significant remaining surface relief ( Fig. 5a and b). Following the switch in drainage location from the basin-edge crevasses in 2012 to moulins within the basin during 2014-2016, the rate of surface meltwater drainage decreased. This is confirmed by the decreasing rate of basin surface elevation uplift during the melt season. From the time of surface meltwater draining into moulins and the observed rapid uplift of the main basin during 2014-2016, we conclude that surface meltwater recharged the subglacial lake every melt season. Moreover, the larger amount of meltwater observed in 2016 corresponded to larger elevation gains. All of these processes indicate that the subglacial lake volume is primarily controlled by supraglacial meltwater filling.
Between 19 April 2019 and 16 January 2020, the surface of the main basin lowered by more than 10 m (Fig. 3). We conclude this surface lowering is most likely due to drainage of the subglacial lake, which is further confirmed by Sentinel-2 images acquired at the end of August 2019 (Fig. 6). Between 24 and 26 August 2019, obvious surface lowering is observed over the main basin and a distinct depression formed at the thumb basin area (Fig. 6a-b), indicating a rapid subglacial lake drainage event occurred during this time. A lack of elevation measurements at the main basin in 2019 prevents us from estimating the exact duration of drainage events. According to the elevation variation at the shallow saddle between the main basin and the thumb basin (crossover F), we speculate the drainage may end in September. The time and duration of this drainage event are consistent with previous large subglacial lake drainage events identified in Greenland Palmer et al., 2015;Willis et al., 2015), which usually starts at a time when the subglacial drainage system becomes efficient and meltwater drains through connected channels . However, the volume of water drained in the 2019 event was likely much less than in 2011, indicating that a large amount of meltwater remained in the subglacial lake. The ubiquity of partial subglacial lake drainage is unknown in Greenland, but similar processes have been observed beneath ice caps in Iceland where the subglacial lakes may become sealed before draining all the water (Björnsson, 2003).
Variations in ice flow speed are consistent with water pressure variations expected during subglacial lake drainage. In August 2019, ice flow immediately downglacier of the basin increased by a factor of 3 over the pre-subsidence values (Fig. 6d) before decreasing back to average values in the following month. We conclude that these abrupt changes resulted from the drainage event as meltwater released from the subglacial lake initially overwhelmed the drainage system, resulting in a larger increase in water pressure and sliding speed. As the subglacial drainage system increased in efficiency and/or the discharge of water decreased as the meltwater was drained, water pressures and sliding speeds declined.
Continued basin surface uplift from 2011 to 2019 suggests that the subglacial lake was not filled by supraglacial meltwater in a single melt season and that water storage persisted after the initial lake collapse basin initially formed. We speculate that the subglacial lake may be located upstream of a topographic ridge that would form a depression in the hydropotential field and, therefore, could store meltwater draining from the surface Palmer et al., 2015). As meltwater is stored, the piezometric head within the lake increases until it exceeds the hydropotential gradient holding it in place, causing discharge. Once discharge begins, melting of channel walls at high water pressures would cause rapid expansion of the drainage system, increasing efficiency of drainage until the piezometric head in the lake lowers and discharge decreases and then ceases. Accumulation of meltwater draining from the surface then begins and continues until another subglacial drainage event occurs.
The elevation profiles through the collapse basin (Fig. 2) indicate that the subglacial lake may have not been fully filled when the drainage event occurred in 2019. This drainage is not associated with high-surface-melt years, and the duration of the drainage event is less than 1 month. We speculate that this subglacial lake exhibits a pattern of slow filling and rapid drainage, similarly to all active lakes beneath the Icelandic ice caps (Livingstone et al., 2022). In contrast, three active lakes beneath the GrIS are characterized by long periods of quiescence . This implies that the timings and behaviors of repeat filling and drainage of Greenland's subglacial lakes are determined not only by water storage volume but also by meltwater input variability (Schoof, 2010) and bedrock relief .
Subglacial lake water in Greenland is sourced from either geothermal and frictional melting or surface meltwater input . The temperature at the bed of the Flade Isblink ice cap is far below the pressure melting temperature, and the ice moves relatively slowly, ruling out the local production of basal meltwater . Moreover, the Flade Isblink ice cap is isolated from the GrIS; hence the subglacial lake is not connected to the subglacial hydrology network beneath the GrIS. Therefore, surface meltwater is likely the only supply for this subglacial lake. Supraglacial meltwater would be routed to the bed through crevasses and moulins and flow toward the ice margin, inducing ice flow variations. A modeling study has estimated that, during an average melt season, about 39 % and 47 % of the surface runoff is drained through crevasses and moulins, respectively, in west Greenland (Koziol et al., 2017). However, only a portion of this surface meltwater would access the ice-bed interface (Nienow et al., 2017). Our results show that 3.0 × 10 6 m 3 of supraglacial water reached the subglacial lake over a 1-month period (7 July to 7 August) during the 2014 melt season. At the same time, total surface runoff produced within the catchment is estimated to be 4.7 × 10 6 m 3 . Thus, only ∼ 64 % of the surface meltwater successfully descended to the bed. The remainder may be locally refrozen in the underlying snowpack (Harper et al., 2012) or firn aquifers that have been detected around the collapse basin area Kuipers Munneke et al., 2014;Miller et al., 2022). Ice slabs are also likely to exist locally (MacFerrin et al., 2019), so meltwater may be restricted or travel via other drainage paths that our study is unable to detect. However, we also cannot rule out the possibility of other drainage paths, subglacial or supraglacial, that we have not resolved.

Conclusion
In the autumn of 2011, a collapse basin about 70 m deep formed in the surface of the Flade Isblink ice cap in northern Greenland due to sudden subglacial lake drainage. Using multi-temporal ArcticDEM and ICESat-2 altimetry data, we document changes in surface elevation of the lake basin and estimate the subglacial lake volume change from 2012 to 2021. The long-term measurements imply that the subglacial lake was most likely recharged by seasonal influx of surface meltwater. The surface of the collapse basin rose by up to 55 m over the 9 years, with 138.2 × 10 6 m 3 of meltwater transported to the subglacial lake during 2012-2017. During our investigation period, a second rapid drainage event occurred in late August 2019, resulting in an abrupt ice velocity change. Compared to the 2011 drainage event, the amount of water drained in 2019 was much smaller and was likely only a portion of the stored water, suggesting partial drainage. In addition, the 2019 drainage was not associated with highsurface-melt years. These findings suggest that the triggering of subglacial lake drainage and subsequent evolution may be controlled by multiple factors and need further investigation. A model of surface melt over the catchment estimates that only ∼ 64 % of the surface meltwater successfully descended to the bed, implying the importance of quantifying the routing of surface meltwater inputs to the ice-bed interface. We have also shown that the new ICESat-2 data have great potential in detecting and monitoring active subglacial lakes beneath the GrIS. Our findings on the Flade Isblink ice cap that the subglacial lake can store meltwater over multiple years and decrease runoff to the ice margin are helpful for better understanding the hydrological processes on the GrIS.  ESA, 2021). The AW3D30 DEM is freely available from the Japan Aerospace Exploration Agency (JAXA) (https: //www.eorc.jaxa.jp/ALOS/en/dataset/aw3d30/aw3d30_e.htm, Tadono et al., 2014;Takaku et al., 2014Takaku et al., , 2020. RACMO2.3p2 Greenland daily runoff data were kindly provided by Brice Noël and can be obtained upon request.
Author contributions. QL and LZ conceived the study. WX and QL processed the data. QL, XC, FH and ZC analyzed the data. IH and MJ provided guidance on the data analysis and result interpretation. QL led the paper writing. QL and IH revised the manuscript.
Competing interests. The contact author has declared that neither they nor their co-authors have any competing interests.
Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.