Altimetry-based ice-marginal lake water level changes in Greenland: Unveiling annual variations in glacial lake outburst �oods linked to runoff

Greenland holds more than 3300 ice-marginal lakes, serving as natural reservoirs for out�ow of meltwater to the ocean. A sudden release of water can largely in�uence ecosystems, landscape morphology, ice dynamics and cause �ood hazards. While large-scale studies of glacial lake outburst �oods (GLOFs) have been conducted in many glaciated regions, Greenland remains understudied. Here we use altimetry data to provide the rst-ever Greenland-wide inventory of ice-marginal lake water level changes, studying over 1100 lakes from 2003–2023, revealing a diverse range of lake behaviors. Around 60% of the lakes exhibit minimal �uctuations, while 326 lakes are actively draining, collectively contributing to 541 observed GLOFs from 2008–2022. These GLOFs vary signi�cantly in magnitude and frequency, with the highest concentration observed in the North and North East regions. Our results show substantial annual differences in the number of GLOFs and the variations are driven by annual difference in meltwater runoff, except for the South West region. Our method detected a 1200% increase in the number of draining lakes compared to existing historical databases. This highlights a signi�cant underreporting of GLOF events and emphasizes the pressing need for a deeper understanding of the mechanisms behind and the consequences of these dramatic events.


Introduction
Globally, proglacial lakes, including ice-marginal lakes, hold approximately 0.43 mm of sea-level equivalent 1 .In Greenland, there are more than 3300 ice-marginal lakes 2 and over the past three decades they have increased number, while changes in their size have remained less clear [1][2][3] .These observed changes have been suggested to be associated with changes in ice sheet surface melt 3 and retreat of the ice margin 2,4 .Covering around 10% of the Greenland Ice Sheet (GrIS) ice margin and 5% of the Peripheral Glaciers and Ice Caps (PGICs) 4 , ice-marginal lakes exert an important in uence on ice dynamics.They have shown to accelerate glacier mass loss and terminus retreat through calving 5,6 and ice velocities at the GrIS margin are roughly 25% higher for glaciers terminating in lakes compared to on land 4 .
The water out ow from ice-marginal lakes can vary from a constant discharge to sudden outburst oods termed jökulhlaups or glacial lake outburst oods (GLOFs).These rapid drainage events have drastic consequences, including alterations in fjord circulations 7 and downstream geomorphology 8 , changes in local ice dynamics 9,10 and bedrock displacements 11 , as well as notable societal impacts 12 .The societal impacts of GLOFs in Greenland is lesser than other regions (such as the Himalayas, the Swiss Alps and Iceland) because of the sparse population density, with few settlements in close proximity to the Ice Sheet and/or surrounding ice caps.Moreover, small scale studies from Greenland have shown that the dynamics of GLOF events have evolved over time, resulting in changes in timing, frequency, release volume, and/or rerouting of the release water 8, 13,14 .Despite recent comprehensive studies of icemarginal lakes across Greenland utilizing optical and Synthetic Aperture Radar (SAR) satellite imagery, and Digital Elevation Models (DEMs) 2,3 , the recorded incidents of GLOFs remain low at 153, with only 25 GLOF locations documented 15 .This count is notably low when considering the substantial number of ice-marginal lakes in Greenland.Moreover, this greatly contrasts the number of reported GLOFs in other glaciated regions globally [15][16][17] in large-scale studies, likely re ecting the differing relevance of GLOFs to societal infrastructure and, in some cases, preservation of life.This suggests a scarcity of documentation rather than a scarcity of GLOF events in Greenland, highlighting the need for further investigation and emphasizing that the phenomenon is largely understudied in this region.Typically, GLOFs in Greenland have been monitored at individual level or focused on small regional areas, employing in situ observations 18 and remote sensing techniques 13,19 .
This study presents the rst-ever, comprehensive large-scale study of water level uctuations in icemarginal lakes across Greenland.Distinguishing itself from previous large-scale studies relying predominantly on optical and SAR satellite imagery 2,3,20 , our approach utilizes satellite and airborne altimetry data acquired from 2003 to present (Data and Methods).We utilize pre-de ned lake outlines 21 to extract altimetry measurement of all ice-marginal lakes larger than 0.2 km 2 (n = 1387).For each lake, we construct a reliable water level time series by applying a simple statistical outlier detection framework, and calculate the largest observed water level difference (dWL) during the observation period (Data and Methods).All lakes with a dWL exceeding 4 m are manually inspected and categorized into one of three general groups: i) lakes with GLOF behavior, ii) lakes without GLOF behavior but with an overall decrease in water level, and iii) lakes without GLOF behavior but with an overall increase in water level (Data and Methods).Finally, we investigate how annual changes in GLOFs relate to variations in runoff and ice dam thickness.

