Time Series Classification of Supraglacial Lakes Evolution over Greenland Ice Sheet

Understanding the hydrologic processes of the GrIS and their implications for its mass balance has been a subject of significant research interest in recent years. This section provides a comprehensive overview of the related works in this domain, encompassing studies on the formation and behavior of supraglacial lakes, their impact on ice dynamics, and the methods used for their detection and monitoring.

Previous studies, such as those by [20] and [21], investigated the evolution of supraglacial lakes and their dependence on topographic features, highlighting the influence of bed topography on their formation and persistence. [22] further explored the seasonal variability of supraglacial lakes, emphasizing their significance in understanding ice sheet dynamics. The formation and drainage of supraglacial lakes have profound implications for ice dynamics and mass loss from the GrIS. [3] demonstrated a shift in the predominant cause of mass loss from ice discharge to meltwater runoff, underscoring the increasing significance of surface melt in the ice sheet’s mass balance. [6] investigated the role of refrozen meltwater in altering firn properties and influencing future meltwater percolation, highlighting the complex feedback mechanisms between surface melt and ice dynamics. Furthermore, studies by [4] and [17] examined supraglacial lakes’ seasonal variability and spatial distribution, providing insights into the factors influencing their formation and evolution. [23] and [24] investigated the role of supraglacial lake drainage in modulating ice sheet dynamics, highlighting the potential for abrupt changes in ice flow and mass loss.

Satellite-based and remote sensing techniques have emerged as valuable tools for detecting and monitoring supraglacial lakes on the GrIS. [25] and [26] utilized multispectral satellite imagery, including data from MODIS (Moderate Resolution Imaging Spectroradiometer) and Landsat satellites, to identify supraglacial lakes and channels, enhancing the spatial coverage of surface meltwater features. [27] explored the use of Sentinel-2 imagery for automated detection of supraglacial lakes, highlighting the potential for improving temporal resolution and monitoring capabilities. In addition to optical imagery, SAR imagery has been increasingly employed to detect supraglacial lakes and buried lakes on the GrIS. [26] demonstrated the utility of SAR for detecting supraglacial lakes, while [18] utilized SAR imagery to identify perennial supraglacial lakes. [17] and [7] further investigated the use of SAR for detecting and monitoring buried lakes, highlighting its capability to penetrate snow and ice.

There are several established machine learning (ML) and deep learning (DL) models widely used for time series classification tasks. Long Short-Term Memory (LSTM) networks and Fully Convolutional Networks (FCN) are among the prominent models used in this domain. Residual Networks (ResNet) and Recurrent Neural Networks (RNNs) are also integral to this field, offering unique time series analysis capabilities. The following defines some popular classification techniques that are selected based on their proven effectiveness and unique methodologies for handling time series data.

LSTMFCNClassifier: Combines LSTM and FCN to leverage both sequential dependencies and local feature extraction.

FCNClassifier: Utilizes FCN to capture local patterns in time series data without fully connected layers, allowing it to handle input sequences of varying lengths.

ResNetClassifier: Employs ResNet with residual connections to facilitate the training of deep networks, effectively capturing complex patterns in the data.

SimpleRNNClassifier: Uses RNNs with recurrent connections to capture dependencies over time, making it suitable for time series classification.

KNeighborsTimeSeriesClassifier: A non-parametric, instance-based learning algorithm that classifies data points based on the majority class among their $k$-nearest neighbors in the feature space.

The seasonal evolution of supraglacial lakes can be viewed as a dynamic system where their state changes over time due to various environmental factors. The Reconstructed Phase Space (RPS) is particularly suitable for this research as it enables the capture of underlying dynamics from time series data, such as remote satellite observations. A dynamical system describes the temporal evolution of a system to capture its dynamics. A phase space represents all possible states of the system that evolve, and the dynamics map describes how the system evolves. The RPS captures the underlying dynamics of a system from time series observations, allowing for the analysis of complex and non-linear behaviors in supraglacial lakes.

According to Takens’ embedding theorem, a time series $x=\{x_{n}\},n=1,\ldots,N$ can be converted into state vectors using time-delay embedding:

where $\tau$ is the time delay and $d$ is the embedding dimension [28]. This embedding reconstructs the state and dynamics of the unknown system from observed measurements. The dimension $d$ should be greater than twice the box-counting dimension of the original system [29]. If $d$ is unknown, it can be estimated using the false nearest-neighbor techniques, and $\tau$ can be determined by finding the first minimum of the automutual information [30].

Gaussian Mixture Models (GMMs) are employed to model the distribution of dynamics represented by the RPS. The seasonal evolution of supraglacial lakes involves complex, non-linear behaviors that are well-captured by GMMs, which can model multiple underlying distributions within the data. A GMM is a weighted sum of $M$ Gaussian distributions:

where $\lambda=\{w_{i},\mu_{i},\Sigma_{i}\}$ are the weights, means, and covariance matrices of the mixture components [31]. These parameters of the GMM are estimated using the Expectation-Maximization (EM) algorithm, which iteratively maximizes the likelihood of the data [32]. The EM algorithm includes two main steps:

1. Expectation (E-step): Calculate the responsibilities for each data point $X_{i}$:

2. Maximization (M-step): Update the parameters using the responsibilities:

The E-step and M-step are repeated until the parameters converge, ensuring the GMM optimally fits the data.

A Bayesian maximum likelihood classifier is used to classify the test data. For each test point $X_{k}$, the likelihoods are computed for each model $a_{i}$:

where $X=\{x_{1},x_{2},\ldots,x_{T}\}$ is the sequence of observations. By leveraging the likelihoods derived from the GMM, the classifier can effectively distinguish between the different types of lakes. After computing all the likelihoods, the class with the maximum likelihood, $\hat{a}$ (i.e., predicted class), is determined using the following equation [33]:

We observe distinct evolutionary changes in supraglacial lakes from satellite data, specifically using the Sentinel-1 (S1, microwave) and Sentinel-2 (S2, optical) satellites. We classify lakes into three distinct categories: refreezing lakes, draining lakes, and buried lakes. From the S1 microwave imagery, we use the horizontally-transmitted, vertically-received (HV) band, previously used for buried lake detection [7, 17]. For each S1 image, we calculate the average HV value both within the lake outline and in the immediate vicinity outside the lake bounds, within a 750$\sim$m buffer of the lake. Figure 2 illustrates an average backscatter timeseries within the lake bounds (blue line, $HV_{lake}$) and from the area surrounding the lake (purple line, $HV_{background}$). By differencing these two backscatter signals using Equation 9, we can identify the backscatter anomaly for a given lake, denoted as $HV_{anom}$.

We also obtain a time series with information from optical imagery for each lake. From each S2 optical image, we determine the water percentage ($p_{water}$) inside the lake by calculating the ratio of water-identified pixels to the total number of pixels within the lake. This calculation is detailed in Equation 10. $N_{water}$ represents the number of pixels identified as water in the S2 image, and $N_{total}$ is the total number of pixels within the lake. Figure 3 illustrates an example time series for all three lake types.

Refreeze lakes form when temperatures decrease towards the end of the melt season, causing the lakes to transition from liquid to frozen. Figure 3(a) illustrates that both $HV_{anom}$ and $p_{water}$ decline at the end of the summer and approach zero in the fall, indicating the lake water has fully refrozen.

The optical timeseries for a buried lake appears similar to that of a refreeze lake, as $p_{water}$ decreases to zero at the end of the melt season. However, microwave images and the backscatter $HV_{anom}$ timeseries in Figure 3(b) indicate that some liquid water remains buried beneath the surface, even at the end of the year. Water strongly absorbs microwave radiation and thus areas with liquid water presence appear relatively dark in the microwave imagery.

Drained lakes undergo a significant reduction in water volume, typically marked by a sharp decline in $p_{water}$. We classify lakes as draining if the lake water is drained to the ice bed and $p_{water}$ approaches zero during the melt season. Figure 3(c) depicts the optical and microwave timeseries of a lake experiencing rapid drainage, with imagery indicating that the lake was drained over the four days between June 01-05.

For this study, we utilize a comprehensive pan-Greenland dataset [7, 34], which includes detailed outlines of supraglacial lakes with a resolution of 30 meters. This dataset consists of 3,846 supraglacial lakes during the 2018 melt season and 6,146 supraglacial lakes during 2019, each with a surface area $>$ 0.05 $km^{2}$. The years 2018 and 2019 are chosen to capture a range of climatic conditions, with 2018 representing a cooler melt season and 2019 representing a warmer melt season [35, 36], providing a comprehensive understanding of supraglacial lake dynamics under different temperature regimes. For this work, we manually label the time series of 777 lakes into three classes: refreeze (189 lakes), drained (392 lakes), and buried (196 lakes). The dataset spans all six subregions of the Greenland Ice Sheet: Northeast (NE), Northwest (NW), North (NO), Central West (CW), Southeast (SE), and Southwest (SW) [37]. This dataset is selected for its ice-sheet-wide coverage and high spatial resolution.

