The study region covers ~42,000 km² of West Greenland, which has an extensive and active supraglacial hydrological system, including a high concentration of supraglacial lakes (Lewis and Smith, 2009; Selmes et al., 2011; Figure 1). Previous remote sensing studies in the area have quantified supraglacial lake area and volume (e.g., Box and Ski, 2007; Georgiou et al., 2009; Williamson et al., 2017) and found that lakes have been forming and draining at higher elevations as melt rates have increased in recent years (Selmes et al., 2011; Howat et al., 2013; Fitzpatrick et al., 2014; Pope et al., 2016), particularly in very warm years, when lakes also drain earlier (Liang et al., 2012). The study region was chosen to match the total extent of the selected Landsat-8 scenes (paths 006, 007, 008, 009, and 010, and rows 011 and 012; Figure 1). Although not all of the images always covered the full study region, all scenes were included to give the best possible temporal resolution (Supplementary Tables 1, 2).

Landsat and Sentinel-1 SAR images were acquired between mid-June and late September 2015 to allow for the detection of summer lake extent and drainage. Additional Sentinel-1 images from October to November 2015 were used to extend the analysis beyond the summer. Figure 2 displays the processing steps for summer lake detection, from the acquisition of the raw imagery to the creation of a binary lake mask. A total of 13 atmospherically corrected Landsat surface reflectance images with <10% cloud cover from 20 June to 25 August were downloaded (Figure 2A; Supplementary Table 1). Scenes from the same day were mosaicked, and the bedrock and ocean were cropped using a mask created from the Greenland Ice Mapping Project (GIMP) Digital Elevation Model (DEM) (Howat et al., 2014).

Sentinel-1 Level-1 images were used since they are georeferenced but minimally processed. The Extra-Wide swath (400 km with 5 sub-swaths) provided the best coverage without compromising the spatial resolution. At the start of the study period, only medium-resolution images were available due to the recent launch of the satellite, so throughout the summer period the highest possible resolution was selected (Supplementary Table 2). Dual-polarization images (HH and HV) were chosen in the Ground-Range Detected mode. A total of 28 summer images spanning 15 June to 30 September was obtained. We used the Sentinel Application Platform (SNAP) toolbox to carry out a radiometric calibration, a single product speckle filter and a topographical correction on each raw Sentinel-1 scene (Copernicus, 2016). This last process used a Range Doppler terrain correction operator, which corrected the terrain and image orientation using bilinear interpolation of the GETASSE30 DEM. The image backscatter values were converted to decibels (dB) using:

(Figure 2D). The two polarizations (HH and HV) of each image were subsequently treated as separate images, rather than as a single scene with two bands. Concurrent images were mosaicked, and the bedrock and ocean were cropped using the GIMP DEM mask (Howat et al., 2014). For some of the Sentinel-1 scenes acquired on the same day, pre-processing resulted in a gap between the two mosaicked edges; the resultant spaces of under a pixel width between the two images were treated as “no data” values, and accounted for in subsequent processing.

Surface lakes were detected in the Landsat imagery using the Normalized Difference Water Index (NDWI), a well-established band-math technique that classically uses a ratio involving the red (B4, 0.64–0.67 μm) and blue (B2, 0.45–0.51 μm) bands (Huggel et al., 2002; Xu, 2006; Doyle et al., 2013; Morriss et al., 2013; Yang and Smith, 2013; Moussavi et al., 2016):

A threshold value of 0.5 was selected empirically as sufficient to define lakes without including surrounding water or slush areas. Sensitivity to this threshold was tested qualitatively by comparison with the original red-green-blue (RGB) false-color images. Higher thresholds produced smaller lake areas, but excluded shallow water that was obviously associated with lakes. Lower thresholds gave larger lake areas but included many pixels that looked like streams or isolated slush patches. The NDWI threshold of 0.5 was applied to all Landsat scenes to create a lake mask for each image (Figure 2B). All of the binary images were then superimposed to create a composite mask of Landsat-derived maximum lake area over the study period. This composite maximum lake area mask was then dilated by two pixels (i.e., 60 m; Figure 2C).

