Daily Snow Water Equivalent Trends over the Great Lakes Basin: Computer Vision and Deep Learning-Based Approach

1. Introduction

The Great Lakes Basin (GLB) is one of the most populated watersheds in North America, with an estimated population of approximately 40 million people. About 63% and 37% of the population live in the U.S. and Canada, respectively [1]. Snow is crucial for sustaining the ecosystems and human communities’ overall well-being in the GLB [2]. Snow covers a pronounced proportion of the Canadian landmass during winter [3], and produces water that recharges groundwater resources, freshwater ecosystems, and contributes to irrigation [4,5].

Snowfall modulates lake ice levels by forming thicker snowpacks [6,7]. Lakes without seasonal ice cover experience high air temperatures [8]. Snow cover also influences soil temperature, soil freezing and thawing processes, and permafrost formation, with a profound impact on carbon and hydrological cycles, especially in the cold regions [9,10]. The timing, duration, accumulation, and melting processes of seasonal snow cover dictate its impacts on the Earth’s thermal regime [9]. The rate and timing of snowmelt affect runoff, which in turn influences plant-available water storage [11]. Rapid fluctuations in snow depth during the winter influence soil temperature and moisture, leading to varying patterns of soil respiration [12]. Increased snowmelt may result in soil attaining field capacity rapidly, stimulating below-root zone percolation and streamflow, whereas slower snowmelt may decrease streamflow [13,14]. For example, using hydrological model simulations and a new snowmelt tracking algorithm, it was shown that although 37% of precipitation falls as snow, 53% of the total runoff in the western United States emanates from snowmelt [15].

The winter climate over North America, driven by both climate-induced events and natural atmospheric cycles (e.g., Pacific Decadal Oscillation and El Niño Southern Oscillation), is changing in unpredictable ways [16]. Diverse variables, including temperature, precipitation, and the snow-albedo feedback, influence snow cover [17,18,19,20,21]. However, climate-induced variables portend the most detrimental effects, causing a shift in storm tracks with predictable consequences on the snow regime [22]. Climate warming is projected to exacerbate changes in snow parameters, snow onset and retreat periods, with unpredictable impacts on the Great Lakes water levels and the viability of its ecosystems [23]. The Great Lakes region experiences local precipitation bands, elevated convection, and lake effect snow bands, which enhance snowfall events [24]. Despite the prevalence of lake-effect snowfall (LES), under high greenhouse gas emission scenarios, snow water equivalent (SWE) is projected to decline in the snowbelt regions to the lee of Lakes Erie and Ontario [25].

Rapid snowmelt, driven by factors such as rain-on-snow, reduces the formation of snowpacks and increases the risk of flooding [26]. Thus, detecting changes in snow parameters is an important piece of information for managing water resources as well as planning and adapting to climate change impacts, yet the relationship between seasonal snow and climate variables is uncertain and non-linear, complicating SWE trend analysis [27]. Additionally, detecting and attributing changes in snow over the GLB is not a straightforward task due to the spatially heterogeneous nature of snow ablation processes [28]. Kunkel et al. [29] further emphasized that non-climatic issues complicate the identification of significant change and trends in snow variables. Therefore, methods that are invariant to stochastic internal variability and seasonal patterns of snow spatial structure are crucial for obtaining accurate estimates of changes in snow parameters.

Under climate change, large-scale losses of snow are projected, raising considerable concerns over future freshwater resources, especially in snow-dominated watersheds, where snowmelt has profound effects on streamflow dynamics [19,30,31]. Therefore, it is imperative to unravel how snow processes have transformed over recent decades and to do so at very high temporal resolutions to inform water resources planning, management, and adaptation to natural hazards such as snowmelt-induced floods. The overarching objectives of this study are to: (a) demonstrate Siamese UNet’s capability to detect basin-scale changes in snow parameters, (b) estimate daily changes in SWE at the watershed scale, and (c) characterize decadal trends and anomalies in SWE distribution over the GLB.

