1. Introduction
Given glacier sensitivity to climate change (Roe and others, Reference Roe, Baker and Herla2017) and associated shifts in meltwater production, glacial lakes are highly dynamic and prone to rapid changes in size (Shugar and others, Reference Shugar2020). Between 1990 and 2018, total glacial lake volume has increased by about 48% worldwide (Shugar and others, Reference Shugar2020). The ongoing melting of glaciers is accelerated by elevation-dependent as well as latitude-dependent warming, which causes mountains and areas at high latitudes and altitudes to warm faster compared to low-lying areas (Mountain Research Initiative EDW Working Group, 2015; Palazzi and others, Reference Palazzi, Mortarini, Terzago and Von Hardenberg2019; Xie and others, Reference Xie, Zhu, Qin, Wang, Xu and Wang2023). Latitude-dependent warming is most pronounced above 45°N and peaks with Arctic amplification above 64°N (Wanishsakpong and others, Reference Wanishsakpong, Luo and Tongkumchum2014). A rapid expansion of glacial lakes at higher latitudes has been observed accordingly (Shugar and others, Reference Shugar2020). However, Zhang and others (Reference Zhang2024) demonstrated that glacial lake changes are highly heterogeneous: detached, glacier-fed lakes expanded with glacier retreat in past decades, while ice-dammed lakes, despite growing in number, often lost volume due to the thinning of ice dams. Although long-term observations of Glacial Lake Outburst Floods (GLOFs) are known to contain problematic reporting bias (Veh and others, Reference Veh, Lützow, Kharlamova, Petrakov, Hugonnet and Korup2022), there is nonetheless, good evidence that the overall expansion of glacial lake volume can contribute to increasing frequency of GLOFs (Nie and others, Reference Nie, Liu and Liu2013; Zheng and others, Reference Zheng2021) with serious impacts on downstream communities (Emmer and others, Reference Emmer2022; Taylor and others, Reference Taylor, Robinson, Dunning, Rachel Carr and Westoby2023).
A glacial lake can be defined as a water body that originates from a glacier (Benn and Evans, Reference Benn and Evans2014). Glacial lakes can occur in various locations relative to the glacier: adjacent (ice-marginal), in front (proglacial), on the surface (supraglacial) or beneath (subglacial) (Yao and others, Reference Yao, Liu, Han, Sun and Zhao2018). There is no standard definition of how to classify or categorize glacial lakes. While some studies define glacial lakes by their location with respect to the glacier that supplies the meltwater (Nie and others, Reference Nie, Liu and Liu2013; Shugar and others, Reference Shugar2020; Wang and others, Reference Wang, Zhong, Zheng, Meng and Liu2023), others are more concerned with the hydrological connection between lake and glacier (Wang and others, Reference Wang, Xiang, Gao, Lu and Yao2015; Zhang and others, Reference Zhang2023). As such, comparing regional inventories can be difficult. Zhang and others (Reference Zhang2024) provide a global overview of regional glacial lake inventories, highlighting significant differences in criteria. For instance, the minimum lake area threshold ranges from 0.002 to 0.5 km2, while the buffer distance from contemporary glacier outlines varies between 0 and 10 km (Zhang and others, Reference Zhang, Wang and An2024a).
Given the difficulty of accessing high-alpine regions, applying remote sensing is fundamental to performing large-scale analysis of glacial lake dynamics (Aggarwal and others, Reference Aggarwal, Jain, Lohani and Jain2016). Satellite data such as Landsat and Sentinel-2 images are commonly used to produce glacial lake inventories (Bhardwaj and others, Reference Bhardwaj2015; Wangchuk and others, Reference Wangchuk, Bolch and Robson2022; Andreassen and others, Reference Andreassen, Nagy, Kjøllmoen and Leigh2022). Aerial images are also occasionally used (Zhang and others, Reference Zhang, Wang and An2024a). Glacial lake mapping is challenging compared to classifying other types of water due to common misclassification of terrain and cloud shadows, seasonal ice cover, diverse levels of suspended sediment, turbidity and periodically outbursting behavior in some cases (Wangchuk and Bolch, Reference Wangchuk and Bolch2020; Abderhalden and others, Reference Abderhalden, Bly, Lappe, Andreassen and Rogozhina2024). For these reasons, many studies still rely on manual or semi-automatic mapping of glacial lakes (Andreassen and others, Reference Andreassen, Nagy, Kjøllmoen and Leigh2022; Zhang, and others, Reference Zhang, Wang and An2024a). For semi-automatic mapping techniques, thresholding of a Normalized Difference Water Index (NDWI) (McFeeters, Reference McFeeters1996) is a common approach (Shugar and others, Reference Shugar2020; Andreassen and others, Reference Andreassen, Nagy, Kjøllmoen and Leigh2022). Shugar and others (Reference Shugar2020) provided the first global attempt to map glacial lake areas and volumes using Landsat images from 1990 to 2018, NDWI thresholding and several post-processing filtering steps. Attempts to automate thresholding approaches usually suffer from high commission and omission errors, mostly related to the misclassification of terrain shadows (Shugar and others, Reference Shugar2020; Kaushik and others, Reference Kaushik, Singh, Joshi and Dietz2022). To address the challenges of high misclassification rates and the extensive need for manual intervention, several studies have employed machine learning or deep learning models that integrate multiple input parameters, including spectral bands, radar data and auxiliary datasets, such as surface elevation and slope (Bhardwaj and others, Reference Bhardwaj2015; Wangchuk and Bolch, Reference Wangchuk and Bolch2020; Kaushik and others, Reference Kaushik, Singh, Joshi and Dietz2022). Such approaches have significantly improved overall accuracy, even at large scales. However, these studies usually rely on the use of digital elevation models (DEMs) to mask steep slopes, which can be problematic for the creation of timely inventories, as DEMs are often outdated and can mistakenly mask areas where glaciers have recently retreated or where new lakes formed.
In the present study, glacial lake outlines for Norway are derived using Sentinel-1 and Sentinel-2 imagery from 2023–24. We apply a random forest machine learning classifier trained on a 10th percentile Sentinel-2/ Sentinel-1 summer composite to classify 'water' versus 'other landcover' without involving DEMs. To reduce misclassification caused by mountain shadows, we test a post-processing step that masks water body detection in areas with high difference values between ascending and descending Sentinel-1 images, based on the assumption that such differences will effectively highlight mountain slopes. We apply a largely automated workflow implemented in Google Earth Engine and Python to enhance the efficiency and reproducibility of the presented method. We use very-high-resolution PlanetScope images to validate the detection accuracy of glacial lakes as well as the delineation accuracy of the lake outlines. Finally, our results are compared with the previous glacial lake inventory of Norway from 2018–19 to analyze lake area changes across different glacier environments. Within selected regions, manual validation and correction were applied to ensure accurate comparison and reliable change detection. The national transferability of the method is assessed by reapplying the workflow to classify glacial lakes in 2018–19 within the selected focus regions and comparing the results to the validated change detection.
2. Study area
The study focuses on all glacier regions in mainland Norway, encompassing a wide latitudinal range from approximately 60°N to 71°N (Fig. 1). Norway hosts over 6000 glaciers and ice patches with a total area of around 2300 km2 (Andreassen, Reference Andreassen2022). These are found at elevations from sea level up to around 2450 m a.s.l. Climatically, the glacier regions are characterized by a cold and humid climate, with mean annual temperatures typically between −2°C and +3°C, and annual precipitation ranging from 500 mm to over 2000 mm depending on the location and topography (Hanssen-Bauer and others, Reference Hanssen-Bauer2017). We further selected six representative focus regions, in which validation, manual correction and regional comparison are carried out: Folgefonna (60.0° N), Hardangerjøkulen (60.5° N), Jostedalsbreen (61.5° N), Svartisen (66.7°), Blåmannsisen (67.3° N) and Lyngen (69.5°) (Fig. 1). These regions capture the diversity of Norway’s glacial environments by spanning from 60° to 70°N from Folgefonna in the South to Lyngen in the North and from maritime to more continental settings. The focus regions include the countries’ six largest ice caps: 1) Jostedalsbreen (458 km2), 2) Vestre Svartisen (190 km2), 3) Søndre Folgefonna (154 km2), 4) Østre Svartisen (125 km2), 5) Blåmannsisen (81 km2) and 6) Hardangerjøkulen (64 km2). The regions include glaciers in a wide range of topographic settings—from plateau glaciers such as Hardangerjøkulen to mountain glaciers in steep alpine terrain in Lyngen (Andreassen, Reference Andreassen2022). Several glacial lake inventories have previously been developed for Norway. The earliest digital outlines were created using manual digitization from Landsat NDWI imagery, supplemented with topographic maps and manual corrections during the 1999–2006 glacier inventory. Additional datasets were produced for earlier (1988–97) and later (2014) periods (see overview in Andreassen and others (Reference Andreassen, Moholdt, Kääb, Messerli, Nagy and Winsvold2021) and Andreassen and Winsvold (Reference Andreassen and Winsvold2022)). The most recent and detailed glacier outlines (GO1819) and glacial lake inventory (GLO1819) employed Sentinel-2 imagery from 2018 (northern Norway) and 2019 (southern Norway). The GLO1819 dataset includes 473 glacial lakes that were likely to be in contact with the glaciers at the time of satellite image acquisition with metadata such as glacial lake ID (‘bresjoID’), glacial lake type (‘innsjoType’) and lake area.
Processing grid for the detection of glacial lakes in mainland Norway with highlighted selected focus regions. Glacier outlines from 2018–19 (Andreassen and others, Reference Andreassen, Nagy, Kjøllmoen and Leigh2022). Coordinates are in EPSG: 25833.
3. Data and methods
3.1. Glacial lake definition
Previous studies have applied a range of criteria to define glacial lakes, most commonly based on lake area, distance to glacier boundaries and minimum glacier size. The choice of threshold values varies with study objectives (Zhang and others, Reference Zhang, Wang and An2024a). In this study, glacial lakes were defined as water bodies located within 200 m of the glacier boundaries from GO1819. This will then include more glacial lakes than GLO1819, but can easily be filtered, e.g. using a 20 m buffer as we do in our change assessment analysis. A minimum glacier area of 0.1 km2 and a minimum lake area of 400 m2 was adopted, consistent with GLO1819. The lake size threshold is smaller than those used in several comparable studies (e.g. 0.01 km2 by Wangchuk and Bolch (Reference Wangchuk and Bolch2020)), though the dataset can easily be filtered to apply larger minimum lake sizes if required.
3.2. Input data
A combination of optical (Sentinel-2) and radar (Sentinel-1) satellite imagery was used, along with auxiliary geospatial datasets, to detect and validate glacial lakes across mainland Norway (Table 1). The primary data source was Harmonized Sentinel-2 Surface Reflectance (SR) imagery, a multispectral optical imaging mission developed by the European Space Agency (ESA) as a part of the Copernicus Programme. It consists of two satellites (Sentinel-2A and 2B) providing global coverage every 5 days and capturing data across 13 spectral bands at 10–60 m spatial resolution (Drusch and others, Reference Drusch2012). To complement the Sentinel-2 images and provide all-weather imaging capability, dual-polarized Sentinel-1 Ground Range Detected (GRD) synthetic aperture radar (SAR) data were also used. Sentinel-1 operates in C-band using the Interferometric Wide (IW) swath mode. The pre-processed GRD imagery provides consistent spatial resolution of 10 m (Torres and others, Reference Torres2012).
Overview of data used for the detection, post-processing and validation of glacial lakes in Norway.
For the six selected focus regions (Fig. 1), very high-resolution PlanetScope and RapidEye composites from August 2024 and Third Quarter (Q3) composites of 2024 were included to assist in visual validation and manual correction of lake outlines. The basemaps provide 2.4–4.8 m resolution multispectral imagery, processed to ensure geometric and radiometric consistency (Planet, 2025). Reference datasets included the 2018–19 Norwegian glacier (GO1819) and glacial lake inventory (GLO1819) (Andreassen and others, Reference Andreassen, Nagy, Kjøllmoen and Leigh2022), which served as the baseline for lake ID tracking and for assessing lake area changes.
Lakes from the Norwegian mapping authority, manually digitized from orthophotos (N50), were additionally used to assist with lake validation (Kartverket, 2023).
3.3. Processing set-up
All processing was carried out using the Google Earth Engine (GEE) Python API, allowing for efficient handling of large datasets and automated analysis. GEE is a cloud-based geospatial analysis platform that enables large-scale processing of satellite imagery and other Earth Observation data (Gorelick and others, Reference Gorelick, Hancher, Dixon, Ilyushchenko, Thau and Moore2017). The detection of glacial lakes was performed within a grid framework consisting of 30 km x 30 km cells, based on GO1819 (Fig. 1).
To minimize boundary effects and ensure full detection of lakes near grid edges, each cell was buffered by 1 km. Any duplicate lake features resulting from overlapping buffer zones were resolved and merged in a later step (Section 3.6). The entire processing workflow for the detection and post-processing of glacial lakes is summarized in Fig. 2.
Processing workflow for the classification of glacial lakes.
3.4. Pre-processing
3.4.1. Sentinel-2 10th-percentile summer composite
All available Sentinel-2 images between 1 June and 15 August for the years 2023 and 2024 with less than 80% scene cloud cover were selected for the glacier regions in mainland Norway. This time window was chosen to maximize the likelihood of snow-free and full-lake conditions, while avoiding long-cast terrain shadows from low solar zenith angles late in the season (Fig. S1-7). Subsequent cloud-masking was applied to all images using the CloudScore+ algorithm (Pasquarella and others, Reference Pasquarella, Brown, Czerwinski and Rucklidge2023). A 10th-percentile composite was then generated from the filtered Sentinel-2 images. This approach emphasizes low-reflectance surfaces like water, bare ground and glacier ice (compared to snow), and is at the same time less affected by cloud shadows compared to a minimum composite. The following bands were selected as predictors for the classification: B2 (blue), B3 (green), B4 (red), B8 (NIR), B11 (SWIR1) and B12 (SWIR2). To enhance water detection, the NDWI was calculated as follows and added as an additional input band:
\begin{equation}NDWI = {\text{ }}\frac{{B3 - B8}}{{B3 + B8}}\end{equation}where 
$B3$ represents Sentinel-2’s green band and 
$B8$ the near-infrared (NIR) band. The NDWI is widely used to classify water bodies as it effectively exploits the difference between the low reflectance of NIR and the relatively high reflectance of Green in water surfaces (McFeeters, Reference McFeeters1996). However, when using the index alone, it is prone to misclassifying ice, shadows and snow as water, especially in mountain environments (Abderhalden and others, Reference Abderhalden, Bly, Lappe, Andreassen and Rogozhina2024), which is why we are opting for a more sophisticated detection approach. To ensure the reliability of the NDWI layer for classification, pixels with a value of exactly 0.0 in either the B3 or B8 band were excluded, as these introduced erroneous NDWI values—likely a result of artifacts from the surface reflectance harmonization process. In addition, a median composite of the ‘Dark Area Pixels’ class from the Sentinel-2 ‘Scene Classification Map’ (SCL) was generated and included as an additional input layer. This class, derived from the cloud and scene classification used in Sentinel-2 Level 2A surface reflectance products (ESA, 2021), provides supplementary information on persistent terrain shadows, thereby improving the classifier’s ability to distinguish water from shadowed terrain.
3.4.2. Sentinel-1 ascending-descending difference
To complement Sentinel-2 optical data and reduce misclassification of terrain shadows, Sentinel-1 VV-polarized imagery was incorporated into the classification process. Because of the side-looking radar geometry, terrain facing toward or away from the sensor experiences foreshortening or radar shadowing, respectively (Richards, Reference Richards2009). These distortions are reversed in ascending and descending passes. By subtracting radar signals (in dB) from these passes, terrain slopes exhibit high difference values, which helps distinguish them from flat surfaces like water:
\begin{equation}S1_{diff}=\vert S1_{asc}-S1_{desc}\vert\end{equation}where 
$S{1_{diff}}$ is the result of subtracting radar signals from the ascending (
${\text{}}S{1_{asc}}$) and descending (
${\text{}}S{1_{desc}}$) images (Fig. 3). The absolut ascending-descending difference was calculated from monthly 10th-percentile composites for July 2024 with vertical-vertical (VV) polarization. VV polarization is sensitive to surface roughness and slope orientation (Richards, Reference Richards2009), which was beneficial to highlight the backscatter difference between terrain and water. The 10th percentile was used to emphasize calm water conditions, which exhibit lower backscatter values compared to wind-affected water that causes higher backscatter values from diffuse reflection. To suppress speckle noise, a focal mean filter with a 100 m kernel was applied. This approach was favored over incorporating a DEM to identify topographic shading, as existing DEMs do not reflect recent glacier retreat, which would limit the algorithm’s ability to detect newly formed lakes at retreating glacier boundaries.
Appearance of slopes vs. water in ascending and descending Sentinel-1 image composites (July 2024). Slopes are causing higher difference values allowing for differentiation between water surfaces and mountain slopes, which otherwise often appear similar in both Sentinel-1 and -2 imagery. Example of Rembesdalskåka, Hardangerjøkulen, Copernicus Sentinel-1 data.
3.5. Classification of glacial lakes with Random Forest
A Random Forest (RF) classifier was used to delineate glacial lakes based on the pre-processed Sentinel-2 and Sentinel-1 input layers. RF is a widely used supervised ensemble learning algorithm that constructs multiple decision trees during training and outputs the mode of their predictions (Breiman, Reference Breiman2001). It is particularly well-suited for remote sensing classification tasks due to its robustness to overfitting, ability to handle high-dimensional data and relatively low sensitivity to noise (Belgiu and Drăguţ, Reference Belgiu and Drăguţ2016; Dirscherl and others, Reference Dirscherl, Dietz, Kneisel and Kuenzer2020). The classifier was trained on a composite image (Section 3.4) within the six focus regions (Fig. 1). A total of 1700 random training points were manually labeled, including 700 for the ‘water’ class and 1000 for ‘other land cover’, representing a range of surface types. The samples were randomly distributed and are imbalanced between the two classes to reflect actual class proportions (Moraes and others, Reference Moraes, Campagnolo and Caetano2024). A binary classification scheme was chosen over a multi-class approach to maximize the sensitivity for the detection of water bodies. The training data were randomly split into 70% for training and 30% for validation to assess classification accuracy. The RF classifier was configured with 800 trees. Following classification, the resulting water mask was cleaned using the ‘connectedPixelCount’-function in GEE with a neighborhood size of 100 pixels and a threshold of 12 connected pixels to filter out noise, which improved computational efficiency during vectorization.
To evaluate shadow-related misclassifications in the post-processing (Section 3.6), a secondary RF classifier was trained using the same input image and configuration, but with a multi-class scheme including the following five classes: ‘water’, ‘glacier ice’, ‘bare ground’, ‘vegetation’ and ‘shadow’. This classification served to generate a likelihood score for each classified water body polygon, indicating where water bodies were likely to be confused with shadows. However, this secondary classification was not used for direct water delineation, as it consistently underperformed at lake margins: In many cases, pixels representing shallow water were misclassified as shadow, reducing the geometric accuracy of the lake outlines compared to the binary classification.
3.6. Post-processing
To ensure that only glacier-associated lakes were included in the dataset, the GO1819 glacier outlines (Table 1) were buffered by 200 meters. Only classified water bodies that intersected with these buffered zones were included, in line with the defined glacial lake criteria (Section 3.1). To reduce the inclusion of false positives, particularly misclassified shadows, two filtering steps were applied based on ancillary raster layers:
1. Radar backscatter filtering: For each classified water body, the mean value of Sentinel-1 ascending-descending backscatter difference was calculated. Water bodies with a mean difference larger than 9 were excluded. To justify the chosen threshold, we identified pixels with slope values >10% (ArcticDEM (Porter and others, Reference Porter2023)) and identified all water pixels (Dynamic World V1 (Brown and others, Reference Brown2022)) within the six focus regions (Fig. 4). The Sentinel-1 difference value was sampled at 4000 random points, and a threshold was set as the maximum of the 85th percentile of water and the 5th percentile of slopes >10%, retaining most water while excluding steep terrain.
Figure 4.Distribution of absolute Sentinel-1 backscatter differences (|S1_asc–S1_desc|) for water and slopes (>10%), displayed as box plots (orange) and violin plots (blue). The horizontal dashed line indicates the chosen threshold, set to retain the majority of water pixels (85th percentile) while excluding steep terrain.
2. Shadow probability filtering: The proportion of each water body’s area classified as ‘shadow’ during the secondary classification was calculated, and only polygons with less than 40% shadow pixels were retained. This threshold was chosen as a balance between removing features likely caused by terrain shadows and preserving genuine water bodies that are partially shadowed, such as those situated in steep valleys. Tests with stricter thresholds removed several known lakes.
Both a filtered and an unfiltered version of the classified water body dataset were exported from GEE, each including metadata on Sentinel-1 backscatter difference, shadow percentage and the processing tile identifier. Following export, the classified water bodies for all processing tiles were merged into a single, national dataset. In the overlapping zones, where polygons from adjacent processing tiles intersected, geometries were dissolved and the metadata from the larger polygon was retained, while areas were re-calculated. To enable comparison with previous glacial lake inventories (Table 1), water bodies absent from the filtered dataset but present in the unfiltered dataset and intersecting with a previously mapped lake were included in the filtered dataset again. To further support temporal analysis, the GLO1819 glacial lake ID was assigned to all new water bodies that spatially intersected a lake. In cases where multiple new polygons overlapped with a single previous lake, all were assigned the same ID and counted as a single lake in the inventory
As a final step, the geometries of the filtered glacial lake dataset were smoothed using the QGIS ‘v.generalize’ function. The ‘Chaiken’ generalization algorithm (Chaikin, Reference Chaikin1974) was selected after several tests with a maximum tolerance of 18, which was applied as an optimal threshold to smooth irregular polygon edges while maintaining the overall shape and area of the lake. A minimum area threshold of 200 m2 was used to remove small artifacts and spurious polygons generated during the generalization process. To further refine the outlines and eliminate remaining angularities, the generalization was applied a second time with a lower tolerance value of 3, which provided visually consistent realistic lake boundaries (Fig. 5). Since ‘v.generalize’ converts multi-part polygons into single-part features during smoothing, a unique polygon ID was assigned prior to this step. After smoothing, the original multipolygon structure was restored using these IDs.
Effect of Chaiken smoothing on glacial lake outlines.
3.7. Validation
To assess the quality of the final glacial lake dataset, a multi-level validation approach was implemented. This included (1) a pixel-based accuracy assessment of the RF classifier, (2) an object-based validation of lake detections across the six focus regions and (3) an evaluation of the delineation accuracy, quantifying the geometric precision of lake outlines.
3.7.1. Detection accuracy
A standard confusion matrix was generated using the 30% holdout of manually collected validation points for ‘water’ and ‘other land cover’ that were not used during training. This matrix provides common accuracy measures such as overall accuracy, producer’s accuracy and user’s accuracy for each class. However, a pixel-based accuracy assessment does not reliably reflect the quality of a classification that is primarily interested in the detection of discrete objects. This is especially true in the presence of large lakes, where correctly labeled pixel abundance may mask omission or commission errors on smaller lakes.
To more accurately reflect detection reliability, we adopted a region-based object validation strategy, as recommended by Radoux and Bogaert (Reference Radoux and Bogaert2017) within the selected focus regions (Fig. 1). Each classified lake was manually compared against Very High Resolution (VHR) PlanetScope basemaps from August and the third quarter (Q3) of 2024 (Table 1), the Sentinel-2 composite used in the classification (Section 3.4), the GLO1819 lake inventory N50 lake layer (Table 1) and occasionally Google Satellite imagery. Each detected object was categorized as either True Positive (TP), where a lake was detected correctly, False Positive (FP), where a misclassified object was labeled as a lake or ‘unclear’, where it was impossible to identify the land cover type from available products (see Table 2). When a lake was visible in the reference imagery, but not classified by our algorithm, it was manually digitized (or derived from previous inventories, where appropriate) and labeled as False Negative (FN). A total of 402 lakes were manually labeled to derive accuracy metrics (Table 2), 33 of which were manually delineated. Precision measures how many of the features identified as positive are actually correct, while recall measures how many of the actual positive items were successfully identified. The F1-score is the harmonic mean of precision and recall, providing a balanced metric (Sokolova and Lapalme, Reference Sokolova and Lapalme2009).
Object-based accuracy metrics to assess the detection performance of glacial lakes using, based on visual classification into true positive (TP), false positive (FP) or ‘unclear’ cases.
3.7.2. Delineation accuracy
In addition to evaluating detection quality, we assessed the delineation accuracy of lake outlines from successfully detected lakes—also referred to as ‘geometric precision’ or ‘positional quality’ (Radoux and Bogaert, Reference Radoux and Bogaert2017). We selected 20 lakes across five of the six focus regions (Fig. 1), to include both, fully ice-free and partially ice-covered examples at different latitudes and in varying topographic settings. Reference outlines were manually digitized from PlanetScope Basemaps (Table 1), offering reliable basis for comparison due to their high spatial resolution (∼3 m) and spatial alignment with Sentinel-2 imagery. The PlanetScope Normalized Surface Reflectance (SR) Basemaps are co-registered to multi-year monthly composite mosaics from Sentinel-2 imagery ensuring high spatial consistency between the two datasets and minimizing systematic geolocation differences (Planet, 2025). To quantify outline differences, we generated evenly spaced perpendicular transects from simplified reference polygons (Fig. 6). For each transect, intersection points with both the classified and reference outlines were identified, and the distance between them was measured. The median of these distances per lake was used to summarize the delineation error. This comparison was performed on both original and smoothed outlines to assess the impact of post-processing on boundary accuracy.
Comparison of classified glacial lake outlines (red) with manually digitized validation outlines (green dashed) for five lakes across different glacier regions in Norway (locations shown in the overview map). Black lines represent validation transects used to compute outline deviation. Background: PlanetScope Basemap August 2024. Lakes shown from regions: (a) Folgefonna, (b) Hardangerjøkulen, (c) Jostedalsbreen, (d) Blåmannsisen, (e) Lyngen.
3.8. Manual correction in focus regions
Within the focus regions, all incorrectly outlined glacial lakes were manually corrected to enable realistic region-based comparisons of glacial lake area changes between GLO1819 and the 2023–24 dataset. Not detected lakes were either manually digitized or copied from previous inventories (GLO1819 or N50), in cases where the lake outline has not changed.
4. Results
4.1. Accuracy of glacial lake detection
4.1.1. Detection accuracy
The performance of the RF classifier, based on pixel-level validation using the validation samples, showed that out of 212 samples, 7 ‘water’ pixels got misclassified as ‘other land cover’, while 8 ‘other land cover’ pixels were erroneously classified as ‘water’ (Table 3). Accordingly, the confusion matrix indicates misclassification rates, yielding a high overall accuracy of 0.97.
Confusion matrix for Random Forest classifier trained on ‘water’ vs. ‘other land cover’ with User’s (UA) and Producer’s Accuracy (PA) and Overall Accuracy in the bottom right.
The object-based validation across the six focus regions (Fig. 1), which accounts for spatial entity detection (Table 4), reveals a high overall recall value of 0.91, meaning that 90% of all actual lakes were correctly identified. However, the overall precision is more moderate at 0.73, indicating that around 30% of the detected lakes were FPs—mostly small, misclassified features with a median area of 1.5 km2. The F1-score, which balances these two metrics, is 0.81 on average, reflecting solid overall performance of the glacial lake detection.
Regional validation of glacial lake detection accuracy. Overview of true positives (TPs), false positives (FPs) and false negatives (FNs) across six glacier regions in Norway. Where it was not possible to identify the absence or presence of a lake, it was labeled ‘unclear’ (UC). The measures accuracy (ACC), precision (PRE), recall (REC) and F1-score (F1) were derived, as well as the median area of false positives (MedFP) and false negatives (MedFN) in m2.
Yet, detection performance varied considerably between regions. The highest overall accuracy, precision and recall were achieved for Hardangerjøkulen, where no FPs or FNs were recorded. Folgefonna also performed well, with a high F1-score of 0.95 driven by both high precision (0.97) and recall (0.93), despite a small number of FNs. Similarly, the Svartisen and Blåmannsisen regions showed strong recall values (0.96 and 1.0), indicating that nearly all true lakes were detected. However, both regions had moderate precision due to a higher number of FPs. In contrast, detection performance in Lyngen was notably lower, with a precision of just 0.36, suggesting a high rate of FPs, although recall remained strong (0.97). Jostedalsbreen exhibited a lower recall (0.67) due to a substantial number of FNs, despite relatively good precision (0.87).
4.1.2. Delineation accuracy
The median distances between the classified and validation outlines for all lakes are 6.2 m after smoothing and 6.5 m for the original outline (Table 5).
Median distance [m] between classified and validation outlines across five glacier regions and two lake conditions (ice-free and partially ice-covered). Values are reported for both the original and smoothed classified outlines. All median errors fall below the 10 m spatial resolution of Sentinel-1 and 2, indicating high geometric accuracy of the (successfully) lake boundaries.
The effect of smoothing the automatically derived glacial lake outlines varied by region and condition (Table 5): Among the glacier regions, Hardangerjøkulen showed the highest delineation accuracy, with a median offset of 4.9 m for smoothed outlines. Folgefonna had the largest deviation with 7.5 m after smoothing, although the difference compared to the original outline (7.6 m) was minimal. For ice-free lakes, the median offset was 5.9 m (smoothed) compared to 6.3 m (original), while for lakes partially covered by ice, the delineation accuracy decreased slightly, with a median of 7.3 m (smoothed) and 7.2 m (original). Across all regions and conditions, the median delineation error remains below the spatial resolution of the Sentinel-1 and 2 sensors (10 m), which were used for classification, indicating a high degree of geometric accuracy in lake boundary detection.
4.1.3. Filtering effects
The initial classification resulted in 6440 potential glacial lakes before any filtering was applied. The filtering step, based on the Sentinel-1 ascending-descending backscatter difference (>9), reduced the number of potential glacial lakes to 1360. A subsequent shadow-based filter using the proportion of shadow-classified pixels further narrowed the dataset to 1246 using an 80% shadow threshold. A stricter threshold of 40% reduced the count to 1023. To support temporal consistency with previous inventories, 359 lakes were re-included into the final dataset through spatial intersection with features from GLO1819 and N50 lake inventories (Table 1), resulting in a total of 1382 glacial lakes after post-processing.
4.2. Glacial lake detections in Norway
As described above, a total of 1382 glacial lakes with a combined area of 124 km2 were automatically detected across mainland Norway prior to any manual correction. To ensure a consistent comparison with the previous Norwegian inventory (GLO1819) and better isolate trends driven by climate and glacier dynamics, we applied a 20 m proximity-to-glacier threshold to both datasets and excluded regulated lakes from the analysis based on the NVE hydropower lake dataset (Table 1). This approach minimizes potential differences arising from threshold variations or artificial lake-level regulation. In 2023–24, 681 glacial lakes with a total area of 36.6 km2 were then identified, representing an increase compared to the 331 lakes with a total area of 30.4 km2 in GLO1819—corresponding to an overall increase in glacial lake area of approximately 21%. A comparison based on the lake identifier (‘bresjøID’) shows that 362 lakes are newly detected in the 2023–24 dataset, while 49 lakes from GLO1819 are no longer present in the updated classification. Previously undetected or new lakes make up for 13% of the total lake area in 2023–24 with a median size of 0.0018 km2. An analysis of lake size distribution reveals a significant difference in scale between the two inventories with a median lake area of 0.004 km2 in 2023–24, compared to 0.015 km2 in GLO1819. This shows that mainly very small lakes were newly detected in 2023–24. It should be noted that the glacial lake detection in 2023–24 reflects a combination of true hydrological changes, such as the emergence of new lakes or lake expansion at retreating glacier fronts, differences in methodology, including the automated nature of the new classification, and the presence of false detections.
A visual comparison between GLO1819 and the 2023–24 dataset reveals that the applied method successfully and automatically detected glacial lakes. We found visible lake expansion between 2018–19 and 2023–24 across different glacial regions in Norway (Fig. 7).
Glacial lake detections showing lake growth and the development/detection of new lakes between 2018–19 (dotted blue outline, labeled with the lake ID) and 2023–24 (pink outline) with a Sentinel-2 false color composite (NIR, Red, Green) from summer images in 2023–24 in the background. Glacier outlines from 2018–19 (GO1819) are displayed in light blue. The example lakes are selected from different glacier regions in Norway: (a) Hardangerjøkulen, (b) Blåmannsisen, (c) Jostedalsbreen, (d), (e) Folgefonna, (f) Lyngen, (g) Folgefonna, (h) Hardangerjøkulen, (i) Jostedalsbreen. The first two rows show examples of glacial lake growth, while the third row shows examples of new lake detections in 2023–24.
In this period, glacier front variation records reveal glacier retreat (Kjøllmoen and others, Reference Kjøllmoen, Andreassen and Elvehøy2024). The lake growth appears to be directly associated with glacier retreat. Our examples show glacier retreat between 50 m (Fig. 7c) and 190 m (Fig. 7a), causing lake area growth between 1% (Fig. 7c) and 200% (Fig. 7f). Seasonal ice cover or artefacts in the 2023–24 Sentinel-2 composite (Fig. 7b) can cause errors in total area and area change.
Some of the newly detected lakes have developed due to glacier retreat (Fig. 7g and h), while others were likely already present during the 2018–19 mapping but not included in GLO1819, probably due to partial or complete ice and snow cover at the time of image acquisition, suboptimal spectral separation or omission during manual correction (Fig. 7i).
FP detections were primarily caused by terrain shadows misclassified as water (Fig. 8d–f). In some cases, adjacent shadows directly connected to the lake surface could not be removed through post-processing filters (Fig. 8e), while incomplete lake outlines were mostly caused by seasonal ice cover (Fig. 8a–c). There are several causes for lakes present in GLO1819 but not detected in 2023–24, which include the presence of ice cover (Fig. 8g), insufficient image quality or small temporary water bodies that were either undetectable or disappeared (Fig. 8h) or several smaller lakes that merged into one expanded lake due to glacier retreat (Fig. 8i).
Misclassification issues in the 2023–24 inventory (pink outline). (a)–(c) Show incomplete lake detections due to seasonal ice cover, while (d)–(f) highlight misclassified terrain shadows. (g)–(i) Show lakes present in GLO1819 and absent in 2023–24 due to misclassification caused by ice-cover (g), being undetected or no longer present (h) or small lakes present in 18/19 that now merged into a larger lake due to glacier retreat (i). Outlines from the previous inventory (GLO1819) are overlaid (blue dotted outline, labeled with lake ID) with a Sentinel-2 composite from summer images in 2023–24 in the background. Glacier outlines (GO1819) are displayed in light blue (labeled with glacier ID). The example lakes are selected from different glacier regions in Norway: (a)–(c)/(g), (h) Folgefonna, (c), (d) Blåmannsisen, (f) Lyngen, (i) Svartisen.
4.3. Lake area changes between 2018–19 and 2023–24 in selected regions
To enable a reliable regional comparison of glacial lake changes between 2018–19 and 2023–24, we used the manually validated and corrected dataset (Sections 3.7 and 3.8), filtered to only TP detections. Out of the 402 manually checked lakes (Section 3.7.1), 24 needed manual correction and 33 lakes were added as they remained undetected by the automated mapping. The change analysis was performed for the six focus regions: Folgefonna, Hardangerjøkulen, Jostedalsbreen, Svartisen, Blåmannsisen and Lyngen (Fig. 1), covering 52% of the glaciers (glacier area) and 74% of the glacial lakes (glacial lake area) referring to GO1819 and GLO1819 (Andreassen and others, Reference Andreassen, Nagy, Kjøllmoen and Leigh2022).
When considering all glacial lakes in the selected focus regions within 20 m from the GO1819 glacier boundaries and excluding regulated lakes, the number of detected glacial lakes increased from 107 in 2018–19 to 162 in 2023–24 (50% increase) (Fig. 9a). All regions consistently saw an increase in lake count, with Svartisen showing the largest gain at 105% and Lyngen the smallest with 17% (Table S1). However, most of this increase is attributable to small lakes (<0.0018 km2).
Changes in glacial lake count (a) and changes in glacial lake area (b) for selected regions sorted from South to North between 2018–19 and 2023–24 with manually corrected outlines excluding regulated lakes within 20 m of the GO1819 glacier boundaries. For each region (left panel), as well as for the total across all focus regions (right panel), lake areas are shown for matched lake IDs present in both inventories (dark bars), and for unmatched IDs representing new or vanished lakes (light bars).
The total glacial lake area in the selected focus regions applying the same criteria was 16 km2 in 2018–19 and 20 km2 in 2023–24, resulting in a total area increase of 22% (Fig. 9b). This reflects both the expansion of existing lakes, the formation of new lakes and the detection of previously unmapped or non-included lakes in GLO1819. Among the regions, Lyngen experienced the most significant glacial lake area growth with 563%, while Svartisen exhibited the smallest change at 4% (Table S1). When considering only glacial lakes that were present in both inventories (i.e. matched lake IDs), the area increase for unregulated lakes was more moderate at 8.6%, suggesting that much of the area gain is due to newly formed or newly detected lakes. Matching lakes increased most prominently in the northernmost region, Lyngen, with a gain of 117%, followed by the southernmost region, Folgefonna, at 57%. Jostedalsbreen’s matching glacial lakes grew by 7%, while Hardangerjøkulen and Blåmannsisen each recorded a modest increase of around 3%. In contrast, Svartisen was the only region showing a decline in matching lake area of −7% (Table S2).
5. Discussion
5.1. Quality and uncertainties of the automated glacial lake detection
5.1.1. Composite creation
The 10th percentile compositing approach proved effective in capturing snow-free and filled lake conditions, providing an advantage over classifications based on median composites or single-scene images, as used in earlier studies (Bhardwaj and others, Reference Bhardwaj2015; Wangchuk and Bolch, Reference Wangchuk and Bolch2020; Andreassen and others, Reference Andreassen, Nagy, Kjøllmoen and Leigh2022; Kaushik and others, Reference Kaushik, Singh, Joshi and Dietz2022) (Fig. 10). This method is particularly well-suited for large-scale automated mapping in regions with persistent cloud cover and short summer seasons.
Sentinel-2 composite from cloud-masked summer images between 1 June and 15 August for 2023–24 in false color (NIR, Green, Blue). Median composite (a) compared to 10th percentile composite (b), which is used for our glacial lake classification, showing the difference in depicting snow- and ice-free as well as filled lake conditions. Example of Rembesdalskåka, Hardangerjøkulen, Copernicus Sentinel-2 data.
However, a key limitation is the loss of accurate acquisition date information and the potential introduction of artefacts such as residual clouds or terrain shadows. The effectiveness of the composite depends strongly on the careful selection of the summer season period. In Norway’s glacier regions, high precipitation and frequent cloud cover (Hanssen-Bauer and others, Reference Hanssen-Bauer2017) often prevent the use of imagery from a single summer season to achieve nationwide coverage (Fig. 11c).
Effect of using one (c), two (b) or three (a) summer seasons (1 June–15 August) on the quality of 10th percentile summer composites from Sentinel-2 Surface Reflectance imagery.
The limited image availability causes increased snow cover, data gaps and residual cloud or shadow (Fig. 11c). Incorporating two summer seasons (Fig. 11b) significantly improves the representation of glacial lakes, capturing mostly ice-free and cloud-free conditions. Adding a third summer season, however, yields only marginal improvements, while increasing the risk of including outdated conditions (Fig. 11a). We found that two years were the best compromise for mapping the glacial lakes in Norway.
Snow depth data from 2022–24 (Figs. S1-S7) shows that the duration of snow- and ice-free conditions generally decreases with both latitude and elevation, with only minor regional differences among the focus regions. Southern glaciers such as Folgefonna (Fig. S2) and Hardangerjøkulen (Fig. S3) typically experience snow-free conditions from June to September below 1000 m a.s.l., whereas further North, at Jostedalsbreen, and in northern Norway at Blåmannsisen, and Lyngen, snow and ice often persist well into late summer, and high-elevation regions may remain snow-covered year-round (Figs. S4-S7). Consequently, the successful detection of glacial lakes at higher elevations can be challenging due to persistent snow cover, especially when combined with unfavorable cloud conditions. Although extending the composite period into September could improve lake visibility in some regions, it would also increase the risk of introducing terrain shadows from low sun angles (Fig. 12b) (Nagy and Andreassen, Reference Nagy and Andreassen2019; Dirscherl and others, Reference Dirscherl, Dietz, Kneisel and Kuenzer2020; Wangchuk and Bolch, Reference Wangchuk and Bolch2020; Andreassen and others, Reference Andreassen, Moholdt, Kääb, Messerli, Nagy and Winsvold2021). Therefore, a conservative summer season period from 1 June until 15 August was selected to balance snow-free conditions with optimal illumination for Sentinel-2 composite generation.
Effect of the summer season length on the quality of 10th percentile summer composites from Sentinel-2 Surface Reflectance imagery. Here, showing increasing cast shadow from low sun angles when including September images (b).
Unlike typical approaches for detecting water in mountain environments (Bhardwaj and others, Reference Bhardwaj2015; Zhao and others, Reference Zhao, Chen and Zhang2018; Wangchuk and Bolch, Reference Wangchuk and Bolch2020; Kaushik and others, Reference Kaushik, Singh, Joshi and Dietz2022), we deliberately avoided using a DEM to mask slopes. The available DEMs are derived from imagery prior to this inventory update and would incorrectly eliminate areas where glacier retreat has created newly exposed lake basins or led to lake expansion. Instead, we applied a novel approach using a Sentinel-1 composite of differences from ascending and descending orbits. While it is contribution to the RF classification was moderate (Fig. 13), the Sentinel-1 difference composite proved effective in removing falsely classified water bodies in the post-processing based on the mean difference value for each classified feature.
Feature importance of input bands to the Random Forest classifier to classify ‘water’ vs. ‘other land cover’. ‘VV_p10’ represents the Sentinel-1 ascending-descending difference composite.
Even though the object-based filtering does not solve the issue of lake-adjacent terrain shadows, a purely pixel-based filtering from SAR imagery would have been unsuitable due to the inherent geometric distortions (e.g. foreshortening, layover), which compromise the spatial fidelity of lake outlines.
5.1.2. Detection and delineation accuracy
While the RF classifier achieved a high overall accuracy, this metric is somewhat misleading as it only reflects the proportion of correctly classified water pixels. As a result, large lakes, containing many more pixels, tend to dominate the accuracy assessment, masking potential misclassifications in smaller lakes. Furthermore, the significant reduction of glacial lake detections following the Sentinel-1 difference post-processing step suggests a lower reliability of the RF classification than suggested by the confusion matrix (Table 3). It is quite clear that the RF classification cannot provide reliable glacial lake classification on its own and heavily depends on the implemented post-processing steps, especially with regard to misclassified terrain shadows. In future work, the potential impact of increasing the training dataset, comparing several machine learning classifiers and further testing of multi-class approaches for the initial lake detection should be considered.
To achieve a more realistic picture of the final dataset’s validity prior to manual editing, we adopted a multi-level validation framework. The detection accuracy assessment with an overall F1-score of 0.8 indicates a generally high reliability of the proposed glacial lake classification with very high recall (0.9) and moderate precision (0.7). In other words, most actual lakes were successfully detected, while some FPs remained. While the relatively high rate of FP detections presents a limitation, the intentional balance between commission and omission errors reflects a study’s design choice: for producing a national-scale glacial lake inventory, over-detection is preferable to omission, as FPs can be more easily removed through manual inspection than missing lakes can be identified and digitized. However, in other applications without manual post-processing, a more balanced setup with e.g. stricter filtering could be advisable. This filtering could for example involve a higher minimum lake area threshold, as FPs were predominantly small features (<0.001 km2), which would fall below thresholds used in similar studies (e.g. 0.01 km2 in Wangchuk and Bolch (Reference Wangchuk and Bolch2020)). Detection performance varied regionally and was likely influenced by topography. The flatter ice caps of Hardangerjøkulen and Folgefonna showed the highest accuracy scores, while regions with steep slopes and narrow valleys, like Lyngen and parts of the Jostedalsbreen area had considerably more FPs and missed detections, largely due to shadow-related misclassifications. Seasonal ice-cover—more prevalent at higher altitudes and latitudes—may also explain some of the missed lakes (Figs. S1-7).
The delineation accuracy assessment that was carried out for successfully detected lakes shows very good agreement between classified and reference outlines, with median deviations below the spatial resolution of the sensor. Even lakes with partial ice cover showed only minimal divergence (∼0.5 m) compared to ice-free lakes. However, our validation method did not account for the misclassification of icebergs as islands inside lake polygons since it only compared the outermost outline intersections at the transects. Smoothing improved the delineation accuracy slightly and was helpful in reducing pixel-level noise at lake edges (Fig. 5). Nonetheless, its value is largely aesthetic, potentially suggesting a level of detail beyond what the data provides.
5.1.3. Filtering effects
The Sentinel-1-based filtering proved effective in eliminating many FPs. However, filtering based on shadow detection was more problematic, as many real lakes—especially those in steep terrain—were excluded due to persistent terrain shadow intersecting with the lake. To address this, lakes filtered out due to shadow overlap were reintroduced if they intersected with outlines from previous inventories. While practical in this context, such a method is dependent on the availability and quality of prior datasets. Transferring this workflow to other regions without high-quality previous inventories would require additional testing and analysis to determine optimal thresholds.
5.1.4. Methodological transferability and limitations
The cloud-based workflow combining Python with GEE and a compositing approach enabled efficient processing at national scale (approx. 6 hours of processing time) without the need for local computational resources. This makes the proposed method both time-efficient and reproducible for future glacial lake inventories in Norway. Manual post-editing to correct the remaining ∼10% of misclassified or undetected lakes within the focus regions required about one additional working day.
To evaluate the transferability of the workflow within Norway, we reran the glacial lake classification for 2018–19 across the focus regions and compared it to the unmodified 2023–24 dataset (both without manual corrections) (Fig. 14).
Glacial lake area changes comparing 2018–19 and 2023–24 lake outlines (before manual correction), both produced with our method. Change is shown in selected regions, sorted from South to North, within 20 m of the GO1819 glacier boundaries, excluding regulated lakes.
Both datasets were filtered to a 20 m buffer around the glacier outlines (GO1819) and regulated lakes were excluded to ensure comparability between the time periods.
The results indicate an overall increase of glacial lake area of 19% (Table S3), consistent with the 9–22% increase observed between the validated dataset GLO1819 and 2023–24 (Fig. 9b). The comparison further shows a significant increase in glacial lake count of 67% (Fig. S3), which is slightly higher compared to the corrected results of the GLO1819 and 2023–24 comparison. However, due to the presence of FP detections (Table 4), the lake count can only reliably be compared after manual correction. All regions show an expansion in lake area, in line with our previous results. Overall, these findings highlight the reproducibility and robustness of our proposed method for the Norwegian context.
Several factors limit the immediate transferability of the workflow to other mountain regions: The high detection accuracy achieved here partially depends on the availability of a detailed and validated previous inventory, which is not always available elsewhere. In such cases, the method remains applicable but would require either more extensive post-processing or a greater degree of manual cleaning, particularly when using small minimum area thresholds. Another constraint is the reliance on Sentinel-1 ascending and descending image pairs for slope masking. This approach is not feasible in regions where only one orbital path is available, e.g. Russia, Antarctica or parts of Canada and Greenland (ESA, 2025b). Moreover, the current classifiers are trained exclusively on Norwegian data. Applying the method in other regions would likely require retraining with locally representative samples, as well as adjustments to the post-classification thresholds S1 ascending-descending backscatter difference and shadow probability (Section 3.5).
That said, Norway presents a particularly challenging test case for automated detection due to persistent cloud cover, short snow-free seasons, low sun angles and the small size of many glacial lakes. In other mountain regions, more favorable imaging conditions may allow for effective compositing, potentially reducing the need for extensive post-processing or allowing for the use of only one summer season.
5.2. Spatial distribution of lake area changes in Norway
Our results reveal a significant increase in both the number and total area of glacial lakes across the six focus regions, both when using the GLO1819 and the newly calculated 2018–19 dataset with our method. The more automated workflow and the use of the dual-season (10th percentile) composite in the 2023–24 mapping enabled better detection of lakes that may have been ice-covered in single-date images used for the earlier inventory. A drawback with our method, however, is that associating the date of mapping to the glacial lake outlines is more difficult as the time span is larger.
To enable meaningful comparisons between the two inventories, the 2023–24 dataset was manually refined and restricted to unregulated lakes within the six selected focus regions (Fig. 1). After harmonizing both datasets, we found that total glacial lake area increased by about 22% since 2018–19, while comparing only lakes present in both inventories yielded a more conservative estimate of roughly 9%. The true rate of change likely falls between these values. This is further supported by the comparison with glacial lake extents from 2018–19 derived using the same classification approach (Fig. 14). Almost all regions show an increase in total lake area, generally associated with glacier retreat (Fig. 7), with the largest relative growth observed in the northernmost region, Lyngen. This pattern aligns with expectations of stronger glacial mass loss at northern latitudes (Mountain Research Initiative EDW Working Group, 2015; Palazzi and others, Reference Palazzi, Mortarini, Terzago and Von Hardenberg2019; Shugar and others, Reference Shugar2020; Xie and others, Reference Xie, Zhu, Qin, Wang, Xu and Wang2023). In particular, Langfjordjøkelen, one of Norway’s northernmost ice caps, recorded the largest cumulative negative mass balance of −31.6 m w.e. and a mean annual balance of −0.9 m w.e. a−1 (Kjøllmoen and others, Reference Kjøllmoen, Andreassen and Elvehøy2024), supporting the observed latitudinal trend in glacier retreat and associated lake expansion.
The observed increase in glacial lake number is consistent with global evidence that accelerated glacier retreat increases meltwater availability and leads to the formation of overdeepened bedrock and moraine-dammed depressions that gradually fill with water and form new lakes (Shugar and others, Reference Shugar2020; Zhang and others, Reference Zhang2024; Zhang and others, Reference Zhang, Wang and An2024a). These global inventories further indicate that regions experiencing the strongest glacial mass loss also exhibit the largest relative increases in both glacial lake number and area, a pattern that is also evident in our results
6. Conclusions
This study presents a largely automated method for mapping glacial lakes across mainland Norway using Sentinel-1 and Sentinel-2 data, a machine learning classifier and a dual-season compositing approach. The method generally proved both efficient and accurate, detecting over 1380 glacial lakes with high detection and delineation accuracy scores. When comparing with Norway’s 2018–19 inventory prior to manual correction, the number of detected lakes has more than doubled and total glacial lake area increased by approximately 21% under consistent criteria. After manual correction in selected focus regions and exclusion of regulated lakes, we estimate that the glacial lake area increased by 9–22% between 2018–19 and 2023–24, where the lower bound represents true lake growth only, while the upper bound also includes newly detected lakes. Regionally, glacial lake growth was most pronounced in northern Norway, consistent with higher glacial mass loss observed in the area. Two key methodological advances contributed to the effectiveness of the workflow: the use of a 10th-percentile Sentinel-2 composite, which reduced the influence of snow and ice cover and the Sentinel-1 ascending-descending difference, which enabled slope masking without relying on DEMs. Despite the improvements, the need for manual correction remains, particularly in steep and high mountain environments with ice-covered lakes year-round, terrain shadows and very small water bodies. These limitations highlight the need for continued refinement of automated approaches and for visual validation to ensure reliable inventories. Further testing of the RF classifier’s robustness by increasing training samples could increase the reliability of the method. Overall, the workflow is time-efficient, reproducible and well-suited for continued monitoring of glacial lakes in Norway. The successful application of our method to classify glacial lakes in 2018–19 shows its transferability within the country. With appropriate retraining and regional calibration, the method holds potential for adaptation to other mountain regions, although its transferability may be constrained by data availability.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2026.10131
Data availability
The 2023–24 glacial lake outlines (with manual corrections) and extent of focus regions described in this paper are available through the Nasjonal Vitenarkiv: https://doi.org/10.58059/4kr8-d860.
Acknowledgements
This paper is a contribution to Copernicus bretjeneste/NVE Copernicus tjenester. We thank two anonymous reviewers as well as the editors H. Jiskoot and G. Zhang for constructive input to the manuscript.
Author contributions
RL and LMA designed the study, RL performed data processing, analysis, validation and led the manuscript writing. LMA supervised the project, contributed to conceptualization, data interpretation and manuscript revision.
Open access
