# Recent Ice Trends in Swiss Mountain Lakes: 20-year Analysis of MODIS Imagery

Manu Tom  Tianyu Wu, Emmanuel Baltsavias, Konrad Schindler  Photogrammetry and
Remote Sensing, ETH Zurich, Switzerland

###### Abstract

Depleting lake ice can serve as an indicator for climate change, just like sea
level rise or glacial retreat. Several Lake Ice Phenological (LIP) events
serve as sentinels to understand the regional and global climate change.
Hence, monitoring the long-term lake freezing and thawing patterns can prove
very useful. In this paper, we focus on observing the LIP events such as
freeze-up, break-up and temporal freeze extent in the Oberengadin region of
Switzerland, where there are several small- and medium-sized mountain lakes,
across two decades (2000-2020) from optical satellite images. We analyse time-
series of MODIS imagery (and additionally cross-check with VIIRS data when
available), by estimating spatially resolved maps of lake ice for these Alpine
lakes with supervised machine learning. To train the classifier we rely on
reference data annotated manually based on publicly available webcam images.
From the ice maps we derive long-term LIP trends. Since the webcam data is
only available for two winters, we also validate our results against the
operational MODIS and VIIRS snow products. We find a change in Complete Freeze
Duration (CFD) of -0.76 and -0.89 days per annum (d/a) for lakes Sils and
Silvaplana respectively. Furthermore, we correlate the lake freezing and
thawing trends with climate data such as temperature, sunshine, precipitation
and wind measured at nearby meteorological stations.

###### keywords:

lake ice detection, climate monitoring, machine learning, semantic
segmentation, support vector machine, random forest, xgboost, satellite image
processing, MODIS, VIIRS, time-series analysis

††journal: RSE

##  1 Introduction

Scientists around the globe are making an effort to understand the changing
climate, to find ways to mitigate the impact of associated extreme weather
conditions, to protect the environment for the future generations. The
repercussions of climate change are foreseen to amplify in the next few
decades. Furthermore, the latest climate models underline the need for urgent
mitigation [Forster et al., 2020]. ”Human activities are estimated to have
caused approximately $1.0^{\circ}C$ of global warming above pre-industrial
levels, with a likely range of $0.8^{\circ}C$ to $1.2^{\circ}C$. Global
warming is likely to reach $1.5^{\circ}C$ between 2030 and 2052 if it
continues to increase at the current rate”, said the IPCC special report on
impacts of global warming [Masson-Delmotte et al., 2018].

Increasing temperature patterns are expected to cause decreasing trends in the
lake ice formation process. Air temperature in the vicinity of a lake affects
the ice build-up process within the lake and vice versa. Moreover, there are
potential positive feedbacks, as frozen lakes (especially when covered with
snow) have higher albedo, and thus lower absorption and evaporation [Slater et
al., 2021, Wang et al., 2018].

Regional weather patterns and lake ice formation processes are inter-
connected. Hence, monitoring the long-term LIP trends can provide integral
cues on the regional and global climate. In addition to its usefulness for
climate studies, lake ice monitoring is also crucial to organise safe
transportation (especially in lakes that freeze only partially), to conserve
freshwater ecology, to trigger warnings against ice shoves caused by wind
during the break-up period, and for winter tourism.

We aim to monitor the main lakes of the Oberengadin region in the Swiss Alps
(which reliably freeze every winter) on a daily basis during the winter
months, to derive the spatio-temporal extent of lake ice. Specifically, we
estimate the four important LIP events: Freeze-Up Start (FUS), Freeze-Up End
(FUE), Break-Up Start (BUS) and Break-Up End (BUE). Using these four dates, we
also estimate the Complete Freeze Duration (CFD) and Ice Coverage Duration
(ICD), refer to Table 1 for definitions.

There are no standard definitions for the ice-on/off dates. Some publications
have termed FUE and BUS as ice-on and ice-off dates, respectively [Hendricks
Franssen and Scherrer, 2008, Tom et al., 2020b]. However, other researchers
(and the NSIDC database, https://nsidc.org/) consider BUE as ice-off [Duguay
et al., 2015]. The Global Climate Observing System (GCOS) pre-requisite for
the accuracy of the ice-on/off dates is $\pm$2 days.

In this work, we focus only on optical satellite data [we have previously also
used webcams [Tom et al., 2020b, Xiao et al., 2018] and Sentinel-1 Synthetic
Aperture Radar [Tom et al., 2020a, SAR] for lake ice monitoring]. Note that,
we do not estimate the ice thickness, only its spatio-temporal extent.

Compared to other sensors, MODIS and VIIRS satellite data have several
advantages such as wide area coverage, good spectral and very good (daily)
temporal resolution, free availability etc. Additionally, compared to other
optical satellites such as Landsat-8, Sentinel-2 and the like, MODIS and VIIRS
offer the best spatio-temporal resolution trade-off for the application of
single-sensor lake ice monitoring, even though the spatial resolution is
moderate (250-1000m Ground Sampling Distance, GSD). In addition, the global
coverage is highly beneficial to eventually scale up to country- or world-wide
monitoring. On the other hand, cloud cover is a bottleneck for optical
satellite data analysis. An important asset is the availability of large time-
series, e.g., MODIS data is available for the entire period since 2000
(contrary to many other sensor data like airborne or terrestrial photography,
webcams etc.). This makes it possible to implement a 20 year analysis and to
derive LIP trends.