1.1. Snowfall Dynamics in the Great Lakes Basin

Synoptic-scale analysis provides evidence that the Great Lakes exert a strong influence on the patterns of snow accumulation and ablation (i.e., Lake-induced snowfall) [32,33,34]. Although Lake-effect snowfall is known to drive an increase in the snowfall/snow water equivalence ratio in the GLB, time series analysis of snowfall within lake-effect areas shows that this pattern likely began in the 1930s and ended in the early 1970s [35]. Lake-effect snowfall is found to increase snowfall in most regions of the Great Lakes. For instance, elevated lake-effect snowfall near the Great Lakes from 1951 to 1980 [36]. In a related study, Burnette et al. [37] identified persistent upward snowfall trends in the lake-effect snowbelt areas of the Great Lakes region until the end of the century. Ellis and Johnson [35] also identified a 40-year upward trend in snowfall from 1930 to 1970, followed by less significant trends after this period.

Contrary to these trends, meteorological variables and weather types are becoming more favourable for snow ablation in the GLB; these include elevated air temperatures, southerly wind flows, less cloud cover, and increased liquid precipitation rates during the winter months [38]. Consequently, several studies have reported declining snowfall events in the GLB. For example, seasonal basin-wide snow depth was found to have declined by approximately 25% from 1960 to 2009, especially in the northern part of Lake Superior [39]. In a related study, analysis of snowfall and total precipitation trends over the Great Lakes from 1980 to 2015 revealed a significant decrease in snowfall at a rate of 40 cm/36 years, along the Ontario snowbelts of Lake Superior and Lake Huron-Georgian Bay [40].

Kunkel et al. presented a 107-year analysis of snow trends over the conterminous United States and found that from 1950 to 2006, snow trends in the central and northwest U.S. regions experienced significant declines in high-extreme snowfall years, while four regions (i.e., Northeast, Southeast, South, and Northwest) showed significant increases in the frequency of low-extreme snowfalls [41]. Snowfall and total precipitation trends were derived for the GLB from 1980 to 2015, revealing a significant decrease in snowfall, at a rate equivalent to 40 cm/36 years, and a considerable reduction of total precipitation of 20 mm/36 years, especially in Lake Superior and Lake Huron-Georgian Bay along the Ontario snowbelts [42]. Furthermore, extreme daily snowfalls have been found shift northward, with large daily changes occurring along the US-Canada border [43]. These findings are corroborated by snow projections in Central and Eastern North America, affirming a decline in extreme snowstorms southward, while the northern regions experience increased snow water equivalent and frequency of snowstorms [44]. In addition, declines in mean annual SWE, high-frequency of low-snowfall years, and less frequent high-snowfall years were projected for the western United States [45]. Furthermore, in western North America, annual snow water storage was reported to have decreased from 1950–2013 as a consequence of early snowmelt, rainfall in spring months, and decreases in winter precipitation [46].

Despite these consistent findings, Kunkel et al. [29] point out temporal inconsistency and inhomogeneity in the data employed to conduct most of these analyses and call for the use of more homogeneous datasets. For instance, using a set of consistent station data, an upward trend in snowfall was identified in Superior and Michigan, largely attributable to warmer surface water and less ice cover [33]. This suggests that gridded snow data hold the potential for a more enhanced and spatial-temporally continuous representation for studying snow. For example, the Snow Data Assimilation System SWE was employed to model snow processes for climate-related water resources studies in the Great Lakes watersheds [47] and stream flow in eastern Canadian basins [48].

2. Materials and Methods

2.1. Study Area and Data