Ice-marginal lake changes across Greenland
Out of the 1387 ice-marginal lakes included in the study more than 80% (n = 1152) had a minimum of two observations after the statistical ltering (Fig. 1), with an observation de ned as the median water level of all altimetry measurements captured within one day.Our results show that 687 ice-marginal lakes have limited surface uctuation between 0 m and 4 m (Fig. 1 and Fig. 2), indicating that these lakes likely have a consistent and continuous water out ow.Some of the lakes may experience larger water level uctuations occurring in between observations; however, as ~ 82% have ve or more observations (Fig. 2), we expect this number to be limited.Ice-marginal lakes with limited water level uctuations are found all around the margin of the GrIS and the PGICs, with the highest relative concentration in the south west (SW) (~ 72%) and south east (SE) (~ 90%) sector.Conversely, the lowest concentration is observed in the central east (CE) sector (~ 47%), whereas the central west (CW), north west (NW), north east (NE) and north (NO) sectors have similar concentrations, ranging from 61 to 63% (Fig. 1).
Our ndings reveal that 465 ice-marginal lakes experience large uctuations exceeding 4m, with an average of 24 observations per lake.We detect 326 lakes exhibiting GLOF behavior (category i), corresponding to more than a quarter of all 1152 ice-marginal lakes with altimetry observations.We detect a total of 541 GLOF events, with 45% of the lakes draining more than once (Fig. 3).As our approach does not capture all occurring GLOFs, this represents a conservative minimum estimate (Data and Methods).Nonetheless, we still nd a notable increase of 301 lakes and 388 GLOF events when compared to existing databases documenting GLOFs in Greenland over the past century 15 .The highest number of GLOF events are detected in 2019 (n = 178), accounting for one-third of all events, including 75 events from lakes that exclusively drained in 2019 (Fig. 2).During the four years with complete ICESat-2 coverage (2019-2022), we identi ed a total of 510 events, among which 170 are one-off events from a single lake.Five lakes drained every year, while 29 lakes experienced three events.Additionally, 81 lakes drained twice, with 55% of them draining every second year.
Lakes demonstrating GLOF behavior are observed across all regions in Greenland, and in many areas we nd them located in close proximity to one another (Fig. 1 and Fig. 3).The NE, NO and south west (SW) sector have the highest number of ice-marginal lakes with GLOF behavior with 101, 78 and 58 lakes respectively.The highest relative concentration of 43% is observed in the CE sector.By calculating the difference between the pre-and post-GLOF water level, we get an estimate of the minimum drainage magnitude of each GLOF event (Data and Method).All sectors, except the SE, contain lakes with a drainage magnitude larger than 50 m, with the largest absolute and relative concentration in the SW with 11 lakes corresponding to 19% of all lakes in this sector.Additionally, the SW sector also have one of the highest concentration of lakes with 25-50 m drainage magnitude along with the CE, both with 21%.
Ice-marginal lakes without GLOF behavior, but with a general decrease in water level during the observational period (category ii), are mainly found in the northern sectors (NE, NW and NO), whereas those with a general increase (category iii) are located in the NE and partly the CE sector (Fig. 1).Additionally, we observe a large number of lakes with an overall water level increase located close to Bredebrae and Storstrømmen in the NE sector (Fig. 1 NE zoom-in).

