Automatic Supraglacial Lake Extraction in Greenland Using Sentinel-1 SAR Images and Attention-Based U-Net

Abstract

With global warming, supraglacial lakes play an important role in ice sheet stability and climate change. They are not only the main factors affecting mass balance and sea-level rise but also the key units of surface runoff storage and mass loss. To automatically map the spatiotemporal distribution of supraglacial lakes in Greenland, this paper proposes an attention-based U-Net model with Sentinel-1 SAR imagery. The extraction results show that compared with the traditional network, this method obtains a higher validation coefficient, with an F1 score of 0.971, and it is spatiotemporally transferable, able to realize the extraction of supraglacial lakes in complex areas without ignoring small lakes. In addition, we conducted a case study in the Jakobshavn region and found that the supraglacial lake area peaked in advance between spring and summer due to extreme melting events from 2017 to 2021. Meanwhile, the supraglacial lakes near the 79°N Glacier tended to expand inland during the melting season.

1. Introduction

The Greenland ice sheet (GrIS) is the second largest continental ice sheet in the world and contains ice equivalent to 7.4 m of sea-level rise [1]. It has been in a state of accelerated melting, with ice loss reaching 3902 ± 342 billion tons from 1972 to 2018 [2,3]. This trend is mainly due to the dynamic thinning of the glaciers and the decrease in surface mass balance, which is dominated by increased melting and surface runoff [4]. With the effect of solar shortwave radiation and precipitation, when firn and ice reach the melting point, the surface runoff begins to be transported in the form of supraglacial streams, which will expand to supraglacial lakes over the low-lying glacial surface. Part of the meltwater flows into the ice sheet along cracks and moulins to form an englacial lake and flows into the bed to form a subglacial lake [5]. During this process, as the main unit of surface water storage and a key factor of mass loss, the supraglacial lake affects the mass balance and dynamics of the GrIS in the following ways. First, compared with the surrounding features, the supraglacial lake has lower albedo, which means it can absorb more energy from the sun, further promoting melting and global warming. Second, supraglacial lake drainage events can increase the flow of ice by hydrofracture. When the water reaches the bed, it will promote the sliding of the base, thus speeding up the movement of the ice sheet. Moreover, infilling and drainage events of the supraglacial lake will lead to bending of the ice, which will trigger ice shelf collapse and cause the glacier to retreat [6,7]. Therefore, to quantify the impact of the supraglacial lake on the mass balance, hydrology, dynamics, and stability of the GrIS, it is necessary to map the continent-wide supraglacial lake system of Greenland and obtain its spatiotemporal variation.

In general, most previous studies on supraglacial lakes evaluated the GrIS on a regional scale, such as by focusing on southwest Greenland, and there is a lack of research on the whole area. The existing supraglacial lake detection methods mainly include field observation, optical remote sensing, and synthetic aperture radar (SAR) remote sensing [8]. So far, field observation is the oldest and most widely used method to detect supraglacial lakes. According to our statistics, more than 100 research stations have been set up in and around the ice sheet, contributing to abundant records on supraglacial lakes. In 2013, Doyle et al. observed the supraglacial lake west of the GrIS and obtained its evolutionary trend [9]. Two years later, Stevens et al. used GPS stations to follow short-term rapid drainage events in Greenland’s supraglacial lakes and discovered a series of drainage networks [6]. Later, in 2019, Chudley et al. described the drainage events of supraglacial lakes around marine-terminating glaciers in the GrIS based on field observation data, suggesting that rapid drainage significantly impacts ice sheet dynamics [10].

Although field observation is the basis for the study of supraglacial lakes, there are still many limitations due to the unique location and environment of polar regions. For example, traditional methods require on-site manual operations, but there are large areas covered by ice and snow that are difficult to reach, so only remote-sensing methods, such as aerial photogrammetry and satellite observation, can be used. Recently, a series of newly launched remote-sensing satellites, including Landsat, Terra, Sentinel, Worldview, Gaofen, etc., have provided massive data for large-scale and systematic research on supraglacial lakes. Currently, the lake extraction method is mainly focused on manual digitization and automatic extraction [11]. Manual digitization, based on visual interpretation, can make full use of prior knowledge and has high detection accuracy, but it is inefficient and easily affected by an expert’s knowledge base. Automatic extraction methods include band threshold segmentation, normalized difference water index (NDWI), supervised and non-supervised learning, and many other methods [12,13,14,15]. As optical images are simple to interpret and contain multiple types of spectral information, they are widely used. However, optical sensors cannot detect polar night conditions and are easily disturbed by clouds, rain, fog, and snow. Therefore, they cannot provide objective and long-term observations of supraglacial lakes [16,17]. In contrast, SAR, an advanced microwave earth observation technology, has the advantages of all-day, all-weather operation, large swath width, and spatiotemporal resolution, making it an important way to observe ice sheets. In particular, the Sentinel-1 SAR satellite has accumulated large-scale, long-period, high-resolution radar observation data since its launch in 2014, providing plenty of images for supraglacial lake research.