Our study area covers the five Great Lakes and their watersheds. Using ArcGIS Pro’s ArcPy Package, the GlobSnow SWE data was clipped to a spatial extent encompassing 72 × 72 pixels. The GlobSnow’s SWE has a spatial resolution of 25 km × 25 km, covering an area of 1,800 km², given the spatial extent of 72 × 72 km. This area encapsulates the GLB and its surrounding ecosystems. We reprojected the data, adopting an Equal Area Projection (e.g., Lambert Equal Area Projection) to maintain the original grid area. Figure 1 depicts the SWE distribution in early, mid, and late February and March over the GLB in 2020. The northwestern part of the GLB, corresponding to the Lake Superior basin, can be seen to have relatively elevated levels of SWE in February. A northward shift in SWE can be observed between February and March. This gradual transition is indicative of the end of the Winter season in the GLB.

2.2. Training Sample Processing and Labeling

We labeled the training data following the methods illustrated in Malik and Colin [49]. Thus, SWE samples at one day interval were labeled as positive (No Change), whereas SWE observations that were two or more days apart were labeled as Change instances. We deployed a spatially-explicit metric – the Structural Similarity (SSIM) index to inspect and exclude ambiguous training examples using 0.98 and 0.90 as thresholds for positive pairs (No Change instances of SWE) and negative pairs (Change instances of SWE), respectively, as shown in [49]. For model development, we used data from 1979 to 2001. The training and validation data splits were 75% (training set) and 25% (validation set). Independent validation consists of data from 2002 to 2018. We separated the data into seasons or months and evaluated the model’s performance across seasons.

2.3. Siamese U-Net Architecture and Training

Our model architecture, with the SSIM index combination, is simplified and depicted in Figure 2. The model receives 72 × 72 GlobSnow SWE data and outputs the same spatial extent at the decoder branch. Using the SSIM index, the resulting outputs from the two-sister decoders are compared, focusing more on changes in the spatial structure SWE. The computed similarity value is passed to a loss function for error estimation and optimization during model training. An exhaustive description of this learning paradigm is outlined in [49,50]. The Model is trained using Keras-Tensorflow framework with NVIDIA GPU memory of 40 GB. Fifty epochs were used with Early Stopping and a patience parameter of 10. Drop-out, batch-normalization, and ReLU activation function were applied.

2.4. Model Accuracy Metrics

To assess the model’s accuracy on the prediction of No Change versus Change instance of SWE, we adopted the metrics: true positive rate (TPR), true negative rate (TNR), false positive rate (FPR), false negative rate (FNR), Precision (PR), F1 score, and overall accuracy (OA) to assess the models’ prediction of SWE similarity and change.

$$Recall\left(TPR\right)={\displaystyle \frac{TP}{TP+FN}}$$

(1)

$$Precision={\displaystyle \frac{TP}{TP+FP}}$$

(2)

$$F1score={\displaystyle \frac{2TP}{2TP+FP+FN}}$$

(3)

$$TNR={\displaystyle \frac{TN}{TN+FP}}$$

(4)

$$OA={\displaystyle \frac{TP+TN}{TP+TN+FP+FN}}$$

(5)

TPR measures the proportion of SWE samples that are correctly predicted as No Change instances; TNR, on the other hand, estimates the fraction of SWE samples that are accurately detected to be Change instances. Thus, FPR = 1 –TNR, and FNR = 1 – TPR. Equations 1 – 5 represent the accuracy metrics. Area Under the Curve of the Receiver Operating Characteristics (AUC-ROC) was also derived to examine TPR and FPR over a range of thresholds.

2.5. Deriving the Daily Snow Water Equivalent Similarity Vector

To estimate daily SWE change events, the pre-trained Si-UNet, after rigorous accuracy assessments, can be employed to make predictions over a series of SWE observations, generating a time-series vector consisting of similarity values between consecutive days. High similarity values are indicative of less pronounced changes; contrarily, low similarity values characterize profound changes in SWE. We formulate the solution as illustrated in [49,50]. The SWEsim vector, which represents daily SWE similarity, is derived using the equation 6 to accumulate the daily SWE similarity values.

$$\overrightarrow{{{SWEsim}_m}_{\left\{t_i,\mathrm{t}\right\}}}=model\left\{{SWE}_m\left({t_i,t}_{i+1}\right)\right\}$$