The dilated mask was applied to each Sentinel-1 scene to produce, for each image, the backscatter values of all pixels that either: (i) appear as surface lakes in the Landsat scenes at some stage of the summer; or (ii) immediately surround those pixels (Figure 2E). The histogram of each masked Sentinel-1 image was then examined and, where a bimodal distribution existed between the lower lake backscatter values and the higher surrounding slush or ice backscatter values, the lowest point between the peaks was selected as a backscatter threshold value to define lake area for that particular image (Figure 2F). The two-pixel buffer value was chosen empirically by application to several images. Sensitivity to the buffer value was tested qualitatively by viewing the resultant backscatter histograms. A smaller buffer obscured the bimodal distribution in some images, whereas a larger buffer did not improve the clarity of the bimodal histogram in any image. Applying a buffer was necessary due to slush surrounding some lakes having similar backscatter values to the lakes, thereby complicating the differentiation of slush and water (Bindschadler and Vornberger, 1992). Threshold values were picked manually from each histogram, and varied from −16 to −29 dB for the HH images, and −22 to −29 dB for the HV images (Supplementary Table 2). The masked Sentinel-1 images were thresholded by their unique saddle value to produce binary images akin to the thresholded Landsat images (Figure 2G). Of the original 28 Sentinel-1 images, 23 were suitable for further analysis, since if no bimodal histogram was present, the image was discarded. Finally, for each time period, the binary images of the two polarizations were superimposed to produce a total Sentinel-1-derived lake area (Figure 2H).

Sensitivity to the Sentinel-1 backscatter threshold values was evaluated for three dates during the summer: days 189, 212, and 237. For each image, its unique threshold value was shifted up and down in 1 dB increments to examine the effect of varying thresholds on total lake area (Supplementary Figures 1, 2). For all images, sensitivity is approximately linear over most of the range but the magnitude varies between the images. On day 212, the sensitivity is ~30 km² dB⁻¹; on day 189 and 237, it is approximately five times greater at 150 km² dB⁻¹ (Supplementary Figure 1). Compared to the chosen threshold value (Supplementary Figure 2B), lowering the threshold diminishes the number of pixels designated as lakes, so some lakes appear patchy (Supplementary Figure 2A). Raising the threshold value increases the number of pixels depicted as lakes, with more lakes identified and each covering a larger area up to the full extent of the mask (Supplementary Figures 2C,D).

On the three dates used for the threshold backscatter sensitivity analysis, there were contemporaneous images from both the Sentinel-1 and Landsat sensors with some areal overlap. For these dates, the lake areas detected from each type of imagery within the overlapping regions were compared, forming the basis of our evaluation of the Sentinel-1 classification scheme with respect to the Landsat NDWI technique, and enabling the most suitable Sentinel-1 polarization to be determined.

Surface lake drainages were easily distinguished from the Landsat imagery: Figures 3A–C show lakes that either disappear rapidly between two images (red and green circles), or more slowly over several consecutive images (yellow circles). These drainages were also detected in Sentinel-1 imagery. Areas of low backscatter created by lakes in one image change to values that are similar to the surrounding ice once the water has drained (Figures 3D–F, yellow circles). An intermediate stage was observed in some instances, where the low lake backscatter flips to very high backscatter values in the next image, before then changing to values that are similar to the surrounding ice (Figures 3D–F, red and green circles). This succession was also noted by Johansson and Brown (2012) using ENVISAT ASAR. We interpreted this sequence as a lake, then a rough heterogeneous surface consisting of patches of slush and ice left behind after the initial drainage producing the high backscatter values, and then a smoother, more homogeneous ice surface once the water within the slush drains and the snow melts. This sequence did not always feature in the imagery, however, so could not be used to inform our drainage detection.

Our drainage detection methodology was used on three datasets: (i) only Landsat data; (ii) only Sentinel-1 data; and (iii) the combined Landsat–Sentinel-1 dataset, using images from 20 June to 30 September. These three separate analyses allowed us to assess the extra information about lake filling and draining that can be obtained using Sentinel-1 data instead of, or in addition to, Landsat data. Our lake drainage detections were derived using Williamson et al.'s (2017) Fully Automated Supraglacial lake Tracking (FAST) algorithm, which was originally developed for MODIS imagery. The FAST algorithm creates a maximum lake extent array from all available summer lake extents from the relevant imagery, and then tracks changes to the individual lake areas within these maximum extents over the sequence of imagery. For each lake, a rapid drainage event was identified if a lake lost >80% of its area within ≤96 h. Although some rapid lake drainage events are known to occur in under 24 h (Das et al., 2008), this four-day threshold is in line with definitions for rapid drainage events used in previous remote sensing studies: two days (Selmes et al., 2011, 2013); four days (Doyle et al., 2014; Fitzpatrick et al., 2014; Williamson et al., 2017); five days (Liang et al., 2012); and six days (Morriss et al., 2013). Slow drainage events were defined as when a lake lost >20% of its area over any time period, after any previously identified rapidly draining lakes were removed from the dataset.