At present, the automatic methods of detecting supraglacial lakes based on SAR mainly include supervised and unsupervised learning, such as random forest, deep learning, Bayesian classification, etc. In particular, the method of deep learning combined with SAR image backscatter coefficients to detect lakes has been developed rapidly. Several studies have used Sentinel-1, Gaofen-3, and other images to carry out the extraction research of supraglacial lakes [18,19,20,21], but most focused only on large lakes and ignored small bodies of water. Compared with traditional methods, such as Canny edge detection and texture monitoring, deep learning is more suitable for batch processing of large amounts of data and high-precision adaptive detection of supraglacial lakes. Because microwaves can penetrate surface features, SAR is able to observe covered supraglacial lakes that optical sensors cannot see. As shown in Figure 1a–d, this refers to supraglacial lakes covered by snow or thin ice (as opposed to those that are deeply buried). The existing studies on covered supraglacial lakes are mainly based on semi-automatic extraction, such as manual digitization, threshold segmentation, polarization classification, etc., which cannot automatically extract a large area [22,23,24,25]. At the same time, due to the complexity of the SAR imaging mechanism and the diversity of surface features on the ice sheets, there are still many problems in supraglacial lake detection that need to be resolved. For example, the backscattering coefficient of SAR is very sensitive to dielectric properties, such as surface humidity, electrical conductivity, geometric properties, surface roughness, particle size distribution, and the internal structure of ice. As seen in Figure 1, the shapes of supraglacial lakes vary greatly (Figure 1e), with some having blurred edges (Figure 1f) and rough surfaces (Figure 1g). Surface features with similar backscattering coefficients to lakes, such as dark ice (Figure 1h), blue ice (Figure 1i), shadows (Figure 1j), ice crevasses (Figure 1k), and wet snow (Figure 1l), will interfere with the detection of supraglacial lakes. Moreover, different radar bands, polarization forms, and incident angles will result in different tones and textures of supraglacial lakes on the images, which makes it difficult to identify them.

Based on Sentinel-1 SAR data, we propose an attention-based U-Net deep-learning framework for Greenland supraglacial lake extraction to overcome most of the described problems. Compared with the traditional U-Net, this method has higher accuracy and can effectively distinguish other surface features similar to lakes. In addition, complete supraglacial lake mapping can be achieved in areas where extraction is difficult, such as slush, blue ice, dark ice, moraine-covered ice sheets, etc. We also performed the extraction and multi-segmentation of lakes covered by snow or thin ice in winter. Subsequently, to verify the method’s spatiotemporal transferability, we used other independent data for testing. The results show that the proposed method achieves higher validation coefficients than other U-Net models, with an average F1 score of 0.971. Furthermore, two specific regions, the Jakobshavn and Nioghalvfjerdsfjorden (79°N) Glaciers, were selected for long-term series analysis. Based on the final mapping results of the extent of the supraglacial lakes, we discuss the seasonal changes of the lakes and the phenomenon of their advancement inland.

2. Materials and Methods

2.1. Study Sites

The GrIS is the second largest ice body in the world, next to the Antarctic, with a maximum length of 2400 km, a maximum width of 1100 km, and an ice area of 1,710,000 km². The study sites used for training and testing the proposed supraglacial lake detection algorithm are evenly distributed across the GrIS (Figure 2).

We selected 10 typical regions as training sites and 6 typical regions as test sites. To make the method universal and convincing, we selected the typical sites based on the following principles. First, the sites had to include all surface features on the ice sheet (e.g., rivers, crevasses, bare ice, wet snow, dry snow, blue ice, dark ice, moulins, bare rock, moraine-covered snow, etc.). Second, the acquisition time of Sentinel-1 images should span different observation times, from the summer melting period to the freezing winter days. Based on these conditions, we selected the training and test sites according to the distinctiveness of hydrological basins. As can be seen in Figure 2, they are mainly divided into north (NO), northeast (NE), southeast (SE), southwest (SW), and central west (CW) [26]. The NE region is mainly around the 79°N Glacier, the fastest outlet glacier in the NE region and the largest remaining ice shelf in the Arctic [27]. Changes in the volume and distribution of supraglacial lakes upstream of the 79°N Glacier have a demonstrable effect on the downstream ice velocity [28]. In the SE region, the tidewater glacier was mainly selected, which is fed by enough snow to flow out of the mountains and down to the sea, including the Helheim Glacier, Heimdal Glacier, Koge Bugt, etc. We selected several land-terminating glaciers in the SW region, including Kangiata Sermia, Russell Glacier, Isunnguata Sermia, etc. In fact, some land-terminating glaciers of the GrIS flow more quickly in the summer, when supraglacial lakes transport the meltwater to subglacial drainage systems [29]. In the CW region, the Jakobshavn Glacier was selected, which is a large outlet glacier; the surrounding area is known as the Jakobshavn Region (JAR), which is an ideal site for studying the distribution changes and storage of supraglacial lakes in fast-flowing glaciers [30]. The NW region includes the Nansen, Sverdrup, and Hayes Glaciers. The NO region mainly covers the Petermann Glacier, one of only seven long-term floating glaciers in Greenland. Large-scale calving events often occur here, which can be helpful in studying the influence of supraglacial lakes on the stability of the ice shelves [31].

2.2. Datasets

2.2.1. Sentinel-1 SAR Data

To achieve continent-wide mapping of Greenland’s supraglacial lakes, we used Sentinel-1 data acquired in interferometric wide swath (IW) mode with dual polarization (HH + HV) provided in the Level-1 Ground Range Detected (GRD) product from ASF Data Search. Sentinel-1 transmits signals with horizontal polarization and simultaneously receives echo signals with horizontal polarization (HH) and vertical polarization (HV) [32]. According to the characteristics of different polarization modes, we chose HH polarization to extract supraglacial lakes in the melting season and added HV polarization in the freezing period, which is explained in detail in Section 2.3.1. Regarding range resolution, the IW mode of S1 has a pixel-level resolution of 10 m, which can achieve more accurate lake monitoring than other modes. In addition, a single Sentinel-1 satellite has a time interval of 12 days, which can be reduced to less than 6 days with the satellite network of Sentinel-1A and Sentinel-1B. Hence, we mapped the supraglacial lake distribution in less than 6 days.

Table 1. Sentinel-1 SAR image information and corresponding basin distribution for selected training and testing datasets.