(6)

Where t and m represent day and month, respectively. Accordingly, $t_i$ denotes SWE sampled on the ith day. The model term in the equation denotes the Si-UNet.

Using the SWEsim vector, we derived long-term median SWEsim from 1981 to 2016 to compare deviations in SWE over 35 years (Equations 7 and 8). We refer to this median as a 35-year baseline (35YB). Deviations from 35YB will be indicative of anomalous SWE trends at basin scales.

$${\mathrm{S}\mathrm{W}\mathrm{E}\mathrm{s}\mathrm{i}\mathrm{m}}_{35YB}=\sum _{i=1}^{N}median\left\{SWEsim\left(1981,1982,1983\dots 2016\right)\right\}$$

(7)

$${\mathrm{S}\mathrm{W}\mathrm{E}\mathrm{s}\mathrm{i}\mathrm{m}}_{div}={SWEsim}_i-{\mathrm{S}\mathrm{W}\mathrm{E}\mathrm{s}\mathrm{i}\mathrm{m}}_{35YB}$$

(8)

3. Results

3.1. Model Accuracy Assessment Report

Figure 3A presents a range of possible thresholds for the model’s predictions of SWE similarity, while the AUC-ROC is depicted in Figure 3B. The F1-score, Precision, and Recall appeared to be constant as the threshold parameter was varied until 0.45, where it began to rise steeply. All the accuracy metrics reached the optimum value at a 0.5 threshold. The F1-score and Recall dropped to 0 when the threshold exceeded 0.5. Contrarily, the false positive predictions reduce to near zero, resulting in Precision scores rising to 1.0 after a 0.5 threshold. The AUC-ROC started at 0.8, increasing gradually to 1.0, and stabilized at this point throughout the range of possible thresholds. These thresholds are values that were greater than or equal to 0.5.

We compared the accuracy of our proposed Si-UNet to the CNN baseline, using the SSIM index and Euclidean Distance (ECD) metric. Further, we investigated the performance of the models given a choice of loss functions – binary cross-entropy loss (BCE) and contrastive loss (CL). Table 1 presents a summary of the models’ accuracy. It can be observed that our model, with the CL and SSIM index, though has a narrow threshold value (i.e., 50%), tends to outcompete the other models, yielding an F1-score of 98.73%. The most challenging instances were the change-pairs (i.e., negative labels), resulting in a TNR of 92.82%. Surprisingly, the base CNN with ECD can be observed to excel at detecting instances of SWE change (TNR) with 99.44% accuracy; the model’s performance at detecting the No Change (TPR) instances tended to be the worst (92.84%), however.

3.2. Empirical Evaluation of Snowfall and Ablation in the Great Lakes Basin

Our model identified three types of sequential snow-related events that create daily variation of SWE trends in the GLB. The first is an alternating increase in snowfall (i.e., snow accumulation) and decreases in snow on the ground (i.e., snow ablation). We termed this sequence a snowfall-ablation (SA) event. This yielded upward- and downward-pointing peaks that characterize sharp transitions between snowfall and ablation events. For example, on the graph in Figure 4, points 2, 11, and 25, corresponding to SWEsim values: 0.523, 0.525, 0.526, represent days with pronounced snowfall, whereas points 7, 9, and 12, aligning with SWEsim: 0.515, 0.516, 0.520, were days that experienced snow ablation. Upward peaks are associated with high SWE similarity, while downward peaks signify low SWE similarity. The second sequence type is snowfall-snowfall (SS), and the third, ablation-ablation (AA). The SS is characterized by two or more consecutive days of snowfall, resulting in pronounced snow accumulation that increases snow depth and, predictably, the amount of SWE (Figure 4, D, E, F). Conversely, AA events are 2 or more consecutive days of snow ablation; hence, resulting in a substantial SWE reduction (Figure 4, A, B, C). Both SS and AA events induce ridges or crests and troughs in daily SWE trends. However, crests are largely generated by SS events in snowfall-dominant months (e.g., January, February, and March), whereas AA events frequently induce troughs. Figure 4 depicts the three types of events, with arrows showing the monotonic pattern and direction of change. The most perceptible pattern (shown in the red rectangles) is SA events in the eastern wing of the Lake Michigan sub-basin. Also noticeable is the SS event in the northwestern part of Lake Ontario sub-basin (shown in the green rectangles). The dashed lines represent simple linear regression (orange) and the 35YB (black). The computed 35YB is ~ 0.5108, aligning closely with the No Change SWE instance (i.e., No Change SWE days for the Si-UNet model prediction). Thus, it can be observed that all the SWEsim values for February 1990 were above the reference 35YB.