Lakes >0.25 km in diameter (areas >0.0491 km² when approximated as perfect circles) have been found to hold a sufficient water volume to hydrofracture through subfreezing ice 1–1.5 km thick on the GrIS (Krawczynski et al., 2009). We therefore only track lake areas that reach 55 pixels (0.0495 km²) in size at least once in the season, since we can be confident that these lakes might be able to initiate hydrofracture (Krawczynski et al., 2009). The smallest lake size typically used in MODIS studies is 0.125 km² (i.e., two 250 m pixels) due to its lower spatial resolution (Selmes et al., 2011, 2013; Williamson et al., 2017). Our algorithm differentiates between “large” and “small” lakes to assess the difference in lake behavior by lake size. “Large” lakes are those ≥0.125 km² (i.e., those resolvable in MODIS); “small” lakes are <0.125 km². The “date of drainage” was defined as the midpoint between the dates of the two images over which the lake disappears (following Doyle et al., 2013). We also calculated an associated error on the lake drainage date as half of the number of days between these two images.

The Landsat and Sentinel-1 imagery was used further to investigate the behavior of surface lakes at the end of the melt season. In the Landsat imagery, freeze-through was detected manually for each lake when it appeared optically similar to the surrounding ice surface (cf. Johansson et al., 2013; Selmes et al., 2013; Luckman et al., 2014). If a lake had not frozen at the end of the available imagery, a freeze-through date was not recorded. When new snowfall obscured freeze-through detection, the day on which the new snow obscured the lake was recorded. This occurred for lakes at higher elevations due to new snowfall observed on the 16 August.

In the Sentinel-1 imagery, freeze-through was also detected manually for each lake. Areas of low backscatter were apparent in the Sentinel-1 images after 16 August (Figure 4); these are interpreted as liquid water beneath a fresh, dry snow cover that is invisible to microwave wavelengths (König et al., 2001). SAR allows the time of freeze-through to be determined directly, as bubbles entrained within frozen lake ice increase the relative backscatter compared to liquid water (Hirose et al., 2008). The backscatter signal of unfrozen and frozen lakes is therefore sufficiently distinct to allow freeze-through identification (Johansson and Brown, 2012). We defined freeze-through as having occurred once the backscatter values within a lake basin were comparable to those of the surrounding glacial ice, similar to our method for detecting lake drainage. As the C-band penetration depth is limited to a few meters, we recognize that a lake may not have frozen through fully, but may have frozen through only to the depth of the SAR wave penetration.

In both the Landsat and Sentinel-1 imagery, freeze-through events were recorded for lakes that never drained, but also for lakes that partially drained earlier in the melt season; there are therefore lakes that appear in both the drainage and the freeze-through categories. The “date of freeze-through” and associated error were calculated in the same ways as for the lake drainage events defined above.

Late in the summer, small dark areas appear in the Sentinel-1 imagery at generally higher elevations than the surface lakes discussed so far (Figure 5B). In the Landsat scenes, these are either never visible (Figure 5A, red circles), or appear as large, ice-covered depressions with a small moat of water around their perimeter (Figure 5A, green circles). These areas of water would have been included in the surface lake detection process described in the Lake Detection section above. The dark areas in the Sentinel-1 imagery have similar backscatter values to the snow-covered unfrozen lakes (see Lake Freeze-Through section above), and are interpreted, therefore, as ponded subsurface liquid water. These features were delineated manually from the Sentinel-1 imagery only, as there were no contemporaneous Landsat images that would have allowed a similar detection procedure as for the surface lakes (see Lake Detection section above). The “date of appearance” and “date of disappearance” of each subsurface lake were determined manually as the midpoint between the images before and after the lake first appeared and was last present, respectively. The associated errors in these dates were calculated in the same way as for the surface lake drainages and freeze-through events defined above.

Figure 6 shows the total lake areas detected in each Sentinel-1 image, depicting separately the areas detected by the HH polarization, the HV polarization, and both polarizations. Some of the image-to-image variation in lake area was due to the varying image coverages and to the error in the total lake area (error bars, Figure 6). The HV polarization typically classified a much greater lake area than the HH polarization. Furthermore, pixels classified using the HH polarization were usually also classified by the HV polarization, so the HH polarization generally added few additional lake pixels to the classification based solely on the HV polarization. There were exceptions; the greatest of these was on day 190, where the HH polarization identified more unique lake pixels than the HV polarization.