The last decades have seen the rise of Machine Learning (ML) as a tool for
data analysis in remote sensing and the Earth sciences. That is, large-scale
statistical data analysis is used to capture the complex input-output
relationships in a data-driven manner. ML is a powerful tool to recognise the
underlying patterns in data where mechanistic models are lacking or too
complicated. We leverage it to create a 20 year time-series of ice cover in
Swiss mountain lakes primarily using the Terra MODIS (https://terra.nasa.gov)
data, and show empirically that the ice formation indeed follows a decreasing
trend. We cast lake ice detection as a 2-class (frozen, non-frozen) supervised
classification problem. As part of our study, we compare the performance of
three popular ML methods: Support Vector Machine [Cortes and Vapnik, 1995,
SVM], Random Forest [Breiman, 2001, RF], and XGBoost [Chen and Guestrin, 2016,
XGB]. Additionally, we assess the sensitivity of these classifiers to the
respective hyper-parameters. We find that a linear SVM offers the best
generalisation across winters and lakes for our data, and derive LIP from the
resulting time-series by fitting a piece-wise linear model per winter.

###  1.1 LIP trend analysis studies

The LIP trends of several lakes with different geographical conditions have
been thoroughly studied and reported in literature. Though most of them use
information from various ice databases (for eg. NSIDC), some studies directly
derived the trends from radar and optical satellite data. In the following
paragraphs, we review the studies that dealt with optical satellite imagery,
especially MODIS data.

Gou et al. [2015] analysed the ice formation trends (2000–2013) in lake Nam Co
(Tibet, area of 1920 $km^{2}$) using MODIS and in-situ data and unravelled the
strong correlations with air temperature and wind speed patterns. This study
found that high wind speeds during winter time could potentially expedite the
freeze-up process. Additionally, this work reported a significant reduction of
the total freeze duration. Another study [Yao et al., 2016] also noted
increasingly shorter freeze duration during the period 2000–2011 when
investigating the lakes in Hoh Xil region (Tibet, 22 lakes with area
$>100~{}km^{2}$), using MODIS, Landsat TM/ETM+, and meteorological data. In
addition, this work estimated late freeze-up and early break-up trends. They
reported that the FUS, FUE, BUS, BUE, CFD and ICD shifted on average by
$0.73$, $0.34$, $-1.66$, $-0.81$, $-1.91$, $-2.21$ d/a respectively. Qi et al.
[2020] used AVHRR, MODIS, and Landsat data to extract the LIP of Qinghai lake
(China, area of 4294 $km^{2}$) for the period 1980–2018. They estimated a
shift of $0.16$, $0.19$, $-0.36$, and $-0.42$ d/a for FUS, FUE, BUS and BUE
respectively, also pointing towards progressively later freeze-up and earlier
break-up. Additionally, they computed the decreasing patterns in ICD ($-0.58$
d/a) and CFD ($-0.52$ d/a). This study also identified the correlations
between the LIP and climate indicators such as the Accumulated Freezing
Degree-Days (AFDD), wind speed, precipitation etc. during the winter season.
To our knowledge, unlike MODIS-based approaches, no work exists in the
literature which performed a multi-winter LIP trend analysis using VIIRS data.

To summarise, most related works reviewed so far have found trends towards
late freeze-up, early break-up and declining freeze duration. To our
knowledge, none of these trend studies applied data-driven ML methods to
identify lake ice, instead they relied on physics-inspired methods based on
empirical indices and thresholds. Beyond analysing the lakes of the
Oberengadin region, in the present work we also show that supervised ML models
are able to detect lake ice with high accuracy, hoping that the outcome may be
useful for future research.

For Swiss lakes, a previous study [Hendricks Franssen and Scherrer, 2008]
verified that the lake ice formation and surrounding air temperature are
heavily correlated. They deducted an empirical relationship between sum of
negative degree days (same as AFDD) and the lake ice build-up process, and
modelled the probability of ice cover via binomial logistic regression. Though
this approach gathered and used the temperature data from 1901–2006 to study
eleven lakes in the lower-lying Swiss plateau, none of the target mountain
lakes were included. Here we focus on the Oberengadin region, with three main
lakes: Sils, Silvaplana and St. Moritz (with the latter very small for the GSD
of MODIS). Moreover, we include another Alpine lake Sihl to check generality.

Our goal in this work is lake ice monitoring using only image data from
optical satellites, which provides a direct, data-driven observation not
influenced by model assumptions about the ice formation process. We see
satellite imagery as an independent information source and consider image
analysis complementary to other methods of lake ice modelling. Furthermore,
compared to our satellite-based approach which can easily analyse the whole
lake area, it is difficult to effectively derive the spatial extent of lake
ice from the point (temperature) observations recorded at one or a few nearby
weather stations, even if they are situated in the immediate vicinity of the
lake (which is not always the case). This is even more true in the Alpine
terrain that we target, due to strong micro-climatic effects.

###  1.2 Lake ice observation with machine and deep learning

ML algorithms have become a standard tool for several environmental remote
sensing research problems, including our earlier works on monitoring lake ice
cover. In Tom et al. [2018] we already investigated pixel-wise classification
of the spatio-temporal extent of lake ice from MODIS and VIIRS imagery with
SVM. Each pixel was classified as either frozen or non-frozen in a supervised
manner. Though this approach achieved strong results (including generalisation
across winters and across lakes with similar geographic conditions), the test
set at the time did not have a complete winter of reference data (including
the critical freeze-up and break-up periods), due to technical problems.
Later, Tom et al. [2020b] did extensive experiments using the data from two
full winters and confirmed the efficacy of SVM for lake ice monitoring with
MODIS and VIIRS. Both these works dealt with small- and mid-sized Swiss
mountain lakes. Xiao et al. [2018] and Prabha et al. [2020] explored the
potential of convolutional neural networks (CNNs) for lake ice detection in
terrestrial webcam images (RGB). They performed a supervised classification of
the lake pixels using the Tiramisu [Jégou et al., 2016], respectively Deeplab
v3+ [Chen et al., 2018] networks, into the four classes: water, ice, snow and
clutter. An integrated approach using both the satellite and webcam
observations was discussed in Tom et al. [2019] to estimate the ice-on and
ice-off dates. Recently, Hoekstra et al. [2020] proposed an automated approach
for ice vs. water classification in RADARSAT-2 data, combining unsupervised
Iterative Region Growing using Semantics (IRGS) and supervised RF labelling. A
deep learning approach to lake ice detection in Sentinel-1 SAR imagery has
been described in Tom et al. [2020a], and achieved promising results,
including generalisation across lakes and winters. Interestingly, even for
Sentinel-1 imagery it proved beneficial to employ a network pre-trained on
close-range RGB data. Very recently Wu et al. [2021] compared the capabilities
of four different ML methodologies: Multinomial Logistic Regression (MLR),
SVM, RF, and Gradient Boosting Trees (GBT) for lake ice observation using
MODIS Top of Atmosphere (TOA) product. They modelled lake ice monitoring as a
3-class (ice, water, cloud) supervised classification problem. The four
classifiers were tested on 17 large lakes from North America and Europe with
areas $>\\!$ 1040 $km^{2}$, and achieved $>\\!$ 94% accuracy. RF and GBT
showed better generalisation performance on this dataset of large lakes.

###  1.3 Operational lake ice / snow products

To our knowledge, the only operational lake ice product at present is the
Climate Change Initiative Lake Ice Cover [Crétaux et al., 2020, CCI LIC]. A
comparison of our results with the CCI LIC product is however not possible,
since none of our target lakes are included in the list of 250 lakes covered
by this product. A second product, Copernicus Lake Ice Extent (LIE,
https://land.copernicus.eu/global/products/lie), is still in pre-operational
stage due to accuracy issues, and coverage only starts in 2017. Though not
designed for lake ice, the MODIS Snow Product [Hall and Riggs, 2016, MSP] and
VIIRS Snow Product [Riggs et al., 2017, VSP] are also possible options for
comparison, since lakes in the Alps are typically snow-covered for most of the
frozen period. We cross check our results with these two snow products, see
Section 4.2.4. More details on all the mentioned products can be found in
Table 7 (A).

###  1.4 Definitions used

We exclude mixed pixels and work only with pixels that lie completely inside a
lake, termed clean pixels. Non-transition dates are the days when a lake is
either completely frozen or completely non-frozen, the remaining days in a
winter season are termed transition dates. By winter season, we denote all
dates from September till May for our purposes.

##  2 Study area and data

###  2.1 Study area

We process four small-to-medium-sized Swiss Alpine lakes: Sihl, Sils,
Silvaplana and St. Moritz, see Fig. 1 and Table 8 (B). While most of the
earlier works [Gou et al., 2015, Qi et al., 2019, 2020, Yao et al., 2016] on
long time-series monitoring of lake ice with MODIS concentrated on larger
lakes, many lakes that freeze are actually small- or medium-sized mountain
lakes, especially outside the (sub-)Arctic regions. The lakes we analyse are
relatively small in area (0.78 - 11.3 $km^{2}$), representative for this
category.

Figure 1: MODIS orthophoto map (RGB composite, red: band 1, green: band 4,
blue: band 3) of Switzerland (left) captured on 7 September 2016. Red and
amber rectangles show the regions Einsiedeln (around lake Sihl) and
Oberengadin (with lakes Sils, Silvaplana and St. Moritz, from left to right)
respectively. Inside each zoomed rectangle on the right, the respective lake
outlines are shown in light green and the nearest meteorological stations
(EIN, SIA and SAM) are marked using pins.

The lakes were chosen according to the needs of a project of the Swiss GCOS
office. For the three small lakes in the region Oberengadin (Sils, Silvaplana,
St. Moritz), located at an altitude $>\\!$ 1750 m, there are long in-situ
observation series (important for climate studies), and they are also included
in the NSIDC lake ice database, although not updated recently. The fourth lake
(Sihl) from the region Einsiedeln has been chosen mainly to test generality,
as it is relatively larger, lies at a lower altitude on the North slope of the
Alps, and has different environmental conditions, see Table 8. The three lakes
in Oberengadin fully freeze every year, whereas lake Sihl does not (but still
freezes in most winters).

For these four lakes there is no reference freeze/thaw data available from the
past two decades. Hence, we study the weather patterns in the regions near the
lakes. For each lake, the temperature and precipitation data recorded at the
nearest meteorological stations are shown in Fig. 2. On the top and bottom
rows, we plot the mean temperature and total precipitation (solid curves)
during the winter months on the $y$-axis, against the winters on the $x$-axis
in a chronological order. Additionally, in both rows, we plot the linearly
fitted trend curve (dotted line) for each station.

Figure 2: For the region near the lakes, the mean winter temperature (row 1)
and total winter precipitation (row 2) are plotted (solid curve) on the y-axis
against the winters shown in a chronological order on the x-axis. The data
from nearest meteorological stations: EIN (Sihl), SIA (Sils and Silvaplana)
and SAM (St. Moritz) are used. The corresponding trends (linear fit, dotted
curve) are also shown with the same colour. Data courtesy of MeteoSwiss.

For lakes Sihl and St. Moritz, the nearest meteorological stations are
Einsiedeln (EIN) and Samedan (SAM) respectively as shown in Fig. 1 (see also
Table 8). Lakes Sils and Silvaplana are located next to each other and hence
share the same meteorological station: Segl Maria (SIA). The station
Einsiedeln is located at a relatively lower altitude and closer to the Swiss
plateau, while the other two stations are in the Engadin valley at higher
altitudes. This explains why the absolute temperature is relatively higher for
Einsiedeln. Additionally, the stations SIA and SAM are located within 20 km
from one another, see Fig 1, and hence have similar temperature and
precipitation values. Exceptionally high winter temperature was recorded at
all three stations in winter 2006–07.

It can be seen from Fig. 2 that during the past 20 winters, at all the three
stations, the mean temperature follows an increasing trend. On the other hand,
precipitation (snow+rain) has a decreasing pattern. While Meteoswiss (MCH) has
reported a significant trend of temperature increase in the Swiss Alps since
1864, they have so far not confirmed a significant precipitation trend
(https://www.meteoswiss.admin.ch/home/climate/climate-change-in-
switzerland/temperature-and-precipitation-trends.html).
Over the shorter period of the past 20 winters, precipitation has been
slightly declining. Warmer winters at higher altitudes in Switzerland could be
linked to a decrease in precipitation, see Rebetez [1996]. The pattern of
precipitation over the 20-year period differs somewhat between Einsiedeln and
the two other (similar) stations, e.g., see the winters 08–09, 12–13.

###  2.2 Data

In our analysis, we use the data from Terra MODIS and Suomi NPP VIIRS
(https://ncc.nesdis.noaa.gov/VIIRS/) satellites downloaded from the LAADS
(https://ladsweb.modaps.eosdis.nasa.gov) and NOAA
(https://www.avl.class.noaa.gov/) databases, respectively. Due to low data
quality, we did not include Aqua MODIS (https://aqua.nasa.gov/modis) imagery.
For MODIS processing, we downloaded the MOD02 (geolocated and calibrated
radiance, level 1b, TOA), MOD03 (geolocation) and MOD35_L2 (cloud mask)
products and pre-processed them using the MRTSWATH
(https://lpdaac.usgs.gov/tools/modis_reprojection_tool_swath/, re-projection
and re-sampling) and LDOPE [Roy et al., 2002, cloud mask] software. For VIIRS,
we downloaded the Scientific Data Record (SDR) data for the imagery bands,
IICMO and VICMO products for the cloud masks, and GITCO (for image bands) and
GMTCO (for cloud masks) for terrain corrected geolocation. VIIRS pre-
processing is done using the following software packages: SatPy
(https://satpy.readthedocs.io/) for assembling the data granules, mapping and
re-sampling, H5py (https://www.h5py.org) for cloud mask extraction, PyResample
(https://resample.readthedocs.io) and
GDAL (https://gdal.org) for re-sampling of cloud masks. As in Tom et al.
[2020b], we use only twelve (five) selected MODIS (VIIRS) bands, chosen with
the supervised XGB feature selection algorithm. The spectral coverage of these
bands is shown in Fig. 14 (C) where MODIS (M) and VIIRS (V) bands are
displayed as blue and red bars, respectively. For each band, the bar width is
proportional to the corresponding bandwidth.

We analyse MODIS data from all winters since 2000–01 (20 winters), and VIIRS
data from all winters since 2012–13 (8 winters). In each winter, we process
all the dates from the beginning of September until the end of May on which at
least 30% of a lake is free of clouds, according to the mask. Fig. 3 displays
more details of the MODIS and VIIRS data that we use as a stacked bar chart
(one colour per lake). For all target lakes, the total number of cloud-free,
clean pixels in each winter is shown on the $y$-axis, against the winters in
chronological order on the $x$-axis.

(a)

(b)

(c)

Figure 3: First row displays the clean, cloud-free pixels (from transition and
non-transition dates) from the four target lakes (Sihl, Sils, Silvaplana, St.
Moritz) used in our experiments. Data from both MODIS (20 winters, 4 lakes)
and VIIRS (8 winters, 3 lakes) is displayed. Second row shows the percentage
of at least 30% Non-Cloudy (NC) days during each winter (derived from the
MODIS cloud mask). Winter 00–01 represents the dates from September 2000 till
May 2001 (and similarly other winters).

In Fig. 3, note that some winters are relatively less cloudy. This explains
why the number of pixels vary across winters, even for the same lake and
sensor. On average there are $16558$, $5899$, $3972$ and $329$ MODIS pixels
per winter for the lakes Sihl, Sils, Silvaplana and St. Moritz respectively.
Similarly, there are $5538$, $1901$ and $1606$ VIIRS pixels for the lakes
Sihl, Sils and Silvaplana respectively. Due to its small size there exist no
clean pixel for St. Moritz in VIIRS imagery bands [Tom et al., 2020b], hence
we exclude it from the VIIRS analysis. The number of pixels per acquisition is
proportional to the lake area and hence varies across lakes even for a given
sensor. Additionally, for a given lake, the number of VIIRS pixels per
acquisition is lower compared to MODIS, due to the higher GSD of VIIRS imagery
bands ($\approx 375m$) compared to MODIS ($250$m). We super-resolved all low
resolution MODIS bands ($500$m, $1000$m) to $250$m using bilinear
interpolation prior to the analysis. This step is not required for VIIRS as
all used bands have the same GSD.

It can be inferred from Fig. 3 that in both VIIRS and MODIS plots, the cloud
patterns of lakes Sils, Silvaplana and St. Moritz are quite similar, due to
geographical proximity (see also Fig. 1). Minor differences exist (in few
winters) between the two very nearby lakes Sils and Silvaplana due to cloud
mask errors (see also Section 3.1). Lake Sihl has a different cloud coverage
pattern compared to other three lakes, due to its lower altitude and different
surrounding topography. Both MODIS and VIIRS have daily temporal resolution,
but the data capture can happen at different times within a day. Consequently,
the cloud pattern (and hence the cloud masks) can differ even between MODIS
and VIIRS acquisitions from the same day. For all four lakes, Fig. 3c shows
the percentage of non-cloudy days (at least 30% cloud-free pixels) in each
winter season. It can be seen that $\approx$40 to $\approx$60% of all days are
not usable due to clouds, significantly reducing the effective temporal
resolution. We notice that the data loss is worse for lake Sihl located near
the plateau.

####  2.2.1 Ground truth

We use the same ground truth as in Tom et al. [2020b], which is based on
visual interpretation of freely available high-resolution webcams monitoring
the target lakes. One label (fully frozen, fully non-frozen, partially frozen)
per day is assigned. Two different operators looked at each image, i.e., a
second expert verified the judgement of the first operator to minimise
interpretation errors. When deciphering a webcam image was difficult,
additional images were used from other webcams viewing the same lake (if
available), images from the same webcam but at other acquisition times on the
same day, and images of the same webcam for the days before and after the
given observation day. We also improved the webcam-based ground truth using
sporadic information available from media reports, and by visually
interpreting Sentinel-2 images, whenever available and cloud-free. No webcam
data is available from the winters before 2016-17. Moreover, the manual
interpretation process is very labour intensive. Thus, ground truth is
available only for two winters (2016–17, 2017–18).

Even though visual interpretation is the standard practice, a certain level of
label noise inevitably remains in the ground truth, due to factors such as
interpretation errors, image compression artefacts, large distance and flat
viewing angle on the lake, etc. Furthermore, the webcams used are not
optimally mounted for lake ice monitoring application and hence do not always
cover the full lake area (or even a major portion of it), even for the
smallest lake St. Moritz. Still, the ground truth serves the purpose, in the
sense that it has significantly fewer wrongly labelled pixels than the
automatic prediction results. We see no possibility to obtain more accurate,
spatially explicit ground truth for our task.

##  3 Methodology

###  3.1 Pre-processing

We perform the same pre-processing steps as in Tom et al. [2020b]. First, the
absolute geolocation error for both sensors (0.75, respectively 0.85 pixels
$x$\- and $y$-shifts for MODIS; 0.0, respectively 0.3 pixels $x$\- and
$y$-shifts for VIIRS) are corrected. The lake outlines are then backprojected
onto the images to extract the clean pixels. Binary cloud masks are derived
from the respective cloud mask products to limit the analysis only to cloud-
free pixels. We noticed some errors in both MODIS and VIIRS cloud masks. The
most critical ones are false negatives, where an actually cloudy pixel goes
undetected. Such cases can corrupt model learning and inference and introduce
errors in the predicted ice maps.

###  3.2 Machine learning for lake ice extraction

We model lake ice detection in optical satellite images as a per pixel 2-class
(frozen, non-frozen) supervised classification problem. For each pixel, the
feature vector is formed by directly stacking the 12 (5) bands of MODIS
(VIIRS), see Fig. 14. We treat snow-on-ice and snow-free-ice as a single
class: frozen. Class non-frozen denotes the water pixels. The class
distributions in winters 2016–17 and 2017–18 are shown in Fig. 15 (D). There
is a significant class imbalance in our dataset, since we include all cloud-
free dates from September till May, of which only a minority is frozen. We
have tested four off-the-shelf ML classifiers: linear SVM [SL], SVM with
Radial Basis Function (RBF) kernel [SR], Random Forest [RF] and eXtreme
Gradient Boosting [XGB] to perform pixel-wise supervised classification in
order to predict the state of a lake.

We recall that SVM is a linear large-margin classifier, which can be extended
to non-linear class boundaries with the kernel trick. The choice of kernel is
critical and depends on the data distribution. In our case we tested both
linear and RBF kernels. RF is an ensemble learning approach which relies on
bagging (bootstrap aggregation) of multiple decision trees constructed from
the data, using randomisation to decorrelate the individual trees. XGB is also
an ensemble method based on (shallower) decision trees, but iteratively learns
the trees with gradient descent, such that each tree corrects error of earlier
ones. XGB is highly scalable and exploits sparsity. We note that several
comparison studies exist in the literature for other applications than ours,
e.g., [Ogutu et al., 2011, Pham et al., 2020, Wainer, 2016].

Note also, while there recently has been a strong interest in deep learning
for remote sensing tasks, it is not suitable for our particular application,
due to the scarcity of pixels with reliable ground truth. The lakes that we
monitor are small and ground truth is available for two winters (see Section
2.2.1 and Fig. 15), which is too little to train data-greedy deep neural
networks. Also, given the large GSD and limited need for spatial context, we
do not expect deep models to greatly outperform our shallower ones, see
Section 4.

###  3.3 LIP estimation

In each winter, using the trained ML model, we process all available
acquisitions that are at least $30\%$ cloud-free and generate pixel-wise
classification maps (one per acquisition). To recover the temporal evolution
(per winter), the percentage of non-frozen pixels is computed from each
classification map and is plotted on the $y$-axis against the acquisition time
on the $x$-axis. Then, as in Tom et al. [2020b], multi-temporal smoothing is
performed as a post-processing step using a Gaussian kernel with standard
deviation $0.6$ days and window width $3$ days. An example MODIS results
timeline for lake Sils from winter 2006–07 is shown in Fig. 4a. Results from
different months are displayed with different colours, see the legend.

(a)

(b)

Figure 4: Piece-wise linear (”U with wings”) curve fitting example. NF
indicates Non-Frozen.

In the post-processed timeline, we find all the potential candidates for the
following four critical dates: FUS, FUE, BUS and BUE, see Table 1 for the
corresponding definitions.

Event | Definition
---|---
FUS | 30% or more of the non-cloudy portion of the lake is frozen
|  and the just previous non-cloudy day should be $<30\%$ frozen
FUE | 70% or more of the non-cloudy portion of the lake is frozen
|  and the just previous non-cloudy day should be $<70\%$ frozen
BUS | 30% or more of the non-cloudy portion of the lake is non-frozen
|  and the just previous non-cloudy day should be $<30\%$ non-frozen
BUE | 70% or more of the non-cloudy portion of the lake is non-frozen
|  and the just previous non-cloudy day should be $<70\%$ non-frozen
ICD | BUE - FUS
CFD | BUS - FUE
Table 1: Key LIP events.

Within a winter, it is possible that $>\\!$ 1 candidates exist per critical
date which all satisfy the respective definition. In order to weed out some
obviously spurious candidates, we enforce the constraint that the four dates
must occur in the following chronological order:
FUS$\rightarrow$FUE$\rightarrow$BUS$\rightarrow$BUE. Then we exhaustively
search for the optimal set of four dates among the remaining candidates. To
that end, we fit a continuous, piece-wise linear ”U with wings” shape to the
per-day values of percentage of non-frozen pixels, such that the fitting
residuals $z$ are minimised (see example fit in Fig. 4b, shown in cyan
colour). In detail, the loss function for the fit is defined as:

$L_{LIP}=\frac{1}{P}\cdot\sum_{i=1}^{N}H_{\phi}(z)$ (1)

where $N$ is the total number of available acquisitions that are at least
$30\%$ cloud-free.

$H_{\phi}(z)=\begin{cases}z^{2}&\mathopen{|}z\mathclose{|}\leq\phi\\\
2\phi\mathopen{|}z\mathclose{|}-\phi^{2}&\mathopen{|}z\mathclose{|}>\phi\end{cases}$
(2)

is the Huber norm of the residual. For the shape parameter $\phi$, we use a
constant value of 1.35 which offers a good trade-off between the robust
$l_{1}$-norm for large residuals and the statistically efficient $l_{2}$-norm
for small residuals [Owen, 2006].

Per lake, we assume that each critical date occurs only once per winter, which
is always true in Oberengadin. Lake Sihl does not always fully freeze. As it
lies outside of the target region and is included mostly to ensure generality
of the ice classifier, we do not extract the LIP events for Sihl. Moreover, we
decide to exclude lake St. Moritz since it is too small for the GSD of MODIS
(only 4 clean pixels), making the fraction of frozen pixels overly susceptible
to noise. We thus prefer to study only the two main lakes in Oberengadin, Sils
and Silvaplana, in terms of long-term lake ice trends. These two lakes fully
freeze every year and typically have a single freeze-up and break-up period.
To further stabilise the LIP estimates we include a weak prior for each
phenological date, in the form of a diffuse Gaussian distribution.

The prior probability ($P$) is given by:

$P=P_{fus}\cdot P_{fue}\cdot P_{bus}\cdot P_{bue}$ (3)

where $P_{fus}$, $P_{fue}$, $P_{bus}$, and $P_{bue}$ are Gaussian normal
distributions for the events FUS, FUE, BUS and BUE, respectively. The prior
formalises the knowledge that freeze-up normally occurs around the end of
December and takes around three days, and break-up occurs around the end of
April over a similar period, for both target lakes. In order not to bias the
estimation, but only to minimise the risk of implausible results, we choose
very wide Gaussians ($\sigma$ = 1 month). Furthermore, we impose a constraint
that the duration of freeze-up (FUE-FUS) and break-up (BUE-BUS) is not more
than two weeks.

Inspired by Qi et al. [2020], we additionally compute further indicators that
can be derived from the four critical dates, namely ICD and CFD, see Table 1
for details. We use $30\%$ as the threshold to estimate the four dates. For
example, a date is considered as FUS candidate if 30% or more of the non-
cloudy portion of the lake is frozen. Some studies based on MODIS [Qi et al.,
2020, Reed et al., 2009, Yao et al., 2016] have used 10% as threshold, while
another approach [Kropáček et al., 2013] even employed 5%. All of them
monitored larger lakes (45 to 4294 $km^{2}$ in area). We empirically found
that for our comparatively tiny lakes the above thresholds are too strict and
a threshold of 30% is needed to ensure reliable decisions. To see why,
consider that in the best case (Sils, cloud-free) a lake has 33 clean pixels,
but that number can go down to as few as 7 (Silvaplana at 70% cloud cover).
Note also that on such small lakes a large portion of all pixels is very close
to the lake’s shoreline, where the residual absolute geolocation error (in the
worst case 0.5 pixel) may have a significant impact.

##  4 Results

In our experiments, we perform a comparison of the performance of various ML
classifiers, derive a 20-year time-series of the key phenological dates, and
perform correlation of our results with the regional weather trends. In
addition to the mean classification accuracy (mAcc), we report a stricter
measure: mean intersection-over-union (mIoU) for fair reporting of the results
on a dataset with imbalanced class distribution. More details can be found in
the following sub-sections.

###  4.1 Choice of machine learning method

We first conduct experiments on the two winters $2016$–$17$ and $2017$–$18$,
for which ground truth is available. These experiments serve to compare the
performance of the SVM, RF, and XGB classifiers and determine which one is
most suitable for our task and dataset. We train the models and report
quantitative results only on non-transition dates, since per-pixel ground
truth is not available for the transition dates.

####  4.1.1 Hyper-parameter tuning

In SVM, the two main hyper-parameters are cost and gamma (for the linear
variant just the cost) which indicate misclassification cost and kernel width
respectively. Higher values of cost means the model chooses more support
vectors which effectively increases the variance and decreases the bias and
results in overfitting and vice versa. When gamma is set to a high value, more
weight is assigned to the points close to the hyperplane compared to the far
away ones. Higher gamma values also result in low bias and high variance thus
causing overfitting. For RF, num_trees represents the number of trees in the
forest while num_variables denotes the number of predictors (features) to
select at random for each split. Note that the minimum_leaf_size is set to $1$
in all our RF experiments. In XGB, the three most important hyper-parameters
are number of trees (num_trees), learning_rate and the tree specific
parameter: colsample_bytree. Learning rate shrinks the contribution of each
new tree to make the boosting procedure more conservative and thus the
resulting model more robust. Colsample_bytree controls the fraction of
features (spectral bands) to be used in each boosting iteration.

For each ML approach, the optimal hyper-parameters are first independently
determined with a grid search. The best-performing hyper-parameters thus
chosen are shown in Table 2. We use these parameters in all further
experiments. However, note that the parameters are dataset-dependent, too.
Overall, RF, XGB and SL exhibited fairly stable results across a range of
hyper-parameters, whereas SR was very sensitive.

Method | Sensor | Hyper-parameters
---|---|---
SL | MODIS |  cost $0.1$
SR | MODIS |  cost $10$, gamma $1$
RF | MODIS |  num trees $500$, num variables $10$
XGB | MODIS |  num trees $1000$, colsample_bytree $1$, learning rate $0.2$
SL | VIIRS |  cost $0.1$
SR | VIIRS |  cost $10$, gamma $1$
RF | VIIRS |  num trees $500$, num variables $3$
XGB | VIIRS |  num trees $500$, colsample_bytree $1$, learning rate $0.3$
Table 2: Optimum hyper-parameters of each classifier and sensor estimated
using grid search.

####  4.1.2 Four-fold cross validation

As the first experiment, we combine the data (independently for MODIS and
VIIRS) of all the available lakes from winters 2016–17 and 2017–18 and perform
$4$-fold cross validation and report the mAcc and mIoU, see Table 3. For both
MODIS and VIIRS, the performance of all classifiers except SL is more or less
the same, with accuracy differences below $1\%$. While SR performs marginally
best on MODIS data, XGB scores well on VIIRS. Though lower than other three
classifiers, the performance of SL is also very good on both sensors.

Sensor | Feature vector | Method | mAcc | mIoU
---|---|---|---|---
MODIS | All 12 bands | SL | $93.4$ | $83.9$
MODIS | All 12 bands | SR | 99.4 | 98.5
MODIS | 10 bands (random) | RF | $98.9$ | $97.2$
MODIS | All 12 bands | XGB | $99.3$ | $98.3$
VIIRS | All 5 bands | SL | $95.1$ | $88.4$
VIIRS | All 5 bands | SR | $97.1$ | $93.1$
VIIRS | 3 bands (random) | RF | $97.6$ | 94.5
VIIRS | All 5 bands | XGB | 97.7 | 94.5
Table 3: Four-fold cross-validation results (in $\%$) on MODIS and VIIRS data.
The data from both winters 2016–17 and 2017–18 are used in this analysis. Mean
classification accuracy (mAcc) and mean intersection-over-union (mIoU) scores
are shown. The best results are shown in bold.

####  4.1.3 Generalisation experiments

In order to study how well the classifiers generalise across space and time,
we train a model on all except one lake (respectively, winter) and test on the
held-out lake (winter).

Fig. 5 displays the results (bar graphs showing mAcc and mIoU) of the four
classifiers for leave-one-lake-out setting on MODIS (top row) and VIIRS
(bottom row) data. It can be seen that the performance varies across lakes and
sensors. For MODIS, XGB has, on average, a narrow advantage over SR and RF,
with SL a bit behind. For VIIRS, on the other hand, RF marked the best
performance closely followed by SL, SR and XGB.

Sihl Sils Silv Moritz Mean
$70$$80$$90$$100$$92.6$$93$$90.5$$82.5$$89.7$$88.8$$93.9$$94.1$$90.2$$91.8$$88.8$$95.4$$94.4$$87.8$$91.6$$88.1$$96.1$$94.4$$88.9$$91.9$mAcc

Sihl Sils Silv Moritz Mean
$70$$80$$90$$100$$78$$86$$82.2$$69$$78.8$$71$$88.3$$88.7$$81.9$$82.5$$70.7$$91.1$$89.2$$77.7$$82.2$$69.5$$92.3$$89.2$$79.6$$82.7$mIoU

Sihl Sils Silv Mean
$70$$80$$90$$100$$94.7$$94.4$$92.8$$94$$91.2$$95.2$$94.2$$93.5$$93.7$$94.8$$93.7$$94.1$$92$$94.6$$93.7$$93.4$mAccSLSRRFXGB

Sihl Sils Silv Mean
$70$$80$$90$$100$$84.8$$89.3$$86.5$$86.9$$77.4$$90.8$$88.9$$85.7$$82.7$$90.1$$88.1$$87$$79.1$$89.6$$88.2$$85.6$mIoUSLSRRFXGB

Figure 5: Generalisation across lakes results on MODIS (top row) and VIIRS
(bottom row) data for the classifiers SL, SR, RF and XGB on lakes Sihl, Sils,
Silvaplana (Silv) and St. Moritz (Moritz). Both mAcc (left column) and mIoU
(right column) are shown.

On both sensors, the best performance (especially in terms of mIoU) is
achieved for the lakes Sils and Silvaplana. This is likely due to them having
the most similar characteristics and imaging conditions, see Table 8. I.e.,
pixels from one of them are representative also of the other one, such that
the classifier trained in one of the two generalises well to the other. Lake
St. Moritz (only for MODIS) has too few clean pixels per acquisition to draw
any conclusions about generalisation. However, we still include it in our
processing to study how far lake ice monitoring with MODIS can be pushed (in
terms of lake area) – indeed, the classification is $>\\!$ 82.5% correct. Lake
Sihl from the region Einsiedeln is different compared to other three lakes
from the region Engadin in terms of area, weather, surrounding topography
etc., c.f. Section 2.1. Hence, the performance on lake Sihl is interesting to
assess geographical generalisation over longer distances. It can be seen in
Fig. 5 that for lake Sihl the SL classifier performs best on both MODIS and
VIIRS data, suggesting that the other (non-linear) models already overfit to
the specific local conditions of Oberengadin.

As a second generalisation experiment, more important for our time-series
analysis, we check how well the trained classifiers can be transferred across
different winters. We train on one winter and test the resulting model on the
held-out winter (leave-one-winter-out), see Fig. 6. We only have data for two
consecutive winters (2016–17, 2017–18) to base the analysis on. Still, we
believe that the experiment is representative for generalisation to unseen
years, since the weather conditions in different years are largely
uncorrelated (c.f. Fig. 2). In particular for the two available winters,
2017–18 was markedly colder than the previous year, see Fig. 2.

Fig. 6 (top row) clearly shows for MODIS that SL adapts best to a new winter,
with significantly higher generalisation losses for XGB, RF, and SR. For VIIRS
(Fig. 6, bottom row), SL also copes best with generalisation across different
years. Here RF can keep up, whereas SR and XGB again suffer higher
generalisation losses.

$2016$-$17$ $2017$-$18$ Mean
$70$$80$$90$$100$$90.6$$93.6$$92.1$$88.8$$85.9$$87.4$$87.1$$89.3$$88.2$$89.5$$89.2$$89.4$mAcc

$2016$-$17$ $2017$-$18$ Mean
$70$$80$$90$$100$$79.6$$82.6$$81.1$$76.1$$67.6$$71.9$$72.9$$74.2$$73.6$$77.7$$74.1$$75.9$mIoU

$2016$-$17$ $2017$-$18$ Mean
$70$$80$$90$$100$$92.2$$92.8$$92.5$$91.3$$91.8$$91.6$$93$$91.9$$92.5$$92.8$$89.1$$91$mAccSLSRRFXGB

$2016$-$17$ $2017$-$18$ Mean
$70$$80$$90$$100$$83.7$$81.7$$82.7$$82.3$$79.4$$80.9$$85.4$$79.8$$82.6$$85.1$$74.8$$80$mIoUSLSRRFXGB

Figure 6: Generalisation across winters results on MODIS (top row) and VIIRS
(bottom row) data for the classifiers SL, SR, RF and XGB. Both mAcc (left
column) and mIoU (right column) are shown.

Overall, all classifiers exhibit a certain performance drop when having to
generalise beyond the exact training conditions. Table 4 shows the detailed
performance drops w.r.t. Table 3. Since mIoU is a stricter measure than mAcc,
the drops are more pronounced. It is interesting to note that in both
generalisation experiments the most robust (least affected) classifier is the
(linear) SL. With all other classifiers the drop is much higher for unseen
winters than for unseen lakes. We conclude that SL is the safest option four
our task, where data from multiple different winters must be processed, and
use it for all further experiments. Recall, though, that we only have a
relatively small dataset at our disposal from few small mountain lakes over
two winters. It is quite possible that the small volume and specific
geographical conditions aggravate the tendency to overfit, and that a higher-
capacity, non-linear classifier will work best if a larger and more diverse
dataset were available, or if the conditions were less variable (large lakes
in smooth terrain).

Sensor | Loss type | SL | SR | RF | XGB
---|---|---|---|---|---
MODIS | across lakes |  3.7/5.1 |  7.6/16.0 |  7.3/15.0 |  7.4/15.8
MODIS | across winters |  1.3/2.8 |  12.0/26.6 |  10.7/23.6 |  9.9/22.4
VIIRS | across lakes |  1.1/1.5 |  3.6/7.4 |  3.5/7.5 |  4.3/8.9
VIIRS | across winters |  2.6/5.7 |  5.5/12.2 |  5.1/11.9 |  6.7/14.5
Table 4: Generalisation loss (across lakes / winters) encountered by each
classifier: SL, SR, RF and XGB. Drop (in % points) for mAcc and mIoU are shown
in black and grey respectively.

###  4.2 Experiments on MODIS data from 20 winters

####  4.2.1 Test-train split

So far, we have used only on parts of the available ground truth for training,
so as to evaluate the method. We now move on to the actual longer-term
analysis, where we process MODIS data from all the 20 winters since 2000–01
(inclusive). Details of the training set for each tested winter is shown in
Table 5. To avoid systematic biases in the estimated ice maps due to
overfitting to a particular year, we proceed as follows: we train the SL model
on all non-transition dates of 2016–2017 and use it to estimate lake ice
coverage for all days in 2017–2018 (including transition dates). We repeat
that procedure in the opposite direction, i.e., we train on all non-transition
days of 2017–2018 and perform inference for all dates of 2016–2017. Then, we
merge all non-transition dates from both winters into a new, larger _two-
winter_ training set, which we further augment with an _auxiliary dataset_.
The latter contains all acquisitions of lakes Sils and Silvaplana captured
during the remaining 18 years in September (when the lakes are never frozen)
and in February (when the lakes are always frozen). The purpose of the
auxiliary dataset is to cover a wider range of weather and lighting conditions
that might not have been encountered in the two winters with annotated ground
truth, for better generalisation. Data of lake Sihl is not included in the
auxiliary set, as it does not freeze reliably, St. Moritz is ignored due to
its negligible number of pixels. The two-winter and auxiliary datasets are
merged and used to train a SL model, which is then used to predict ice cover
maps for the 18 remaining winters.

Test set | Training set
---|---
winter 2016–17 | winter 2017–18
winter 2017–18 | winter 2016–17
winters till 2015–16, from 2018–19 | winters 2016–17 and 2017–18, auxiliary set
Table 5: Test-train split for the MODIS data from 20 winters.

####  4.2.2 Qualitative results

Exemplary qualitative results on some selected days (one per winter) of lake
Sihl are shown in Fig. 7. The respective dates are displayed below each sub-
figure. The lake outline overlaid on the MODIS band $B_{1}$ is shown in green.
Pixels detected as frozen and non-frozen are shown as blue and red squares,
respectively. The results include fully-frozen, fully non-frozen and partially
frozen days.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

(o)

(p)

(q)

(r)

(s)

(t)

Figure 7: MODIS qualitative classification results (overlaid on band $B_{1}$
from the respective day) for lake Sihl on selected dates from the past 20
winters using the SL classifier. Blue and red squares are overlaid on the
pixels detected as frozen and non-frozen respectively.

####  4.2.3 Additional check using VIIRS data from 8 winters

A direct quantitative analysis is not possible, since no ground truth is
available for 18 out of the 20 winters. In order to validate our MODIS
results, we additionally process the VIIRS data from 8 winters (since winter
2012–13, inclusive) and compare the results to MODIS. Since a pixel-to-pixel
comparison is not straight-forward due to different GSDs, we fit the timelines
per winter for each lake as described before (Fig. 4a) and compute absolute
differences (AD) between the daily estimates for the percentage of frozen
pixels. The AD is computed only on dates when both MODIS and VIIRS
acquisitions are present, and when the lake is at least 30% cloud-free. The
ADs are then further aggregated to obtain a Mean Absolute Difference (MAD) for
each winter. Fig. 8a shows, for each lake, the mean and standard deviation of
the MAD across the 8 common winters. The low mean values (3.5, 5.8 and 4.3
percent respectively for Sihl, Sils and Silvaplana) show that our MODIS and
VIIRS results are in good agreement, especially considering that a part of the
MAD is due to the difference in GSD between MODIS (250m GSD) and VIIRS
($\approx$375m GSD). Note also that, the acquisition times during the day (and
hence the cloud masks) can differ; and that, although the absolute geolocation
has been corrected for both sensors, errors up to 0.5 pixel can still remain
[Aksakal, 2013] and affect the selection of clean pixels near the lake shore.

(a)

(b)

(c)

(d)

Figure 8: Row 1 shows the comparison of our MODIS and VIIRS results for the 8
common winters (left sub-figure) and comparison of our MODIS (20 winters) and
VIIRS (8 winters) results with the respective snow products (right sub-
figure). Row 2 shows the deviations between the two snow products and our
webcam-based ground truth. Row 3 displays per-winter MAD (MODIS) for each
lake.

####  4.2.4 Comparison with MODIS and VIIRS snow products

We compare our MODIS (20 winters) and VIIRS (8 winters) results to the
respective operational snow products: MSP (collection 6, MOD10A1), VSP
(collection 1, VNP10A1F). For the regions of interest, VSP has some data gaps,
hence the comparison is done whenever VSP data is available. For actual snow
cover mapping, errors of 7-13% have been reported for MSP [Hall and Riggs,
2016]. Our findings are in line with this: for the two 2016–17 and 2017–18
(non-transition days only) we observe an error of 14% w.r.t. our ground truth,
see Fig. 8c.

Similar to the evaluation methodology explained in Section 4.2.3, for each
lake, we first estimate the percentage of frozen pixels per day using our
MODIS and VIIRS results. Since a pixel-to-pixel registration is difficult in
the presence of absolute geolocation shifts and/or GSD differences, the daily
percentage of frozen pixels is also computed from the snow products and the
MAD is estimated for each winter. See Fig. 8d for our MODIS results vs MSP
comparison. For the three lakes, the per-winter MAD is shown on $y$-axis
against the winters on $x$-axis. We again exclude lake St. Moritz because of
its minuscule area.

Overall, the 20-year time-series inter-comparison (per-lake mean and standard
deviation of MAD, Fig. 8b) does not suggest large, systematic inconsistencies.
On average, our MODIS and VIIRS results deviate by mean MAD values of 14-18%
and 12-19% respectively. These deviations are only a little higher than the
estimated error of the snow products, and relatively stable across different
years.

It is important to point out that the snow products are an imperfect proxy for
lake ice, because a lake can be frozen but not snow-covered, especially near
freeze-up when it has not yet snowed onto the ice. Also, mixed ice and water
cases go undetected in MSP [Hall and Riggs, 2016]. Fig. 8c shows that the snow
products are less consistent with the manually annotated ground truth than our
ice maps. In fact, most deviations between our estimates and the snow products
occur around the transition dates, mostly freeze-up. Additionally, MSP and VSP
use a less conservative cloud mask than we do (accepting not only confident
clear and probably clear, but also uncertain clear as cloud-free). Despite
these issues, the inter-comparison provides a sanity check for our results.
For completeness, we note that our algorithm has similar issue and thin ice is
sometimes confused with water: firstly, snow-free ice is rare and
underrepresented in the training set. Secondly, it appears predominantly near
the transition dates (especially freeze-up) when we do not have pixel-accurate
ground truth. Thirdly, thin ice and water are difficult to distinguish, we
observed that even human interpreters at times confused them when interpreting
webcam images.

It is interesting to note that, for both sensors, the mean MAD is inversely
proportional to the lake area (see Fig. 8b). This hints at residual errors in
the products’ geolocation, which would affect smaller lakes more due to the
larger fraction of pixels near the lake outline. Besides the $<\\!$ 0.5 pixel
inaccuracy of our maps, inaccurate geolocation of the snow products has been
reported (more for MODIS, less for VIIRS) especially for freshwater bodies,
due to uncertainties in gridding, reprojection etc. [Hall and Riggs, 2016].

####  4.2.5 LIP trends using MODIS data

As discussed in Section 3.3, we fit ”U with wings” polygon to each winter to
estimate the four critical dates: FUS, FUE, BUS and BUE. Sometimes, these
phenological dates are defined such that a second, consecutive day with
similar ice conditions is required to confirm the event. We do not enforce
this constraint, because, quite often, the days after a potential freeze-up or
break-up date are cloudy, and looking further ahead runs the risk of pruning
the correct candidates.

Using the estimated LIP dates from 20 winters, we plot their temporal
evolution for lakes Sils (top) and Silvaplana (bottom) in Fig. 9. On the
$y$-axis, all the dates from 1 December to 1 June (we skip September till
November since no LIP events were detected during these months), while on the
$x$-axis we show the winters in a chronological order. In each winter, the
non-frozen, freeze-up, frozen and break-up periods are displayed in cyan, red,
yellow and dark green colours, respectively.

(a)

(b)

Figure 9: 20 winter temporal LIP characteristics of lakes Sils and Silvaplana
estimated from MODIS using SL classifier.

It can be seen from Fig. 9 that the freeze-thaw patterns of both lakes vary
considerably across winters. For lake Sils (Silvaplana), on average, the FUS
occurred on 3 January (5 January) followed by a freeze-up period of 3 (3) days
until FUE on 6 January (8 January). Additionally, on average, the lake
remained fully frozen (CFD) for 113 (108) days until BUS on 29 April (26
April) and the break-up period lasted 1 (1) day until BUE on 30 April (27
April). The average number of days from FUS to BUE (ICD) is 117 (112).

The Oberengadin region with lakes Sils and Silvaplana is a single valley (Fig.
1 and Table 8) and hence the two have similar weather conditions. Silvaplana
is relatively deeper but has smaller area than Sils, making them comparable in
terms of volume, too. So similar LIP patterns can be expected. However, the
clouds above the lakes (especially on partly cloudy days), and the associated
cloud mask errors, can cause small differences. In winter 2016-17, the ice-on
date of the two lakes, confirmed by visual interpretation of webcams, differ
by 7 (low confidence) to 10 (medium confidence) days, see also Tom et al.
[2020b].

In most of the winters, the LIP characteristics of these lakes derived using
our approach are in agreement, see Fig. 9. However, there are some outliers
too ($>\\!$ 10 days deviation).

A notable outlier is the break-up period in winter 2009-10. For Sils
(Silvaplana), BUS and BUE were both estimated as 19 May (28 April). This drift
primarily happened because of a huge data gap due to clouds and cloud mask
errors. During the period from 28 April till 20 May, Silvaplana had $>\\!$ 30%
cloud-free MODIS acquisitions only on 28 April, 29 April, 8 May and 20 May,
and the lake was detected as non-frozen on all these dates. However, Sils had
MODIS acquisitions on 28 April, 29 April, 5 May and 19 May. On 5 May the lake
was detected as 100% frozen due to a false negative cloud mask, although
break-up had started on the two earlier dates (75%, respectively 60% frozen)
and the lake was ice-free on May 19. We also checked the Landsat-7
acquisitions on 20 April 2010 and 22 May 2010 and found that both lakes were
fully covered by snow on the former date and fully non-frozen on the latter
date. No cloud-free Landsat-7 data is available in between these two dates.
For Sils, the actual BUS probably happened on 29 April ($>\\!$ 30% non-frozen)
and BUE soon after (likely on 30 April, since the BUE of Silvaplana was
detected as 28 April and Sils was detected $<\\!$ 70% non-frozen on 29 April).
However, both dates went undetected until 19 May, because of the clouds in
combination with the maximum allowed duration of 2 weeks for the break-up.

In winter 03-04, the freeze-up periods of Sils (FUS on 1 January, FUE on 2
January) and Silvaplana (FUS and FUE on 14 January) were also detected far
apart, again due to a data gap because of clouds. Sils was estimated 68% and
90% frozen on 1 and 2 January, respectively, so they were chosen as FUS and
FUE. On lake Silvaplana, the sequence for 1-5 January was
4%$\rightarrow$13%$\rightarrow$0%$\rightarrow$21%$\rightarrow$0% frozen. Then
14 January and 21 January were both found 100% frozen, so the fitting chose 14
January as both FUS and FUE. No cloud-free MODIS data exist on the
intermediate dates 6-13 January and 15-20 January, and we could also not find
any cloud-free Landsat-7 images between 21 December 2003 and 29 January 2004
(both inclusive) to check, but could confirm 0% ice cover on 20 December and
100% cover on 30 January. Connecting all the dots, we speculate that FUS and
FUE of Silvaplana occurred soon after 5 January.

In winter 2013-14, our method asserts FUE of Sils on 26 December and of
Silvaplana on 6 January. Between those dates there were a number of partially
frozen dates, but with more ice cover for Sils than Silvaplana. Additionally,
2-5 January were cloudy, leading the fitting to chose the earlier date for the
former, but the later one for the latter. We again checked with Landsat-7 that
on 15 December both lakes were fully non-frozen, whereas on and 25 January
both lakes were fully snow-covered. There exist no cloud-free Landsat-7 image
in between these two dates to pin down the dates more accurately.

In some winters, there is almost no freeze-up and/or break-up period detected
by our algorithm. This is partly a byproduct of the relatively loose threshold
needed to estimate the initial candidates for our small lakes (see Section
3.3), bringing the start and end dates of the transition closer together; and
also influenced by frequent cloud cover during the critical transition dates
(often more than half of all days c.f. Section 3c). For instance, if a couple
of adjacent dates are cloudy during break-up (and the real BUS occurred during
one of these dates) and on the next non-cloudy day, the lake is estimated
70.1% non-frozen, then our fitting will choose this date as both BUS and BUE.

We go on to analyse the freeze-up (FUS, FUE) and break-up (BUS, BUE) patterns,
by plotting time-series of the four critical dates over the past 20 winters
for lakes Sils and Silvaplana, see Fig. 10. Here also, we plot all the dates
from 1 December till 1 June in each winter on the $y$-axis and the winters in
chronological order on the $x$-axis. Additionally, per phenological date, we
fit a linear trend. Progressively later freeze-up and earlier break-up is
apparent for both lakes.

(a)

(b)

Figure 10: 20 winter ice freeze-up (FUS, FUE) and break-up (BUS, BUE) trends
for lakes Sils and Silvaplana from MODIS estimated using SL classifier.

In each winter, we also derive the remaining LIP events (ICD, CFD) listed in
Table 1. Their trends are shown in Fig. 11, with the duration in days on the
$y$-axis and the winters in chronological order on the $x$-axis. Obviously,
ICD and CFD are decreasing for both lakes.

(a)

(b)

Figure 11: 20 winter ice duration (ICD, CFD) trends for lake Sils and
Silvaplana from MODIS estimated using SL classifier.

Quantitative trend values that we estimated for all the LIP events of Sils and
Silvaplana are shown in Table 6. As explained above, we correct obvious
failures of the automatic analysis, and set the following corrections for lake
Sils: BUS and BUE occurred on 29 April and 30 April respectively in winter
2009–10. Similarly, for Silvaplana: FUS and FUE occurred on 6 January in
winter 2003–04 and FUE occurred on 26 December in winter 2013–14. For
completeness we also fit trends without the correction – these differ only
slightly and confirm that the corrections hardly impact the overall picture.
The trend towards earlier break-up (BUS, BUE) is more pronounced than the one
towards later freeze-up (FUS, FUE), for both Sils and Silvaplana. It is
interesting to note that the decrease in freeze duration (ICD and CFD) is
stronger for the slightly smaller lake Silvaplana.

Lake | FUS | FUE | BUS | BUE | ICD | CFD
---|---|---|---|---|---|---
Sils | 0.23 | 0.31 |  -0.46/-0.47 |  -0.32/-0.34 |  -0.55/-0.57 |  -0.76/-0.78
Silvaplana |  0.45/0.37 |  0.38/0.36 | -0.51 | -0.45 |  -0.9/-0.82 |  -0.89/-0.87
Table 6: Estimated LIP trends (black). Results before manual correction of the
automatic results are shown in grey.

####  4.2.6 Correlation of LIP events with meteorological data

We have also studied the (centred and normalised) cross correlation
$\in[-1,1]$ between the LIP events (corrected version) and climate variables
such as temperature, sunshine, precipitation and wind during the 20 winters.
The results are shown in Fig. 12 for the lakes Sils (top) and Silvaplana
(bottom). Temperature (measured at 2m above ground) and precipitation data
were collected from the nearest meteorological station SIA. However, we used
the sunshine and wind measurements at station SAM, since these were not
available for the complete 20 winters time span at SIA. We did not use the
cloud information (number of non-cloudy pixels) from MODIS data as a measure
of sunshine duration, since that would ignore the evolution throughout the
day, and suffers from a non-negligible amount of cloud mask errors.

FUS FUE BUS BUE ICD CFD
$-1$$-0.75$$-0.5$$-0.25$$0$$0.25$$0.5$$0.75$$1$$0.5$$0.51$$-0.55$$-0.5$$-0.64$$-0.66$$-0.55$$-0.55$$0.41$$0.36$$0.6$$0.6$$0.2$$0.17$$-0.6$$-0.6$$-0.6$$-0.64$$-5\cdot
10^{-2}$$-0.12$$0.56$$0.58$$0.32$$0.39$$-0.21$$-0.14$$-0.1$$-5\cdot
10^{-2}$$-0.21$$-0.23$Correlation

FUS FUE BUS BUE ICD CFD
$-1$$-0.75$$-0.5$$-0.25$$0$$0.25$$0.5$$0.75$$1$$0.47$$0.43$$-0.56$$-0.51$$-0.62$$-0.63$$-0.53$$-0.52$$0.44$$0.38$$0.58$$0.6$$0.3$$0.2$$-0.55$$-0.53$$-0.57$$-0.62$$-0.22$$-0.19$$0.51$$0.53$$0.4$$0.41$$-0.17$$-0.13$$-8.15\cdot
10^{-2}$$-9.1\cdot
10^{-2}$$-0.26$$-0.22$CorrelationMWTAFDDSunshinePrecipitationWind

Figure 12: Bar graphs showing the 20 winter correlation of the LIP events
(Sils on top row, Silvaplana on bottom) with climate variables.

Mean Winter Temperature (MWT) corresponds to the air temperature (in ∘C)
averaged over the whole winter season
(September till May). As expected, Fig. 12 shows that MWT has strong negative
correlation with the freeze durations CFD and ICD, negative correlation with
the break-up events BUS and BUE, and positive correlation with the freeze-up
events FUS and FUE. We conclude that, indeed, as winters got warmer over the
past 20 years the lakes froze later and broke up earlier. For both lakes, the
relationship of MWT with CFD is shown in Fig. 13 (row 1). Further significant
correlations are displayed in E.

Figure 13: Correlation of LIP events and weather variables: MWT (∘C) and CFD
(days) are shown in first row, AFDD (∘C) and FUS in second row, total winter
sunshine (hours) and CFD in third row, total precipitation (mm) in the months
January to May (J2M) and BUE in last row. Results for lakes Sils and
Silvaplana are displayed in left and right columns respectively.

AFDD represents the cumulative sum (of daily mean temperature) on the days
with average air temperature below the freezing point ($0^{\circ}$C) in a
winter season. AFDD is a popular proxy for the thickness of ice cover [Beyene
and Jain, 2018, Qi et al., 2020]. For both Sils and Silvaplana, AFDD has
strong positive correlations with ICD and CFD, strong negative correlation
with the freeze-up events FUS and FUE, and moderate positive correlation with
ice break-up events BUS and BUE, see Fig. 12, again indicating that in colder
winters (higher AFDD) the freeze-up occurs earlier and the break-up later,
leading to longer freeze duration. The relatively weaker correlation for the
break-up indicates that freeze-up play a larger role for that event. As an
example, the correlation with FUS is shown Fig. 13, for further plots see E.

To study the effect of sunshine on LIP events, we correlate the total winter
sunshine (hours) with the freeze length events ICD and CFD, total sunshine in
the months September to December (S2D) with the freeze-up events FUS and FUE,
and the total sunshine from January to May (J2M) with the break-up events BUS
and BUE. Here, we assume that the sunshine in the months after freeze-up has
no connection with freeze-up events. Similarly, we assume that the sunshine in
the early winter months (September till December) does not affect the break-up
events. We notice strong negative correlation of the total winter sunshine
with ICD, CFD and break-up events. The more sunshine in the months near break-
up, the earlier the ice/snow melts, which also reduces the total freeze
duration. An example correlation with CFD is visualised in Fig. 13, further
significant correlations are displayed in E.

We also check the relationship between the LIP events and the total
precipitation (snow and rain) during the winter months. Similar to sunshine
analysis, we correlate the total precipitation during the months from
September till December, January till May and September till May to the
freeze-up (FUS, FUE), break-up (BUS, BUE) and freeze duration events (CFD,
ICD) respectively, see Fig. 12. Notable are the break-up events with good
positive correlation. More precipitation in the months January to May (likely
to be predominantly snow), favours later break-up, and vice-versa. Correlation
with BUE is shown in Fig. 13 and with BUS in E.

Finally, inspired by Gou et al. [2015], we also looked at the effect of wind
on the LIP events, which may also influence lake freezing. We correlated the
mean winter wind speed (km/h) with CFD and ICD, mean wind speed from September
to December with FUE and FUS, and mean wind speed from January to May with BUS
and BUE. However, we did not find any significant correlations, see Fig. 12.

###  4.3 Discussion

In any ML-based system, the variety in the training dataset has a critical
influence on the model being learnt. Our dataset consists of small lakes and
has significant class imbalance. This is a biased, but realistic scenario,
representative of mountain lakes in sub-Arctic and temperate climate zones.
For supervised classification, proper ground truth information is available
only from the winters 2016–17 and 2017–18. Using the data from these two
winters, we first performed thorough experimentation with different
classifiers. Additionally, we did an inter-comparison of the individual
performance of these ML classifiers. For the four-fold cross validation
experiment, the highest performance on MODIS and VIIRS data were reported for
SR and XGB classifiers respectively, see Table 3. However, SR, RF and XGB
suffered a significant generalisation loss compared to the SL counterpart.
I.e., non-linear classifiers tended to overfit on our relatively small
dataset. We emphasise that it does not contradict the findings of Wu et al.
[2021] for a dataset of larger lakes with many more training pixels. There, RF
and GBT performed very well. We empirically selected SL as the most suitable
classifier for our lakes of interest, and we recommend to repeat the empirical
exercise when moving to different geographical conditions. Our MODIS and VIIRS
results validate each other in a relative sense ($<5.6~{}\%$ MAD in the worst
case, lake Sils), but could be subject to a common bias. In the absence of
ground truth there is no way to assess our absolute accuracy, but as an
external check against a methodologically different mapping scheme we inter-
compared our results against the respective operational snow products. The
deviations were in the expected range (MAD $<20~{}\%$).

Besides its advantages, like dense coverage, regular time-series and
homogeneous observation conditions, satellite data also has disadvantages. The
trade-off between spatial and temporal resolution makes it difficult to
monitor smaller lakes – with 21 MODIS (9 VIIRS pixels) for lake Silvaplana and
only 4 MODIS pixels for lake St. Moritz, our study goes to the limit in that
respect. A further, often-named obstacle for optical satellite observation are
occlusions due to clouds, which significantly reduce the effective temporal
resolution and also cause irregular gaps in the time-series. These
unpredictable data gaps are particularly troublesome for ice phenology,
because the critical events occur over a short time and at times of the year
where cloud cover is frequent in sub-Arctic and mid latitudes.

Data gaps due to clouds are the main source of error in our LIP estimation,
besides cloud mask errors, confusion between water and thin/floating snow-free
ice, and quantisation effects around hard thresholds. This makes phenological
observations challenging – in particular the uncertainties of our predictions
are largest during freeze-up, because of the frequent, but short-lived
presence of snow-free-ice. Still, it appears that our classifier copes better
with the reflectance of ice than simple index-based snow products, and it is
likely that more training data and, if available, additional spectral bands
could further improve its detection.

As said before, a limiting factor for small lakes is the GSD, as decision
based on very few pixels become overly unreliable and prone to statistical
fluctuations, and even small geolocalisation errors have a large effect. Our
work is also on the challenging ends of the spectrum in terms of local weather
conditions: in a drier climate the observations would be less affected by
cloud cover (we process lakes with as little as 30% cloud-free area to obtain
sufficient temporal coverage), and fewer clouds also means fewer cloud-mask
errors.

##  5 Conclusion

In this paper, we have described time-series monitoring of ice in selected
Swiss lakes, based on MODIS optical satellite images. On the one hand, we have
tried to push the limits of MODIS data for the analysis of small-to-medium
sized lakes, and have shown that even for such high-Alpine lakes it is
possible to derive meaningful correlations between the 20-winter LIP trends
and climate data. On the other hand, we have confirmed that a dedicated ML
scheme maps lake ice more accurately than classical index- and threshold-based
approaches.

As expected, our results point towards later freeze-up (FUS at a rate of 0.23
d/a for lake Sils, respectively 0.45 d/a for Silvaplana and FUE at a rate of
0.31 d/a for lake Sils, respectively 0.38 d/a for Silvaplana), earlier break-
up (BUS: -0.46 d/a for lake Sils, respectively -0.51 d/a for Silvaplana and
BUE: -0.32 d/a for lake Sils, respectively -0.45 d/a for Silvaplana) and
decreasing freeze duration (ICD: -0.55 d/a for lake Sils, respectively -0.9
d/a for Silvaplana and CFD: -0.76 d/a for lake Sils, respectively -0.89 d/a
for Silvaplana). We also observed significant (but not surprising)
correlations with climate indicators such as temperature, sunshine and
precipitation.

Our ML approach is generic and easy to apply to other sensors beyond MODIS and
VIIRS (given training data). Importantly, the VIIRS is projected to ensure
continuity well into the future, opening up the possibility to establish an
even longer time-series. One solution for the cloud issues of optical
satellites is to complement/replace them with radar observations, e.g.
Sentinel-1 SAR. We have done preliminary research in this direction [Tom et
al., 2020a]. SAR-optical data fusion holds great promise, particularly in view
of the GCOS requirement to monitor lake ice at $\pm$2 days temporal
resolution.

In the future, it will be interesting to cross-check and validate existing
temperature-based models like Hendricks Franssen and Scherrer [2008] against
our data-driven results. Beyond a simple inter-comparison, temperature
measurements could be used to eliminate gross prediction errors, and to bridge
temporal gaps in optical satellite image-based predictions that occur due to
clouds or fog.

We expect that ML-based ice detection itself could still be further improved
with additional training data. I.e., pixel-accurate annotations during
transition dates, as well as for more winters and a wider variety of lakes.
Unfortunately, gathering such data is not only a considerable, tedious effort,
but also poses its own challenges. In most locations and for older data, no
corresponding webcam data (or similar regular photography) is available; even
when available, its coverage is almost invariably incomplete; and even with
usable webcam and satellite imagery, manual annotation is not trivial and
prone to mistakes exactly in the situations that are most critical also for
computational analysis (such as thin, black ice). We speculate that, given the
enormous archive of unlabelled satellite data, approaches such as
unsupervised, semi-supervised or active ML may be applicable and could improve
the detector.

Yet another interesting research direction is to close the gap between
knowledge-driven, model-based top-down models of lake ice formation and data-
driven, bottom-up machine intelligence. Introducing expert knowledge about ice
formation and associated physical constraints into ML models could also reduce
the need for training data, and get the best of both worlds in terms of
accuracy as well as interpretability of the model. How to best bridge the gap
between statistical ML models and physical process models is an open question
and an active research direction in the Earth sciences and beyond.

##  6 Acknowledgements

This work is part of the project Integrated lake ice monitoring and generation
of sustainable, reliable, long time-series funded by the Swiss Federal Office
of Meteorology and Climatology MeteoSwiss in the framework of GCOS
Switzerland. We are grateful to Damien Bouffard (EAWAG, Switzerland) for
providing advice regarding the correlation of meteorological variables and LIP
events.

## References

  * Aksakal [2013] Aksakal, S.K., 2013.  Geometric Accuracy Investigations of SEVIRI High Resolution Visible (HRV) Level 1.5 Imagery.  Remote Sens. 5, 2475–2491.
  * Beyene and Jain [2018] Beyene, M.T., Jain, S., 2018\.  Freezing degree-day thresholds and Lake ice-out dates: Understanding the role of El Niño conditions.  Int. J. Climatol. 38, 4335–4344.
  * Breiman [2001] Breiman, L., 2001.  Random forests.  Mach. Learn. 45, 5–32.
  * Chen et al. [2018] Chen, L.C., Zhu, Y., Papandreou, G., Schroff, F., Adam, H., 2018.  Encoder-decoder with atrous separable convolution for semantic image segmentation, in: European Conference on Computer Vision,, Munich, Germany, 8–14 September.
  * Chen and Guestrin [2016] Chen, T., Guestrin, C., 2016\.  XGBoost: A Scalable Tree Boosting System, in: International Conference on Knowledge Discovery and Data Mining, San Francisco, USA, 13–17 August.
  * Cortes and Vapnik [1995] Cortes, C., Vapnik, V., 1995\.  Support-Vector Networks.  Mach. Learn. 20, 273–297.
  * Crétaux et al. [2020] Crétaux, J.F., Merchant, C., Duguay, C., Simis, S., Calmettes, B., Bergé-Nguyen, M., Wu, Y., Zhang, D., Carrea, L., Liu, X., Selmes, N., Warren, M., 2020\.  ESA Lakes Climate Change Initiative (Lakes_cci): Lake products, Version 1.0. Centre for Environmental Data Analysis.  Available online: http://dx.doi.org/10.5285/3c324bb4ee394d0d876fe2e1db217378 (accessed 16 February 2021).
  * Duguay et al. [2015] Duguay, C., Bernier, M., Gauthier, Y., Kouraev, A., 2015\.  Remote sensing of lake and river ice, in: Tedesco, M. (Ed.), Remote Sensing of the Cryosphere. Wiley-Blackwell, Oxford, UK, pp. 273–306.
  * Forster et al. [2020] Forster, P.M., Maycock, A.C., McKenna, C.M., Smith, C.J., 2020\.  Latest climate models confirm need for urgent mitigation.  Nat. Clim. Chang. 10, 7–10.
  * Gou et al. [2015] Gou, P., Ye, Q., Wei, Q., 2015\.  Lake ice change at the Nam Co Lake on the Tibetan Plateau during 2000-2013 and influencing factors.  Progress in Geography 34, 1241–1249.
  * Hall and Riggs [2016] Hall, D.K., Riggs, G.A., 2016\.  MODIS/Terra Snow Cover Daily L3 Global 500m SIN Grid, Version 6, Boulder, Colorado, USA, NASA National Snow and Ice Data Center Distributed Active Archive Center.  Available online: https://doi.org/10.5067/MODIS/MOD10A1.006 (accessed 16 February 2021).
  * Hendricks Franssen and Scherrer [2008] Hendricks Franssen, H.J., Scherrer, S.C., 2008\.  Freezing of lakes on the Swiss plateau in the period 1901–2006.  Int. J. Climatol. 28, 421–433.
  * Hoekstra et al. [2020] Hoekstra, M., Jiang, M., Clausi, D.A., Duguay, C., 2020\.  Lake Ice-Water Classification of RADARSAT-2 Images by Integrating IRGS Segmentation with Pixel-Based Random Forest Labeling.  Remote Sens. 12, 1425\.
  * Jégou et al. [2016] Jégou, S., Drozdzal, M., Vázquez, D., Romero, A., Bengio, Y., 2016.  The One Hundred Layers Tiramisu: Fully Convolutional DenseNets for Semantic Segmentation, in: International Conference on Computer Vision and Pattern Recognition Workshops, Las Vegas, USA, 26 June – 1 July.
  * Kropáček et al. [2013] Kropáček, J., Maussion, F., Chen, F., Hoerz, S., Hochschild, V., 2013.  Analysis of ice phenology of lakes on the tibetan plateau from modis data.  Cryosphere 7, 287–301.
  * Masson-Delmotte et al. [2018] Masson-Delmotte, V., Zhai, P., Pörtner, H.O., Roberts, D., Skea, J., Shukla, P.R., Pirani, A., Moufouma-Okia, W., Péan, C., Pidcock, R., Connors, S., Matthews, J.B.R., Chen, Y., Zhou, X., Gomis, M., I., Lonnoy, E., Maycock, T., Tignor, M., Waterfield, T., 2018.  IPCC: Summary for Policymakers, in: Global Warming of 1.5°C. An IPCC Special Report on the impacts of global warming of 1.5°C above pre-industrial levels and related global greenhouse gas emission pathways, in the context of strengthening the global response to the threat of climate change, sustainable development, and efforts to eradicate poverty, World Meteorological Organization, Geneva, Switzerland, 32 pp, Available online: https://www.ipcc.ch/sr15/ (accessed 16 February 2021).
  * Ogutu et al. [2011] Ogutu, J.O., Piepho, H., Schulz-Streeck, T., 2011.  A comparison of random forests, boosting and support vector machines for genomic selection.  BMC Proc. 5, S11.
  * Owen [2006] Owen, A.B., 2006.  A robust hybrid of lasso and ridge regression, Technical Report, Stanford University.  Available online: https://statweb.stanford.edu/~owen/reports/hhu.pdf (accessed 16 February 2021).
  * Pham et al. [2020] Pham, T., Yokoya, N., Xia, J., Ha, N., Le, N., Nguyen, T., Dao, T., Vu, T., Pham, T., Takeuchi, W., 2020\.  Comparison of Machine Learning Methods for Estimating Mangrove Above-Ground Biomass Using Multiple Source Remote Sensing Data in the Red River Delta Biosphere Reserve, Vietnam.  Remote Sens. 12, 1334\.
  * Prabha et al. [2020] Prabha, R., Tom, M., Rothermel, M., Baltsavias, E., Leal-Taixe, L., Schindler, K., 2020\.  Lake Ice Monitoring with Webcams and Crowd-Sourced Images.  ISPRS Ann. Photogramm. Remote Sens. Spatial Inf. Sci. V-2-2020, 549–556.
  * Qi et al. [2020] Qi, M., Liu, S., Yao, X., Xie, F., Gao, Y., 2020.  Monitoring the Ice Phenology of Qinghai Lake from 1980 to 2018 Using Multisource Remote Sensing Data and Google Earth Engine.  Remote Sens. 12, 2217\.
  * Qi et al. [2019] Qi, M., Yao, X., Li, X., Duan, H., Gao, Y., Liu, J., 2019.  Spatiotemporal characteristics of Qinghai Lake ice phenology between 2000 and 2016.  J. Geogr. Sci. 29, 115–130.
  * Rebetez [1996] Rebetez, M., 1996.  Seasonal relationship between temperature, precipitation and snow cover in a mountainous region.  Theor. Appl. Climatol. 54, 99–106.
  * Reed et al. [2009] Reed, B., Budde, M., Spencer, P., Miller, A.E., 2009\.  Integration of MODIS-derived metrics to assess interannual variability in snowpack, lake ice, and NDVI in southwest Alaska.  Remote Sens. Environ. 113, 1443–1452.
  * Riggs et al. [2017] Riggs, G.A., Hall, D.K., Román, M.O., 2017.  NASA S-NPP VIIRS Snow Products Collection 1 (C1) User Guide.  Available online: https://nsidc.org/sites/nsidc.org/files/technical-references/VIIRS-snow-products-user-guide-final.pdf (accessed 16 February 2021).
  * Roy et al. [2002] Roy, D., Borak, J., Devadiga, S., Wolfe, R., Zheng, M., Descloitres, J., 2002\.  The MODIS Land Product Quality Assessment Approach.  Remote Sens. Environ. 83, 62–76.
  * Slater et al. [2021] Slater, T., Lawrence, I., R., Otosaka, I., N., Shepherd, A., Gourmelen, N., Jakob, L., Tepes, P., Gilbert, L., Nienow, P., 2021.  Review article: Earth’s ice imbalance.  Cryosphere 15, 233–246.
  * Tom et al. [2020a] Tom, M., Aguilar, R., Imhof, P., Leinss, S., Baltsavias, E., Schindler, K., 2020a.  Lake Ice Detection from Sentinel-1 SAR with Deep Learning.  ISPRS Ann. Photogramm. Remote Sens. Spatial Inf. Sci. V-3-2020, 409–416.
  * Tom et al. [2018] Tom, M., Kälin, U., Sütterlin, M., Baltsavias, E., Schindler, K., 2018.  Lake ice detection in low-resolution optical satellite images.  ISPRS Ann. Photogramm. Remote Sens. Spatial Inf. Sci. IV-2, 279–286.
  * Tom et al. [2020b] Tom, M., Prabha, R., Wu, T., Baltsavias, E., Leal-Taixe, L., Schindler, K., 2020b.  Ice Monitoring in Swiss Lakes from Optical Satellites and Webcams using Machine Learning.  Remote Sens. 12, 3555\.
  * Tom et al. [2019] Tom, M., Suetterlin, M., Bouffard, D., Rothermel, M., Wunderle, S., Baltsavias, E., 2019\.  Integrated monitoring of ice in selected Swiss lakes, Final Project Report.  Available online: https://arxiv.org/abs/2008.00512 (accessed 16 February 2021).
  * Wainer [2016] Wainer, J., 2016.  Comparison of 14 different families of classification algorithms on 115 binary datasets.  arXiv pre-print, arXiv:1606.00930v1 .
  * Wang et al. [2018] Wang, W., Lee, X., Xiao, W., Liu, S., Schultz, N., Wang, Y., Zhang, M., Zhao, L., 2018.  Global lake evaporation accelerated by changes in surface energy allocation in a warmer climate.  Nature Geosci. 11, 410–414.
  * Wu et al. [2021] Wu, Y., Duguay, C.R., Xu, L., 2021.  Assessment of machine learning classifiers for global lake ice cover mapping from MODIS TOA reflectance data.  Remote Sens. Environ. 253, 112206\.
  * Xiao et al. [2018] Xiao, M., Rothermel, M., Tom, M., Galliani, S., Baltsavias, E., Schindler, K., 2018\.  Lake ice monitoring with webcams.  ISPRS Ann. Photogramm. Remote Sens. Spatial Inf. Sci. IV-2, 311–317.
  * Yao et al. [2016] Yao, X., Li, L., Zhao, J., Sun, M., Li, J., Gong, P., An, L., 2016.  Spatial-temporal variations of lake ice phenology in the Hoh Xil region from 2000 to 2011.  J. Geogr. Sci. 26, 70–82.

##  Appendix A Operational lake ice products

Details of operational lake ice products are shown in Table 7.

Product | Availability | Spatial | Temporal | Input
---|---|---|---|---
|  | res | res | sensor(s)
CCI Lake Ice Cover | from 2000 | 250m | daily | MODIS, VIIRS
|  |  |  | Sentinel-1,-3
MODIS Snow Product | from 2000 | 500m | daily | MODIS
VIIRS Snow Product | from 2012 | 375m | daily | VIIRS
Lake Ice Extent | from 2017 | 250m | daily | MODIS
Table 7: Comparison of the operational lake ice / snow products. Note that
lake ice extent is still a pre-operational product. Res indicates resolution.

##  Appendix B Characteristics of the target lakes

Details of the target lakes are shown in Table 8

| Sihl | Sils | Silvaplana | St. Moritz
---|---|---|---|---
Lat (∘N), Long (∘E) | 47.14, 8.78 | 46.42, 9.74 | 46.45. 9.79 | 46.49, 9.85
Altitude (m) | 889 | 1797 | 1791 | 1768
Depth (Max, Avg) | 23, 17 | 71, 35 | 77, 48 | 42, 26
Area (km${}^{\text{2}}$) | 11.3 | 4.1 | 2.7 | 0.78
Volume (Mm${}^{\text{3}}$) | 96 | 137 | 140 | 20
Meteo station | EIN | SIA | SIA | SAM
Lat (∘N), Long (∘E) | 47.13, 8.75 | 46.43, 9.77 | 46.43, 9.77 | 46.53, 9.88
Altitude (m) | 910 | 1804 | 1804 | 1708
Table 8: Details of the lakes. (primary source: Wikipedia). Maximum and
average depths are shown in m. Last three rows display information about the
nearest meteorological stations.

##  Appendix C MODIS and VIIRS bands

MODIS and VIIRS band spectrum are displayed in Fig. 14

Figure 14: Spectrum of the MODIS (M, blue) and VIIRS (V, red) bands used in
our analysis.

##  Appendix D Class imbalance in our dataset

Details on class imbalance in our dataset are shown in Fig. 15

Sihl (M) Sihl (V) Sils (M) Sils (V) Silv (M) Silv (V) Moritz (M)
$0$$5{,}000$$10{,}000$$15{,}000$$20{,}000$$4{,}137$$1{,}919$$2{,}345$$894$$1{,}736$$739$$157$$13{,}568$$4{,}598$$3{,}019$$1{,}051$$1{,}965$$765$$191$Number
of pixelsFrozenNon-Frozen

Sihl (M) Sihl (V) Sils (M) Sils (V) Silv (M) Silv (V) Moritz (M)
$0$$5{,}000$$10{,}000$$15{,}000$$20{,}000$$1{,}005$$198$$1{,}858$$722$$1{,}169$$591$$124$$11{,}804$$4{,}311$$2{,}435$$784$$1{,}574$$621$$140$Number
of pixels

Figure 15: Bar graphs showing the class distribution in our dataset from the
winters 2016–17 (top) and 2017–18 (bottom). The total number of clean, cloud-
free pixels from the non-transition dates that are at least 30% cloud-free are
shown. M and V denote MODIS and VIIRS respectively. Silv and Moritz represent
lakes Silvaplana and St. Moritz respectively.

##  Appendix E Correlation of weather data and LIP events

Correlation of LIP events with temperature variables are shown in. Fig. 16.
Similarly, correlations with sunshine and precipitation are displayed in Fig.
17.

Figure 16: Correlation of MWT (∘C) with ICD (days) and BUS are shown in rows 1
and 2 respectively, AFDD (∘C) with ICD, CFD (days) and FUE are displayed in
rows 3, 4 and 5 respectively. Results for lakes Sils and Silvaplana are
displayed in left and right columns respectively.

Figure 17: Correlation of total winter sunshine (hours) with ICD (days) is
shown in row 1, total sunshine from January to May (J2M) with BUS and BUE are
displayed in rows 2 and 3 respectively. Correlation of total precipitation
(mm) from January to May (J2M) and BUS is shown in last row. Results for lakes
Sils and Silvaplana are displayed in left and right columns respectively.