3.3. Great Lakes Basin Daily SWE Trends

Using the ${\mathrm{S}\mathrm{W}\mathrm{E}\mathrm{s}\mathrm{i}\mathrm{m}}_{div}$ vector, we present 30-year daily SWE similarity trends from 1989 to 2018 (Figure 5, Figure 6 and Figure 7). The y-axis and x-axis represent SWE deviations from 35YB and days, respectively. The blue trend lines denote daily variability in SWE, and the red lines are linear regression conditioned on ${\mathrm{S}\mathrm{W}\mathrm{E}\mathrm{s}\mathrm{i}\mathrm{m}}_{div}$. Positive (i.e., positive SWE anomaly) and negative (i.e., negative SWE anomaly) deviations are values above and below the 35YB, respectively. Zero values, on the other hand, signify no deviation from the 35YB. Positive and negative deviations are largely associated with high and low SWE scenarios, respectively. Both deviation types are a function of snowfall and SWE reduction induced by ablation events.

In April (Figure 5), strong statistically significant increases in SWE similarity were detected in 1991 and 1993 (R² = 0.28 and 0.3, respectively). Contrarily, 1989, 1992, 1996, and 1997 all exhibited a significant decline in SWE similarity (R² ranged from 0.22 – 0.46). In March, less substantial increases and decreases in SWE trends dominated. It can be observed that the declining trends occurred when the SWE similarity decreased in mid-March (e.g., 13–17 March). These trends are more obvious in March 1990 and 1995. Two statistically significant but opposing trends were observed in February 1994 and 1995. While SWE declined significantly in 1994 (R² = 0.39), it increased significantly in 1995 (R² = 0.57). January and February depicted a predominantly increasing SWE similarity and trend. While the SWE trends in 1996 and 1997 were not statistically significant, the trends for most years were characterized by substantial positive SWE trends. One notable exception was February 1994, in which the trend was significantly declining (R² = 0.39).

Figure 6 shows daily changes in SWE from 1999 to 2008. Significant positive SWE trends can be observed in April, with R² ranging between 0.41 – 0.72. The strongest positive trends were observed in 2000 (R² = 0.47) and 2005 (R² = 0.72). The 2000 trend was driven largely by snow accumulation (i.e., SS events) occurring after April 5th, whereas the 2005 trend can be attributed to AA events starting on April 5th and progressing gradually until the end of April. A moderately strong negative trend was observed in 2002 (R² = 0.21). SWE trends in March were predominantly weakly positive. The most significant positive trend was observed in 2000 (R² = 0.32). February and January were characterized by positively significant trends in SWE. February 1999 and 2005 recorded high SWEsim values with strong positive trends (R² = 0.3 and 0.60, respectively), whereas January 2001, 2003, and 2008 portrayed competing positive trends (R² = 0.43 – 0.44).