Discussion
Water level time series, trends and simultaneous GLOFs Our analysis of water level time series for more than 1150 ice-marginal lakes reveals a diverse range of lake behaviors, along with large differences within lakes of the same category.This variability is in uenced by several factors, such as lake size, shape and location, catchment area, runoff, damming glacier size, and thickness, and, notably, the density and timing of observations.The latter has proved hugely important for properly gauging the water level variations for a large quantity of lakes on a Greenland-wide scale.
Within the category of lakes exhibiting GLOF behavior, we nd that those experiencing recurrent events tend to drain at decreasing water levels over time.While in some instances, this may be linked to the timing of the observations taken prior to the event, we observe a similar pattern for lakes with dense observation coverage, as exempli ed by Iluliallup Tasia and Lake Isvand (Fig. 4.).Furthermore, upon examining optical satellite images of the selected lakes, we also identify a reduction in the pre-GLOF lake area (Fig. 4.).Similar patterns have been observed in studies of individual lakes in Greenland and attributed to a thinning of the damming glacier 8,10,13 .Understanding whether this constitutes a general trend across Greenland requires a longer and more consistent dataset, which will become available as more data is continuously collected.However, large-scale studies of GLOFs from other glaciated regions detect only a moderate correlation between magnitude and glacier thinning 17 .The thinning of the damming glacier also in uence those lakes, which show no GLOFs and a general decrease in water level during the observational period (category II) (Fig S6).This can be induced by the lowering altitude of a potential spillway.However, in cases where a substantial decline in water level is observed from the early to the more recent observations, there is a potential risk of overlooking GLOFs that occur between observations.Lakes displaying a general increase in water level (category iii) are likely in uenced by the thickening or advance of the damming glacier (Fig. S7), potentially coupled with changes in runoff from the catchment area into the lake.Alternatively, the rise in water level could be due to GLOFs occurring before the acquisition of the rst altimetry observations, indicating that we are only obtaining measurements during the subsequent lling period.This was con rmed using optical satellite images for a small cluster of lakes dammed by Budol Isstrøm (Fig. 5).Two of the lakes drained simultaneously in late August 2017, while the third drained in the spring or summer of 2018, all prior to the rst ICESat-2 observation.
Given the observed close proximity of various draining lakes (Figs. 1 and 3), we hypothesize that simultaneous GLOFs are relatively common and may occur at multiple locations across Greenland.A similar case of simultaneous GLOFs is observed at A.B. Drachmann Glacier, the neighboring glacier of Budol Isstrøm, where two lakes drained simultaneously in 2022, following more than 3 years of water level increase (Fig. S2).The simultaneous occurrences of these GLOFs might be initiated by changes in the damming glacier, such as an acceleration in ice velocity leading to the opening of cavities 23 or a decrease in (sub)glacial meltwater 24 , possibly induced by a sudden drop in air temperature 8 .Alternatively, it may be a cascading effect, where the sudden release of meltwater from one lake triggers an immediate response in others, as it has been observed for supraglacial lakes 25 .However, the mechanism responsible for the simultaneous GLOFs is believed to be highly localized, as none of the other lakes dammed by Budol Isstrøm or A.B. Drachmann Glacier drained during the same period.

Annual GLOF response to changes in meltwater runoff
Throughout the observational period, we nd marked annual uctuations in the number of observed GLOFs (Fig. 6).From 2008 to 2018, the GLOF count is limited by the availability of observations.However, following the introduction of ICESat-2 in October 2018, there is a notable increase in both the annual number of observations and GLOFs (Fig. 6).Between 2019 and 2022, the number of observations remains almost constant, and thus cannot explain the substantial differences in the annual number of observed GLOFs during this period.Concurrently, there has been continuous thinning of the damming glaciers across all regions (Fig. 6 and S3).Ice dam thinning has been attributed to drainages 8, 10,23 , which would imply a consistent increase in the annual number of GLOFs from 2019 to 2022.
We nd that the annual number of GLOFs follows the annual variability in RACMO-derived runoff (Fig. 6 and Fig. 7), indicating that the number of GLOFs is higher in years with high runoff, and vice versa when runoff is low.The correlation is evident both at a Greenland-wide scale and at a regional scale.Previous work 3 has suggested increased runoff as a factor driving regional changes in total lake volume along the CW and SW margins from 1987-2010, but it has never before been linked to large-scale variations in GLOFs.Larger volumes of runoff accelerate the lling of ice-marginal lakes, causing them to reach their maximum threshold more rapidly and indicates that system responds almost instantaneously to these changes.However, regional differences in the covariance and the negative trend in the SW region imply that uctuations in runoff have varying in uence on the GLOF cycle in the different regions.In the SW region, the close clustering of three out of four years (2020, 2021, and 2022) in Fig. 7 indicates that the orientation of the trend is sensitive to the data from 2019.The SW has the highest runoff, receiving more than two times the amount of meltwater during an average year (2020) compared with some of the northern regions during extreme melt years (2019) (Fig S4).The extreme runoff in the SW in 2019 does not result in more GLOFs, pointing to potential alternative controls.The additional meltwater received during extreme years could possibly result in earlier ll-ups and, consequently, earlier drainage dates in all regions.Determining whether the lakes experience earlier outbursts would require longer and more temporally consistent time-series data.However, recent large-scale studies of GLOF changes during the past decade from other glaciated regions report earlier outbursts 17 .
Given the increase in runoff over the past 40 years (Fig. S5) and the anticipated acceleration in the coming decades in Greenland 26 , we hypothesize that the number of annual GLOFs has grown during the past decades, and will persist as an upward trend in the future.This may be further accelerated with persistent future thinning of damming ice.However, as runoff from the GrIS is currently 60% more variable form year-to-year compared to the recent three decades 27 , we would also expect large interannual variations.This contrasts with recent large-scale observations from Alaska, which show an unchanged frequency of GLOFs 16 , and points to potential differences in GLOF behavior between glaciated regions in the Arctic.