Table 1 compares the lake areas detected by Sentinel-1 and by Landsat for the three days when there were overlapping images. For each day, the total lake area classified by the Landsat NDWI scheme was closer to the area produced using the Sentinel-1 HV polarization than the HH polarization, with the best match given when the lake areas produced from both the HV and HH polarizations were combined (Table 1). This validates the final stage in our Sentinel-1 lake detection methodology: the use of both the HV and the HH polarization detections in combination for defining lakes from Sentinel-1 imagery (Figure 2H).

Examples of the Landsat and Sentinel-1 lake detections on the three days with contemporaneous images are shown in Figure 7. The Sentinel-1-detected lake shapes, sizes and areas were generally similar to those detected in the Landsat scenes: for example, on day 189 (Figures 7A–C) and day 237 (Figures 7G–I). However, the larger lakes depicted by our Sentinel-1 scheme on day 212 were more poorly defined than those detected by Landsat (Figures 7D–F), while the Sentinel-1-derived lakes on day 237 (Figure 7I) showed pixels designated as water that were not shown as water in the Landsat images (Figures 7G,H). Overall, the Sentinel-1 lake detection calculated a greater total lake area than the Landsat NDWI classification within the same region, and the percentage differences between the lake detections varied enormously from 13.6 to 110.9% (Table 1). The total root mean square error (0.17 km²) was comparable to values found in previous studies, but the total root square error per km² (0.43, as defined by Liang et al., 2012) was typically an order of magnitude larger, although the previous studies were comparing two types of optical imagery (Sundal et al., 2009; Selmes et al., 2011; Liang et al., 2012; Williamson et al., 2017). The Pearson's r-value in Table 1 (0.92) shows that the lake areas determined by each method were highly correlated. The mean individual lake areas across all lakes detected from both Landsat and Sentinel-1 images were 0.235 and 0.115 km², respectively.

An assessment of lake areas with elevation was made as a further component of our Landsat–Sentinel-1 comparison. We found no systematic bias in the lake area detections with respect to elevation for either dataset. Early in the melt season, Sentinel-1 tended to overestimate lake area at higher elevations and marginally underestimate areas at mid-elevations compared to Landsat, although the differences were small.

There were three instances where the Sentinel-1 scheme had difficulty depicting lake areas; examples are shown in Figure 8. First, there were times when lakes were visible in Landsat, appearing, as expected, with low backscatter values in the HV polarization, but appearing with very high backscatter values in the HH polarization (Figures 8A–C). Second, there were times when lakes visible in Landsat were invisible in the HH polarization, but were seen in the HV polarization (Figures 8D–F, red circles). Third, there were times when lakes absent in the Landsat imagery appeared as low backscatter in the Sentinel-1 images, comparable to values representing water, in both the HH and HV polarizations (Figures 8D–F, green circles).

The combined Landsat–Sentinel-1 image dataset had a mean temporal resolution of 3.1 days, compared to 4.7 and 8.6 days for the Sentinel-1 and Landsat datasets, respectively. Consequently, the combined image analysis detected a greater number of drainage events with smaller date uncertainties than the single image analyses (Table 2; Supplementary Figure 3). The overall uncertainty in the mean drainage date for the combined Landsat–Sentinel-1 analysis (±4.0 days) was half that for the Landsat analysis (±8.0 days; Table 2). Compared to the Sentinel-1 analysis, the Landsat analysis detected more drainage events (Table 2), despite the smaller number of scenes (12 vs. 22) and the typically smaller coverage of each Landsat mosaic compared to each Sentinel-1 mosaic: the mean area compared to the full study region was 46 vs. 74%, respectively (Supplementary Tables 1, 2). This highlights the potential bias in the single imagery datasets that have a poorer temporal resolution compared to the combined dataset.

The combined dataset allowed large and small lake areas to be tracked; the mean lake size was 0.089 km². Compared to the total lake population (1,749 lakes), 11% of lakes were large and drained rapidly, 11% were small and drained rapidly, 31% were large and drained slowly, and 28% were small and drained slowly. The mean drainage date from the combined image analysis and the Sentinel-1 analysis alone is 16 July (±4.0 days). The date for the Landsat analysis is slightly earlier (13 July; ±8.0 days) and is affected by the dataset containing fewer images and ending earlier in the season (Table 2, Supplementary Tables 1, 2). The mean drainage date of the combined analysis shows a noticeable difference according to the lake size: compared to large lake drainages, small lake drainages occur an average of six days earlier for rapid drainage events, and three days earlier for slower events (with an overall mean of 4.5 days earlier for both types of event).

