Perennial Supraglacial Lakes in Northeast Greenland Observed by Polarimetric SAR

: Supraglacial liquid water at the margins of ice sheets has an important impact on the surface energy balance and can also inﬂuence the ice ﬂow when supraglacial lakes drain to the bed. Optical imagery is able to monitor supraglacial lakes during the summer season. Here we developed an alternative method using polarimetric SAR from Sentinel-1 during 2017–2020 to distinguish between liquid water and other surface types at the margin of the Northeast Greenland Ice Stream. This allows the supraglacial hydrology to be monitored during the winter months too. We found that the majority of supraglacial lakes persist over winter. When comparing our results to optical data, we found signiﬁcantly more water. Even during summer, many lakes are partly or fully covered by a lid of ice and snow. We used our classiﬁcation results to automatically map the outlines of supraglacial lakes, create time series of water area for each lake, and hence detect drainage events. We even found several winter time drainages, which might have an important effect on ice ﬂow. Our method has problems during the peak of the melt season, but for the rest of the year it provides crucial information for better understanding the component of supraglacial hydrology in the glaciological system.


Introduction
Each summer, numerous supraglacial lakes (SGLs) appear in the ablation zone of the Greenland Ice Sheet (GrIS). These lakes form from accumulated meltwater in surface depressions [1] and can influence the mass balance of the ice sheet in several ways. The lake water has a lower albedo than the surrounding glacier, which increases the absorption of shortwave radiation and causes additional melting [2,3]. Furthermore, the large amounts of liquid water stored in SGLs can suddenly drain via crevasses or hydrofracturing, reach the base of the ice sheet and locally increase ice velocities due to reduced basal traction [4][5][6]. However, not all SGLs drain each summer. As observed by optical imagery, several lakes freeze over at the end of the ablation season [7] and get covered by the first snowfall.
Most studies investigating SGLs on the GrIS use optical satellite images, e.g., from Landsat [8,9], ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer [10,11]), MODIS (Moderate-resolution Imaging Spectroradiometer [9,[11][12][13]), Sentinel-2 [14,15] or WorldView [16][17][18]. An automated detection of SGLs from optical data is usually performed via straightforward thresholding on a band ratio as the Normalized Difference Water Index (NDWI). However, there are several disadvantages of the detection of SGLs from optical imagery. These observations are limited by daylight conditions; hence, no information exists during the polar night. Furthermore, cloud-free conditions are necessary, which is especially problematic during spring and autumn.
Finally, the classification from optical data detects open water surfaces. However, when a lake freezes over and is covered by fresh snow, it cannot be identified any longer, but this does not necessarily mean that the lake has disappeared. Energy-balance and phase transition modeling demonstrated that under the climatic conditions on the GrIS, ice lids form up to a thickness of~1-3 m, but below that depth, the water stays liquid [19].
Synthetic Aperture Radar (SAR) is able to overcome the disadvantages of optical data and monitor SGLs all over the year. As an active sensor, SAR does not rely on sunlight. Furthermore, depending on the frequency, the signals can penetrate clouds and even snow cover. Johansson and Brown [20] and Miles et al. [14] used SAR to complement optical data in observing SGLs. However, both studies show that SAR and optical imagery often disagree. They found that SAR data alone are not suitable for an automated classification of SGLs and lake outlines. Therefore, they drew lake outlines manually in the SAR images or used polygons from optical imagery.
The interpretation of a SAR image is significantly more complex compared to optical imagery. The backscattered power from an ice sheet is sensitive to the dielectric properties of the material, such as wetness or changes in conductivity, but also depends on variations in geometric properties such as surface roughness, grain size distribution and the internal structure [21]. The C-band radar signals penetrate solid ice by 1-2 m and fresh snow up to 10 m [22]. As all these effects cumulate in the backscattered power measurement, a simple backscatter threshold is not suitable for the discrimination of SGLs. Additional observational information is provided by polarimetric SAR (PolSAR; see Figure 1b). In our region of interest, Sentinel-1 (S1) transmits the signal in horizontal polarization but records the return signal in horizontal polarization (HH) and in vertical polarization (HV). According to Freeman and Durden [23], measurements of a full polarized SAR system (HH, HV and VV) can be decomposed into the simple scattering mechanisms of single bounce, double bounce and volume scattering. In the absence of a VV channel for our S1 data, single bounce and double bounce scattering cannot be separated. However, double bounce scattering is usually attributed to targets such as trunks or buildings, which are not expected in our area of interest. Hence, we consider the HH signal as a measure for surface scattering and the HV signal for volume scattering (compare Figure 1c). With the help of this additional information, we developed an automatic processing which is able to detect SGLs and other surface types on the GrIS from S1 PolSAR alone. This allows results which are independent of optical data, and can, hence, be interpreted and compared without any bias. We focus on a region of north-eastern Greenland covering the major outlet glaciers of the North East Greenland Ice Stream (NEGIS, Figure 1a). This region is of particular interest as it drains about 12% of the GrIS [24] and is expected to host the greatest inland expansion of SGLs by the end of the 21st century [25]. At its main outlet glaciers Nioghalvfjerdsbrae (also known as 79°North Glacier, 79NG) and Zachariae Isstrøm (ZI), a general increase in ice flow velocity has been observed in recent years [26]. For 79NG it has been shown that ice flow velocities are also varying seasonally and can locally be increased by the drainage of SGLs [27]. Therefore, our study demonstrates not only the potential of PolSAR to monitor supraglacial hydrology; it also provides crucial information to better understand the processes in this key region of the GrIS.
The following section about data and methods gives details about the Sentinel data used, our newly developed classification scheme and the definition of SGLs, based on this classification. The results section starts with some examples and insights into the outcome of our classification, followed by two different kinds of validation. After that, we use the water class to define outlines of SGLs in this region and analyze several time series to detect lake drainage events. The discussion puts these results into context and discusses the strengths and remaining issues of our method. The final conclusion summarizes this manuscript.

Sentinel-1 Backscatter
We used the backscattered power of the Sentinel-1 (S1) SAR satellites to distinguish between different surface classes. The combination of the 12-day repeat cycle of Sentinel-1A (launched in April 2014) and Sentinel-1B (launched in April 2016) provides a constellation with a 6-day repeat cycle [28]. In this study we used repeat cycle 74 (see Figure 1), which observes the margin of the NEGIS in Interferometic Wide (IW) swath mode with a total swath width of 250 km and a nominal resolution of 5 × 20 m [28]. Since 21 May 2017, both satellites have provided observations in dual polarization mode for this region.
We processed the Single Look Complex (SLC) S1 observations using the GAMMA Software [29] by calibrating for instrumental noise, mosaicing the burst IW SLC data, multilooking the images to reduce the speckle and georeferencing the scenes to polar stereographic coordinates (EPSG:3413), resulting in a ground pixel resolution of 50 × 50 m. In order to correct for topographic effects, we transferred the radar backscatter normalization convention from the application of the ground area (σ 0 ) to a normalization by the area perpendicular to the the line-of-sight (γ 0 ) [30] using the ArcticDEM [31].

Classification
In contrast to the supraglacial lake classification from optical images, a single threshold is not suitable for SAR data [14,20]. Many different dielectric and geometric surface and subsurface parameters have an impact on the backscattered power γ 0 . Therefore, we developed our PolSAR-based Bayesian classification scheme (henceforth called PBC) using the additional information provided by polarimetric SAR. γ 0 HH and γ 0 HV of each pixel in a SAR scene (see Figure 2 left) could be interpreted as a 2D-domain. In this domain, pixels of a specific surface type are lying in a specific region in the HH-HV plane. This allows a significant reduction of ambiguity, compared to the 1D domain of HH alone. As the pixel values of HH and HV still showed a high correlation, we reduce the second dimension to the "independent" additional information by subtracting both polarizations (HH-HV). However, we still found some ambiguity between different surface types in this 2D domain, and hence, introduced an additional geometric dimension. Parameters such as snow wetness, density and grain size depend on large scale topographic and climatic conditions, changing rather smoothly over scales of 10 to 100 km. In contrast, SGLs form at local depressions, which have a well defined size of hardly more than 2 km in diameter [13]. In order to quantify the relationship of a pixel value to its neighborhood, we developed a geometric "anomaly index." Therefore, we first appl a 25 km floating median filter, which is significantly larger than the expected SGLs ( γ 0 (x, y) = MED{γ 0 (x−r... x+r, y−r... y+r)} with r = 12.5 km). In order to avoid aliasing effects from non-glaciated areas, we limit this filter to glaciated areas only by an ice sheet mask [32]. Then, we calculate the absolute anomaly A abs of a pixel value γ 0 x,y as its difference with respect to the filtered image (i.e., the high-pass filtered version of the image): A abs (x, y) = γ 0 (x, y) − γ 0 (x, y). To account for regions with different variabilities within the filter window, we normalize A abs by the median absolute deviation within the filter window ( γ 0 (x, y) = MED{|γ 0 (x−r... x+r, y−r... y+ r) − γ 0 (x, y)|} ) to a relative anomaly A: In order to obtain a combined anomaly index from both polarizations, we use the root sum square according to Examples of this anomaly index are shown under classification dimensions in Figure 2 and in Figure 7.
To identify pixels of liquid water in these data, we developed an algorithm based on the machine learning approach of naïve Bayesian classification (e.g., [33]). According to its value in the dimensions HH, HH-HV and the anomaly, we determine the probability that a pixel belongs to a specific radar reflection type. The approach is referred to as "naïve" as it assumes that all dimensions are independent of each other. As we use the difference HH-HV instead of HV and we assume that our spatial anomaly index is also highly independent of the backscattered power itself, we presume that this precondition is met very well.
The conditional probability for this classification is defined by a selection of representative training areas (polygons) for characteristic radar reflection classes. Fahnestock et al. [21] identified four different facies of the GrIS in SAR data. As the main focus of our classification lies in water bodies, we use a slightly different, simplified version of these facies. Besides a class for liquid water, we merged the dry snow zone and the percolation zone into a class called dry. Additionally, we merged the wet snow zone and the bare ice zone into our class called wet/icy. Furthermore, we found a well pronounced signature of crevassed areas; hence, we added another class called crevassed. Based on visual inspection of the three classification dimensions from S1 and by additionally checking quasi-simultaneous optical imagery from S2, we selected four representative training polygons for dry, five polygons for wet/icy, four crevassed areas and 18 water bodies (i.e., SGLs, some clearly visible in optical images, some buried under the surface and only visible by radar; see [14,34]). In order to reflect the variability of the backscatter during a seasonal cycle, we analyzed the time series of S1 and S2 during 2018 and defined a valid epoch range for each polygon. This means that we manually checked each epoch for each polygon in the optical and SAR data and used this epoch for training only if we could clearly see the respective surface type. In consequence, to account for the large extent of wet snow during the peak of the melt season, we also defined one polygon of wet/icy at higher elevations, which is only valid for a few scenes with extreme melt. Due to the lack of optical data in winter to cross check, we limited the epoch range for training to March-September. However, the SAR time series shows that the very early and very late epochs are also representative for the rest of the year.
With the help of these training polygons, we defined a characteristic spectral signature of each class over all three dimensions (Figure 2 top right). Separately for each class, we located each pixel of the training polygons in the 3D classification space, obtaining a class pixel point cloud. This point cloud was then used to obtain a conditional likelihood of this class, which was later used for classification. Therefore, we created a discrete 3D raster with the same dimensions as the point cloud of the training data and a spacing of 0.5 dB in the backscatter dimensions and 1.0 for the anomaly. We binned the point cloud to this raster by setting raster cells which contained data of this class to value of 1; cells without training data were set to 0. This binning discretisized the point cloud of each class to binary values (0/1), keeping the spectral signature but ignoring the data density in each cell. Any 3D bin which contained any data of this class had a conditional probability for the class of 100%; all bins without any data points had a probability of 0%. Then we applied a highpass boxcar filter with a diameter of 5 bins in each dimension to this binary binned space. This transformation from binary to fuzzy probabilities accounts for some fuzziness in the real spatial signature and also reduces the influence of noise. This resulted in a significantly smoother fuzzy probability space (as displayed on the sidewalls of Figure 2 top right). A bin with all similar bins in a 2 bin neighborhood remained at a probability of 100%, if only 50% of the cells in the filter window contained data, the fuzzy space probability was 50% and an isolated classified bin (which is quite likely an outlier) obtained a probability of 1/5 3 = 0.8% in the fuzzy space. We created such a fuzzy probability space for each class. These grids define the characteristic spectral signature of each surface type and were used as conditional probability in the Bayesian classification.
Then, we classified all the SAR scenes using the Bayesian probability grids. According to the values of a pixel in the three classification dimensions, we obtained a probability for each class from the respective fuzzy probability space (i.e., the spectral signature of this class) and attributed this pixel to the class with the highest probability. However, if no class had a probability of more than 5%, we marked this pixel as unclassified. Furthermore, if the most likely class and the second most likely class were closer than 5%, we also marked this pixel as unclassified, as the classification was ambiguous. Using this approach we were able to classify~85% of the grounded ice pixels in each scene during the winter season. During the melt season, this fraction decreased up to a minimum of 45% of the pixels as the separation between water, wet/icy and crevassed was often ambiguous.

Supraglacial Lake Definition
Based on the pixel classification results, we extracted the lake outlines automatically. In order to distinguish lakes from thin surface melt water layers or classification outliers, we used the following criteria: On a grounded ice sheet, SGLs form in local depressions, which are controlled by bedrock topography [1]. Hence, we assume that the location of the lakes is stable over the years. Furthermore, optical imagery shows that the vast majority of SGLs exist for at least a few weeks. Radar profiles [34] give evidence that the water in several lakes even stays liquid after the end of the melt season. In contrast, water covered surfaces after melt events drain or freeze relatively fast.
In consequence, we considered lakes as stationary and marked a pixel as lake if it was classified as water in a sufficient number of scenes. We used a threshold of 13 of the 159 scenes, which is about 1 /12 of the scenes covering 3 years. This means that if a pixel was classified as water for at least one month per year (on average), it was assumed to be a lake. We derived the outline polygons of each connected set of marked pixels and denoted them as individual lakes. This enabled us to analyze the time series of water area for each lake separately.

Complementary Optical Data
To validate our PBC results, we created a similar time series of grids from optical imagery (S2) and applied a very basic classification scheme for SGLs. We created 50 × 50 m resolution mosaics for the whole area of interest. Due to the different orbit period of S2, we created a quasi-simultaneous image for each epoch of S1 by composing the S2 scenes with the least cloud cover within ±1 day (see some examples in Figures 1 and 3). From the red and blue band of these images we calculated the Normalized Difference Water Index for ice (NDW I ice [16], henceforth just called NDW I here) according to NDW I = (blue − red)/(blue + red) and classified water using a threshold of 0.25, (as, e.g., [15,18]). This threshold is able to identify very shallow water also. However, visual inspection showed that sometimes other dark objects such as cloud shadows were also classified as water. Therefore, we applied an additional cloud shadow filter using the band difference thresholds according to Moussavi et al. [35] (green − red > 0.09).   Figure 1) in S1 polarimetric SAR (HH and HH-HV), our classification result and corresponding S2 optical images. Figure 3 shows a sequence of six examples of PolSAR data from S1 and our automated classification results. A comparison with the quasi-simultaneous cloud-free optical S2 images revealed that large parts of many SGLs are still lake-ice covered, even in the peak of the melt season. At the end of the melt season, a lid of lake-ice forms and the first snowfall covers the lakes, so that they cannot be identified in optical data any longer. However, the radar signals penetrate the ice cover and show that the liquid water is still there. This demonstrates that a large fraction of the supraglacial hydrological system is not visible for optical data.

