Introduction
The Greenland Ice Sheet has exhibited increasingly negative mass loss since the early 2000s, driven primarily by increased surface melting and changes in ice dynamics (Meredith and others, Reference Meredith, Pörtner, Roberts, Masson-Delmotte, Zhai, Tignor, Poloczanska, Mintenbeck, Alegría, Nicolai, Okem, Petzold, Rama and Weyer2019; Otosaka and others, Reference Otosaka2023; The GlaMBIE Team and others, 2025). Around the margins of the Greenland Ice Sheet, meltwater accumulates in ice surface depressions during summer, forming tens of thousands of supraglacial lakes (Hu and others, Reference Hu2021; Dunmire and others, Reference Dunmire2025). In a warming climate, supraglacial lakes are increasing in number, area and elevation range due to enhanced surface melting, a trend projected to continue under future climate scenarios (e.g., Leeson and others, Reference Leeson2015; Igneczi and others, Reference Igneczi2016; Otto and others, Reference Otto, Holmes and Kirchner2022; Fan and others, Reference Fan, Ke, Luo, Shen, Livingstone and Lea2025).
Supraglacial lakes are most abundant during the melt season, when meltwater ponding lowers albedo and can locally enhance ice temperature and increase firn density (e.g., Lüthje and others, Reference Lüthje, Pedersen, Reeh and Greuell2006; Law and others, Reference Law, Arnold, Benedek, Tedesco, Banwell and Willis2020; Jullien and others, Reference Jullien, Tedstone, Machguth, Karlsson and Helm2023). Throughout summer, these lakes may drain slowly through overspill, partially or rapidly through hydrofracture (e.g., Das and others, Reference Das2008; Tedesco and others, Reference Tedesco, Willis, Hoffman, Banwell, Alexander and Arnold2013; Dunmire and others, Reference Dunmire2025). Lake drainage events during the melt season are well documented and impact ice sheet dynamics through short-term local influences on ice motion (e.g., Das and others, Reference Das2008; Tedesco and others, Reference Tedesco, Willis, Hoffman, Banwell, Alexander and Arnold2013; Doyle and others, Reference Doyle2014; Stevens and others, Reference Stevens2015). Lakes that do not drain during summer freeze over at the end of the melt season. Of all supraglacial lakes present across the Greenland Ice Sheet during the melt seasons of 2018 and 2019, 
$\sim$35% and 
$\sim$45%, respectively, persisted beyond summer (Dunmire and others, Reference Dunmire2025).
Over winter, supraglacial lakes can be categorised into four types based on the behaviour they exhibit, each having a different impact on ice dynamics (Fig. 1). First, lakes that drained during summer leave an empty lake basin behind and are therefore not a winter lake. Second, lakes that freeze over at the end of the melt season may sporadically drain over winter, generating a delayed, sudden release of stored meltwater into the drainage system (Schröder and others, Reference Schröder, Neckel, Zindler and Humbert2020b; Benedek and Willis, Reference Benedek and Willis2021; Maier and others, Reference Maier, Andersen, Mouginot, Gimbert and Gagliardini2023; Dean and others, Reference Dean, Scharien, Willis and McDougall2026). Third, lakes may freeze-through entirely, forming solid ice blocks that contribute to an expanding area of impermeable ice slabs that ultimately enhance runoff potential (Dunmire and others, Reference Dunmire, Banwell, Wever, Lenaerts and Datta2021; Reference Dunmire2025; Jullien and others, Reference Jullien, Tedstone, Machguth, Karlsson and Helm2023). Fourth, lakes may retain liquid water during winter, insulated by a cover of ice and snow, referred to as buried lakes (Koenig and others, Reference Koenig2015). Of all meltwater lakes observed across the Greenland Ice Sheet over summer and winter in 2018 and 2019, buried lakes comprised 21% and 19%, respectively (Dunmire and others, Reference Dunmire2025). Such buried lakes present a form of perennial meltwater storage beyond the melt season and act as a buffer against meltwater discharge from the ice sheet, temporarily mitigating the ice sheet contribution to sea-level rise (e.g., Koenig and others, Reference Koenig2015; Dunmire and others, Reference Dunmire, Banwell, Wever, Lenaerts and Datta2021). In the following melt season, buried lakes may then overflow early in the season due to limited capacity to store additional meltwater, triggering the formation of meltwater channels (Koenig and others, Reference Koenig2015; Lampkin and others, Reference Lampkin, Koenig, Joseph and Box2020).
Overview of the types of supraglacial lake behaviours over the course of a hydrological (May–April) year and the resulting glaciological implications. Lakes form and may drain during the summer months (blue section). During the winter months (red section), lakes that drained in summer expose the empty lake basin. Those lakes that persist into winter will either drain, fully freeze through or remain liquid (i.e., buried lakes). Each winter lake type results in different glaciological implications, with examples listed here (grey section). Lake images are from Copernicus Sentinel-2 scenes acquired in 2020.
Figure 1 Long description
A flow diagram shows supraglacial lake processes over a May–April hydrological year. In summer, a lake either drains or does not drain. Summer drainage produces an empty lake basin with no lake in winter, resulting in no winter lake behaviour and greater surface cooling. A lake that persists after summer can follow one of three winter pathways: freeze-through, burial beneath ice ans snow, or drainage. Freeze-through leads to densification of ice or firn, latent heat release, and increased runoff potential. Burial indicates liquid water presence and meltwater retention, leaving less capacity for future storage. Winter drainage modifies ice flow and the subglacial drainage system. Optical imagery is available for observing the summer processes, while radar imagery is available throughout summer and winter.
While supraglacial lake monitoring has, to date, largely focused on the melt season due to the long-term availability of optical satellite imagery, there has been a growing effort to monitor lake activity year-round using radar satellite data (e.g., Johansson and Brown, Reference Johansson and Brown2012; Miles and others, Reference Miles, Willis, Benedek, Williamson and Tedesco2017; Zheng and others, Reference Zheng2023; Dunmire and others, Reference Dunmire2025). Radar satellites provide high spatio-temporal coverage and are largely unaffected by atmospheric conditions or seasonal darkness, both of which impact optical imagery. Supraglacial lakes can be identified in radar imagery based on their backscatter signature; however, the signal is highly complex and thus presents major challenges for automated methods of classification and interpretation (e.g., Fahnestock and others, Reference Fahnestock, Bindschadler, Kwok and Jezek1993; Rott and Davis, Reference Rott and Davis1993; Parrella and others, Reference Parrella, Hajnsek and Papathanassiou2016; Chen and others, Reference Chen, Li, Chen, Ju and Cheng2022; Li and others, Reference Li, Chen, Lin, Hooper, Chen and Cheng2023).
Sentinel-1 dual-polarised C-band Synthetic Aperture Radar (SAR) imagery has been used to capture lake behaviours based on changes in backscatter intensity over winter, with specific temporal backscatter patterns linked to winter drainage, refreezing and buried lake presence. However, research on winter lake behaviours based on SAR backscatter to date has focused either on locally or regionally isolated lakes (Koenig and others, Reference Koenig2015; Miles and others, Reference Miles, Willis, Benedek, Williamson and Tedesco2017; Dunmire and others, Reference Dunmire2020; Schröder and others, Reference Schröder, Neckel, Zindler and Humbert2020b; Benedek and Willis, Reference Benedek and Willis2021; Li and others, Reference Li, Lhermitte and López-Dekker2021; Deakin and others, Reference Deakin, Willis and Dell2025; Dean and others, Reference Dean, Scharien, Willis and McDougall2026) or larger-scale mapping of winter lake persistence at specific points in time without considering the range of lake behaviours over the entire winter season (Johansson and Brown, Reference Johansson and Brown2012; Dunmire and others, Reference Dunmire, Banwell, Wever, Lenaerts and Datta2021; Zheng and others, Reference Zheng2023). As a result, a consistent framework for monitoring winter lake processes across space and time is still lacking, and a systematic assessment is needed to understand the role of lakes in winter-time meltwater routing, drainage and storage, and their implications for the ice sheet dynamics.
In this study, we investigate whether temporal patterns in backscatter intensity reflect lake behaviour across space and time, and whether these patterns can be generalised to support SAR-based classification of winter lake behaviour at larger scales. We apply a combined method of: 1) winter lake drainage detection using the approach of Benedek and Willis (Reference Benedek and Willis2021); and 2) k-means time series clustering as an unsupervised approach for classifying all other winter lake behaviours without the need for training data. To assess the generalisability of our method for winter lake behaviour classification, our analysis focuses on two study sites on the Greenland Ice Sheet, Russell Glacier and the Northeast Greenland Ice Stream (NEGIS), across two climatically contrasting years, 2018 and 2019. The results provide the basis for developing a framework for automated, year-round monitoring of supraglacial lake activity, complementing well-established melt season monitoring and contributing to a more comprehensive understanding of supraglacial lake dynamics in a changing climate.
Data and methods
Study sites
Two study sites representing contrasting glaciological regimes were selected (Fig. 2). Russell Glacier in Southwest Greenland (4667 km
$^2$) is land-terminating and characterised by relatively slow ice flow (Moon and others, Reference Moon, Gardner, Csatho, Parmuzin and Fahnestock2020) (2018–19 mean annual ice velocity over the study site extent = 94 m a
$^{-1}$; based on NASA MEaSUREs ITS_LIVE data (Gardner and others, Reference Gardner2025)). In contrast, the two main outlet glaciers of the North East Greenland Ice Stream (NEGIS; Nioghalvfjerdsbræ and Zachariæ Isstrøm, 18589 km
$^2$) are marine-terminating and among the fastest-flowing regions of the Greenland Ice Sheet (Mouginot and others, Reference Mouginot2015) (2018–19 mean annual ice velocity over the study site extent = 224 m a
$^{-1}$; based on NASA MEaSUREs ITS_LIVE data (Gardner and others, Reference Gardner2025)). Higher ice velocities are associated with increased ice deformation and surface roughness, which can affect the radar backscatter signal (Fahnestock and others, Reference Fahnestock, Bindschadler, Kwok and Jezek1993). Assessing radar backscatter across sites with different glaciological conditions may therefore reveal the extent to which findings can be generalised across the ice sheet.
Location of the two study sites in Greenland: (a) Russell Glacier; (b) Northeast Greenland Ice Stream (NEGIS). Black polygons represent outlines of lakes investigated in this study from Glen (Reference Glen2025) for (a) (derived from optical imagery over the 2018 and 2019 melt seasons) and Schröder and others (Reference Schröder, Neckel and Humbert2020a) for (b) (derived from polarimetric SAR using data from 2017 to 2020). Background images: World Imagery; Esri, Vantor, Earthstar Geographics and the GIS User Community.
Figure 2 Long description
The figure consists of two detailed satellite maps separated by a locator map of Greenland. The locator map marks panel a at Russell Glacier on the western ice-sheet margin and panel b at the Northeast Greenland Ice Stream, or NEGIS, on the northeastern margin. Panel a shows a single study-area boundary extending inland from the western margin, with investigated lake outlines scattered across the area. Panel b shows a larger study-area boundary near the northeastern coast, containing a denser distribution of investigated lake outlines.
The two study sites are also subject to contrasting regional climatic conditions, resulting in differences in supraglacial lake formation and distribution. Southwest Greenland records the regionally highest mean annual air temperatures across the ice sheet, whereas Northeast Greenland experiences consistently lower mean annual air temperatures (2018–19) (Zhang and others, Reference Zhang2019). Despite these differences, multi-year, ice-sheet-wide lake mapping (2017–22) has revealed the highest concentration of buried lakes along the southwest and northeast margins of the Greenland Ice Sheet (Zheng and others, Reference Zheng2023), indicating favourable conditions for winter lake persistence and possibly freeze-through and drainage in both study areas.
Supraglacial lake formation exhibits strong inter-annual variability linked to climatic forcing (Zhang and others, Reference Zhang2019; Dunmire and others, Reference Dunmire2025; Glen and others, Reference Glen2025). To test whether our workflow captures inter-annual variability in winter lake formation, two consecutive years with contrasting climatic conditions in both areas were selected for analysis: 2018–19 winter following a low summer temperature and low melt year, and 2019–20 winter following a high summer temperature and high melt year (Zhang and others, Reference Zhang2019).
Sentinel-1 pre-processing
Radiometrically terrain-corrected Sentinel-1 Ground Range Detected products were acquired through the Alaska Satellite Facility Vertex data portal. Level-1 Sentinel-1A/-1B acquisitions in interferometric wide (IW) swath mode, with a 6 day repeat cycle, were selected from September 1st to April 30th for 2018–19 and 2019–20 (Supplementary Material S1). To minimise any signal variability due to ‘look angle’ differences, we only use imagery from ascending flight directions and from relative orbits 90 (incidence angle range = 30
$^{\circ}$–49
$^{\circ}$ for image, 38.6
$^{\circ}$–44.2
$^{\circ}$ for study site extent) and 74 (incidence angle range = 30.4
$^{\circ}$–45.6
$^{\circ}$ for image, 35.0
$^{\circ}$–44.9
$^{\circ}$ for study site extent) for Russell Glacier and NEGIS, respectively (Benedek and Willis, Reference Benedek and Willis2021; Chen and others, Reference Chen, Li, Chen, Ju and Cheng2022). We did not normalise backscatter for different look angles, as has been recently conducted (e.g., Li and others, Reference Li, Chen, Lin, Hooper, Chen and Cheng2023; Dean and others, Reference Dean, Scharien, Willis and McDougall2026). The defined time period begins with the onset of refreezing in the autumn and concludes prior to the onset of melting the following spring (Brangers and others, Reference Brangers, Lievens, Miège, Demuzere, Brucker and De Lannoy2020), hereafter referred to as ‘winter’. We refer to the remainder of the year prior and after winter as ‘summer’, representing the melt season. Although the duration of the melt season can vary across years and sites, we defined the winter period conservatively to exclude any end of summer backscatter signal and to maintain consistency and comparability between study sites and years.
The following pre-processing parameters were selected for the Sentinel-1 data product: backscatter coefficient normalisation to sigma-nought (sigma0) radiometry, calibration to decibel (dB) scale, application of an Enhanced Lee speckle filter (7 x 7), and set output pixel spacing to 20 m. The HV cross-polarised channel was selected for radar backscatter analysis as its greater signal-to-noise ratio and enhanced ice-water contrast improve detection of subsurface water layers compared to the co-polarised HH channel (Miles and others, Reference Miles, Willis, Benedek, Williamson and Tedesco2017).
Extraction of backscatter time series
For both study sites, we used existing supraglacial lake outlines as an entry point to our radar-based supraglacial lake assessment, consistent with approaches adopted in earlier studies (Miles and others, Reference Miles, Willis, Benedek, Williamson and Tedesco2017; Benedek and Willis, Reference Benedek and Willis2021; Zheng and others, Reference Zheng2023). Supraglacial lake outlines were taken from Schröder and others (Reference Schröder, Neckel and Humbert2020a) at NEGIS and Glen (Reference Glen2025) at Russell Glacier. NEGIS lake outlines were created using a polarimetric SAR approach with year-round Sentinel-1 data from 2017 to 2020, providing a more comprehensive lake dataset than can be generated from optical imagery, which may miss perennially buried lakes (Schröder and others, Reference Schröder, Neckel, Zindler and Humbert2020b). Russell Glacier lake outlines were produced by applying the Normalised Difference Water Index (NDWI) to optical imagery over the 2018 and 2019 melt seasons (Glen and others, Reference Glen2025). For both sites, we use the combined maximum summer lake extent as our lake outlines. We filled ring-shaped lake outlines in the optically derived Russell dataset, on the assumption that these included lake areas covered by floating ice, for consistency with the NEGIS data, as Sentinel-1 can detect liquid water beneath floating ice covers, whereas NDWI can only detect open water surfaces. Overall, a total of 430 lakes at NEGIS (lake areas ranging from 0.1 to 25.9 km
$^2$) and 77 lakes at Russell Glacier (lake areas ranging from 1.1 to 6.1 km
$^2$) were included in the final dataset.
To examine the temporal evolution of supraglacial lake backscatter, we extracted the mean backscatter signal within each individual lake outline across all Sentinel-1 images to create the time series later used in the lake drainage detection and clustering analysis. We normalised the lake backscatter signal with respect to the surrounding ice surface, following previous studies (e.g., Benedek and Willis, Reference Benedek and Willis2021; Li and others, Reference Li, Lhermitte and López-Dekker2021; Dunmire and others, Reference Dunmire2025), by subtracting the mean backscatter from a 1 km buffer outside each lake outline from the respective mean lake backscatter. Finally, to further examine misclassifications post-clustering, we also extracted backscatter intensities along a transect corresponding to the longest profile across each lake (Li and others, Reference Li, Lhermitte and López-Dekker2021). Lake transects were extended by 1 km at either end to include the surrounding ice backscatter.
Detection of lake drainage
The first step of our method was to detect lakes that drain during winter, following Benedek and Willis (Reference Benedek and Willis2021). Drainage events are temporally isolated and therefore inherently difficult to capture using clustering approaches (Alqahtani and others, Reference Alqahtani, Ali, Xie and Jones2021). However, the backscatter signature of winter drainage events in SAR time series is well established Schröder and others (Reference Schröder, Neckel, Zindler and Humbert2020b); Benedek and Willis (Reference Benedek and Willis2021); Maier and others (Reference Maier, Andersen, Mouginot, Gimbert and Gagliardini2023); Dean and others (Reference Dean, Scharien, Willis and McDougall2026), allowing us to identify such lake behaviour. Specifically, a sudden increase in backscatter between consecutive acquisitions, anomalous relative to the rest of the time series (
$ \gt 1.5$ than average backscatter increase between consecutive images (Benedek and Willis, Reference Benedek and Willis2021)) and followed by a sustained, stable backscatter level for the remainder of the season, is interpreted as indicative of lake drainage (Schröder and others, Reference Schröder, Neckel, Zindler and Humbert2020b; Benedek and Willis, Reference Benedek and Willis2021; Dean and others, Reference Dean, Scharien, Willis and McDougall2026)(Fig. 3a). We used the script from Benedek and Willis (Reference Benedek and Willis2021) to identify any winter-draining lakes prior to clustering.
Example of a winter lake drainage event detected on Russell Glacier in 2018. (a) Change in mean backscatter of the lake compared to the surrounding ice background before, during and after the drainage event in November. (b) Pre- and post-winter ArcticDEM (2 m) differencing (10 Oct 2018 and 10 Apr 2019, respectively, Supplementary Material S4, Porter and others (Reference Porter2022)) over the same lake (A), showing a decrease in elevation over the lake body and thus supporting the detection of a winter drainage event. For comparison, lakes labelled (B) and (C) are neighbouring winter lakes that did not drain in winter 2018 and do not show any elevation change over this period.
Figure 3 Long description
The figure combines a backscatter time series with an elevation-difference map for a winter drainage event at Russell Glacier. Panel a shows that the lake mean backscatter rises gradually from about minus 20 decibels in September to about minus 16 decibels in October, then increases sharply to approximately minus 10 decibels in early November. It subsequently remains increases from one timestep to the next to the surrounding ice background mean. Thereafter both the lake and background mean remains comparatively stable at approximately minus 12 to minus 11 decibels for the remainder of the time series. Panel b maps elevation change from October 2018 to April 2019 using a scale from minus 3.5 to plus 3.5 metres. Lake A contains a strong red area indicating surface lowering of at least 3.5 metres, whereas neighbouring lakes B and C show no significant elevation change.
Clustering of mean backscatter time series
In the second step of our method, we classified all other winter lake behaviours using a time series clustering method. We selected a k-means algorithm for clustering as it is fast, computationally efficient at handling large amounts of data, and unsupervised (Vlachos and others, Reference Vlachos, Lin, Keogh and Gunopulos2003; Aghabozorgi and others, Reference Aghabozorgi, Seyed Shirkhorshidi and Ying Wah2015). Prior to clustering, all time series data were normalised to a [0, 1] range to ensure comparability across samples. All analyses were performed in Python (version 3.11.9) using the Scikit-learn (Pedregosa and others, Reference Pedregosa2011) and sktime (Löning and others, Reference Löning, Bagnall, Ganesh, Kazakov, Lines and Király2019; Király and others, Reference Király2025) libraries for time series normalisation and clustering, respectively. Clustering was then performed using the TimeSeriesKMeans algorithm from the sktime library, employing Euclidean distance as the similarity metric for time series of equal length (Alqahtani and others, Reference Alqahtani, Ali, Xie and Jones2021).
The k-means time series clustering approach requires the number of clusters k as a user-defined input (Aghabozorgi and others, Reference Aghabozorgi, Seyed Shirkhorshidi and Ying Wah2015; Alqahtani and others, Reference Alqahtani, Ali, Xie and Jones2021). To objectively determine the optimal number of clusters, a combined Elbow and Silhouette analysis was performed on the time series data for cluster numbers from 2 to 10. The Elbow method was used to identify the point of diminishing returns in inertia reduction, while the Silhouette analysis evaluated cluster compactness and separation (Rousseeuw, Reference Rousseeuw1987; Mihevc and Podržaj, Reference Mihevc and Podržaj2025). Based on the balance of these criteria, k = 4 was selected as the optimal number of clusters across all years and study sites (Figs S1–S4). Euclidean distance was selected as the metric to assess similarity between time series. Time series clustering was performed separately for each study year and site to ensure independent and context-specific clustering outputs, capturing both inter-annual and regional variability.
Validation of winter lake classification
Lake behaviour associated with a distinct temporal backscatter pattern was classified into four different winter lake types: 1) freeze-through; 2) buried; 3) drainage; or 4) an empty lake basin with no winter lake presence (Fig. 1). This classification was validated using two independent datasets that indicate some, or all, of these lake processes: 1) start and end of summer lake labels derived from optical imagery; and 2) radargrams from NASA’s Operation IceBridge.
Sentinel-2 optical satellite imagery was manually investigated to derive categorical labels characterising lake conditions at the end of the melt seasons preceding and the onset of the melt season following each winter (following previous observations, e.g., Koenig and others (Reference Koenig2015); Miles and others (Reference Miles, Willis, Benedek, Williamson and Tedesco2017); Otto and others (Reference Otto, Holmes and Kirchner2022); Dunmire and others (Reference Dunmire2025)). Pre-winter labels defined whether the lake had either: i) drained during the melt season and left an empty lake basin (‘drained‘); or ii) retained water that froze over and was covered with an ice lid at the end of the melt season (‘refreeze’; Table 1). All lakes optically labelled ‘refreeze’ therefore represent and are considered to be ‘winter lakes’, whereas ‘drained’ optical labels represent the absence of winter lakes. Post-winter labels captured lake onset at the start of the following melt season, with meltwater either: i) accumulating as an open water surface with no ice lid in an empty lake basin (‘refill‘); or ii) already present as liquid water beneath an ice lid that melted gradually from the margins of the lake inwards over consecutive timestamps (‘liquid’; Table 1). Lakes optically labelled ‘liquid’ typically became visible before surrounding empty lake basins began to ‘refill’, and retained liquid water through winter (buried lakes). Lakes optically labelled ‘refill’ either indicate no winter lake being present, or a winter lake that was frozen through and thus began to accumulate water on the surface of the previously frozen lake (Table 1).
Overview of optical imagery labels from the preceding and following summer, as well as the presence of subsurface attenuation in Operation IceBridge (OIB) echograms from April/May, with the presented combination of labels being indicative of the four winter lake types shown in Fig. 1 and used to validate the Sentinel-1 winter lake classification results.
Table 1 Long description
The table has four columns: optical imagery label in the preceding summer, optical imagery label in the following summer, whether subsurface attenuation is present in an Operation IceBridge echogram from the following April or May, and the resulting winter lake validation label. Four combinations are listed. A lake labelled drained in the preceding summer and refill in the following summer, with no subsurface attenuation, is labelled no lake. A lake labelled refreeze in the preceding summer and refill in the following summer, with no subsurface attenuation, is labelled freeze-through. A lake labelled refreeze in the preceding summer and liquid in the following summer, with subsurface attenuation present, is labelled buried. A lake labelled refreeze in the preceding summer and refill in the following summer, with no subsurface attenuation, is labelled drainage.
Radar echograms from the US National Aeronautics and Space Administration’s Operation IceBridge mission were available for both study sites during spring 2019 and were acquired wherever the flight track crossed a supraglacial lake basin (Paden and others, Reference Paden, Li, Leuschen, Rodriguez-Morales and Hale2014) (Figs S5 and S6). Each echogram was manually inspected for subsurface attenuations of radar backscatter signal, indicating the presence of internal layers of buried water (Koenig and others, Reference Koenig2015; Chen and others, Reference Chen, Howat and De La Peña2017; Lampkin and others, Reference Lampkin, Koenig, Joseph and Box2020; Zhu and others, Reference Zhu, Zhou, Zhu, Wang and Zhang2023). At Russell Glacier, a total of 17 lakes intersected with an IceBridge profile acquired between 5 and 16 May 2019. At NEGIS, a total of 27 lakes intersected with an IceBridge profile acquired between 5 and 23 April 2019. Since these echograms were available at the very end of, or just after, our winter time series, any detected subsurface water representing buried lakes is interpreted to have persisted throughout winter 2018 and is used as additional validation to optical labels.
Results
Detection of lake drainage
Only five lakes showed evidence of winter lake drainage: a single lake at Russell Glacier in 2018 (image pair dates: 7–13 November 2018); and four lakes at NEGIS, two in 2018 and two 2019 (image pair dates: 25 September–1 October 2018, 30 November–6 December 2018, 8–14 October 2019, 26 October–1 November 2019). Of these five lakes, an ArcticDEM image pair from the onset and end of the winter period was available for the one drained lake at Russell Glacier in 2018, and showed a negative elevation change within the lake basin consistent with the identified drainage event (Fig. 3b). Additionally, all five drained lakes were labelled as ‘refreeze’ in the preceding summer and ‘refill’ in the summer following the winter when drainage occurred, corroborating the drainage behaviour interpretation (Table 1).
Time series clustering
The k-means clustering approach captured four distinct temporal patterns that were mostly consistent across both study areas and both winter periods (Fig. 4). Cluster A captures a rapid, coherent rise in backscatter within the first 1–2 months, followed by a stable high plateau for the remaining winter. Cluster B exhibits a steady, gradual increase in signal across the entire winter that does not reach a plateau. Cluster C represents more subtle patterns, with only a slight increase between the start and end of winter, with the exception of Russell 2019, which showed a very slight decrease. The individual time series in cluster C are punctuated by large backscatter fluctuations, spanning the full range of normalised values [0–1]. Cluster D is characterised by a gradual decrease in backscatter across winter, also punctuated by large fluctuations, with the exception of Russell 2019, which showed less pronounced fluctuations across all time series around the cluster mean.
Time series clustering output for Russell Glacier and NEGIS for 2018 and 2019. All individual normalised lake time series are shown in grey (n = number of lake time series per cluster), with the cluster mean plotted in red. The lake behaviours, interpreted post-clustering, are labelled beneath the cluster title.
Figure 4 Long description
A four-by-four matrix compares clusters A to D across Russell Glacier in 2018 and 2019 and NEGIS in 2018 and 2019. Each plot spans September to April and shows individual normalised lake backscatter time series in grey and their mean in red. Cluster A, interpreted as freeze-through, rises rapidly from near zero during September and October before remaining on a high plateau. Cluster B, interpreted as buried lakes, increases more gradually throughout winter and approaches its maximum by March or April. Clusters C and D are both interpreted as no lake. Cluster C is variable and fluctuates around intermediate normalised values, whereas Cluster D generally decreases through winter. Number of time-series included per clusters A, B, C and D are respectively 20, 21, 10 and 25 for Russell 2018; 11, 16, 21 and 29 for Russell 2019; 65, 60, 124 and 179 for NEGIS 2018; and 113, 126, 82 and 107 for NEGIS 2019.
Cluster interpretation: Linking patterns to lake processes
We interpreted the temporal patterns in clusters A–D (Fig. 4) by: i) comparing individual mean lake backscatter time series to the mean background backscatter time series of the surrounding ice (Fig. 5, 2.1); and ii) inspecting all lake backscatter profiles along transects across each lake basin (examples for each cluster shown in Fig. 5, 2.). Lakes in clusters A and B exhibit distinctly lower backscatter signal than the surrounding ice at the start of the winter time series (Fig. 5, 2.1). This behaviour is consistent with observations of both terrestrial and supraglacial lakes that contain liquid water beneath an ice lid at the end of summer, which typically display lower backscatter relative to adjacent water-free ice surfaces due to: a) the strong dielectric contrast at the ice-water interface, causing specular reflection of backscatter away from the sensor, and b) enhanced microwave absorption by liquid water (Morris and others, Reference Morris, Jeffries and Weeks1995; Hirose and others, Reference Hirose, Kapfer, Bennett, Cott, Manson and Solomon2008; Antonova and others, Reference Antonova2016).
Example backscatter time series information for lakes in clusters A–D from Russell Glacier in 2018. (1) Lake cessation in summer preceding the investigated winter, captured with Sentinel-2 optical imagery. (2) C-band SAR backscatter from the lakes over winter, with (2.1) showing the mean lake signal (blue) against mean background signal (orange); and (2.2) showing backscatter along the lake transects (lake outlines are marked by vertical dotted lines, with an additional 1 km showing the surrounding ice surface on either side). (3) Lake onset in summer following the investigated winter, captured with Sentinel-2 optical imagery. Lake images are from Copernicus Sentinel-2 scenes acquired in 2018 and 2019.
Figure 5 Long description
A four-row matrix presents Russell Glacier lakes from 2018. Rows show Cluster A, freeze-through; Cluster B, buried; and Clusters C and D, no lake. Columns show optical images before winter, winter mean lake and background backscatter, lake-transect backscatter, and optical images after winter. In Cluster A, the lake refreezes before winter and refills afterwards. Mean backscatter rises from about minus 18 to minus 12 decibels in September, converging with the background, while the initially low transect signal increases towards surrounding-ice values. In Cluster B, the lake refreezes before winter and is liquid afterwards. Mean backscatter rises gradually but remains below the background, with low backscatter persisting within the lake boundaries. In Clusters C and D, the lakes drain before winter and refill afterwards. Their lake and background means remain similar and nearly constant, and their transects lack the persistent low-backscatter signal seen in Cluster B.
For lakes in cluster A, the backscatter time series show a pronounced increase in backscatter across the lake basin within the first 1–2 months of winter, after which the signal converges with the higher backscatter of the surrounding ice. This early-winter backscatter increase is typically attributed to progressive ice-cover thickening as the internal structure of the congealed ice changes and an increased amount of microwave energy is scattered off entrained bubbles and reflected back to the sensor (Morris and others, Reference Morris, Jeffries and Weeks1995; Hirose and others, Reference Hirose, Kapfer, Bennett, Cott, Manson and Solomon2008). This lake freeze-through transition occurs over a timescale of weeks to months, and once lake backscatter is comparable to the surrounding ice, the basin is interpreted as fully frozen through (e.g., Miles and others, Reference Miles, Willis, Benedek, Williamson and Tedesco2017). In contrast, lakes in cluster B display a gradual increase in mean backscatter that consistently remains below the background ice signal throughout winter, a pattern evident across most of the basin transect. Such behaviour has previously been linked to the persistence of liquid water beneath an insulating ice cover, representing overwinter water storage and is thus characteristic of buried lakes (Koenig and others, Reference Koenig2015; Deakin and others, Reference Deakin, Willis and Dell2025). We thus interpret lakes in cluster A as those that completely freeze-through and lakes in cluster B as buried lakes that retain liquid water throughout winter.
For both clusters C and D, the mean backscatter from the lake basins is barely distinguishable from that of the surrounding ice (Fig. 5, 2.1). Although the transects exhibit substantial fluctuations both within and outside the basins, there is no systematic difference between the signal from within and outside the lake basin (Fig. 5, 2.2). The backscatter trends over winter differ between clusters, with generally positive slopes in cluster C and negative slopes in cluster D (Fig. 4), likely reflecting minor variations in surface or near-surface ice properties. However, as there is no consistent offset between the lake and surrounding ice, we interpret these basins as empty. Therefore, both clusters C and D represent the same winter lake type (No lake) despite being separated into two clusters due to slightly different temporal backscatter patterns. We suggest reasons for the separation into two clusters in the Discussion.
Validation of clustering-based winter lake classification
Each cluster A–D consistently matched one combination of optical labels from the preceding and following summer across both years and both sites (Table 1). For the majority of lakes in clusters A and B, the optical label from the preceding summer was ‘refreeze’, i.e., lakes covered by an ice and snow lid at the end of summer (Fig. 5). In the following summer, lakes in cluster A were largely classified as ‘refill’ as they appeared empty post-winter and refilled later in summer. Of all lakes in cluster A, 57.8% were labelled ‘refreeze’ & ‘refill’, corroborating their classification as lakes that completely freeze-through (Figs S7 and S8). Lakes in cluster B were primarily labelled ‘liquid’ in the following summer and appear as already water-filled beneath an ice lid after retaining water through winter (Fig. 5). Of all lakes in cluster B, 62.8% were labelled ‘refreeze’ & ‘liquid’ (Figs S7 and S8), supporting their interpretation as buried lakes. A total of 86.8% of all lakes in clusters C and D were labelled ‘drained’, as empty basins at the end of the preceding summer, and ‘refill’ in the following summer (Figs S7 and S8), supporting their interpretation as empty lake basins with no winter activity (Fig. 5).
For two lakes at Russell Glacier and nine lakes at NEGIS, all of which were labelled as ‘liquid’ post-winter in the optical data, radar echograms confirmed a subsurface horizon interpreted as liquid water buried underneath an ice lid, e.g., Fig. 6. For 11 lakes at Russell Glacier and 18 lakes at NEGIS, all of which were labelled ‘refill’ post-winter in the optical data, no such horizon was found in corresponding radar echograms, suggesting no subsurface water here, in agreement with the optical labels. Agreement between optical labels and radar echograms therefore provides independent support for the manual labelling of buried lake presence derived from optical data (Table 1).
(a) Operation IceBridge radio echogram from 14 May 2019 (Paden and others, Reference Paden, Li, Leuschen, Rodriguez-Morales and Hale2014) over the example lake from cluster B in Figure 5, with visible subsurface reflectance signature. (b) Copernicus Sentinel-2 image of the example lake from 13 September 2018. The red line indicates the profile across the example lake in Panels A and B.
Figure 6 Long description
The figure pairs an Operation IceBridge radio echogram with an optical image of the same buried supraglacial lake. Panel a shows a five-kilometre flight profile across an undulating ice surface. Between approximately 3.3 and 4.7 kilometres along the profile, the surface forms a shallow depression and a dark, near-horizontal subsurface reflection occurs several metres below it, indicating subsurface attenuation associated with liquid water. A red horizontal line with a diamond at its left end and a circle at its right end marks the section crossing the lake. Panel b shows the corresponding elongated, ring-shaped lake in a Sentinel-2 image. The same red transect crosses approximately east-west through the upper part of the lake, with north towards the top of the image.
Clustering performance
The performance of the time series clustering was evaluated by assessing how many lakes in each cluster correctly matched the pre- and post-winter optical labels (Table 1). Overall classification accuracy was 75.0% in 2018 and 74.0% in 2019 for Russell Glacier, and 79.4% in 2018 and 70.8% in 2019 for NEGIS (Fig. 7). Across both years and both sites, around half the misclassified lakes (12.4%) were either empty lake basins that should have been winter lakes (freeze-through or buried), or vice versa, i.e., winter lake presence/absence was incorrectly identified (Figs S9–S13). The other half of misclassifications (12.6%) were winter lakes that were correctly distinguished from empty lake basins, but were incorrectly classified as freeze-through instead of buried, or vice versa. The clustering approach therefore achieved an overall accuracy of 87.6% in separating winter lakes from empty lake basins (up to 88.5% at NEGIS in 2018), whereas further subdivision of winter lake behaviours into freeze-through and buried types reduced the overall accuracy to 75.0% (Fig. S13).
Cluster performance based on comparison between cluster interpretations and optical lake labels. Categories are freeze-through lakes (cluster A), buried lakes (cluster B), empty lake basins (clusters C and D) and misclassified lakes (from any cluster). Labels contain the percentage of lakes and the count in brackets.
Figure 7 Long description
A 100-percent stacked bar chart compares classification outcomes for Russell Glacier and NEGIS in 2018 and 2019, plus the combined results. For Russell 2018, 36.8 percent, or 28 lakes, are no winter lake; 22.4 percent, or 17, are freeze-through; 15.8 percent, or 12, are buried; and 25.0 percent, or 19, are misclassified. For Russell 2019, the corresponding values are 50.6 percent, or 39; 11.7 percent, or 9; 11.7 percent, or 9; and 26.0 percent, or 20. For NEGIS 2018, they are 64.7 percent, or 277; 5.4 percent, or 23; 9.3 percent, or 40; and 20.6 percent, or 88. For NEGIS 2019, they are 35.0 percent, or 150; 18.0 percent, or 77; 17.8 percent, or 76; and 29.2 percent, or 125. Overall, 48.9 percent, or 494 lakes, are no lake; 12.5 percent, or 126, are freeze-through; 13.6 percent, or 137, are buried; and 25.0 percent, or 252, are misclassified. No lake is the largest category in every dataset, and misclassification ranges from 20.6 to 29.2 percent.
Investigation of misclassified lakes
Of all lakes analysed across both years and both sites, a total of 25.0% were misclassified. All misclassified lakes were further inspected based on their backscatter transects across lake basins. While an explanation for the misclassification of 10.7% of lakes could not be clearly determined and remains unknown, 89.3% of the misclassifications can be explained by four types of error sources (Tables S1 and S2).
Mixed-type backscatter signal accounts for 19.1% of misclassified lakes. These lakes exhibit multiple behaviours within a single basin, with some parts freezing through (Fig. 8b) and other parts remaining liquid (Fig. 8a and c). The clustering operates on the basis of a single time series averaged across each lake basin; therefore, such basins are assigned to a single class, even when multiple physical processes are occurring simultaneously at sub-basin scale. As the mean time series is derived by averaging backscatter across all sub-basin pixels, behaviour covering the largest area of the lake basin exerts the greatest influence on the resulting signal.
Change in backscatter signal over time across selected lake transects, representing examples of the four misclassification categories. Lake basin outlines are marked by vertical black dashed lines, with surrounding ice backscatter shown for an extra 1 km at either end of the transect.
Figure 8 Long description
Four vertically arranged transect plots show mixed, masked, ambiguous, and fluctuating misclassification categories. Each shows backscatter in decibels along a lake transect. Curves progress from purple in September 2018 to yellow in April 2019. Dashed lines mark lake-basin boundaries, with one kilometre of surrounding ice at each end. In the mixed example, sections show contrasting temporal patterns: regions a and c have low early-winter backscatter that gradually increases to surrounding ice levels, while region b increases more quickly. In the masked example, adjacent regions d and e show contrasting high and low backscatter compared to surrounding ice. In the ambiguous example, region f shows low early-winter backscatter followed by an increase, although it is unclear whether it reaches surrounding ice values. In the fluctuating example, backscatter varies sharply and irregularly within the basin, with repeated peaks and troughs but no coherent temporal trend.
Masked backscatter signal causes 12.7% of misclassified lakes. Spatial averaging of backscatter across the basin masked local variability, producing a mean signal characteristic of an empty lake. Lake transects reveal that these basins contain areas of strongly increased backscatter, likely driven by surface roughness (Fig. 8d), alongside areas with signatures of refreezing or buried water (Fig. 8e).
Ambiguous backscatter signal explains 21.8% of misclassified lakes. Their temporal evolution of backscatter is not sufficiently distinct across the lake basin to clearly distinguish between buried and freeze-through behaviour (Fig. 8f).
Fluctuating backscatter signal characterises the remaining 35.7% of misclassified lakes, making them unclassifiable by either the clustering approach or manual inspection (Fig. 8, bottom panel).
Spatio-temporal comparison of lake types
Based on the results from our combined lake drainage extraction and clustering approach, we can compare the spatial distribution of the different winter lake types between both sites and both years (Figs 9 and 10). At Russell Glacier, a larger proportion of lakes persisted into winter in 2018 (winter lakes = 39.0%, Table 2), following a relatively cold melt season, compared to 2019 (winter lakes = 23.4%, Table 2), following a relatively warm melt season. In both years, the majority of these ‘winter lakes’, which include drained, freeze-through and buried types, are concentrated at higher elevations, with comparatively fewer at lower elevations (Fig. 9).
Spatial distribution of different winter lake types at Russell Glacier during winters of 2018 and 2019, derived from our combined winter drainage extraction and time series clustering approach (Table 2). Black dashed lines represent 100 m ice surface elevation isolines (NASA MEaSUREs GIMP, (Howat and others, Reference Howat, Negrete and Smith2015)).
Figure 9 Long description
Two side-by-side maps show winter lake classifications across Russell Glacier in 2018 and 2019. North is upwards, and dashed 100-metre elevation contours generally increase from about 1,000 metres in the west, on the left margin of the map extent, to 1,800 metres and higher towards the east, on the right margin of the map extent. Lakes are classified as buried, freeze-through, drainage, no lake, or misclassified. The classifications are spatially intermixed, although buried and freeze-through lakes occur more frequently in the higher-elevation part of the mapped region. Misclassified lakes are scattered throughout, and drainage events are rare with only 1 in 2018 near the spatial cluster of winter lakes. Compared with 2018, the 2019 map contains more no-lake classifications and fewer freeze-through and buried lakes. Winter lakes can be seen only above 1,400 metres, clustering at higher elevations.
Spatial distribution of different winter lake types at NEGIS during winters of 2018 and 2019, derived from our combined winter drainage extraction and time series clustering approach (Table 2). Black dashed lines represent 100 m ice surface elevation isolines (NASA MEaSUREs GIMP, (Howat and others, Reference Howat, Negrete and Smith2015)).
Figure 10 Long description
Two side-by-side maps show winter lake classifications across NEGIS in 2018 and 2019. North is upwards, and dashed 100-metre elevation contours decrease from approximately 1,200 metres in the west, on the left, to 300 metres towards the lower-elevation coastal margin in the east, on the right. Lakes are classified as buried, freeze-through, drainage, no lake, or misclassified. In 2018, no-lake classifications are more widespread, particularly across the western and southern parts of the lake region, while buried and freeze-through lakes are scattered mainly across the central and eastern areas. In 2019, buried lakes are more numerous and occur throughout the region, including a prominent group of elongated lakes in the south. Freeze-through lakes are also more numerous than in 2018. Misclassified lakes remain scattered across both maps, and drainage events are rare.
Count and percentage of lakes by winter lake type at Russell Glacier and NEGIS during winters of 2018 and 2019, derived from our combined winter drainage extraction and time series clustering approach.
Table 2 Long description
The table reports winter lake classifications for 77 lakes at Russell Glacier in each year and 430 lakes at NEGIS in each year. At Russell Glacier in 2018, 28 lakes, or 36.3 percent, were no lake; 1, or 1.3 percent, drainage; 17, or 22.1 percent, freeze-through; 12, or 15.6 percent, buried; and 19, or 24.7 percent, misclassified. In 2019, 39 lakes, or 50.6 percent, were no lake; none were drainage; 9, or 11.7 percent, freeze-through; 9, or 11.7 percent, buried; and 20, or 26.0 percent, misclassified. At NEGIS in 2018, 277 lakes, or 64.4 percent, were no lake; 2, or 0.5 percent, drainage; 23, or 5.3 percent, freeze-through; 40, or 9.3 percent, buried; and 88, or 20.5 percent, misclassified. In 2019, 150 lakes, or 34.9 percent, were no lake; 2, or 0.5 percent, drainage; 77, or 17.9 percent, freeze-through; 76, or 17.7 percent, buried; and 125, or 29.0 percent, misclassified.
In contrast, fewer winter lakes were detected at NEGIS in 2018 (15.1%) than in 2019 (36.1%; Table 2). Their spatial distribution differs between years, with winter lakes concentrated at lower elevations in 2018 but distributed more broadly across the study area in 2019 (Fig. 10).
Discussion
Lake behaviour classification from C-band backscatter time series
Our combined method of lake drainage extraction and unsupervised time series clustering from C-band backscatter was successful at identifying lake behaviours over winter, with 75.0% of all analysed lakes correctly classified across both sites and both years (Fig. 7). Winter lakes were distinguished from empty lake basins that had drained during the summer with 87.6% accuracy, and could be further partitioned into fully refreezing and buried lakes (Fig. S13). The temporal backscatter signatures of different winter lake behaviours were reproduced across both study sites (Fig. 4) and are in agreement with previous studies for Greenland (Miles and others, Reference Miles, Willis, Benedek, Williamson and Tedesco2017; Benedek and Willis, Reference Benedek and Willis2021; Dunmire and others, Reference Dunmire2025) and Antarctica (Li and others, Reference Li, Chen, Lin, Hooper, Chen and Cheng2023; Deakin and others, Reference Deakin, Willis and Dell2025). This demonstrates the potential for using temporal C-band SAR signatures to classify lakes during winter, complementing optical observations acquired in summer and enabling the development of year-round lake inventories.
Furthermore, our winter lake classification approach demonstrates two clear advantages over previous studies in accurately detecting and classifying winter lake behaviours. First, it does not rely on any pre-defined training data (Dunmire and others, Reference Dunmire2025). Backscatter signatures from C-band SAR provide a useful basis for distinguishing winter lake behaviour with an unsupervised k-means clustering approach, successfully capturing the majority of cases despite some misclassification. Second, our clustering approach uses the complete winter time series of backscatter for lake classification, therefore capturing their full evolution over the course of winter. As has been previously found (e.g., Miles and others, Reference Miles, Willis, Benedek, Williamson and Tedesco2017; Benedek and Willis, Reference Benedek and Willis2021; Deakin and others, Reference Deakin, Willis and Dell2025), our observations show high temporal variability in lake evolution throughout winter, with: i) drainage events occurring infrequently and sporadically, and ii) duration and timing of lake freeze-through varying widely. However, previous winter lake mapping approaches have classified winter lakes based on single-date SAR acquisitions or temporally limited SAR mosaics around mid-winter (Dunmire and others, Reference Dunmire, Banwell, Wever, Lenaerts and Datta2021; Zheng and others, Reference Zheng2023), therefore capturing only a snapshot of winter lake evolution. As a result, winter lake behaviour may be misclassified or overlooked. For example, a lake identified as buried may still undergo freeze-through after the selected SAR acquisition date, while drainage or complete freeze-through occurring prior to acquisition may go undetected. Using the full winter time series preserves temporal information, enabling more accurate classification of winter lake behaviour.
Spatio-temporal variability and glaciological implications of winter lake behaviour
Our approach successfully captures spatio-temporal variability in winter lake formation, distribution and behaviour across sites and years (Figs 9 and 10). Our findings agree with previous studies showing that: (a) winter lakes generally occur at higher elevations in Southwest Greenland than in Northeast Greenland (Zheng and others, Reference Zheng2023; Dunmire and others, Reference Dunmire2025); (b) in the northeast, buried lakes formed more frequently and at significantly higher elevations in 2019 compared to 2018 (Dunmire and others, Reference Dunmire, Banwell, Wever, Lenaerts and Datta2021; Zheng and others, Reference Zheng2023); and (c) at Russell Glacier, more supraglacial lakes refroze to become winter lakes in 2018 than in 2019 (Glen and others, Reference Glen2025). This pronounced inter-annual and spatial variability has been attributed to differences in summer and autumn air temperatures, snow cover extent and firn properties, which together create conditions that either promote or inhibit winter lake persistence (Dunmire and others, Reference Dunmire, Banwell, Wever, Lenaerts and Datta2021; Reference Dunmire2025; Zheng and others, Reference Zheng2023).
The ability to reliably locate these winter lakes and identify their winter behaviour is crucial for assessing their impact within the ice sheet system (Fig. 1). Our results reveal that up to 
$\sim$ 50% of all lakes can persist into the winter season (Figs S9 and S12), with between 5.3% (NEGIS, 2018) and 22.1% (Russell Glacier, 2018) of all lakes identified to freeze through fully over the course of a winter (Table 2). In-situ measurements on the Larsen C Ice Shelf, Antarctic Peninsula, revealed that the refreezing of ponded meltwater led to the formation of subsurface ice layers with increased density and reduced permeability compared to the surrounding firn, and produced englacial heating upon refreezing (Hubbard and others, Reference Hubbard2016). Similarly, ice lenses—formed by locally refrozen meltwater that has percolated through the firn—can connect to larger ice layers and contribute to the expansion and thickening of ice slabs observed around the entire margin of the Greenland Ice Sheet, leading to significantly reduced local permeability, increased firn densification and resulting in increased runoff potential (e.g., Machguth and others, Reference Machguth2016; Jullien and others, Reference Jullien, Tedstone, Machguth, Karlsson and Helm2023; Culberg and others, Reference Culberg, Michaelides and Miller2024; Tedstone and others, Reference Tedstone2025).
Furthermore, we identified that between 9.3 and 17.7% (NEGIS 2018 and 2019, respectively) of lakes do not freeze through completely over a winter but retain buried water (Table 2). These buried lakes represent a form of meltwater storage and delayed release, which can occur either during winter or in the following melt season by reducing the lake’s capacity, leading to overflow and formation of meltwater channels (Koenig and others, Reference Koenig2015; Lampkin and others, Reference Lampkin, Koenig, Joseph and Box2020). We detected only a small number of drainage events in winter (Table 2), and previous research suggests such isolated winter drainage events have only limited impact on ice dynamics (Schröder and others, Reference Schröder, Neckel, Zindler and Humbert2020b; Benedek and Willis, Reference Benedek and Willis2021). However, recent studies observed cascading drainage events of a cluster of winter lakes, which resulted in accelerated ice flow over large areas of the basin and modifications to the drainage system (Maier and others, Reference Maier, Andersen, Mouginot, Gimbert and Gagliardini2023; Dean and others, Reference Dean, Scharien, Willis and McDougall2026). Thus, robust identification of buried lake locations is a prerequisite for reliably projecting where clustered winter drainage events may occur and potentially affect ice-flow dynamics.
Sub-basin complexity of winter lake behaviours
Although averaging backscatter across entire lake basins offers a practical means of reducing two-dimensional, pixel-scale SAR data and is widely adopted in lake studies (e.g., Miles and others, Reference Miles, Willis, Benedek, Williamson and Tedesco2017; Benedek and Willis, Reference Benedek and Willis2021; Dunmire and others, Reference Dunmire2025), it inevitably obscures spatial variability within individual basins. Our investigation of misclassified lake behaviours shows that three out of four identified misclassification sources (panels ‘Mixed’, ‘Masked’, ‘Ambiguous’; Fig. 8), comprising 53.6% of all misclassified lake behaviours (Tables S1 and S2), arise from the inability of the basin-averaged time series to capture the sub-basin spatial and temporal complexity of the backscatter signal.
The backscatter signal covering the largest proportion of a lake basin is most strongly represented in the mean time series due to its weighting during spatial averaging of pixel-scale values. This mean backscatter time series is then assigned to a single cluster representing a specific winter lake behaviour. Consequently, the classification may not fully capture lake dynamics where multiple physical processes occur simultaneously within the lake basin. Although misclassification of lake behaviour was largely attributed to non-uniform backscatter, this limitation applies more broadly to all lakes exhibiting spatial variability in backscatter across their basins.
Our results reveal the non-uniform temporal evolution of winter lakes across their basin and highlight the spatial complexity of winter lake behaviour on a sub-basin scale. Li and others (Reference Li, Lhermitte and López-Dekker2021) identified mixed-type lake behaviour on Amery Ice Shelf, Antarctica, using InSAR-derived coherence, indicating spatial heterogeneity within individual lake basins, with signals consistent with both refreezing and drainage processes. Based on our findings, mixed-type lake behaviour might be a widespread occurrence in winter. Our observations of spatially variable backscatter patterns across a lake basin (Fig. 8a–c and f), interpreted as representing different degrees of refreezing, raise important questions regarding the factors determining which parts of a lake basin refreeze more quickly and fully as opposed to those that retain liquid water due to, e.g., lake depth or ice lid thickness. Future work should investigate more sophisticated methods that incorporate both temporal and spatial dimensions of backscatter variability (i.e., functional or pixel-based time series analysis) to fully capture the complexity of winter lake behaviours and accurately assess their glaciological implications.
Method and C-band data limitations
While the time series clustering approach presented here offers new insight into winter lake monitoring, several methodological limitations must be acknowledged. First, like previous studies, our approach relies on pre-existing lake outlines, e.g., derived from SAR or summer optical imagery (Miles and others, Reference Miles, Willis, Benedek, Williamson and Tedesco2017; Benedek and Willis, Reference Benedek and Willis2021; Zheng and others, Reference Zheng2023). These datasets provide a simple entry point to year-round lake monitoring, but the optical imagery in particular may miss small or fully buried winter lakes that do not present as open water during the melt season. Future studies should derive lake outlines from SAR imagery on a larger scale to address this limitation (Schröder and others, Reference Schröder, Neckel, Zindler and Humbert2020b).
Second, the choice of k-means clustering with a fixed number of clusters and a selected Euclidean distance metric imposes methodological constraints. Pre-selection of the cluster number determines centroid positions and directly impacts the separation of time series data into clusters (Aghabozorgi and others, Reference Aghabozorgi, Seyed Shirkhorshidi and Ying Wah2015; Alqahtani and others, Reference Alqahtani, Ali, Xie and Jones2021). Moreover, k-means struggles to handle noise or outliers in data (Aghabozorgi and others, Reference Aghabozorgi, Seyed Shirkhorshidi and Ying Wah2015; Alqahtani and others, Reference Alqahtani, Ali, Xie and Jones2021). Winter drainage events, which manifest as abrupt, isolated jumps in backscatter, therefore need to be treated separately rather than captured within the clustering framework (Alqahtani and others, Reference Alqahtani, Ali, Xie and Jones2021). Additionally, the winter drainage detection method applied in this study may underestimate the number of winter drainage events due to conservative thresholds used for filtering time series (Benedek and Willis, Reference Benedek and Willis2021). Although the four drainage events identified at NEGIS are consistent with the findings of Dean and others (Reference Dean, Scharien, Willis and McDougall2026), they reported several additional drainage events in both years, possibly due to their use of higher-temporal-resolution time series.
Third, there are uncertainties associated with the supporting datasets used for validation. ArcticDEMs and Operation IceBridge radargrams offer limited spatial and temporal coverage, providing independent verification of winter drainage and buried lakes, respectively. Optical labels derived manually from Sentinel-2 imagery are subjective and affected by cloud cover and coarse temporal resolution. Determining whether a lake basin contains liquid water beneath an ice lid—particularly in the post-winter period when assessing the persistence of buried lakes—is challenging using optical imagery alone, as cloud cover can introduce multi-day data gaps and ice lids may obscure subsurface conditions.
Fourth, interpretation of SAR backscatter remains challenging because the signal is influenced not only by dielectric properties related to water content but also by geometric properties of the remotely sensed target surface (Fahnestock and others, Reference Fahnestock, Bindschadler, Kwok and Jezek1993; Rott and Davis, Reference Rott and Davis1993; Li and others, Reference Li, Lhermitte and López-Dekker2021). In this study of winter lakes, backscatter is predominantly controlled by dielectric contrasts between ice, snow and water. However, several observed features could not be explained by dielectric effects alone. These include strongly fluctuating backscatter signal (Fig. 8, bottom panel) and localised areas of increased backscatter (Fig. 8d) within a subset of misclassified lake basins, as well as minor but consistent differences between the backscatter time series of clusters C and D (Fig. 4) even though both correspond to empty lake basins. We suggest that these complex backscatter responses are caused by geometric effects, such as variations in surface roughness or internal ice structure (Rott and Davis, Reference Rott and Davis1993; Li and others, Reference Li, Lhermitte and López-Dekker2021), for which no supporting ground-truth data are currently available. Additionally, the incidence angle affects backscatter signal across an image that may interfere with classification. While using imagery from a single orbit can help minimise incidence angle-related effects, future work could incorporate normalisation methods and thereby improve temporal resolution by including all available SAR imagery from different orbits (Dean and others, Reference Dean, Scharien, Willis and McDougall2026).
Finally, unknown and spatially variable C-band penetration depth (between 
$\sim$1 and 2 m on bare ice and up to 
$\sim$10 m on dry snow/firn) imposes physical limits on the detectability of buried lakes, which have been found at depths of 
$\sim$10 m (Rignot and others, Reference Rignot, Gogineni, Joughin and Krabill2001; Koenig and others, Reference Koenig2015). Lakes classified as freeze-through may, in some cases, retain water below the C-band penetration depth, causing deep buried lakes to be missed or misinterpreted. The use of lower-frequency radar such as L-band, which can penetrate deeper into ice surfaces, could help overcome this limitation. L-band SAR sensors can penetrate 
$\sim$5–10 m deeper into firn and ice than C-band systems (Rignot and others, Reference Rignot, Gogineni, Joughin and Krabill2001), offering enhanced sensitivity to subsurface meltwater and refreezing processes. Upcoming missions such as ESA’s ROSE-L and NASA–ISRO’s NISAR, alongside existing ALOS L-band data, provide promising opportunities to evaluate the complementary information of different SAR wavelengths (Li and others, Reference Li, Lhermitte and López-Dekker2021; Shukla and others, Reference Shukla2024).
Despite these limitations, the strong agreement between clustered winter lake types and independent optical labels and radar data, and the internal consistency of the derived temporal signatures across sites and years, indicate that the approach is robust for identifying dominant winter lake behaviours. The growing Sentinel-1 archive, now a decade in length with a 6–12 day repeat cycle, provides extended multi-year time series data to examine inter-annual variability and potential trends in winter lake persistence, refreezing and drainage frequency across the ice sheet. Ultimately, future research should move towards a year-round classification of supraglacial lake behaviours that integrates both summer and winter observations. Quantifying the volume of drained and refrozen water in summer, together with the winter distribution and storage of buried and freeze-through lakes, would allow assessment of the net hydrological contribution of supraglacial lakes to basin-scale and ice-sheet-wide dynamics. Such an integrated framework would represent a major step towards understanding the role of supraglacial lakes in the evolving Greenland Ice Sheet hydrological system.
Conclusion
This study demonstrates that C-band backscatter time series from Sentinel-1 SAR provides a robust basis for identifying and characterising winter supraglacial lake behaviour on the Greenland Ice Sheet. Unsupervised time series clustering achieves 87.6% accuracy in distinguishing winter lakes from empty lake basins, and further captures distinct temporal signatures associated with winter lake behaviour of freeze-through and buried lakes (75.0% accuracy) across differing study sites and contrasting climatic years. Our assessment demonstrates that temporal evolution of C-band backscatter contains the critical information needed to classify winter lake processes accurately without training data. We conclude the following key findings:
(1) Our results demonstrate that winter lake behaviours can be reliably identified from temporal C-band backscatter without training data. Characteristic time series signatures directly reflect physical processes within lake basins that freeze, drain and store subsurface water.
(2) The characteristic signatures of the C-band time series lake behaviours are consistent across contrasting climatic regions and between climatically contrasting years, indicating robustness for spatial and temporal transferability of our method.
(3) Substantial sub-basin variability shows that winter lake behaviour is not uniform within individual basins. Mixed-type lakes highlight the need for further development that integrates spatial as well as temporal information to fully capture the sub-basin complexity in winter lake behaviours.
Together, these results provide new insights into how C-band backscatter can be utilised to complement existing supraglacial lake studies and contribute towards continuous, year-round monitoring of surface hydrology across the Greenland Ice Sheet. Extending such approaches across longer time periods and incorporating multi-frequency radar data will further enhance understanding of evolving winter hydrological processes and their impact on ice sheet dynamics under a warming climate.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2026.10191.
Data availability statement
The data presented in this study and code required to reproduce the time series clustering workflow can be accessed here: https://doi.org/10.5281/zenodo.18786243. Lake outline polygons for NEGIS were accessed through Pangaea Data Center: https://doi.pangaea.de/10.1594/PANGAEA.918707 (Schröder and others, Reference Schröder, Neckel and Humbert2020a; Reference Schröder, Neckel, Zindler and Humbert2020b). Lake outline polygons for Russell Glacier were accessed on Zenodo: https://doi.org/10.5281/zenodo.16577959 (Glen and others, Reference Glen2025; Glen, Reference Glen2025). Sentinel-1 RTC products were obtained from the Alaska Satellite Facility’s Vertex data portal: https://search.asf.alaska.edu/#/. Sentinel-2 data were accessed through the Copernicus Browser: https://browser.dataspace.copernicus.eu/. Arctic Digital elevation models were obtained from the Polar Geospatial Center:https://www.pgc.umn.edu/data/arcticdem/. MEaSUREs Greenland Ice Mapping Project (GrIMP) Digital Elevation Model data (Howat and others, Reference Howat, Negrete and Smith2015) and Operation IceBridge radar echograms (Paden and others, Reference Paden, Li, Leuschen, Rodriguez-Morales and Hale2014) were accessed through the NASA National Snow and Ice Data Center: https://nsidc.org/data.
Acknowledgements
This work was supported by the Natural Environment Research Council (NERC grant reference number NE/S007423/1), with J.O.’s studentship through the ENVISION Doctoral Training Partnership. M.M. was supported by the UK NERC Centre for Polar Observation and Modelling (NE/Y006178/1), and the UKRI Frontier Research grant GLOBE: The Greenland Subglacial Lake Observatory (EP/Y02642X/1). H.M. was supported by Research England under the Expanding Excellence in England (E3) funding stream, awarded to MARS: Mathematics for AI in Real-world Systems in the School of Mathematical Sciences at Lancaster University.
The ArcticDEM data were accessed through the Polar Geospatial Center, University of Minnesota. DEMs provided by the Polar Geospatial Center under NSF-OPP awards 1043681, 1559691, 1542736, 1810976 and 2129685. Geospatial support for this work provided by the Polar Geospatial Center under NSF-OPP awards 1043681, 1559691 and 2129685.
We would like to thank Ian Willis, an anonymous reviewer and the Scientific Editor Hester Jiskoot for their detailed and valuable comments.
Author contributions
J.O., K.E.M. and A.A.L. conceptualised the research. J.O. carried out the main body of work, performed all the data analysis and wrote the manuscript. K.E.M., A.A.L., J.M. and M.M. contributed to the scientific content, technical details and overall structure of this paper. H.M. provided technical input on the time series clustering workflow. All co-authors contributed to the discussion of results and editing of the manuscript.
Open access