The mean elevation of all draining lakes is comparable for the three analyses: 1,100 m for the combined analysis, 1,170 m for the Landsat analysis, and 1,059 m for the Sentinel-1 analysis (Table 2). Overall, drainage tends to occur later in the summer with increasing elevation (Supplementary Figure 3). The combined analysis also shows that small lake drainage events occur at lower elevations than large lake drainages (overall mean 159 m difference). Specifically, small rapid events were observed to occur an average of 148 m lower than large rapid events, and small slow events were 169 m lower than large slow events.

A total of 621 lake freeze-through events was identified from the Landsat and Sentinel-1 imagery (Table 3). The Landsat-detected dates (Table 3; Supplementary Figure 4A) have a smaller range and larger errors (Supplementary Figure 4C) than the Sentinel-1-detected dates (Supplementary Figures 4B,D). The mean freeze-through date detected by Sentinel-1 (30 August) is over 2 weeks later than that seen in Landsat (12 August; Table 3). The mean elevation of freezing lakes is ~100 m higher than the mean elevation of lakes that drained (cf. Tables 2, 3), although there is no statistically significant correlation between lake freeze-through date and elevation for either set of imagery (Table 3).

The mapped subsurface lakes show a similar distribution and density to the buried lakes identified by Koenig et al. (2015) using IceBridge Ice Sounding Radar (Figure 9). We interpret our subsurface lakes to be similar features to the Icebridge-detected buried lakes. There are differences between the two datasets, which is to be expected as the data were acquired in different years. Koenig et al. (2015) note that buried lakes often appear in particular locations for just a single year, though can recur, and they map buried lakes for 2009, 2010, and 2012 (Figure 9A). In contrast, the subsurface lakes mapped in this study are from 2015 only (Figure 9B).

A total of 249 subsurface lakes was identified from the Sentinel-1 imagery (Supplementary Figure 5). The lakes tend to be at much higher elevations than the surface lakes described above, although some extend to the lower elevations where draining and freezing surface lakes occur. The mean elevation of the subsurface lakes (1,593 m; Table 4) is higher than that for the surface lakes that freeze through (1,270 m; Table 3). Eighteen of the subsurface lakes at the highest elevations (all above 1,725 m) were never visible in the Landsat imagery. The rest were partially visible in the Landsat data because of melting around their perimeters (e.g., Figure 5), which was detected as water by the NDWI classification.

The mean appearance date of the subsurface lakes is late in the ablation season (9 August), while the mean disappearance date is approaching the winter (7 October). As with surface lake freeze-through, there is no statistically significant correlation between the date of either subsurface lake appearance or disappearance with elevation (Table 4). The errors associated with the dates of subsurface lake appearance and disappearance increase at the highest elevations in the south of the study region (Supplementary Figures 5C,D) due to a restricted number of Sentinel-1 scenes over this area. This is also reflected in the mean error for the subsurface lake appearance and disappearance dates (±6.0 and ±5.5 days, respectively; Table 4), which is higher than the mean error of surface lake drainage dates (±4.0 days; Table 2) and freeze-through dates (±3.9 days; Table 3) from the Sentinel-1 dataset. The largest subsurface lake detected (15.9 km²) is more than double the size of the largest surface lake. Figure 10 shows the comparative sizes of surface and subsurface lakes for all of the identified lakes over the full study region.

Sentinel-1's ability to image through clouds and in the dark, together with its higher temporal resolution than Landsat and higher spatial resolution than MODIS make it a valuable tool in the mapping of surface lakes on the GrIS. In our study, the mean lake size identified from Sentinel-1 (0.115 km²) was smaller than the minimum lake size used in most MODIS studies (0.125 km²). We were able to detect 2,297 lakes within the Sentinel-1 dataset over a study region of ~42,000 km². This compares to the detection of 1,120 lakes across the whole of west-southwest Greenland (an area approximately three times greater in size) by Selmes et al. (2011) using MODIS.