Time Series
On the other hand, this sequence of examples in Figure 3 also demonstrates the limitations of polarimetric SAR in monitoring SGLs. During the peak of the melt season, large amounts of liquid meltwater are stored in the porous firn matrix (up to saturation). This increases the dielectric contrast everywhere on the surface and hampers the separation between the lakes and their surroundings by SAR. Our PBC scheme is trained to separate different surface zones but during the peak of the melt season, water, wet/icy and crevassed in the lower lying parts of our region of interest cannot be separated anymore by PolSAR.
A time series of the area, classified as supraglacial water, is shown in Figure 4a. In order to allow a glaciological interpretation of this time series, Figure 4 refers to the grounded parts of the drainage basins of 79NG and ZI only. Compared to the full extent of the classified S1 observations (see Figure 1), this excludes significant parts of the ice sheet margin north and south of the large ice streams. The outline of these two drainage basins was defined by Krieger et al. [36]. Furthermore, we limit these basins to the grounded part of the ice sheet according to the outline polygons from IMBIE2 [37,38]. Figure 4a shows that within these two basins we observe a relatively stable area of liquid water of~250 km 2 persistently over the three years. As explained above, very wet surface conditions during the melt season cause the classification to fail, which appears as mimima in the time series where maxima would be expected.
An additional binning of the time series into classes of 100 m elevation steps (Figure 4b) shows that the majority of the water is located in the range of 700-1200 m above sea level in this region. Compared to the total area of the respective zones, the water covers about 2-3% of the surface area between 500 m and 1400 m. Even if more water would typically be expected at lower altitudes, it is not surprising that only very few lakes were detected there, as the lowest 500 m of altitude are relatively steep and crevassed in this region. Elevation bin medians of the backscattered power and the difference HH-HV clearly show the melt periods. In regions with an overall low backscatter in HH and a high difference in HH-HV, different surface types cannot be separated any longer by PolSAR. Similarly binned median ERA-5 reanalysis data [39] confirm temperatures above melting during these periods.
In Figure 4a, we also compare the results from PBC to classification results from quasisimultaneous S2 optical imagery (see Section 2.4). The resulting area from the image pixels, classified as water from optical data, is shown by the black bars. In the PolSAR time series we see that the classification of the subsequent scenes provides very consistent results, except during the peak of the melt season, where the PBC gets less reliable. In contrast, this is the only time of the year where the optical data observes any significant amount of water at all. Nevertheless, this comparison shows that the optically visible water is only a small fraction of the total water area.
In order to quantify these differences, we compare the classified water area right after the temperatures dropped below zero again in late summer. On 19 August 2017 we identify a total water area of 465 km 2 at 79NG and ZI from PolSAR but only 60 km 2 from optical data. At the end of the relatively cold summer of 2018 (26 August) we observe 303 km 2 of water by PolSAR and 36 km 2 from S2. The summer of 2019 caused extreme melting. On 21 August 2019 our approach classified 1088 km 2 of liquid water with a significant fraction at higher elevations of up to 1500 m and more. However, many of these lake-like features already disappeared in the next scene six days later and several more disappeared after 12 days (two scenes later), leaving a classified area of 520 km 2 on 2 September. Hence, we assume that many of the classified features were very shallow accumulations of meltwater, refreezing within a few days up to a week, which we do not consider as real, significantly deep supraglacial lakes. During this peak at 21 August 2019 we observed 166 km 2 of supraglacial water from optical imagery, which decreases to 129 km 2 on 27 August and only 32 km 2 of optically visible water on 8 September 2019. In contrast, PolSAR still classified 432 km 2 as water on 8 September 2019. In conclusion, this comparison shows that even during the melt season, the water area detected by PolSAR is about 6-8 times larger than what is visible in optical data. Only a few weeks after each melt season, no water is visible in optical S2 data anymore, but the SAR data shows that the lakes are still there. Shallower lakes freeze through during the winter, as can be seen at the lower part of the region in Figure 3. When all the water in the lake is frozen, the dielectric contrast at the water/lake-ice-interface at the top of the water body vanishes. This causes a significant reduction in surface scattering, and hence, the HH-HV difference can no longer distinguish between the frozen lake water and the surrounding glacier ice. However, lakes with a water depth from several meters up to tens of meters (which is the case for many SGLs; see [17,27,[40][41][42]) have a heat capacity which allows them to stay liquid over the whole winter if they are well isolated by a significant layer of ice and snow [19].

Validation
In order to validate the results, a comparison between an optical and a PolSAR-based classification can be done. However, during most of the year, the water is buried under a snow/ice layer and optical data cannot observe any water. During the peak of the melt season the optical data can detect liquid water on the surface, but during this time, the results of our PolSAR-based method is quite unreliable. Furthermore, a significant amount of liquid water is still buried which also cannot be detected by optical data. Therefore, we used two different approaches to validate our results.
In the first validation we used the same polygons as we used to train the classification scheme, but we checked how the algorithm performed at different epochs. The training polygons cover the period March to September 2018; hence, we chose three epochs of 2017 and three epochs of 2019 for validation.
For each year we selected one epoch in spring, one in the peak of the melt season in summer and one epoch in autumn. As in Section 2.2, each polygon was verified by visual inspection of S1 and S2. If any of the polygons did not apply for this period, it was slightly modified or withdrawn. Then, we compared the total area of the test polygons of each class to the classification results, which is displayed in Figure 5. This validation shows that 75-90% of the water polygons were classified correctly, except during the melt season, where the classification agrees with only~50%. Besides the classification of water, we note that dry shows a nearly 100% match with the test regions, except in August 2019, where we hardly see any dry snow in this region at all. In contrast, the polygons for wet/icy and crevassed yield a relatively large fraction of unclassified and even misclassified cells. This indicates that the separation between these two classes (and water as well) is often difficult. Better classification results for crevasses might be obtained by another spatial parameter as large crevasses often show typical stripe patterns. However, areas with smaller crevasses do not show this pattern as the pixel size of 50 × 50 m is too large to resolve these structures. Nevertheless, our classification is focused on classifying water. The other classes are just used to distinguish water from other surface types, which works quite reliable, except during the melt season. In a second validation we compare our results to snow-penetrating radar profiles from an airborne mission (Figure 6). During a campaign in April 2018 an Ultra Wide Band Microwave (UWBM) snow radar system [43] with a frequency range of 2-18 GHz was operated from AWI's Basler BT-67 aircraft in this region. We compared the UWBM-profiles, observed between 14 April and 18 April 2018 with PolSAR from 16 April 2018, as shown in Figure 6. Reflections from SGLs were picked manually from the radar profiles. We checked for an agreement between the airborne detection and our classification results within a maximum distance of 200 m to allow for some uncertainties of the positioning and the orientation of the airborne radar, the UWBM footprint size, uncertainties in the geolocation of the SAR images or observational noise, which could cause gaps in the classification. For 90.4% of these lake picks our PolSAR-based classification detected water as well. Furthermore, the unmatched UWBM picks (cyan in Figure 6) are all located downstream of supraglacial lakes. We interpret these features as remainders of former lake ice, which followed the ice flow but still shows a significant contrast to the layers below. Some examples for the agreement of both data sets are shown in detail in Figure 6. In HH-HV, lake A shows a strong contrast against its surroundings in April 2018. Comparing with quasi-simultaneous optical imagery shows that everything is covered by ice and snow at this time and no water is visible. When comparing with S2 scenes from the previous summer, the very bright part in the upslope direction was clearly visible as open water, while the darker part of the lake in the downslope direction was ice covered during the whole summer season of 2017 as well. The two features B and C appear significantly more fuzzy against their surroundings in HH-HV, making their classification as water from PolSAR more questionable. Comparing with optical imagery from the previous summer, the two features were hardly identifiable as water. Only with the help of radar signals can these buried water bodies be identified and the UWBM profiles confirm the classification results of PBC.  Figure 3) and the corresponding radar profiles (d).
The white arrows in (b) denote the flight direction, i.e., the orientation of the radar profiles from left to right.

Lake Drainage
In order to analyze the state of individual SGLs, we have to define lake outline polygons. As described in Section 2.3, we use the criterion of a persistent classification as water over a sufficient number of scenes to distinguish lakes from short term water covered surfaces or classification outliers. We join such connected cells and if the area of these cells exceeds 0.1 km 2 we define them as an individual lake of significant size. According to this definition we find a total number of 1055 lakes in the area covered by the PolSAR observations between 2017 and 2019. 433 of these lakes are located within the basins of 79NG and ZI.
Based on these lake outline definitions, we can analyze the water area of each lake over time, and hence, detect rapid drainage events. Therefore, we create lake area time series by summing up the classified water pixels of each epoch. In order to reduce the noise, we smooth these series using a floating median over three epochs (covering a time interval of 12 days). Now, we scan these time series for rapid lake drainages, which we define as a rapid change in lake area water fraction (compared to the total area within the outline) from more than 30% to less than 10%. Previous studies with optical data used significantly higher thresholds for filled lakes (e.g., 80 % in [15]). However, several lakes vary in area significantly between different years, so that in cold years they do not reach levels as 80% of the maximum lake area at all. As we also observed drainages of such partly filled lakes we chose a lower threshold to detect drainage events. All such rapid declines of water area are displayed as circles in Figure 8. The vast majority of such events are located at the very beginning of the melt season. To illustrate this temporal coincidence, we plotted the average backscatter HH over the drainage basin area in the background of Figure 8 (as shown before in Figure 3d). After the onset of melt or even rain events, the wet ice sheet surface looks very similar to the lakes themselves in the radar images, and as a consequence, the water within the lakes can no longer be detected reliably from PolSAR, causing these apparent drainages. Visual inspection of simultaneous optical imagery prove that these lakes are not draining. Therefore, another criterion is necessary to distinguish between real and false drainages.
To justify that a disappearance of classified water is really related to a lake drainage, we calculate similar time series for the average HH and HH-HV backscatter within the lake outline (as in [44,45], but we use both polarizations). As it can be seen in Figure 7a, an open water lake surface (as observed during the summer) causes a specular reflection, which means that the side-looking radar system hardly receives any return power at all. After a lake drainage, the signal is reflected from the former lake bottom, which is significantly rougher, resulting in an increase in HH. Hence, a drainage of an open water lake is accompanied by an increase in HH. During the winter time (Figure 7b), the HH-channel is less suitable to detect drainages because the signal has to penetrate a layer of ice and snow first to reach the buried lake. Hence, the HH-channel is already strongly influenced by the interactions of the radar signal with the top layers and the presence or absence of a flat water table is less prominent in HH. In contrast, the HV-channel is significantly more affected by a winter time drainage. After the strong dielectric contrast of the water body is gone, the signal can penetrate even deeper below the former lake bottom, increasing the HV-reflection significantly more, compared to HH. Hence, the relatively high HH-HV difference of a buried lake decreases rapidly and both channels become almost balanced after a lake drainage. Both kinds of rapid changes in backscatter (HH increase and HH-HV decrease), however, can also have other sources, e.g., larger scale weather effects. To justify that these changes affect the lake only, but not its surroundings, we also check the absolute anomaly A abs (see Equation (1)) for such significant jumps.
In consequence, we add two alternative criteria using the radar reflection time series to exclude outliers in the drainage event detection from the lake area time series alone. If a significant increase in HH and A abs HH of more than 4 dB occurs simultaneously with the vanishing water area, the drainage is confirmed as a typical summer time drainage. If, in contrast, the drainage event is accompanied by a drop in the difference HH-HV and A abs HH−HV of more than −2 dB, the event is confirmed as a typical winter time drainage.
All drainage events, detected from these time series, are displayed in Figure 8. The majority of apparent drainage events occur at the onset of the melt season, which results in an overall increase in HH-HV, displayed in the background. The backscatter time series within the lakes, however, show that these events are false drainages, caused by a disrupted classification due to surface melting. Based on the above mentioned criteria, we were able to separate between real and false drainages. Even during the summer season, where melt hampers the classification, some drainages are identified in the SAR data and confirmed by optical images. We find that drainages occur during the whole year, not only in summer. We observe drainages of a significant area in autumn, many weeks after the end of the

Discussion
Our PolSAR-based Bayesian classification is based on the characteristic spectral signature of different types of glacier surface (dry, wet/icy, crevassed and water). Therefore, the selected set of representative training areas is crucial for successful results. A large number of training examples helps to cover the full range of surface characteristics within one class. We used training data which covers all seasons of 2018. For the water class, we chose some open water meltwater lakes, which were clearly visible in optical data during the peak of the melt season, but also some regions which were not so obvious from optical data. Some of theses regions are ice covered parts of only partially open lakes, others do not show any open water during the summer of 2018. Nevertheless, the HH-HV PolSAR signature, the UWBM radar and also the unusually flat surface, visible in optical data, indicate buried SGLs, so we selected some of them as training areas as well. While it is important to cover the full range of surface characteristics for one class, it is also important to visually check each epoch of each training polygon in all available datasets (HH, HH-HV, optical) to assure that the algorithm is trained correctly.
We did explicitly not use any additional information for classification. While information about, e.g., the day of year or the surface temperature could help to interpret a pixel's spectral signature correctly, we aim at doing this classification from PolSAR data alone. Especially during the summer melt onset or rain events, our algorithm has a large problem in distinguishing the classes water, wet/icy and crevassed. However, this is not a problem with the algorithm but of the SAR data itself. If the whole surface gets covered by a meltwater layer, volume scattering decreases dramatically everywhere and the remaining spatial variation in HH is more dominated by parameters which are not directly related to our classes, as e.g., small scale surface roughness. One could argue that a combination of SAR and optical data would help to mitigate this problem. However, we see that most of the lakes are at least partly ice-covered, even during the melt season. Therefore, water classification based on a combination of optical and SAR imagery needs some additional assumptions on the persistence or disappearance of supraglacial water and cannot be just performed on a per-pixel basis.
Our results show that a significant amount of water is stored on the ice sheet surface at 79NG and ZI. While optical data would indicate the presence of water only during a few weeks in summer, airborne radar (e.g., [34] or the profiles presented here) has shown that the lakes persist also during winter. Beyond such campaign-style information, SAR imagery is able to observe surface and subsurface characteristics with a significantly better spatio-temporal coverage [20,44,46]. Our method allows to use this SAR data to map supraglacial water in high resolution in space and time automatically even in winter. Our data show that the water covered area remains relatively constant over the 3 years of observation. This includes the relatively cold year 2018 but also the extreme summer melt year 2019. At the end of the melt season, some of the shallower water bodies freeze through, but after October the area of water remains relatively constant.
After identifying SGLs from our classification results, we created and analyzed time series for each lake. We identified several drainage events, distributed all over the year, not only in summer. Many of these events do not appear isolated but cascade simultaneously over several lakes within a local vicinity. This indicates that the events that are responsible for a drainage often have a larger scale effect (as e.g., the opening of a long crevasse) or that one drainage triggers other drainages in the surroundings as well [47,48].
Since 2017 Sentinel-1 observed the region of 79NG and ZI regularly in dual-pol mode, which is crucial especially for the classification of water during the wintertime. After 3 years of data, we have observations during very different characteristic years. According to Ignéczi et al. [25] this region will host the largest spread of SGLs over the next decades, hence, observations of the current state are essential to quantify the effect of these future changes. Therefore, consistent time series over several years are necessary to distinguish between an extreme year and a long term trend.
Beyond this specific region, many previous studies focused on SGLs on the western part of the ice sheet in the vicinity of Jakobshavn Isbrae [8][9][10]15,18]. Johansson and Brown [20], Miles et al. [14] and Benedek and Willis [44] demonstrated that also SAR imagery is applicable there. In contrast to the manual mapping of SGLs, our PBC algorithm, once trained adequately, is able to map and monitor many SGLs over large regions automatically. The margin of the GrIS is covered practically completely by S1 in dual polarization mode. Some first tests with other orbits show that this classification, based on the spectral signature, provides comparable results also at different orbits. However, a Greenland-wide inventory would require further investigations towards the adaptation of this algorithm on an ice sheet scale, which is beyond the scope of this paper.