Lake outlines and identi cation of ice-marginal lakes
To accurately estimate changes in water levels of lakes using altimetry observations, precise and contemporary lake outlines are essential.We use lake outlines from the recently released mapping of Greenland conducted by The Danish Agency for Data Supply and Infrastructure (SDFI) 21 .The dataset has a high accuracy and internal consistency and is based on high-resolution satellite images (0.5 x 0.5 m) captured between 2018 to 2022, during a time when more than 96% of the altimetry data used in our analysis is captured.The dataset contains more than 150.000 lakes, each assigned a unique lake ID, with 10,073 identi ed as ice-marginal.This identi cation is based on the lakes proximity to the GrIS and peripheral glaciers, with the criterion of being situated within a 100-meter buffer from the corresponding SDFI ice mask.We remove all lakes with an area smaller than 0.2 km 2 , as smaller lakes have limited altimetry coverage and consequently contain a high fraction of potential inaccurate observations, leaving a total of 1387 ice-marginal lakes to be included the study.

Altimetry data
Lake water levels are determined by combining altimetry observations from four different datasets (Table 1).We used altimetry data products, such as ATL06 and GLAH06, as they are pre-processed and have been previously used for studying uctuations in lake water levels 10,13,28 .All ICESat-2 data was downloaded from the NSIDC homepage 29 between July 1st and July 10th, 2023, with the most recent observations captured on April 16, 2023 (Table 1).The ATL06 is chosen over the ATL13 dataset, which is speci c for inland surface water as it has extremely sparse coverage compared to the ATL06 dataset.
Contrary, the ATL03 photons dataset exhibits an exceptionally high density of observations, necessitating a substantial computational capacity, likely with minimal additional information.

Estimation of water level time series
All altimetry datasets are merged and spatially ltered using the lakes outlines, resulting in a comprehensive dataset of 15,745,605 altimetry measurements spanning all 1387 ice-marginal lakes.We apply a simple outlier detection framework to remove erroneous measurements located inside the lake area and calculate lake-speci c water level time series.The framework is based on the statistical variability of our observations and is designed to effectively handle lakes with varying numbers of observations: 1.For each lake, we calculate the median elevation of each unique day (Observation median) based on all measurements acquired on that day.Subsequently, we determine the median elevation across all unique days (lake median WL) (Fig. 8) to prevent bias due to variations in the number measurements across different days.2. Next, we calculate the absolute difference from all Observation medians to the lake median WL and calculate the median absolute deviation (MAD) (Fig. 8).We then de ne an upper and lower threshold using Eqs.( 1) and ( 2), which we use to identify potential outliers.
Upper threshold = lake median WL + C * MAD Eq. ( 1) Lower threshold = lake median WL -C * MAD Eq. ( 2) This detection process employs a conservative outlier detection approach with C set to 3. We opt for the MAD method as it offers greater robustness in handling large outliers observed in our dataset (Fig. 8).
3. If any Observation median falls outside of the upper or lower bounds, we mark them as potential outlier.However, we do not lter the observations immediately, but perform additional test to check if they are actual outliers: 4. If over 30% of measurements during a day are within the threshold, we remove measurements outside and recalculate a new Observation median (Fig. 8) 5.If over 70% fall outside, we test the validity of the measurements.We lter out all measurements if they fail to meet the following criteria: i) total measurements > 3, ii) Observation STD < 20m, and iii) difference between Observation median and lake median WL < 100m (Fig. 8). .Lastly, we perform a nal check on all Observation medians varying more than 50m from the lake median WL, to capture potentially undetected outliers.If the Observation median varies more than 10m from both the previous and following observation (date) within a 200-day window, the Observation median is removed (Fig. 8).This is based on the assumption that observations obtained within close proximity in time are likely to be more similar.
Following the outlier removal, we recalculate the Observation median water level of all remaining lake measurements and subsequently generate time-series illustrating the variations in water level for each individual lake.Next, we calculate the largest observed water level difference (dWL) (max.water levelmin.water level) of each lake to identify lakes with notable water level changes.