Sentinel-1 has a dual-polarization capability and it is important to understand how lakes appear in the two polarizations, and the implications of this for mapping lake areas. We identified three instances (the examples shown in Figure 8) where particular care is needed in detecting lakes from Sentinel-1 imagery. It is known that several factors affect radar backscatter, including seasonal timing of acquisition, look direction, incidence angle, surface and near-surface moisture, dielectric constant, and the structure/composition of the surface features (White et al., 2015). Furthermore, these properties can affect different radar polarizations in different ways. First, we interpret instances of widespread and consistently high and uniform backscatter across lakes in the HH polarization but not in the HV polarization (Figures 8A–C) as a product of wind-induced surface waves on the lake. In these circumstances, the HV polarization can better map open water because the backscatter is less dependent on surface roughness (Scheuchl et al., 2004) and is largely independent of incidence angle and wind direction (Vachon and Wolfe, 2011).

Second, we interpret instances of invisible lakes in the HH polarization but visible in the HV polarization (Figures 8D–F, red circles) in terms of a particular structure/composition of the ice-sheet surface. These instances typically occur in the upper ablation area, where wet snow and firn give a low dielectric contrast between the lakes and surrounding areas. Furthermore, topography is undulating in these areas, suggesting that the backscatter from the surface roughness features masks that associated with the presence of water (Töyrä and Pietroniro, 2005). Finally, this effect is more common in images with a low incidence angle, suggesting that these causes are exacerbated at shallow look angles. The HV backscatter response is largely unaffected due to its greater signal-to-noise ratio (Partington et al., 2010; Nagler et al., 2015, 2016). Third, we interpret instances of dispersed areas of very low backscatter, typical of that associated with lakes but where no lakes are visible (Figures 8D–F, green circles), in terms of subsurface water. These instances typically appear in the upper ablation area and lower accumulation area, where contemporaneous Landsat imagery suggests that there may be deep snow-filled surface depressions, and where the later imagery often depicts the presence of surface water.

Despite the instances where Sentinel-1 has difficulty detecting lakes, particularly in the HH polarization, our surface lake detection method performs well at determining total lake areas that are comparable to those derived by the Landsat NDWI classification (Figure 7). The overestimation of the total area derived by Sentinel-1 compared to Landsat (Table 1) is likely due to Sentinel-1 detecting shallower lake water in addition to the deeper water identified by the NDWI classification with its threshold value of 0.5. A lower NDWI threshold would have delineated more water and may have reduced this lake area difference.

The Sentinel-1 lake detection method could be improved in several ways. First, the application of the Landsat-derived composite lake mask with a fixed dilation value may have obscured the bimodal histogram distribution between water and ice backscatter, which led to the discarding of five of the Sentinel-1 images. This effect was greatest at the start and end of the melt season, where smaller lake areas were surrounded by a greater number of ice pixels within each Landsat-detected maximum lake area. Our small buffer of two pixels minimized this although, for future work, the use of a dynamic dilation proportional to lake size through the melt season would reduce the number of ice pixels and may perform better. Alternatively, it might be possible to determine a threshold value directly from the Sentinel-1 images that distinguishes water and slush, suppressing the need for the Landsat mask. Second, we were unable to independently validate the Sentinel-1 lake detection method against the Landsat NDWI due to the use of the Landsat lake detections as a mask in the Sentinel-1 procedure (Figure 2E). Including additional high-resolution imagery, such as Sentinel-2 or WorldView, or lake perimeters from field data, would allow an entirely independent validation of the Sentinel-1 lake detection method. This would subsequently provide greater confidence in the Sentinel-1-derived lake areas. Third, higher effective spatial resolution could be obtained for the Sentinel-1 imagery if a finer-resolution DEM were used in the pre-processing stages, potentially resulting in the Sentinel-1 imagery having superior spatial resolution to Landsat.

The very small number of rapid drainage events, and the large errors associated with the timing of the slow drainage dates for the Landsat analyses highlights the much lower temporal resolution of the Landsat images compared to the Sentinel-1 and combination datasets (Table 2). The very small error associated with the rapid drainage dates (±1.0 days; Table 2) is therefore an artifact of the small number of events that could be detected. Many of the slow drainages detected by the Landsat analysis may have been included erroneously, as snowfall events would have prevented the NDWI from detecting water; these may consequently have been marked as slow drainages, with this effect potentially also affecting the combination analysis due to its inclusion of Landsat data. By comparison, the Sentinel-1 dataset detected many more rapid drainage events, as well as identifying their drainage dates with higher certainty (Table 2). Sentinel-1 can image through fresh snow, providing a further advantage over Landsat in allowing accurate determination of slow lake drainage events at any point during the melt season. The number of rapid drainage events identified from the combined dataset was greater still, and the error associated with defining the drainage dates was consequently small. The low mean errors produced from the Sentinel-1 and combination analyses reflect their superior temporal resolution compared to the Landsat analysis, and therefore their greater suitability for assessing rapid changes to surface lakes.