Daily changes in the GLB SWE are depicted in Figure 7. A predominant negative trend characterized the SWE pattern from April 2008 to 2018. While a significant positive SWE trend was observed in 2010 (R² = 0.39), the remaining years exhibited negative SWE trends (R² = 0.2-0.59). The decline frequently occurred between April 13th and 17th for most of the years. March was characterized by frequent positive trends, with highly significant SWE increase occurring in 2014 and 2015 (R² = 0.37 and 0.27, respectively). Contrary to March trends, February and January depicted a substantial prevalence of non-significant trends. In February 2009, 2016, and 2017, SWE exhibited increasing trends (R² = 0.22, 0.20, and 0.36, respectively), while the rest of the years showed no statistically significant change. For example, January largely portrayed no significant trends, and where trends were observed (e.g., 2010 and 2018), they tended to be weakly significant.

3.4. SWE Anomaly Frequency in the Great Lakes Basin

The frequency of SWE deviation (in days) from the 35-YB throughout January, February, March, and April (JFMA) is depicted in Figure 8, Figure 9 and Figure 10. The y-axis denotes the frequency of days having values above or below the median, whereas the x-axis represents deviation from the 35-YB. The red dotted vertical line demarcates deviations below and above the median (i.e., values below and above zero).

Figure 8 depicts SWE deviations from the 35YB SWEsim between 1989 and 1998, through JFMA. The histogram bins are more clustered and continuous for most of the months and years. Aside from January 1994 and 1995, when SWE was largely below the baseline, the rest of the years had either high SWEsim values or a balanced frequency above and below the baseline. Contrary to this trend, from February 1989 to 1998, SWE exhibited consistently higher and clustered SWEsim values. The deviations were between -0.005 and +0.01 with frequencies starting at ~2 days. SWE anomaly patterns in March also align with February, as positive deviations (i.e., SWEsim values above the 35-YB) were prevalent. The deviations for February, March, and April exhibited balanced frequency from 1993 to 1995. Similarly, April was characterized by symmetrical frequency of values below and above the baseline. The only notable exceptions were 1990 and 1996. These years incurred two contrasting SWE deviations: pronounced daily changes in SWE, resulting in high frequencies of deviations below the 35YB (April 1990), and fewer daily changes, yielding high frequencies of deviations above the 35YB (April 1996).

From 1999 to 2008 (Figure 9), it can be observed that January exhibited the lowest SWE similarity values (i.e., values below 0.5108), while there was a balance between daily occurrences of low and high values in April. Contrarily, the frequency of daily occurrences of high SWEsim values appeared to be elevated in February and March. A notable exception was in 2007, wherein JFMA exhibited values below the 35YB, with deviations ranging from -0.01 to -0.03.

Figure 10 presents SWE deviations between 2009 and 2018 compared to the 35-YB. It can be deduced that January and February tended to have anomalously low SWEsim values, ranging from -0.01 to -0.05. Except for 2009 and 2018, where higher SWEsim values were recorded in March, the remaining years elicit a balanced frequency of daily variability in SWE. Similarly, April 2010, 2012,2013, and 2014 were characterized by balanced frequencies of daily changes in SWE relative to the 35-YB.

4. Discussion

In this study, we developed a U-Net model inspired by the Siamese model architecture for detecting anomalies in classification tasks [51]. Our model detected daily changes in SWE with 98.73% accuracy (F1-score), outcompeting the baseline CNN and U-Net models without spatially-explicit metrics (i.e., the SSIM index). The AUC-ROC of 0.997 further accentuates the model’s robustness in detecting changes in SWE. Guided by the SSIM index, the model evaluates SWE similarity based on structural differences between pairs of consecutive SWE samples. This reduces the models’ tendency to capture spurious structural differences that are not related to the underlying spatial processes driving snow accumulation and ablation [52,53].