ID	Acquisition Date	Path	Study Region	Flight Direction	Type
S1-01	10 August 2020	207	NE	Descending	Training/Test
S1-02	27 August 2021	141	SE	Descending	Training
S1-03	18 October 2021	112	SE	Descending	Test
S1-04	29 September 2021	10	SE	Descending	Training
S1-05	7 July 2017	54	SW	Descending	Training/Test
S1-06	12 February 2021	90	SW	Ascending	Training/Test
S1-07	17 July 2020	90	SW	Ascending	Training
S1-08	21 February 2019	54	CW	Descending	Test
S1-09	7 August 2020	127	CW	Descending	Training
S1-10	28 August 2020	90	CW	Ascending	Training
S1-11	11 February 2019	90	NW	Descending	Training
S1-12	2 January 2019	26	NO	Descending	Test
S1-13	27 September 2021	74	NO	Ascending	Training

shows the specific information of Sentinel-1 images selected for training and test datasets in this paper. S1-02, S1-04, S1-07, S1-09, S1-10, S1-11, and S1-13 are training datasets, and S1-03, S1-08, and S1-12 are test datasets. S1-01, S1-05, and S1-06 are part training and part test datasets. All images were acquired in 2017, 2019, 2020, and 2021 from the melt-rich summer in July and August to the freezing period in January and February. Using these images, we created a dataset containing 19,850 patches, which were applied to train the attention-based U-Net model.

2.2.2. Auxiliary Data

• Sentinel-2 data

To verify the difference between SAR and optical image extraction results, we selected Sentinel-2 multispectral imagery (MSI) with a cloud cover of less than 15% during the austral summer (July–August). We used the Level-1C product, which was orthorectified and is a map-projected image containing top-of-atmosphere (TOA) reflectance data from the USGS EarthExplorer. Considering the limitations of optical images, we only used the Sentinel-2 images and binary classification results extracted by Sentinel-1 images for comparison and supplementation.

2.

DEM data

Although most of the Greenland ice sheet is flat, the periphery, especially the ablation zone, is rugged and produces numerous shadows. In the post-processing phase, we used ArcticDEM to reduce the shadow confusion due to excess terrain relief. The ArcticDEM data were constructed from in-track and cross-track high-resolution imagery acquired by optical imaging satellites covering Greenland, such as WorldView-1, WorldView-2, WorldView-3, and GeoEye-1. Here, we selected 10 m resolution DEM products to assist with shadow removal [33].

3.

Coastline data

The areas outside the ice sheet default to non-supraglacial lake areas. In addition, the division of hydrological basins is helpful to determine the regional statistics of supraglacial lakes and understand their development process and their influence on the GrIS. Since supraglacial lakes in coastal Greenland have similar backscattering coefficients to seawater, ice-marginal lakes, and lakes on floating tongues or ice shelves, we integrated the ice sheet boundary and its hydrological basin products from previous studies [26].

2.3. Method

2.3.1. Data Preprocessing

For the preprocessing of Sentinel-1 SAR images, we used the open-source Python package called Snappy. The steps included orbit correction, thermal noise removal, radiometric calibration, filtering, terrain correction, and linear to dB conversion (Figure 3). First, we used orbit correction for image registration and then reduced the noise inherent in the SAR satellite system by removing thermal noise. Radiometric calibration converted the radar’s backscattered signal into a unitary physical quantity, sigma0. Next, coherent speckle noise was filtered by the refined Lee filter. Geocoding and terrain radiation calibration were carried out according to the terrain correction. Finally, to enhance the image contrast, the unit of the backscattering coefficient was converted from linear to dB.

After the above steps, subsections of the preprocessed images were normalized into 0,1. Min–max normalization was used to standardize all data into the same range and accelerate the convergence of the attention-based U-Net model. This was performed by subtracting the minimum of the image and then dividing it by the difference between the maximum and minimum values. Then, we manually labeled the normalized images in the ArcGIS software to create a dataset for training. The labeled images were cropped into patches of 256 × 256 pixels, with an overlap rate of 0.5. Then, 3970 samples were augmented to 49,886 (positive: negative samples = 1:2) through rotation, flipping, brightness stretching, adding Gaussian noise, and manual screening. Then, the samples were divided into training and test datasets at a ratio of 8:2. In addition, for the problem of terrain overlay, we calculated the slope according to ArcticDEM, and removed the area with a slope > 15° to reduce the interference of shadows in the extraction of lakes. We used the Greenland ice sheet boundary product to solve the interference problem of ocean pixels, ice-marginal lakes, and lakes on floating tongues or ice shelves. The result was projected in the WGS 1984 Stereographic North Pole coordinate system to provide 10 m resolution supraglacial lake distribution products.

2.3.2. Supraglacial Lake in SAR

To understand the characteristics of supraglacial lakes in SAR images, we performed statistical analysis of the backscattering coefficient on the preprocessed Sentinel-1 SAR images and curve fitting. As seen in Figure 4, the backscattering coefficients of the supraglacial lake present a Gaussian distribution, mainly distributed in the −30 to −15 dB range. However, compared with the HV polarization (Figure 4b), the difference in backscattering coefficient distribution between supraglacial and non-supraglacial lakes is more obvious in the HH polarization (Figure 4a). Moreover, in the HV polarization, the backscattering coefficient of the covered supraglacial lake is in the range of −25 to −5 dB, and that of the non-supraglacial lake is mainly distributed in the range of −30 to 0 dB. Based on the scattering properties of the ground object, surface scattering mainly occurs in the melting season, when the supraglacial lake is full of water and has a smooth surface in the absence of wind, and most of the energy is reflected away from the radar and generally appears as a small backscattering coefficient on the SAR images. When the temperature decreases, the lake’s surface gradually forms thin ice with abundant sharp edges, which strongly contributes to the volume scattering, so the covered supraglacial lake is brighter on the SAR images, and the backscattering coefficient is increased [26]. Overall, the backscattering properties of lakes vary with the season; lakes in the melt-rich season are more clearly defined in HH polarization, and lakes covered with thin ice are more pronounced in HV polarization. Therefore, when creating the training dataset, we selected the HH polarization images in summer and the HV polarization images in other seasons.