The temporal resolution of the combined dataset (3.1 days as a seasonal average) is approaching that of MODIS's daily resolution and may even exceed MODIS's effective temporal resolution due to MODIS's inability to image through cloud. Our rapid lake drainage results are comparable to those from previous work based on MODIS. Of our total lake population, 22% were defined as draining rapidly (11% large and 11% small lakes). The 11% large rapid drainage events is similar to Selmes et al.'s (2011) 13% for the whole of southwest Greenland over 2005–2009, Fitzpatrick et al.'s (2014) 28% for Russell Glacier over eleven melt seasons, and Williamson et al.'s (2017) 21% for the Paakitsoq region and 15% for the Store Glacier region for the 2014 ablation season, all based on MODIS. However, we observed an additional 11% rapid drainage events for small lakes (<0.125 km²), which would have been missed in these MODIS-based studies.

The mean rapid drainage dates obtained from our combined Landsat–Sentinel-1 analysis further highlights the advantages of this dataset over MODIS data. We observe that the mean drainage date of small, rapidly draining lakes occurs earlier (27 July ±1.3 days) than for larger lakes (2 August ±1.6 days; Table 2). Morriss et al. (2013) used MODIS imagery to study rapid drainage events from 2002 to 2011 over West Greenland, generally finding that many rapid drainage events occurred in early to mid-July. However, between 2004 and 2008 they observed a number of events clustering in early to mid-August, and in 2007 they even observed a cluster in early September. Furthermore, Williamson et al. (2017) obtain mean drainage dates during 2014 of 12 July and 7 July for the Paakitsoq and Store Glacier regions, respectively, from MODIS data. However, they note that few drainages were observed beyond 18 July due to widespread cloud cover in these regions. Thus, our combined Landsat–Sentinel-1-derived mean drainage date for large lakes is likely more accurate than MODIS-derived dates, which are affected by cloud.

The combined dataset detected a greater number of slow lake drainage events (57%) at a lower mean elevation when compared to the rapid drainage events, similar to the results of Selmes et al. (2013). Slow lake drainages also occurred across a wider elevation range than rapid events (Figure 10). This might point to a glaciological control preventing rapid drainage events at low elevations. Lake area and volume is primarily controlled by the ice-surface topography, which varies with ice-sheet elevation (Chu, 2014). Towards the ice-sheet edge, the surface is steeper than the interior, and large depressions are rare (Sundal et al., 2009). Lakes forming in these locations are more likely to overflow and drain slowly, perhaps because the depressions are too small to reach the critical water volume necessary to initiate hydrofracture (Krawczynski et al., 2009) or because water inputs via nearby crevasses/moulins are incapable of causing the basal slip or uplift necessary to generate tensile fracture across the lake basins (Stevens et al., 2015).

The combined Landsat–Sentinel-1 dataset allowed small surface lake areas to be tracked at the high spatial resolution of both sensors. Thus, lake areas were tracked that were smaller (mean 0.089 km²) than those resolvable in MODIS (minimum size 0.125 km²). The mean elevation of small draining lakes is 159 m lower than for large lake drainages, for both rapid and slow draining lakes (Table 2). Therefore, previous work using MODIS may not have considered the full span of GrIS lake drainages due to its low spatial resolution, limiting the minimum trackable lake size to the size of the “large lakes” in this study. Models of the basal hydrology or basal drag of the GrIS require accurate calculations of the locations, timings and magnitudes of surface water inputs (Schoof, 2010; Hewitt, 2013; Werder et al., 2013; Banwell et al., 2016; Koziol et al., 2017). Our work holds promise in providing a methodology for deriving a dataset that can be used to either prescribe accurate surface water inputs to ice-sheet models, or for testing algorithms of surface melt, routing and lake drainage that can provide such inputs.