Determining lake characteristics
All lakes with a water level difference larger than 4 m were selected (n = 503) for subsequent analysis as we are interested in lakes with substantial uctuations.Within this subset, we examine the changes in water level and correct undetected outliers by manually removing them or by minimizing the lake area used for the spatial ltering using a negative buffer on the original lake area.In rare instances, we refer to optical imagery to con rm whether speci c altimetry observations are indeed outliers.Following this additional ltering, we are left with 465 lakes with a water level difference exceeding 4 m.
We classify these 465 lakes based on their water level changes and divide them into three general groups: i) lakes with GLOF behavior, ii) lakes without any apparent GLOF behavior but with an overall falling water level during the observational period, and iii) comparable to (ii) but with an overall rise in water level during the time series.There may be overlaps between categories as lakes experiencing GLOFs might also exhibit rising or falling water levels.However, we assign each lake to only one category due to the relatively short observational time span of the majority of our observations.For each of the lakes exhibiting GLOF behavior, we manually determine signi cant key parameters such as the number of drainages, drainage year, and the pre-and post-drainage water level.For some lakes, such as Lake Hullet (Fig. 4), limited observations make it challenging to discern seasonal variations and precisely estimate the magnitude and timing of GLOFs.Nevertheless, GLOFs occurring in 2020, 2021, and 2022 can be identi ed.Conversely, for lakes with more frequent observations, such as Iluliallup Tasia and Lake Isvand (Fig. 4), we can more accurately determine the seasonal evolution and GLOF characteristics.We de ne the magnitude of each GLOF event as the difference between the pre-and post-drainage water level.The drainage magnitude serves as a valuable indicator for assessing the scale of these events.However, as it depends on the timing of the observations it should be interpreted as a minimum magnitude.Additionally, it does not provide a direct measurement of the volume of released water, as the metric is dependent on the lake's bathymetry.Consequently, ice-marginal lakes with extensive surface areas and drainage magnitudes of less than 10 meters can still discharge signi cant volumes of water, while smaller lakes with larger drainage magnitudes may release less volume.