2.3.3. Attention-Based U-Net Architecture

We proposed an attention-based U-Net deep-learning model with a classic encoder–decoder structure to extract supraglacial lakes, shown in Figure 5. The architecture consists of eight residual blocks, four downsampling blocks, four squeeze and excitation (SE) blocks, one atrous spatial pyramid pooling (ASPP) module, and four upsampling blocks. In the network, residual structures replace ordinary units of the traditional U-Net [34]. In contrast, residual structures add fast connections to single forward propagation, allowing deeper networks to be trained without degradation while extracting more supraglacial features. In the residual unit, batch standardization and Leaky-ReLU activation are performed after each convolution. This not only reduces the model’s sensitivity to initialization parameters but also has a regularization effect, to a certain extent. Each downsampling block consists of a 3 × 3 convolution block with BatchNorm (BN), an activation function Leaky-ReLU, and a 2 × 2 pooling layer. Symmetrical to it, the upsampling stage is performed by stacking transposed 2 × 2 and 3 × 3 convolutions. In addition, to extract more information in the feature map after encoder convolution, a lightweight attention module SE is introduced. The SE block is used to adaptively extract image features after each residual unit, which enables the network to focus on the specific segmentation task, thus enhancing detail segmentation [35]. Here, ASPP can extract multi-scale contextual information and multiple effective fields of view [36], ensuring the integrity of lake extraction and reducing discrete blocky segmentation. This is important, considering that the appearance of Greenland’s lakes varies widely. Furthermore, the attention-based U-Net uses skip connection to perform feature mapping between the corresponding coding and decoding blocks to supplement part of the lost information, ensuring that the finally restored feature map incorporates more low-level features, so that multi-scale prediction and deep supervision are possible.

In the process of supraglacial lake training, the advantages of the proposed U-Net are as follows. First, it can realize end-to-end training with a small amount of data, which is suitable for processing SAR data with small sample size due to manual annotation. Second, since the semantic information of supraglacial lakes in SAR images is relatively simple, such as lake and non-lake areas, high-level semantic information and low-level features are very important, making use of skip connection, U-shaped structure, SE block, and ASPP block. This gives the model a strong ability to distinguish between rivers, dark ice, blue ice, wet snow, etc., and lakes, to effectively reduce misidentification and extraction. Finally, compared with other traditional U-Net networks, the proposed model pays more attention to the details and completeness of lake segmentation, so it has higher accuracy and will not ignore small lakes.

To mitigate the resolution degradation caused by multiple instances of downsampling, ASPP is used as the transition layer of the network (Figure 6). The ASPP module can capture image context information of various scales, thus facilitating the inclusion of semantic knowledge of multiple scales in the extracted feature map. The feature map will pass through the ASPP module after the encoder module and then be used as the input for the decoder module. Specifically, the ASPP module consists of one 1 × 1 convolution and four parallel 3 × 3 convolutions with dilation rates of 2, 4, 8, and 12, respectively (Figure 6a). At the same time, an image-level feature is produced by global average pooling. The resulting features from the two branches are concatenated and passed through another 1 × 1 convolution.

The structure of the SE module, as shown in Figure 7, is as follows: squeeze and excitation are both included, W and H represent the width and height of the feature, C represents the number of channels, and the size of the input feature is H × W × C. First, in the squeeze module, the features U_c of each image channel are compressed into a 1 × 1 × C vector through global average pooling. The $z\in R^c$ statistic is generated by shrinking U through its spatial dimensions H × W in Equation (1):

$$z_c=F_{sq}\left(U_c\right)=\frac{1}{W\times H}{\sum }_{i=1}^W{\sum }_{j=1}^H{u}_c\left(i,j\right)$$

(1)

where $i\in W,j\in H$. Then, in the excitation module, to prevent the model from becoming complex and to account for generalization, two fully connected layers are set around the nonlinearity to act as a bottleneck to parameterize the paired gate mechanism. The specific operation is as follows:

$$s=F_{ex}\left(z,W\right)=\sigma \left(g\left(z,W\right)\right)=\sigma \left(W_2\delta \left(W_{1}z\right)\right)$$

(2)

where $\delta$ refers to the ReLU function, and $W_1ϵR^{\frac{C}{r}\times C}$ and $W_2ϵR^{C\times \frac{C}{r}}$. Finally, the normalized weights are weighted to achieve the purpose of extracting specific information, and the output of the entire network module is subjected to a rescaling operation.

2.3.4. Loss Function

In image segmentation, we use the binary cross-entropy (BCE) loss function and cross-entropy loss function to approximate the value between the target and the expected output. The weight update rate of the loss function is proportional to the size of the error. The BCE loss function calculates the cross-entropy of binary classification tasks, as shown in Formula (3):

$$L=-\frac{1}{N}{\sum }_{i=1}^N\left(y_i\log \left(\widehat{y_i}\right)+\left(1-y_i\right)log\left(1-\widehat{y_i}\right)\right)$$

(3)

where N is the total number of samples, $y_i$ represents the true value, and $\widehat{y_i}$ represents the predicted value. The cross-entropy loss function is used to calculate the cross-entropy of multi-classification tasks, as shown in Formula (4):

$$L=-\frac{1}{N}{\sum }_{i=1}^N{\sum }_{c=1}^K{y}_{ic}\log \left(p_{ic}\right)$$

(4)

where K represents the number of categories, $p_{ic}$ represents the probability that the ith sample belongs to category C, and ${\displaystyle \sum }_{c=1}^K{p}_{ic}=1, i=1,2, \cdots , N$. $y_{ic}\in \left\{0,1\right\}$ and can be regarded as a one-hot code; that is, if the ith sample belongs to category C, the corresponding value is 1; otherwise, it is 0. Finally, the attention-based U-Net model, consisting of 11.3 million trainable parameters, is optimized with a learning rate of 0.0001 and a batch size of 16.