Surface lake freeze-through dates could be manually identified from the Sentinel-1 imagery. This has been limited in previous studies that used optical sensors due to their inability to image late in the summer when cloud cover is high, solar illumination is low, and snowfall is common. Previous freeze-through studies have used proxies, such as the MODIS Land Surface Temperature product, to distinguish between freeze-through and slow lake drainage on the assumption that freeze-through occurs once the temperature drops below a defined threshold (Selmes et al., 2013). These approaches are therefore likely introduce errors. Accurate detection of freeze-through dates from the Landsat imagery in this study was impossible after the snowfall on day 228 (Figure 4). By contrast, Sentinel-1 can penetrate the shallow snow cover to return the backscatter signal of liquid water beneath, allowing more accurate freeze-through dates to be determined, even during periods of heavy cloud cover (Supplementary Figure 4).

The mean elevation of lake freeze-through was 1,270 m, similar to the ~1,500 m reported in 2005–2009 using MODIS across the whole of west-southwest Greenland (Selmes et al., 2013). Additionally, the mean elevation of freezing lakes is ~100 m higher than that for draining lakes (Tables 2, 3); in the MODIS study, this difference was ~200 m (Selmes et al., 2013). There was no statistically significant trend with elevation (Table 3; Supplementary Figure 4), suggesting that factors other than, or in addition to, the lapse rate control lake freeze-through.

The highest freezing lake elevation observed in this study was 1,712 m and, in general, lakes froze at higher elevations than those that drained (Figure 10). Sundal et al. (2009) found virtually no surface lakes above 1,600 m during the 2003 melt season, yet our study found a number of lakes above 1,600 m, some of which drained and the rest of which froze-through during 2015. The higher surface lake elevations observed in our study compared to Sundal et al. (2009) is in line with other work that indicates an increasing presence of water at higher elevations with increasing air temperatures (Bartholomew et al., 2010; Howat et al., 2013; Morriss et al., 2013; Fitzpatrick et al., 2014). If these higher elevation lakes increasingly drain rapidly rather than freeze-through, they could increase ice velocities in the short term if the water is accommodated in an enlarging subglacial cavity network (Zwally et al., 2002; Doyle et al., 2014), or decrease velocities in the long term if discharges are sufficient to cause channelization (Bartholomew et al., 2012; Sole et al., 2013). An improved understanding of the precise controls on hydrofracture would help to estimate whether surface lakes will be able to drain at these higher elevations in the future.

The only previous investigation of subsurface lakes used the IceBridge Snow Radar to assess subsurface lake depths, but it was unable to calculate their areal extents (Koenig et al., 2015). Sentinel-1 imagery has not only allowed subsurface lake areas to be measured, but has also allowed their seasonal lifespan to be calculated. Compared to the surface lakes that either drain or freeze, subsurface lakes appear quite late in the melt season, after the majority of drainage events have occurred (Tables 2–4; Supplementary Figure 5). Subsurface lakes freeze in the early winter, much later than the freeze-through of surface lakes, and are present at generally higher elevations than the surface lakes (Figure 10). The subsurface lakes extend up to nearly 1,900 m, with an elevation range of 1,179–1,882 m (Table 4); Koenig et al. (2015) also mapped the majority of their buried lakes between 1,000 and 2,000 m, and did not observe any above 2,000 m.

Although, the Sentinel-1-detected subsurface lakes show a similar distribution and density to those detected by IceBridge (Figure 9), the C-band SAR of Sentinel-1 is biased toward shallower subsurface features, while the L-band radar of IceBridge images deeper features. The similar distribution and density of these features for three years of IceBridge data and one year of Sentinel-1 data suggests that IceBridge only detected up to a third of the features due to its restricted flight line coverage, despite its greater penetration depth. Koenig et al. (2015) concluded that liquid water buried in lakes is present throughout the winter. Our work has allowed an analysis of lake behavior into the winter, and due to the late mean date of lake disappearance (7 October), we assume these lakes freeze-through. The discrepancy in the results may be due to the different penetration depths of the C-band and L-band radars, and/or the data extending to different dates in the autumn and early winter. The later freeze-through dates of subsurface lakes (7 October) compared to surface lakes (30 August) is likely due to the greater insulating effect of snow cover over the former compared to the latter. As also noted for the draining and freezing lakes, there is no statistically significant relationship between the date of subsurface lake appearance or disappearance with elevation, suggesting that factors other than or additional to the lapse rate control freeze-through. Thus, Sentinel-1 imagery allows the lifespan of the buried lakes to be determined, enhancing knowledge about the year-round distribution of liquid water within the surface layers of the GrIS.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.