We utilize satellite imagery from two sources: microwave imagery (S1 satellite) and optical imagery (S2 satellite). Data is acquired using the Google Earth Engine (GEE) [38]. S1 imagery is already preprocessed within GEE by applying thermal noise removal, radiometric calibration, terrain correction, and conversion to decibels via log scaling. For each supraglacial lake identified in 2018 and 2019, the dataset includes S1 imagery from January 1 to December 31 of the respective year. For optical imagery, we use S2 Level-1C orthorectified top-of-atmosphere reflectance, specifically Bands 2 (Blue, 20 m), 3 (Green, 20 m), 4 (Red, 20 m), 10 (Cirrus, 60 m), and 11 (SWIR 1, 20 m). To generate a complete annual time series of $HV_{anom}$ and $p_{water}$ for each lake, we linearly interpolated between S1 and S2 observations and applied a 12-day smoothing filter to $HV_{anom}$ to reduce variability across different S1 orbits. Our analysis focuses on time series data from May 1 to December 31 to capture the seasonal dynamics of supraglacial lakes during the melt season.

In this study, we classify supraglacial lakes into three categories: refreeze, drain, and buried using Reconstructed Phase Space (RPS) and Gaussian Mixture Models (GMM). The methodology is divided into two main phases: training and testing. The training phase begins by selecting one representative sample for each class. We then construct the RPS for these samples using time lag ($\tau$) and embedding dimension ($d$). Instead of using false nearest neighbor or automutual information methods, we employ a grid search technique to select $\tau$ and $d$ from a range of [2, 30]. For each combination of $\tau$ and $d$, we construct the RPS for the representative samples of each class. We then initialize GMMs for each class with 10 mixtures. For each GMM, we use $k$-means clustering to generate 10 different sets of initial parameters and select the best-performing set. Subsequently, we train the GMM models $M=\{M_{1},M_{2},M_{3}\}$ for each class—refreeze, drain, and buried—using the EM algorithm. The trained GMMs serve as probabilistic models for each class of supraglacial lakes. This approach is computationally efficient, requiring only one representative sample per class, which reduces the computational burden compared to other machine learning models that need larger training datasets.

During the testing phase, we create the RPS for each instance in the dataset using the current combination of $\tau$ and $d$ in the grid search. Each instance’s RPS is evaluated against the trained GMMs, and the likelihood score is computed for each GMM. Using the maximum likelihood classifier, the class with the highest likelihood score is assigned to the test instance, classifying each supraglacial lake as either refreeze, drain, or buried.

To evaluate the performance of our model, we calculate accuracy for each combination of $\tau$ and $d$ during the grid search. The optimal values $\tau^{*}$ and $d^{*}$ are chosen based on the highest accuracy achieved. A step-by-step outline of the proposed methodology is provided in Algorithm 1. The dataset and the code are available on GitHub¹¹1https://github.com/ehfahad/TSC-of-Supraglacial-Lakes-Evolution-over-GrIS.

The RPS-GMM model is trained with two sets of features: $HV_{anom}$ alone and a combination of $HV_{anom}$ and $p_{water}$. The reason for having two sets of features is that we aim to evaluate the impact of including the water percentage feature on the model’s performance. Our dataset includes time series data from 777 lakes, categorized into 189 refreeze, 392 drained, and 196 buried lakes, across all six subregions of the GrIS for the years 2018 and 2019.

As shown in Figure 4, the model trained with only the backscatter difference achieves an accuracy of 85.46%. Incorporating the water percentage substantially improves the accuracy to 89.70%. This improvement underscores the value of adding water percentage with the backscatter signal. Despite an imbalanced dataset, the weighted averages of precision, recall, and F1 score also reflect better performance across all classes when incorporated $p_{water}$ with $HV_{anom}$.

To further validate the RPS-GMM model, we compare its performance with several established ML and DL models using the sktime time series package [39]. Each model is evaluated through $5$-fold cross-validation, with 80% of the data ($\sim$622 samples) used for training and 20% for testing in each fold.

The results, detailed in Table I, show that the RPS-GMM model outperforms all compared ML and DL models in terms of accuracy. For example, the highest accuracy among existing models using only the backscatter difference is 84%, by the LSTMFCNClassifier. When incorporating both features, most models experience decreased performance due to the added complexity and limited training sample size. The LSTMFCNClassifier shows a slight improvement, possibly due to its inherent ability to capture long-term temporal dependencies, achieving 87% accuracy compared to the RPS-GMM model’s 89.70%.

The superior performance of the RPS-GMM model, especially when including the water percentage feature, highlights the effectiveness of a comprehensive feature set for accurately classifying supraglacial lakes. Notably, the RPS-GMM model, trained on only one representative sample per class, significantly outperforms models trained on 80% of the data in $5$-fold cross-validation. This efficiency demonstrates the model’s robustness in capturing supraglacial lake dynamics with minimal training data. In contrast, existing ML and DL models showed varying effectiveness, with deep learning models generally underperforming, mainly when both features were included. This underperformance may be due to the complexity of the time series data and the relatively small dataset size, which might not be sufficient for effective training of deep learning models.