1. Introduction
The volume, distribution, and movement of subglacial water under the Greenland Ice Sheet is difficult to quantify and, as a result, is poorly documented. Nevertheless, one component of subglacial hydrology, the study of subglacial lakes, has recently received attention. After the initial discovery of two subglacial lakes in northwest Greenland (Palmer and others, Reference Palmer2013), recent summaries (Bowling and others, Reference Bowling, Livingstone, Sole and Chu2019; Livingstone and others, Reference Livingstone2022; Fan and others, Reference Fan2023) have shown that more exist and that often the presence of subglacial lakes can be inferred from the measurement of surface elevation change when other possible causes of the change, e.g., ice dynamics, can be discounted. The lakes with assumed episodic water in- and outflow are termed ‘active’, but other subglacial lakes have been identified using radio-echo sounding without any apparent periodic change in ice elevation above the lake (Palmer and others, Reference Palmer2013; Livingstone and others, Reference Livingstone2022). Although the discovery of subglacial lakes is relatively recent in Greenland they have been known and studied in Iceland for many decades. For example, geothermal energy creates a subglacial lake in the Grímsvötn caldera, Vatnajökull, leading to periodic outburst floods, or jökulhlaups (Björnsson, Reference Björnsson2003; Guðmundsson and others, Reference Guðmundsson, Sigmundsson, Björnsson and Högnadóttir2004).
One of the first investigations linking vertical ice surface displacements to episodic subglacial water movement was achieved in West Antarctica with interferometric SAR data from overlapping ascending and descending passes coupled with shifts estimated from speckle tracking (Gray and others, Reference Gray2005). This link to episodic movement of subglacial water to- and from subglacial lakes was further proven in Antarctica by subsequent work using ICESat-1 (Fricker and others, Reference Fricker, Scambos, Bindschadler and Padman2007; Smith and others, Reference Smith, Fricker, Joughin and Tulaczyk2009). Sudden drainage of a subglacial lake in western Greenland (Howat and others, Reference Howat, Porter, Noh, Smith and Jeong2015; Palmer and others, Reference Palmer, McMillan and Morlighem2015) was documented by using a combination of satellite imagery (Landsat and Worldview), laser altimetry (ICESat-1) and NASA IceBridge flights. More recently Andersen and others (Reference Andersen2023) also used repeat-pass interferometry with Sentinel-1 data to identify what they describe as episodic drainage cascades beneath the Northeast Greenland Ice Stream (NEGIS; Fahnestock and others, Reference Fahnestock, Abdalati, Joughin, Brozena and Gogineni2001).
While these studies have been revealing and helpful in visualizing the behavior of subglacial water, we should acknowledge that the picture may be complicated by ice dynamic effects. For example, nonsteady state subglacial water movement may change the basal drag and, as a result, the along-flow strain rate, which can then affect the vertical component of the ice surface displacement (Sergienko and others, Reference Sergienko, MacAyeal and Bindschadler2007). However, this process is relatively slow, potentially many times longer than the times associated with the significant height changes observed with active subglacial lakes or with the repeat-cycle of satellite repeat-pass SAR interferometry. Further, we cannot assume that the water volume associated with the episodic vertical movement represents the total volume of subglacial water moving under the ice or contained in a subglacial lake (Li and others, Reference Li, Lu and Siegert2020).
Improved radar altimeter data from CryoSat, the advent of ICESat-2, and the availability of time-stamped 2 m resolution ArcticDEM strips have enhanced our ability to detect and study active subglacial lakes in Greenland. Fan and others (Reference Fan2023) have now identified 18 active subglacial lakes under the Greenland Ice Sheet based primarily on ICESat-2 data. DEMs from the TanDEM-X program have also been used with CryoSat and ArcticDEM data to monitor active subglacial activity in Greenland (Sørensen and others, Reference Sørensen2024).
In this study we report on the discovery of three additional active subglacial lakes, based on the unusual history of height change in part of the ablation area of the northeastern Greenland Ice Sheet. We document what appears to be the largest active subglacial lake found to date under the Greenland Ice Sheet and show that the temporal and spatial resolution of CryoSat swath mode outputs (Gray and others, Reference Gray2013; Gourmelen and others, Reference Gourmelen2018) can provide a new window on the temporal flux of water into and out of active subglacial lakes. The temporal height change of the three lakes suggests that they are weakly hydrologically connected but the exact subglacial water pathways remain elusive.
2. Study area
The study area (Fig. 1) is in northeastern Greenland close to three large outlet glaciers which originate from the NEGIS and drain ~12% of the Greenland Ice Sheet (Khan and others, Reference Khan2022). The location of the study area containing the subglacial lakes is in the ablation zone and is indicated by the arrow in Figure 1a.
Figure 1. (a) Illustration of the location of the study area in the northeast of the Greenland Ice Sheet. (b) Surface topography is illustrated by color superimposed on a Landsat 9 image from August 2022 resampled to a latitude longitude grid. (c) The box in (b) is shown in (c) but now color is used to indicate the yearly ice movement in m a−1. The outlines of the lakes L1 and L2 are shown by the white lines and indicate the lake extent after outflow. The Lake L3 outline is dotted and indicates the approximate lake extent when water is present.
The area containing the subglacial lakes is in the peripheral region of the Greenland Ice Sheet between two large outlet glaciers (the Zachariae Isstrøm and Storstrømmen Glacier) which drain part of the NEGIS. Figure 1b shows that particularly lakes L1 and L2 are situated in areas where there are dips in the surface topography. Also, the ice movement (Fig. 1c) is very slow at the lake sites (Joughin, Reference Joughin2022). The bed in this area is part of the mountain range that rings this part of Greenland and to our knowledge there have been no remote sensing flights or fieldwork undertaken to measure ice thickness or bed elevation directly over the central L1/L2 area shown in Figure 1c. Inevitably this reduces the accuracy of the Bedmachine dataset (Morlighem and others, Reference Morlighem2017) and limits our ability to model subglacial water flow.
3. Methodology
3.1 Cryosat
CryoSat was developed by the European Space Agency (ESA) and launched in 2010 to provide radar altimetry data over the cryosphere, including the large ice sheets and smaller ice caps. The experimental SARIn (synthetic aperture radar interferometric) mode of the radar altimeter allowed geocoded footprints that were smaller than those possible with conventional radar altimetry. This new mode continues to provide key results for the periphery of the ice sheets and for smaller ice caps and glaciers around the world. The excellent satellite longevity has transformed this ‘proof-of-concept’ altimeter into a key contributor to the estimation of worldwide glacial ice loss (Otosaka and others, Reference Otosaka2023; Tepes and others, Reference Tepes2021).
We derived the CryoSat heights from ESA Level 1b data from fall 2010 to the end of 2022 using the processing described in Gray and others (Reference Gray2019) and Gray (Reference Gray2021). Swath mode results (Gray and others, Reference Gray2013; Gourmelen and others, Reference Gourmelen2018) have been used because they can provide height data in local surface depressions, which are difficult to access with the conventional ‘point-of-closest-approach’ (POCA) processing. The delay-Doppler algorithm (Raney, Reference Raney1998) used to produce the L1b data provides waveform data every 300 m along-track. Each waveform has 1024 values for the power, differential phase and coherence as a function of delay time. The swath mode heights are estimated from the waveform values delayed after the waveform leading edge characteristic of the POCA.
Swath mode processing provides many more estimates of terrain height than with POCA processing, but care needs to be taken to minimize the contribution from range-ambiguous regions. Initially waveform values with power less than −165 dB, and coherence less than 0.85 are removed. Further editing is undertaken based on a comparison of the GIMP DEM (Howat and others, Reference Howat, Negrete and Smith2014) with the derived surface heights. While the off nadir look angle for swath mode heights is normally derived directly from the interferometric phase angle, occasionally one needs to add or subtract 2π radians to the phase to achieve the most likely solution. This is also based on the comparison of the three possible solutions with the GIMP DEM. In the current work special care was taken to ensure that likely real values were not eliminated if they differed significantly from the interpolated GIMP height. This was necessary because heights over the subglacial lakes could vary by many tens of meters over relatively short time periods.
The temporal height change over two of the lakes was calculated using the point-to-point comparison method described in Gray and others (Reference Gray2015 and, Reference Gray2019). In the current context this assumes that one can define an area on the ice surface above the subglacial lake for which the relative height does not change. The surface topography above the lake need not be flat but should rise and fall without a change in the surface shape. For Lake 1, the largest subglacial lake, the height difference between DEMs obtained in 2007, 2012 and 2021 is used to define the area which is assumed to be floating and does not change shape with time. There is a major outflow early in 2012 which reduces the area that is supported by water. All the CryoSat swath mode data are collected for this reduced area and separated into time bins spanning the period from the fall of 2010 to the end of 2022. Typically, the best temporal resolution is achieved by binning the data into a sequence of 28-day periods, which is close to the satellite orbital subcycle. Each time bin has coverage on 2 or 3 days which are normally separated by 2 days. Then each height estimate for the selected area in one time bin is compared with all the height estimates in another time bin. If the separation between the centers of the two estimates is less than a preset distance, here set as 200 m, the height change could be included in the average and std dev. of the height change between those two time periods. If the magnitude of the height change was outside the range of the mean plus or minus two standard deviations, then that value was rejected. The average and std dev. of the height differences for all possible combinations of time periods is calculated and used as described in Gray and others (Reference Gray2015). Errors in the CryoSat heights are dominated by the unknown contribution from range ambiguities, which are hard to quantify (Gray, Reference Gray2019; Haacker and others, Reference Haacker, Wouters and Slobbe2023). By comparing the CryoSat results with more accurate height change values from ICESat-2 and the smoothed 2 m ArcticDEMs, we can estimate that the potential error in the points in the CryoSat height change temporal histories. This is done on a case-by-case basis.
A variation of this method was also used to estimate the height and volume change associated with the February 2012 outflow from the largest lake. While the DEM difference approach provides excellent height difference at high spatial resolution, the dates at which the ASTER and ArcticDEMs are available may not allow height change estimation between the desired times. By identifying time periods before and after the first large outflow, we can create a low-resolution height change image from CryoSat data by repeating the point-to-point comparison method for a sequence of relatively small footprints across the test area. We use the 28-day height change plot from CryoSat to identify the periods with relatively small height change before and after the outflow and use footprints of size 0.01° in latitude (~1.1 km) and in 0.05° longitude (~1.15 km), sampled at one half these values. For the 2012 outflow event we used all the CryoSat data from 455 days prior to the outflow (9-Nov-2010 to 07-Feb-2012) and all the data from 94 days after the outflow (31 March 2012 to 3 July 2012). These data were then binned in the above footprints and the height change between the two periods estimated. The problem with this approach is that there are height changes within these time windows which can lead to potential errors in the result. Nevertheless, if the height change associated with the outflow is much larger than any height change within the before and after windows, the results will still be useful. If there are less than 10 height change values in the 28-daytime period for any footprint the average is not calculated leading to no data values in the associated plot or image. All the heights from the CryoSat processing are with respect to the WGS 84 ellipsoid, as are all the other heights quoted in the paper.
3.2 DEMs
Time specific DEMs were downloaded from the US Polar Geospatial Centre in the form of 2 m resolution strips (Porter and others, Reference Porter2022). For the study site the DEMs are available from 2011 to 2022 and have been generated primarily from the WorldView series of satellites using the Surface Extraction with TIN-based SETSM algorithm (Noh and Howat, Reference Noh and Howat2014). The DEMs have been smoothed, resampled to a 20 m grid and co-registered. In comparing the surface height change above the subglacial lake between different dates, we look for stable terrain on rock outcrops adjacent to the subglacial lake, pick a reference DEM and then adjust the other DEM to match the height at the reference area. It is not necessary to have complete coverage of both lakes in the two DEMs. If there is coverage of a common floating area, then we can estimate the height change between the DEM dates. These height changes are unlikely to be in error by more than 1 m.
The ASTER DEMs were generated using MicMac ASTER (Girod and others, Reference Girod, Nuth, Kääb, McNabb and Galland2017) using ASTER Level 1A reconstructed unprocessed instrument data downloaded from NASA Earth Data (www.earthdata.nasa.gov/). Again, the DEMs were resampled to a common grid in the same polar stereographic projection as the 2 m ArcticDEMs. While the ASTER DEMs are reported to have a vertical uncertainty of 10 m and a horizontal uncertainty of 2 m (Girod and others, Reference Girod, Nuth, Kääb, McNabb and Galland2017), the possibility of a bias in ice surface height change between an ArcticDEM and an ASTER DEM has been minimized by using a reference area over stable terrain and adjusting the ASTER DEM to match the ArcticDEM over that area. However, the poorer quality and larger height noise of the ASTER DEMs is clear when we difference a resampled ASTER DEM with an ArcticDEM.
3.3 ICESat-2 data
We use the ATL06 dataset (Smith and others, Reference Smith2023) to provide geolocated surface heights above the WGS 84 ellipsoid. The data were acquired by the Advanced Topographic Laser Altimeter System (ATLAS) instrument on board the Ice, Cloud and land Elevation Satellite-2 (ICESat-2) observatory. Data are available from October 2018, but we have concentrated on results from 2019. As ICESat-2 is in a 91-day repeat orbit, the coverage is limited both spatially and temporally. However, there are eleven passes over the central part of the subglacial lake L1 in 2019 and, if the shape of the ice surface does not change with time, then it is possible to compare lake heights at time intervals less than the repeat period. The average height difference between the ICESat-2 data and the height obtained by interpolating a reference ArcticDEM to the position of the ICESat-2 data are calculated for the 2019 ICESat-2 passes over the central part of the subglacial lake. This provides a relative surface height change at the times of the ICESat-2 passes which can then be matched with the CryoSat height changes and referenced to the winter of 2010/11.
4. Results
Our analyses provide evidence for the presence of three subglacial lakes in the study area, primarily based on changes of ice surface elevation over time. The location of the lakes is provided in Figure 1, and we describe the results for each lake below beginning with Lake 1, the largest lake. Lake 2 (78.13N, 23.64W, 600–650 m) is higher in elevation than Lake 1 (78.11N, 23.02W, 540–600 m), with Lake 3 (78.04N, 22.65W, 390–410 m) the lowest. All the heights are with respect to the WGS 84 ellipsoid.
4.1 Lake 1: height and volume changes
The difference between an ASTER DEM from 26 April 2007, and a smoothed ArcticDEM from 15 July 2021 shows a large drop in ice surface elevation (Fig. 2a) between 2007 and 2021. The central region delineated by the green line encloses an area of 42.3 km2 with an average height loss of 40 m indicating a volume change of ~1.7 km3. A larger height loss of ~54 m was obtained for the same area using a 2009 ASTER DEM suggesting that the surface height above the lake had increased by ~14 m between spring 2007 and summer 2009. The apparent noise in the DEM difference is partially due to the data limitations, particularly for the ASTER DEM which has larger height noise, and poorer spatial resolution and geocoding than the smoothed and resampled 2 m ArcticDEMs. When we compare an ArcticDEM from May 2012 with the July 2021 ArcticDEM the opposite trend is observed (Fig. 2b). Now the central area delineated by the black line in Figure 2b is lower in May 2012 than in July 2021. Also, the light blue area with a small change in height difference is smaller (~22.6 km2) than the red area in Figure 2a. This suggests that sometime after summer 2009 but prior to May 2012 there was a major outflow and ice that was supported by water prior to the outflow around the outer edges of the subglacial lake are now grounded and not supported by water.
Figure 2. (a) Height change between an ASTER DEM from 26 April 2007 minus the DEM from 15 July 2021. (b) The elevation difference between the ArcticDEM from 19 May 2012 and the ArcticDEM from 15 July 2021. The green line in both (a) and (b) is an estimate of the ice surface area that was supported by the subglacial lake prior to February 2012. The black line encloses an estimate of the surface ice area that appears to have been supported by the subglacial lake throughout the study period and encloses the same area as the white line around L1 in Figure 1c.
The Sentinel-2 color image in Figure 3a shows that there is a significant amount of glacial debris on and entrained in the ice surface which is moving up and down with the periodic subglacial water in- and outflow. The ArcticDEM from 15 July 2021 has been used as the reference DEM and the heights for the other DEMs are adjusted such that the average height difference of a common stable area between that DEM and the reference DEM is zero. The profile height plots in Figure 4 are interpolated from the corrected DEMs. Figure 3b illustrates the ArcticDEM strip from 15 July 2021 and the positions of the profiles from two ASTER and 14 ArcticDEMs. Height profiles from various dates are shown in Figure 4.
Figure 3. (a) Sentinel-2 color image from 25 August 2022 showing the subglacial Lake 1 post-outflow boundary in black and the pre-outflow boundary in green. (b) The ArcticDEM from 15 July 2021 of the same area, where the surface height color bar varies from 350 m (blue) to 850 m (yellow). The positions of the two outlines in Figure 2 are repeated here in black and green together with the positions of the six height profiles shown in Figure 4. Red dots mark every km along each profile and correspond to the plot km markings in Figure 4.
Figure 4. (a) to (f) Surface heights for different acquisition dates along each of the six profiles in Figure 3 for the 14 ArcticDEMs (2012–2021) and the 2 ASTER DEMs (2007 and 2009; dotted red and orange lines). The color for each height profile indicates the acquisition date.
The temporal change in the surface height above the subglacial lake can be estimated using CryoSat data. To get the best temporal resolution in height change we need to identify as large an area as possible which rises and falls as a unit and does not change in shape or tilt with time. The area identified with the black outline in Figure 2 has been used and includes a total of 16 513 swath mode height estimates spread over the CryoSat lifetime. These data have been divided into 156 time periods with an average length of 28.4 days (Fig. 5). Typically, the coverage in each time period occurs on 2 or 3 occasions separated by 2 days. On average there were 104 height samples in each time bin, but the number varied from 2 to 493. Results from 3 of the 156 bins were removed from Figure 5 as potentially unreliable due to the uneven distribution of samples.
Figure 5. Average surface height change of the area over Lake 1 enclosed by the black outline in Figure 2b with respect to the winter of 2010/11. The connected blue dots are the CryoSat results while the red dots are from analysis of the ArcticDEMs. The CryoSat and ArcticDEM results have been matched at the 19 May 2012 position as described in the text.
We use 14 time stamped ArcticDEMs to validate the CryoSat derived average height change plot. The average height of the same area as that used for the CryoSat analysis was calculated for each DEM even if it did not cover all the area. Then the height difference of this area between the 19 May 2012 DEM and the other 13 DEMs was calculated using a reference area on adjacent solid rock to remove any biases. An offset was added to all the differences such that the 2012 DEM difference value matched the interpolated CryoSat value on 19 May 2012. The agreement between the other 13 DEM differences, marked as red dots in Figure 5, and the interpolated CryoSat data in blue is good; the mean difference is 0.8 m with a std dev. of 0.4 m. This result supports our contention that the error in the CryoSat height change estimates in Figure 5 is unlikely to be greater than ±2 m.
The plot of the average height change of this area (Fig. 5) shows 2 large subglacial jökulhlaup-like events. The first began after 7 February 2012: by March 6 the average ice level had gone down by ~36 m, and by 4 May the level had gone down by ~49 m. Using the area defined by the black outline in Figure 2b these height losses correspond to a volume loss of over 1.1 km3 of water. However, the total volume lost will be larger because with a sloping bed at the lake edge there will be some water lost which was not accounted for in this calculation.
Figure 6 illustrates surface height change derived using CryoSat data between two time periods spanning the first major water outflow early in 2012. All the CryoSat data from 9 November 2010 to 7 February 2012 was binned into the 1.1 by 1.15 km footprints, as described in section 3.1, to create a pre-outflow pseudo DEM. This was compared with a similar post-outflow DEM created using CryoSat data from 31 March 2012 to 3 July 2012. From Figure 5 we see that the surface height is not constant particularly over the later time interval. However, the major height change occurs prior to the start of the second period after which the height begins to increase. This implies that the estimate of water outflow may be noisy as the samples from the second interval can occur anywhere in the 94 days which spans the dip in the height change plot centered around 15 May 2012 (Fig. 5). Some of the footprints had insufficient CryoSat data for a credible solution and these are not colored in Figure 6. The area enclosed by the white line is 86 km2 and the total volume loss is 3.1 ± 0.3 km3. This area is larger than that identified as floating prior to the 2012 outflow in Figures 2 and 3 based on the difference of two time-specific DEMs. Because the pre- and post-ouflow time periods are large, 455 and 94 days respectively, we can't assume that all this volume change happened in the short time period shown in Figure 5. Noting that the surface height for the area enclosed by the black outline in Figure 6 dropped by 37 ± 2 m in 28 days this is equivalent to an average outflow of 325 to 364 m3 per second. But from Figure 6 we know that this outflow estimate must be a fraction of the total outflow and therefore the peak outflow rate could be many times this estimate and be comparable to some of the peak discharges from the Grimsvötn jökulhlaups in Iceland (Björnsson, Reference Björnsson2003).
Figure 6. The surface ice height change spanning the first observed large subglacial outflow from Lake L1 is illustrated in color superimposed on a Landsat image. The white line is the boundary for the estimates of area and water volume change quoted in the main text and the black line outlines the area which we believe is floating throughout the studied time period.
Lake 1 refills at a rate approximately one third to one half of that at which it drained, but the filling is also episodic and occurs during most summers since 2012, likely triggered by the summer melt. Summer 2018 is the exception to this trend as there is no height increase at that time, instead the lake filling started in early April 2019 but experienced a large outflow towards the end of July 2019, possibly triggered by the large melt in the summer of 2019 (Fig. 5).
To further investigate the unusual filling and drainage patterns in 2019 we compared each of 11 2019 ICESat-2 passes over the central part of Lake 1 with the ArcticDEM from March 2020. By interpolating the DEM data to the position of the ICESat-2 data, the height differences between the two dates were determined and averaged. These differences are plotted in Figure 7 and again compare favorably with the CryoSat results. As the CryoSat data reflects the height change from the winter of 2011–2012 we need to add an offset to the ICESat height change data in Figure 7 to match the CryoSat data in Figure 5. The ICESat-2 results suggest that the filling began after 11 March and stopped by 9 April and occurred at an average of 0.34 m day−1. The maximum rate of the 2019 height decrease is 0.91 m day−1 from the CryoSat data (18 August to 15 September 2019) but slightly larger at 1.02 m day−1 with the ICESat-2 results (22 August to 9 September 2019).
Figure 7. Comparison of the 2019 height change of the central region of Lake 1 using CryoSat and ICESat-2 data.
After the 2012 outflow event the lake elevation had not achieved its pre-2012 level. The second large drainage event occurred in 2019, but in the interim, there had been drainage events in 2014 and 2016. The refill began in the late summer for most of these years.
4.2 Lake 2: height and volume changes
There are other areas in Figure 6 which are also undergoing periodic height changes characteristic of episodic subglacial water movement. The area labeled ‘L2’ in Figure 6 is at the center of active subglacial Lake 2. Figure 8 illustrates the difference in height between two ArcticDEMs, the 19 May 2012 DEM minus one from 14 April 2021, and shows the shape and outline of Lake 2. There are some other areas with height changes that appear also to be related to nonsteady state water movement, for example the blue area at 78N, 23.6W marked with an arrow.
Figure 8. The color overlay on the Landsat image illustrates height difference using the 9 May 2012 DEM minus the 14 July 2021 DEM. The white outline marks the area used for the CryoSat temporal Lake 2 height change in Figure 9 and the black outline marks the region used as the DEM reference area. The black arrow shows the position of another area with apparent episodic height change.
All the CryoSat swath mode data inside the white line in Figure 8 have been used to document the temporal change in surface height above Lake 2. The area is relatively small (7.4 km2) and the numbers of CryoSat height samples in some of the 28-day windows are too small to provide reliable results leading to the gaps in the surface height history (Fig. 9). There are sufficient ArcticDEMs covering this area that we can match 15 different DEMs to the CryoSat results and show that the height change between the various dates over Lake 2 are consistent with the CryoSat results in Figure 9. However, in this case the uncertainty in individual CryoSat height change estimates is larger than for Lake 1 but is unlikely to be greater than ±5 m.
Figure 9. Average temporal change in the surface height above Lake 2 with respect to winter 2010/11. The dashed blue line connects the CryoSat data points, and the red dots are derived from the height difference of pairs of ArcticDEM strips.
Figure 9 shows that Lake 2 experienced sudden height loss and outflow on two occasions, in late fall 2012 and in early 2019. The onset of these outflows corresponds in time to periods in which the height of Lake 1 is increasing in elevation (Figs 5, 7). The height changes derived from the DEMs and marked as red dots in Figure 9 have been registered to the CryoSat data by matching the June 2013 value to the CryoSat results.
4.3 Lake 3: height and volume changes
Figure 10 illustrates the area including Lake 3 using the height difference between the ArcticDEM from 2 April 2013 and the ArcticDEM from 15 July 2021 to illustrate the position and extent of the subglacial lake. The surface of Lake 3 is the lowest of the three and it is not as clearly defined as Lake 1 and 2. Lake 3 occupies the predominantly blue area in Figure 10 which has profile P2 running through it. Unfortunately, this area is too small to provide adequate CryoSat data to generate a temporal height change plot like Figures 5 or 9. However, surface heights were interpolated from the available time referenced ArcticDEMs. By inspection of the height profiles in Figure 11, particularly profile P2, we see an elevated surface on two occasions from the DEMs of 29 July 2020 and 15 July 2021. At the other times covered by these ArcticDEMs the surface heights are relatively stable. Consequently, we suspect that this area may be different from Lakes 1 and 2 in that there may not be significant water under the ice at the lowest elevations to designate it as a persistent subglacial lake, but rather as an intermittent subglacial lake.
Figure 10. The color overlay on the Landsat image illustrates the height difference between the ArcticDEM from 2 April 2013, minus the ArcticDEM from 15 July 2021 when the Lake 3 level was higher than in 2013. The white lines indicate the positions of the four height profiles in Figure 11.
Figure 11. a to d: Historical surface height from the 4 profiles, P1 to P4, over Lake 3 shown as white lines in Figure 10. Profile 1 (a) runs from the land out onto the ice, profiles 2 (b) and 3 (c) are north to south and profile 4 (d) runs from the north and ends on the glacial debris feature also visible in Figure 1c.
The lowest surface elevation in this area is ~330 m and is interpolated from the earliest of the ArcticDEMs (2 April 2013) along profile P1 (Fig. 11a), close to the land boundary. This profile crosses what appears to be an exposed channel close to the land. From profiles 1 and 4 we see that the surface height of the part of the channel adjacent to the land increases from 2 April 2013 to 29 July 2020, then it decreases on 3 May 2021, and increases again to a high level on 15 July 2021. We have included in Figure 11a a profile from the July 2007 ASTER DEM which also shows a relatively high elevation.
Profile 3 (Fig. 11c) crosses the glacial debris feature adjacent to the lake L3 boundary shown in Figure 1c. This feature is also covered by profile 2 (Fig. 11b) at distance 3.7 km. The debris feature appears to be entrained in the ice surface and to move up and down with the ice, and even move along profile 3 (1.25 to 1.3 km in Fig. 11c). Profile P4 (Fig. 11d) initially follows along the channel adjacent to the land and then along the debris feature. The height variations adjacent to the land mirror the variations seen in Figure 11a, but the surface height does increase to ~430 m along the debris feature at the end of profile 4.
5. Discussion
The nature of the height change history of the central region of Lake 1 (Fig. 5) implies both episodic filling and outflow. The central region of Lake 1 is at local elevation minimum (Fig. 4, profile E), and in the absence of inflow or outflow there would not be a significant surface elevation change in any summer even with significant melt. Indeed, as the same mass of water occupies less volume than that of ice, in this situation there would be a small decrease in surface elevation. However, Figure 5 shows clear height increases, usually greater than 10 m, in almost all the late summers except for 2018. These height increases must be due to significant water inflow from higher elevations. As the area is in the ablation zone, we expect that the cumulative yearly height loss due to runoff and accumulation will be negative but offset by some height gain due to ice dynamics. Using the regional climate model MAR (Fettweis and others, Reference Fettweis2013), local cumulative yearly accumulation and meltwater runoff in this region varies from a low in 2018 (−0.36 m) to highs in 2012 (−2.0 m) and 2019 (−2.56 m).
The surface height history over Lake 2 (Fig. 9) also provides clear evidence that most of the water must be coming from regions with a larger hydraulic pressure, presumably higher elevations. The lake filling between the two sudden outflows is continuous, although there appeared to be some upward jumps in Figure 9 possibly related to enhanced input due to summer melt. It took 6 years to rise nearly 50 m i.e., on average ~8 m per year. Using the MAR model, the cumulative yearly surface accumulation and melt for the Lake 2 area varies from around −0.5 m to at most −2 m water equivalent, so there must be significant water coming in from higher elevations, either from the west or north-east. It is possible that the lake is getting water from a region of high geothermal heat flux, although a recent evaluation of seven heat flow maps does not indicate this (Zhang and others, Reference Zhang2024).
While Lake 2 outflows do correspond to times at which Lake 1 is filling, there is not a strong link between the in- and outflows of Lakes 1 and 2, as shown by the surface height time histories in Figures 5, 9. The large 2019 outflow from Lake 1 may have lead to the height increase in Lake 3 shown by the P2 profile from the 29 July 2020 ArcticDEM (Fig. 11b).
The dearth of reliable ice thickness and bed elevation data in this area precludes quantitative modeling of subglacial water movement. However, given the importance of surface slope in dictating water driving potential at the glacier bed (Shreve, Reference Shreve1972), surface elevation and elevation change data allow us to speculate on possible water pathways through this region. It does appear possible that water at the bed may flow from the region of higher speed close to the NEGIS towards Lakes 1 and 2 simply because the surface does tend to slope towards the subglacial lakes. Further upstream in the NEGIS there are troughs at the edges of the ice stream (Franke and others, Reference Franke2020), which would tend to contain any subglacial water under the ice stream, but that may not be the case for all the boundary in this region. In the recent study of the episodic water flow under the NEGIS (Andersen and others, Reference Andersen2023), the estimated water volume for the drainage cascades was in the range 106 to 107 m3, much less than our estimates of water traversing the Lake 1, Lake 2 region. For example, the spring 2019 Lake 1 filling of just the central region of Lake 1 (11 March to 9 April, Fig. 7) involved the input of ~2.2 × 108 m3 of water and the subsequent outflow (24 July to 9 September 2019) from the same region was ~8.5 × 108 m3 of water. Unless the water volume for the drainage cascades is only a fraction of the total subglacial water flow under the NEGIS, it appears unlikely that we can ascribe the NEGIS as a significant source for additional water input to our area. However, Mouginot and others (Reference Mouginot, Bjørk, Millan, Scheuchl and Rignot2018) have suggested that some of the water associated with the slow surge behavior of Storstrømmen and L. Bistrup may come from significant ice melt at the base of the NEGIS.
Figure 11a shows that in profile P1 there is a large variation in surface elevation in what appears to be a trough cut adjacent to the land. Note that the surface elevation in this trough is relatively high in 2007 (~395 m) but at a minimum (~325 m) in April 2013, and that it increases with time until 2020. We suggest that this was part of the pathway for the large 2012 outflow from lake L1 that cut through ice and increased the depth of this trough, which then began to refill by the freezing of water. While the subglacial pathways remain elusive, we suggest that the temporal height change of this trough adjacent to the land could be related to the large 2012 outflow.
Figure 11b, profile P2, shows that the input subglacial water flow to lake L3 was greater than the outflow sometime prior to July 2020 and again after 3 May 2021 but before 15 July 2021 and that this produced a rise in surface elevation over a quite wide area as shown in Figure 10. Figure 11c illustrates the variation in height across profile P3, which includes the debris arm visible in the Landsat image in Figure 10. Note that although the ice is essentially stagnant here the debris feature is supported by ice and may have moved slightly during the 8 years, as implied by the movement of the height modulation at 1.25 to 1.3 km in profile P3. This is probably the reason that the difference DEM image (Fig. 10) shows height modulation as surface features move slightly between the times of the DEMs.
6. Conclusions
Height change data derived from swath-mode processed CryoSat data have lead to the discovery of what appears to be the largest active subglacial lake found to date under the Greenland Ice Sheet. The reliability of the all-weather, day-night CryoSat data has been central in creating the temporal surface elevation record in a way that is impossible with other satellite sensors. Lake 1 shows periodic in- and outflows, and the rate of the outflows implies that they could be characterized as subglacial jökulhlaups. Subglacial Lake 2 is higher in elevation, ~650 m in comparison to ~580 m for Lake 1, and smaller in size but has also shown two jökulhlaup-like outflows separated by over 6 years. The Lake 2 filling between the outflows appears to increase steadily with time and did not show strong summer meltwater inputs although there are some summer height increase which may reflect local summer melt. As the yearly Lake 2 surface height increase is larger than we would expect from local summer melt, the input water source remains uncertain. Lake 3, the lowest and smallest subglacial lake, is different from the other persistent subglacial lakes in that there may not be water under the area after the episodic outflow, hence usage of the term ‘intermittent’.
There is not a strong hydrological link between the subglacial lakes although we did note that at the time of the Lake 2 outflows Lake 1 surface elevation was increasing. Also, the increase in surface elevation of Lake 3 shown by the profile P2 from the 29 July 2020 ArcticDEM could reflect water input from the large Lake L1 2019 outflow.
The 2 m resolution time specific ArcticDEMs have provided not only validation of the CryoSat height history but also high-spatial resolution and spatially extensive height change imagery between the acquisition dates. These DEMs can reveal small active subglacial lakes not possible with altimeters like CryoSat or ICESat-2. While the study of subglacial hydrology remains a challenging area, we hope that this work has shown the potential for complementary use of altimeters like ICESat-2 and CryoSat-2, coupled with high quality time specific DEMs like those available from the Polar Geospatial Center of the University of Minnesota.
Acknowledgements
This work was supported by the European Space Agency through the provision of CryoSat-2 data. NSIDC facilitated the provision of ICESat-2 data. DEMs were provided by the Polar Geospatial Center under NSF-OPP awards 1043681, 1559691 and 1542736. C. Dow was supported by the Canada Research Chairs Program (950-231237) and D. Burgess was supported by the Climate Change Geoscience Program, Lands and Minerals Sector, Natural Resources Canada. Support for L. Copland was provided by a Discovery award from the Natural Sciences and Engineering Research Council of Canada. We thank Dr Guðmundsson, an anonymous reviewer, the editor, Dr William Colgan, and associate editor Dr Ralf Greve for their constructive help in improving this paper.
Open access