We computed a 35-year median baseline SWE similarity (i.e., 35YB) of ~ 0.5108 and identified three types of alternating snow-related events using a sample SWE times series in February, which is the SWE month dominant month in the GLB: (1) is snowfall-snowfall (SS), (2) ablation-ablation (AA), and (3) snowfall-ablation (SA). These events characterize transitions between snowfall and ablation processes (Figure 4). The points 2, 11, and 25, in Figure 4, correspond to SWEsim values: 0.523, 0.525, and 0.526. These were days with pronounced snowfall. Conversely, in Figure 4, points 7, 9, and 12 were days experiencing high ablation. Thus, the SWEsim values tended to be relatively low, 0.515, 0.516, 0.520, respectively. The most perceptible pattern is SA events in the eastern wing of the Lake Michigan watershed. Also noticeable is the SS event in the northwestern part of Lake Ontario (Figure 4). While the northern and southern parts of the GLB, respectively, experience peak ablation events in April and February, most of the regions in the basin experience peak ablation in March [54]. This pattern may have an impact on March SWE anomalies by increasing the frequency of days incurring SWEsim below the baseline. Ablation events have been found to decrease within the northeastern and northwestern Lake Superior basins, while occurring frequently in the eastern Lake Huron and Georgian Bay basins [54]. This pattern is obvious in Figure 4, where SWE persists around the Lake Superior sub-basin throughout February.

Alternating positive and negative deviations are signs of snowfall-ablation, which creates sharp peaks in the time series (i.e., trend line). The trend line becomes virtually flat (e.g., R² < 0.001), with less pronounced crest and trough regions when the magnitude and frequency of alternating snowfall and ablation events are equivalent (e.g., Figure 5, February 1989). Figure 5, Figure 6 and Figure 7 present 35-year daily SWE trends over the GLB. Linear regression conditioned on the days of SS, AA, and AS revealed complex patterns of the daily changes. Significant positive daily SWE trends followed elevated frequency of daily SWEsim values (e.g., SWEsim values equivalent to or above 35-YB). Conversely, negative trends are largely driven by a higher frequency of anomalously low daily SWEsim values. The trends were frequently above the 35-YB for most of the days in the early 2 decades (i.e., from 1989 to 2008). In April 1989, 1992, 1996, and 1997, SWE trends significantly declined (R² range from 0.22 – 0.46), suggesting that April fundamentally marks the end of the winter season in the GLB. The remaining years were characterized by substantial positive SWE trends for January and February. In the western United States, a decline in SWE has been observed, with an average reduction of 15 – 30% (25 – 50 km³) of April 1st SWE since the mid-century [55]. Increased snowfall days have been reported in the core winter months, December to February, due to lake-effect snowfall within the lake belt regions relative to inland areas [56]. Monthly snow trends in snowfall-to-precipitation ratio were also analyzed for northern New England, showing a significant decline in this ratio for March and December [57]. In March, declining trends occurred when the SWE similarity decreased in mid-March (e.g., 13–17 March). This pattern is indicative of a shortening winter season length, characterized by limited snowfall in November and March [56]. While SWE declined significantly in February 1994 (R² = 0.39), it increased significantly in 1995 (R² = 0.57). This interannual variability may be linked to large-scale climate processes [16,58,59]. Annual snow variability is also attributable to the changing nature of rain-snow transition, driven largely by near-0 °C temperatures that vary across climatological zones and regions, impacting snowstorm line locations [60,61].

Figure 8, Figure 9 and Figure 10 show daily SWE anomalies relative to the 35-YB. A significant shift to SWE similarity values below the 35-YRB occurred after 2009, with January and February being the months with high frequency of days when SWE fell below the baseline (–0.01 to –0.04). The late arrival of snow in parts of southern Canada likely drives large SWE dissimilarity, resulting in SWEsim anomalies below the 35-YB. Suriano and Leathers [62] study partly corroborates the decreasing pattern of daily SWE over the GLB in February 2009. The researchers concluded that the ablation event on 10 February 2009 culminated in extensive snow cover reduction across the entire basin and a basin-scale average ablation of 5.34 cm. The ablation event was attributed to advective surface air flow, and dewpoint temperatures above 0 degrees, coupled with elevated wind speeds in eastern Lake Erie, and west-north of Lakes Superior and Michigan. This probably drove -0.007 shift in SWE below the 35-YB between 9 – 10th February 2009. Likewise, a decrease in snowfall east of Lake Michigan has been identified and linked to both reductions in lake-effect and system snowfall, with system snowfall being dominant in the western region of the lake. The ablation events were primarily driven by snowfall sensitivity to regional temperature, and rain-on-snow events [63]. Despite the reported reduction in snowfall, SWE remained above the 35-YB for most of the days in February, suggesting that February is potentially a high SWE month in the GLB. This is substantiated by the findings that high SWE years often experience greater snowfall, less accumulated melting degree, and SWE above normal snowfall in February [64].