2.3.5. Post-Processing

After the training and prediction of the attention-based U-Net model, we obtained the supraglacial lake extent map extracted by SAR. Then, the post-processing steps were carried out for further refinement. First, there is the problem that shadows caused by the radar look angle and terrain relief have similar scattering coefficients to water. We utilized the 10 m ArcticDEM and its derived slope map to reduce the misclassification of lakes by removing areas with slopes greater than 15° [37]. Second, for coastline masking, we used the coastline data of the Greenland ice sheet to better solve the interference problem of ocean pixels and ice-marginal lakes. Third, we projected the result in the WGS 1984 Stereographic North Pole coordinate system to provide 10 m resolution supraglacial lake distribution products.

Furthermore, in some regions, such as the 79°N Glacier, we used the NDWI and normalized difference snow index (NDSI) to extract supraglacial lakes in Sentinel-2 optical images. The results were compared and supplemented with the SAR results. Specifically, we first obtained the mask of rock and seawater by calculating NDSI > 0.85 and B2 < 0.4. Then, the areas with B11 > 0.1 and B10 > 0.01 were excluded as cloud cover. Finally, the pixels with NDWI > 0.18 and B3-B4 > 0.09 were selected to map the lake extent. The threshold selection was based on the results of Moussavi et al. [16].

$$\mathrm{NDWI}=\frac{\mathrm{Blue}-\mathrm{Red}}{\mathrm{Blue}+\mathrm{Red}}$$

(5)

The Blue and Red bands in Equation (5) correspond to B2 and B4 in S2.

$$\mathrm{NDSI}=\frac{\mathrm{Green}-\mathrm{SWIR}}{\mathrm{Green}+\mathrm{SWIR}}$$

(6)

The Green and SWIR bands in Equation (6) correspond to B3 and B11.

2.3.6. Accuracy Assessment

As a data product for supraglacial lake extracted from SAR is not yet available, and there are some differences between optical and SAR results, in this paper, we used artificially labeled supraglacial lake data to randomly select 10 images for accuracy assessment. At present, the evaluation standards that are often used for image segmentation include precision (P), recall (R), accuracy (A), and F1 score (F1). Accuracy means the percentage of correctly predicted samples divided by total samples. However, in the case of unbalanced samples, it cannot measure the quality of the results, so we calculate the F1 score at the same time. F1 considers both R and P and ranges from 0 to 1; the higher the value, the better the accuracy. In more detail, F1 is defined as the harmonic mean of R and P, where R is the fraction of examples classified as positive, and P is the fraction of true positive examples among the examples that the model classified as positive. Following the evaluation, indices are used, where TP indicates the number of positive classes predicted as positive classes, TN indicates the number of negative classes predicted as negative classes, FP indicates the number of negative classes predicted as positive classes, and FN indicates the number of positive classes predicted as negative classes.

$$\mathrm{A}=\frac{\mathrm{TP}+\mathrm{TN}}{\mathrm{TP}+\mathrm{TN}+\mathrm{FP}+\mathrm{FN}}$$

(7)

$$\mathrm{R}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}$$

(8)

$$\mathrm{P}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FP}}$$

(9)

$${\mathrm{F}}_1=\frac{2\times (\mathrm{R}\times \mathrm{P})}{(\mathrm{R}+\mathrm{P})}$$

(10)

3. Results

We used several methods to test the applicability of our attention-based U-Net deep-learning model to supraglacial lake features in different environments, which is a necessary condition for our approach to be widely used in Greenland. The success of the model and the remaining challenges are described in detail below, and the results of the supraglacial lake extraction are discussed.

3.1. Automatic Extraction of Supraglacial Lakes

To make the results more convincing, we selected SAR images from different regions and at different times for testing. They were distributed in the vicinity of typical glaciers and ice shelves, such as the Petermann, 79°N, Koge Bugt, Kangiata, Russell, and Jakobshavn Glaciers. The meltwater in these regions is more abundant and diverse, making it easier to verify the results.

As can be seen, the surface characteristics of supraglacial lakes in different regions vary significantly in shape, size, depth, and location. However, whatever the changes, the proposed method can extract jagged, arc-shaped, or funnel-shaped lakes with diameters ranging from several kilometers to tens of meters, and the distribution is intricate (Figure 8a–h). Specifically, precise edges can be detected for lakes with a minimum diameter of 0.1 km (Figure 8c), blurred edges, or low visual contrast (Figure 8d–f). In addition, in the ablation zone of Greenland, the segmentation of the lake is affected by slush areas (Figure 8d). As slush areas contain large amounts of wet snow, with the increased water content, the backscattering gradually becomes similar to the lakes. Likewise, topographic shadows with similar backscatter coefficients will be easily confused with lakes (Figure 8b). Although the terrain is flat in most parts of the GrIS, some areas still have large topographical undulations at the edge of the ice sheet, forming mountain shadows. However, the proposed attention-based U-Net model can effectively separate supraglacial lakes from slush areas as well as shadows, solving the problem of difficult lake extraction in these areas.

Greenland’s melting typically runs from May to September. It continues at a moderate pace throughout the summer, with a surge in mid-July. Here, we present the supraglacial lake extraction result near the Jakobshavn Glacier in mid-July 2020 (Figure 8e,f). During this period, under the promotion of temperature, rainfall, dark ice, and other factors, the backscattering coefficient of the glacier surface is very small, and there exists a lot of wet snow. This poses a great challenge for lake extraction, and some snow with high water content may be classified as lakes. However, the attention-based U-Net model can solve this problem well, effectively distinguishing lakes from wet snow. Figure 8g,h show a partially snow-covered supraglacial lake extracted at the end of August, which cannot be seen by optical images. We find that even in summer, covered supraglacial lakes, as an important part of summer meltwater transport, are abundant. In the ablation zone, the lake’s appearance is round and larger (Figure 8g).