Accuracy assessment the of water level time series
To validate our ndings, we conduct a comparison between the water level uctuations derived from altimetry data and observed changes in lake area in optical satellite images for 110 lakes that drain.As some 110 lakes drained more than once during the observations period our validation include a total of 184 GLOF events.For each of these events, we manually examined optical satellite images to con rm whether a drainage in the altimetry-based water level coincided with alterations in the lake area.Out of these, we successfully con rm 95% of the GLOFs.For lakes with more detailed mapping of changes in the lake area (Figs. 4 and 5), we observed a high agreement with the observed uctuations in water level and the timing of the GLOFs.Nevertheless, there may be instances where we erroneously document a drainage event using our altimetry record.This could be attributed to various factors, such as erroneous altimetry observations, frontal uctuations of the damming glaciers, or the presence of large oating icebergs in the lakes.
When comparing the generated water-level time series with ndings from comprehensive, lake-speci c studies 8,10 , it becomes apparent that our method may not detect all drainage events, nor the speci c timing of event (Fig. 9).In some cases, this is attributed to GLOFs occurring during periods when altimetry observations are limited, such as before the launch of ICESat-2.Additionally, drainage events may be missed due to a lack of observations, particularly with smaller lakes that receive only sparse altimetry coverage throughout the year (Fig. 9, Russel Glacier Ice-dammed lake).Expectedly, our largescale approach does not capture all drainage events or their precise timing, as in detailed lake-speci c studies.This emphasizes the trade-off between detail and analyzing a broader dataset spanning numerous lakes across Greenland.Consequently, our results on GLOF occurrences represent a conservative minimum estimate.
For some lakes, we observe time-series with unusual and abrupt water-level uctuations, such as the sudden increase observed in 2011 at Lake Tininnilik, which drained in 2010 (Fig. 9).These instances primarily occur following drainage events when the lake area is reduced.During lower water levels, we include numerous altimetry observations from outside the post drained lake area, i.e. from the surrounding bathymetry/bedrock.Consequently, observations from the same over ight often exhibit a much larger internal variance, whereas observations from different days may signi cantly differ due to variations in the drained lake bathymetry.In some cases, re ning the lake outline for a more precise spatial ltering of the observations considerably improves the water level time series.Lakes showing a general decrease in water level, which typically coincides with a decrease in lake area, may also be in uenced by spatial ltering issues.In addition, if the lake-terminating glacier front has undergone substantial retreat or advance during the observational period, there is a risk that measurements might include elevation observations of the ice rather than the lake water level.However, decreasing the lake size too drastically, may end up in excluding observations with potential important information.Striking the right balance in selecting the most appropriate criterion for outlier detection is crucial.Being overly conservative in the statistical outlier removal may also result in the inadvertent exclusion of abrupt water level changes as outliers.
Automated water detection frameworks applied to optical images, commonly used in studies of icemarginal lake 2,3 , rely heavily on cloud-free images with a distinct water body re ectance.Thus, these approaches would likely not have been able to capture the observed changes in Fig. 5 as the lakes are often covered by ice or snow, which is not uncommon for ice-marginal lakes located in the northern regions.This limitation emphasizes one of the primary advantages of using altimetry data for largescale, Greenland-wide studies of lake changes, as it is largely independent of clouds, snow, ice, and lake water turbidity.

Annual runoff
The mean annual runoff are determined at Greenland wide and at regional scale using the RACMO2.3p2with a 1km resolution 30 .

Ice dam elevation
The elevation of the ice dams is calculated using the PRODEM 31 dataset, which contains annual summer elevations of the ice sheet marginal zone between 2019 and 2022 at a 500 m resolution 32 .For all lake exhibiting GLOF behavior, we identify the centroid point and extract the annual elevation value from the closest PRODEM pixel, located 1500 m from the lake centroid to avoid erroneous observations from the frontal, uctuating part of the ice margin.For each lake, we rank the ice dam elevations from one to four, assigning the highest rank (1) to the year when the ice dam had its maximum elevation.Finally, we calculate the mean annual rank of all ice dams.This implies that a higher mean annual rank corresponds to thinner glaciers.Time-series of water level and lake area changes for luliallup Tasia, Lake Isvand and Lake Hullet.These lakes have previously been found to produce substantial GLOFs 22 .Errors bars indicate the STD of all measurements.

Declarations Data and code availability
Evolution of lakes dammed by Budol Isstrøm.A) series of water level and lake area changes.Note the concurrent GLOF of lake 29377 and in 2017.B) Optical Sentinel-2 images the observed max.and min.lake area of each lake from 2017-2023.Due to ice and snow cover the lakes, the minimum and maximum lake extent may be di cult to determine.observed lake changes would not be detected using automated water detection frameworks applied to optical images, which highlights one of the primary advantages of altimetry data.Note the general agreement between uctuations in lake area water level.Variations in runoff and annual number of GLOFs.A) Annual changes in GLOFs and runoff.Runoff anomaly is relative to a 2000-2022 baseline.B) and C) show the covariance between the annual number of GLOFs and annual mean runoff at B) a Greenland-wide level and C) within and across regions.For cross-regional comparison, we the annual GLOFs and runoffs by calculating the ratio of the total amount across all years (2019-2022) within a region.Thus, the sum each region to one on both axis.Note that the SE region is not included in the correlation plots (B and due to the amount annual GLOFs.

Figures
Figures

Figure 1 Overview
Figure 1

Figure 2 Number
Figure 2

Figure 3 Overview
Figure 3

Table 1
Overview of the types of altimetry data included