Kluver and Leathers [65] identified snowfall regions in the Conterminous United States, pointing out two dominant snowstorm track regions in the GLB. The northern part of the Upper-Midwest was found to experience increasing trends in snowfall frequencies, whereas areas (mainly around Lake Michigan) showed no significant frequency of snowfall events. These findings reveal the complex nature of daily SWE trends in the GLB. However, Musselman et al. [66] emphasize the need to examine snowmelt and SWE trends in tandem. They compared long-term changes in snowmelt and SWE in western North America and found that winter snowmelt trends increased about three times the decline in SWE trend, suggesting that snow water resources may be declining more rapidly than portrayed in SWE analysis alone. Snowmelt trends are highly sensitive to ambient temperatures, while SWE trends are modulated by precipitation variability, both of which are linked to climate change signals [67,68,69,70]. Thus, the prediction, planning, and management of future water resources and infrastructure may become more complicated due to changes in SWE and snowmelt in the changing climate [71].

The limitations of our study can be distilled into the modelling framework (i.e., Siamese U-Net, SSIM index, and threshold parameters for determining change in SWE), SWE datasets, and the spatial scale of the analysis. Although we compared the U-Net to the baseline CNN, exploring the capabilities of other deep learning models to characterize changes in SWE could provide valuable insights into how model sensitivity to daily changes in SWE translates into the decadal trends. For example, as pointed out in previous studies, the Siamese Attention model may be more sensitive to daily SWE variations [50]. Additionally, the threshold parameter adopted is model-dependent and could affect the magnitude of SWE trends if different models were used to derive the SWEsim vector. In the context of spatial scale, the heterogeneous nature of snow cover and the relatively large geographic extent of the GLB imply that spatial processes such as snow accumulation and ablation likely exist at sub-basins [72]; therefore, examining SWE variability for each of the Great Lakes’ sub-basins could help disentangle the complexities of SWE trends and variability over time. Furthermore, climatological factors can complicate the identification of significant change and trends in snow variables [29]. Therefore, synoptic-scale analysis of changes in snow and the underlying SWE could provide more insight into localized-scale influences of climatic drivers and lake-effect snowfall on daily SWE trends.

5. Conclusions

We developed a convolutional Siamese U-Net to detect and characterize changes and trends in SWE over the GLB. Our model effectively detected daily changes in SWE with 98.73% accuracy and an AUC-ROC of 0.997, outperforming the baseline CNN model. The model detected decadal trends in daily SWE, allowing a fine-temporal resolution glimpse into the changes in snow parameters across the Great Lakes regions. Except for a few years (1994, 1996, and 1997), daily changes in SWE largely followed positive trends from 1989 to 2008, especially in January and February. Most of the trends in March were non-significant, suggesting a possible contraction of the Winter season to a few snow-dominant months (i.e., December, January, and February). We found minimal deviations from a 35-YB SWE similarity between 1989 and 2008. While SWE similarity tended to be frequently high and above the 35-YB in the first 2 decades (i.e., 1989 – 2008), it declined significantly in the latter decade (i.e., 2009 – 2018). By examining previous studies on snow in the GLB, we linked these trends to characteristic patterns of snowfall and snow ablation processes, driven by natural climate variability, global climate warming, and lake-effect snowfall, the latter two being the most significant drivers of daily SWE trends in the GLB.