Further, to enhance the understanding of winter lakes, we implemented a multi-classification of lakes according to the freezing conditions. Specifically, they were divided into supraglacial lakes (unfrozen) and frozen supraglacial lakes (Figure 9a–d). We can find that, in winter, there exists a lot of former lake ground, which used to be a lake but gradually dried up over time and became a sink on the glacier surface (Figure 9a,b). In SAR images, former lake ground and lake ice can be well distinguished. Because the former lake ground is multi-year ice, and its surface is rough, when there is no water or little water, the backscattering is mainly dominated by volume scattering, showing high brightness on the SAR image. However, the SAR signal could detect the water body beneath the lake ice due to its penetration. Therefore, it is not as bright as the former lake ground in the SAR image. Only after 8–10 m of snow is accumulated, or the surface has undergone more melting, is the former lake ground no longer clearly visible in SAR imagery. However, if the former lake ground is filled with water, we consider the area to be a lake.

The overall F1 score across all test scenes was calculated as an average of 0.971, with R = 0.988, P = 0.956, and A = 0.994. The non-supraglacial lake class yielded high values: P = 0.999, R = 0.995, and F1 = 0.997 compared to the supraglacial lake class with P = 0.913, R = 0.981, and F1 = 0.946. Furthermore, to prove the validity and scientificity of the proposed algorithm, in this paper, we used the traditional U-Net, U-Net + ASPP, and attention-based U-Net models. The network segmentation results comparing the same segmentation task are shown in

Table 2. Segmentation results of U-Net, U-Net + ASPP, and attention-based U-Net.

Method	A	R	P	F1
U-Net	0.979	0.845	0.849	0.839
U-Net + ASPP	0.993	0.979	0.954	0.966
Attention-based U-Net	0.994	0.988	0.956	0.971

. The results show that the F1 score of this method is higher than that of the other two methods, with an average F1 of 0.971. It can be seen from Figure 10 that attention-based U-Net performed better on the details of edge segmentation and overall accuracy. The experiments were implemented in Python 3.8 and PyTorch1.5 with an NVIDIA RTX 3080 and an Intel i9 GPU. Predicting a single Sentinel-1 scene with a size of 1.08 GB (116,171 × 117,852 pixels) takes about 42.35 s compared to 1591.78 s for a CPU.

To compare the similarities and differences of supraglacial lake extraction between optical and SAR images, we selected Sentinel-1 and Sentinel-2 images covering part of 79°N in the same period around 30 August 2021, for supraglacial lake detection, and obtained the maximum lake range (Figure 11a). According to the statistics and comparison of the results, the total supraglacial lake area is 215 km², of which optical images can extract 71.9%, and SAR images can extract 57.2%, with an overlap of 62.8 km². SAR images can show more of the covered lakes close to the ice sheet’s interior, and rivers are not misclassified as lakes. In Greenland, especially during the melt season, rapid drainage events of supraglacial lakes will occur under the influence of temperature, precipitation, or topography, resulting in transient supraglacial streams or rivers on the glacier surface (Figure 11b–d). There are many lakes that are too shallow or changing too fast to be captured by the SAR. Therefore, a supraglacial lake product with higher temporal and spatial resolution can be produced by combining the results of Sentinel-1 SAR and Sentinel-2 MSI images.

3.2. Spatiotemporal Change Analysis of Supraglacial Lakes

Using the proposed attention-based U-Net deep-learning method, we extracted the contours of supraglacial lakes in different periods near the Jakobshavn and 79°N Glaciers. These regions have multiple supraglacial hydrological features, including all the areas (dark ice, blue ice, slush, shadow, and wet snow) that are difficult to extract. Figure 11 shows the seasonal changes of supraglacial lakes around these two glaciers from 2017 to 2021, with a time interval of less than 6 days. The statistics of the Jakobshavn region were acquired at the beginning of spring and summer (June and July; Figure 12a). In most cases, when the peak values are excluded, the area gradually increased from June to July and fluctuated with the seasons. During the melting season, including July and August, the total lake area in the 79°N region gradually changes, with a peak mainly occurring around mid-August and then decreasing at the end of the season (Figure 12b). In general, the area of a supraglacial lake should develop steadily with the change of seasons. However, according to the final statistics, both the Jakobshavn and 79°N Glaciers experienced dramatic changes in 2019, 2020, and 2021. Jakobshavn entered a higher stage ahead of schedule in mid-June 2021, and the peaks in 2019 and 2021 were much higher than those during the same period in other years. At the same time, the 79°N Glacier also showed that the peaks in 2019 and 2021 were much larger than those in other years, and in 2020, the supraglacial lake entered the melting period ahead of schedule. According to the analysis, the Greenland region experienced extreme melting events in 2019 and 2021 [38], promoting the formation of supraglacial lakes. As a result, in those two years, the area was much higher than in other years. In addition, because the temperature of 79°N continued to rise in 2020, the lake area entered the melting period early. The specific temperature change is shown in Figure 12c,d.

Figure 13 shows the inland migration of the supraglacial lakes near the 79°N Glacier between July and August in 2017–2021 under the influence of the topography; in the figure, blue and red represent the lakes in July and August, respectively. The scale of lakes was much larger in 2019 and 2021 than in the other three years. The total peak area of supraglacial lakes in these two years was more significant than in 2017, 2018, and 2020, especially in regions with high altitude. However, in years with lower total lake area, such as 2018, the distribution of lakes was more skewed toward lower elevation. Overall, the distribution of supraglacial lakes from 2017 to 2021 expanded inland at higher altitude. If the meltwater is abundant in the Greenland ice sheet, the location of supraglacial lakes is mainly controlled by topographic factors. However, the expansion of supraglacial lakes at high altitude is likely controlled by the spike in surface mass balance, especially near the end of the melt season [28]. Second, significant reductions in the lake area can be attributed to abrupt climate changes or successive drainage events, from which the meltwater structure, in this case, can be inferred.