Conclusions
Previous studies used optical satellite images to map and monitor SGLs on the GrIS. Optically derived results suffer from the fact that liquid water is often not visible due to snow cover, missing daylight or bad weather conditions. Some of these studies used SAR to overcome these limitations [14,20,44]. However, using only one polarization, they were not able to map lake outlines automatically but used results from optical classifications or manually picked outlines instead.
In this paper we present an algorithm which is solely based on PolSAR data and is able to map and monitor SGLs automatically. We demonstrated that this is especially important, as manual outline detection is not suitable on a larger scale and outlines from optical data significantly underestimate the real extent of SGLs. At the end of the melt season, when we can expect reliable results from both techniques, the area of detected supraglacial water from SAR is about 6-8 times larger than the water area from optical data. A more detailed look revealed that many lakes, which are clearly visible in PolSAR and also in airborne radar profiles, do not show an open water surface in optical images. These lakes might have formed during an extreme melt year but are later covered by an ice lid all over the year and hence are not detectable in optical imagery. Other lakes show only a partly persistent ice cover. Here, the ice lid moves on with the ice flow, which results in an open upstream part and a buried downstream part of these lakes.
This demonstrates that satellite-based radar imagery is crucial to study the full extent of the supraglacial hydrological system of Greenland. The additional information provided by polarimetric SAR allows to separate surface from volume scattering, which is fundamental to automatically mapping supraglacial water. With an expected increasing role of melt water on the GrIS in a changing climate, this information provides an important data basis for future time series of changes.