Lake outlines and identification 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 identified as ice-marginal. This identification 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 km2, 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 fluctuations in lake water levels10,13,28. All ICESat-2 data was downloaded from the NSIDC homepage29 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 specific 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.
Table 1
Overview of the types of altimetry data included
Mission | Dataset | Temporal coverage |
ICESat/GLAS | GLAH06 | 20 February 2003 to 11 October 2009 |
IceBridge ATM | ILATM1B V1 | 31 March 2009 to 8 November 2012 |
ILATM1B V2 | 20 March 2013 to 20 November 2019 |
ICESat-2 | ATL06 | 14 October 2018 to 16 April 2023 |
Estimation of water level time series
All altimetry datasets are merged and spatially filtered 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-specific 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:
-
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.
-
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 define 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).
-
If any Observation median falls outside of the upper or lower bounds, we mark them as potential outlier. However, we do not filter the observations immediately, but perform additional test to check if they are actual outliers:
-
If over 30% of measurements during a day are within the threshold, we remove measurements outside and recalculate a new Observation median (Fig. 8)
-
If over 70% fall outside, we test the validity of the measurements. We filter 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 final 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 level – min. 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 fluctuations. 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 filtering using a negative buffer on the original lake area. In rare instances, we refer to optical imagery to confirm whether specific altimetry observations are indeed outliers. Following this additional filtering, 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 significant 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 identified. 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 define 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 significant 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 findings, we conduct a comparison between the water level fluctuations 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 confirm whether a drainage in the altimetry-based water level coincided with alterations in the lake area. Out of these, we successfully confirm 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 fluctuations 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 fluctuations of the damming glaciers, or the presence of large floating icebergs in the lakes.
When comparing the generated water-level time series with findings from comprehensive, lake-specific studies8,10, it becomes apparent that our method may not detect all drainage events, nor the specific 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 large-scale approach does not capture all drainage events or their precise timing, as in detailed lake-specific 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 fluctuations, 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 overflight often exhibit a much larger internal variance, whereas observations from different days may significantly differ due to variations in the drained lake bathymetry. In some cases, refining the lake outline for a more precise spatial filtering 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 influenced by spatial filtering 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 ice-marginal lake2,3, rely heavily on cloud-free images with a distinct water body reflectance. 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 large-scale, 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.3p2 with a 1km resolution30.
Ice dam elevation
The elevation of the ice dams is calculated using the PRODEM31 dataset, which contains annual summer elevations of the ice sheet marginal zone between 2019 and 2022 at a 500 m resolution32. 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, fluctuating 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.