4. Discussion

4.1. Extraction of Supraglacial Lakes

Although there are many studies on supraglacial lakes, the extraction of lakes on the GrIS is still challenging due to the variety of surface features. Compared with the traditional U-Net deep-learning algorithm, the attention-based U-Net showed improved accuracy, precision, and recall rates, and finally obtained an F1 score of 0.971 (Table 2). Further, it can be clearly seen from Figure 9 that the proposed algorithm is more accurate in mapping the outlines of lakes, with a lower probability of missed or false detection. Specifically, for lakes with different shapes, low contrast, and blurred edges, a clear extent can be extracted adaptively according to the context information (Figure 8a–h). Even lakes with an area less than 0.1 km² can be accurately monitored. When the diameter of the lake is close to the spatial resolution, the boundary may be blurred. Therefore, with higher-resolution SAR images, smaller lakes can be detected. Moreover, several studies have shown that it is challenging to extract supraglacial lakes in slush regions [16,18]. However, the proposed method solved this problem and correctly distinguished features similar to supraglacial lakes, such as dark ice, blue ice, shadows, wet snow, and crevasses (Figure 8d–f and Figure 9b).

Eventually, the whole model can run on a GPU. Compared with using CPU only, the calculation time for a 1.03 GB image is reduced from 1591.78 s to 42.35 s. This shows that the proposed method can process large-scale SAR image data very quickly and is suitable for batch extraction of supraglacial lakes. In terms of overall accuracy, the F1 score is lower for the supraglacial than the non-supraglacial lake class. On the one hand, this is because the training dataset was not large enough, and the proportion was unbalanced, so it is necessary to increase the data volume on this basis. On the other hand, because the applied ice sheet boundary products are not precise enough, some seawater and lake pixels were still confused.

Around the Petermann and Kangiata Sermia Glaciers, we found many covered supraglacial lakes during the 2020 melt season (Figure 8c–h), which is consistent with events reported in other independent studies [31]. Although none of these lakes can be seen in optical images, they contribute a significant amount of meltwater during the summer and therefore cannot be ignored. As can be seen from the comparison of the Sentinel-1 and Sentinel-2 extraction results of supraglacial lakes (Figure 11), in a total area of 215 km², the Sentinel-1 mapping extent includes the supraglacial and covered supraglacial lakes, and the lake range moves toward the center of the ice sheet. However, only using SAR cannot capture the very shallow and rapidly changing lakes in time. When considering the integrity of the hydrological process on the ice sheet surface, it is necessary to overcome the limitations of a single sensor and fuse data from multiple sources, such as optical and SAR data, to realize the mapping of supraglacial lakes and calculate the maximum range mapping product. It should be pointed out here that the problem of overestimating lake area due to rapid lake flow should also be considered.

When calculating the lake area, there is an objective possibility of classifying supraglacial river into the supraglacial lake, although the water body of the river is very small. This is because the proposed attention-based U-Net model is better at extracting polygon features than linear features, such as rivers. It is evident in Figure 11 that the results extracted with NDWI confuse lakes and rivers, but in our method, rivers are not classified as lakes. For some water bodies that exist on the ice shelves or floating ice tongues and ice-marginal lakes, we used Greenland ice sheet boundary shapefiles for cropping. That is, only lakes in the GrIS are considered, and other water bodies are not counted in the lake area. Furthermore, we used the Greenland-wide ice-marginal lake inventory [39] to solve the interference problem of the ice-marginal lake. For supraglacial lakes in winter, according to the backscattering characteristics of water, snow, and thin ice, we added HV polarization data because they contribute more to volume scattering, which is the main scattering mode on the surface of supraglacial lakes in winter. In addition, the existing model was further subdivided to achieve multi-class detection, including supraglacial, frozen supraglacial, and non-supraglacial lakes (Figure 8h). Here, supraglacial lakes refer to bare water bodies and water covered by snow, and frozen supraglacial lakes refer to water bodies with thin ice on the surface. In fact, as the temperature decreases, the lake surface gradually freezes, and part of the lake has water inside, which can remain throughout the whole winter and become an essential part of the englacial hydrological system. Nevertheless, it is difficult to rely on Sentinel-1 imagery to monitor frozen supraglacial lakes at present because the penetration depth of the C-band radar is determined by the physical properties of the medium through which it passes and varies with water content, but it can reach a maximum depth of several meters [26], so it can only achieve superficial layer detection. Therefore, more accurate observations of frozen supraglacial lakes should be combined with the results of ground-penetrating radar. In addition, when the lake surface is windy, or the lake begins to drain, the backscattering coefficient will increase, interfering with the extraction.

4.2. Spatiotemporal Change Analysis

In Section 3.2, we performed statistical analysis of seasonal variation of less than 6 days from 2017 to 2021 in the region near 79°N and the Jakobshavn Glacier. The results show that our method is well suited for multitemporal analysis of GrIS surface hydrological networks. Due to the incomplete coverage of Sentinel-1 data from the summer of 2014 to 2016, we did not cover the entire Sentinel-1 time series but selected 2017 to 2021 for research. In that 5-year period, 2019 and 2021 were characterized by high temperature and extreme melting, and we found that the total lake area in these 2 years had a clear peak effect compared with other years. Especially near the Jakobshavn Glacier, the lake area should regularly increase with increasing temperature in June and July, but due to extreme melt events, the area rose sharply from 9 to 12 July in 2019 and 2021. Since we only calculated the statistics for the Jakobshavn lake area changes in June and July, if the changes occur on a longer time scale, the impact of extreme melting may be more significant, so we will continue to carry out long-term series analysis there. For the 79°N Glacier, the temperature exceeded the melting point in July 2021, the highest on record, but the rapid increase in the lake area occurred in August. This delay may have been caused by the slow melting of thick snow on the glacier’s surface. Moreover, the changes of glacial lakes in 2020 were somewhat abnormal compared with other years and peaked earlier, in early July, which may have been caused by the combined effect of temperature and precipitation before July. All of this highlights the relationship between the development of supraglacial lakes and climatic factors. In the early stage of lake development, even small changes in the temperature and precipitation will cause the lake area to increase rapidly. Taken together, extreme climate events, higher temperatures, days beyond the melting point, and snow thickness at the start of the melting season are all important factors that influence changes in the area of supraglacial lakes. However, the specific factors affecting the distribution of lakes still need further clarification [5,40,41,42].

We also found a clear trend of lake expansion from lower to higher elevation in the 79°N Glacier in the NE region during the 2017–2021 melt season. Although we only have 5 years of research, we believe that the change in supraglacial lakes in the 79°N Glacier region is a persistent phenomenon; in particular, the distribution of lakes in 2018 was different than in the other 4 years. They were rarely distributed at an altitude greater than 900 m, and the area of individual lakes was small. In 2021, more lakes were distributed above 900 m, and the individual lake area was larger. In theory, under the combined action of topography and gravity, supraglacial lakes should be distributed from high to low altitude, but the opposite was found. This may be due to the fact that, at the beginning of the melt, high-altitude lakes are closer to the interior of the ice sheet, and the bottom ice sheet and surface snow are thicker and less affected by warm sea water and temperature, which slows down the development of lakes in high-altitude areas to a certain extent. Therefore, in July, low-altitude glacial lakes develop rapidly. As the temperature increases, the surface of the high-altitude area begins to melt, and lakes begin to form, with a peak in mid- to late August. At the same time, a surface hydrological system develops in the low-elevation area. Many tributaries and moulins allow the lakes to drain from the surface or inside into the ocean or the bottom of the ice sheet, so they gradually disappear, leading to migration of lakes to higher altitudes, but in this process, we did not consider other meltwater channels in our analysis and only focused on supraglacial lakes. However, river-like linear features are seen in the extraction results, suggesting that there may have been more liquid water on and within the glacier than discussed here. How supraglacial lakes transport meltwater between glaciers and subglacial features is unclear. Therefore, it is necessary to use automatic extraction methods to produce large-scale, long-period, high-spatial-resolution, and short-time-interval polar lake products. Based on continuous spatiotemporal change products, research on the Greenland meltwater network can be carried out.

5. Conclusions

Under the influence of global warming, ice sheet surface melting has intensified, further promoting the development of supraglacial lakes, which has a non-negligible impact on the stability of ice sheets and climate change. Compared with optical remote sensing, SAR has more advantages in all-day, all-weather detection of supraglacial lakes, identification of covered lakes, and long-term sequence analysis. This study used an attention-based U-Net deep-learning network to map supraglacial lakes in the GrIS using Sentinel-1 SAR imagery. Based on traditional convolutional neural networks, we improved the U-Net, including the SE attention module, residual connections, and ASPP. The primary purpose was to improve the accuracy of segmentation and the universality of the method. The attention-based U-Net was trained on data covering 49,886 patches of 10 regions and was evaluated with other temporally and spatially independent regions. Ice sheet boundaries and ArcticDEM topographic products were combined to further reduce false lake classification. Compared to traditional convolutional neural network models, our method had the highest accuracy, returning an average F1 score of 0.971 and accuracy of 0.994. The evaluation indices of the supraglacial lake class are P = 0.913, R = 0.981, and F1 = 0.946. The results of all test datasets show that the proposed method achieves good performance on complex regions (blue ice, slush, dark ice, moraine-covered snow) and various morphologies (small size, blurred edges, low contrast, rough surface, ice and snow cover), and it correctly extracts similar features (moulins, ice crevasses, shadows). For winter SAR images, combined with dual polarization data, multi-class segmentation is achieved according to the frozen state of the lakes, including unfrozen, frozen, and non-supraglacial lakes. Additionally, this method utilizes GPU acceleration, making the time to process a single image 30 times shorter, which is suitable for rapid processing of massive data.

Furthermore, we carried out multi-source data fusion for the surface of the 79°N Glacier in August 2020 and automatically generated the maximum extent of the glacial lake corresponding to Sentinel-1 and Sentinel-2 images. The final area is 215 square kilometers, which illustrates the necessity of using both SAR and optical data to capture rapidly changing supraglacial lakes. In addition, according to the 2017–2021 supraglacial lake dynamics of the Jakobshavn area and 79°N and the supraglacial lake products of less than 6 days, we found that extreme melt events in 2019 and 2021 primarily affected the peak lake areas in both regions, making them higher than in other years. In particular, the Jakobshavn Glacier area showed dramatically changed lake peaks between June and July. During the melting season, the supraglacial lakes of the 79°N Glacier have a clear tendency to expand inland; that is, they migrate from low to high altitude. At the same time, there is a certain lag between the maximum lake area and maximum temperature. All of this indicates the influence of climate, topography, surrounding features, and other factors on the development of supraglacial lakes.

Overall, the above results highlight the superiority of the attention-based U-Net model for supraglacial lake segmentation in Sentinel-1 SAR images. In the future, we will continue to integrate multi-source data, including the available optical, SAR, and elevation data, to carry out refined hydrological mapping of the entire Greenland ice sheet. At the same time, a follow-up study will further investigate the polarization information to achieve a physical interpretation of the ice sheet’s surface features, enhancing the interpretability of the deep-learning model. Furthermore, for the analysis of supraglacial lake products, the main factors affecting the distribution of lakes and their controlling effects on englacial and subglacial hydrology should be studied in combination with various other factors, such as climate and topography, in the future. In this way, a complete description of the development and impact of hydrological processes on the ice sheet’s surface can be achieved.