Methodology:

We classified Landsat-8 OLI images acquired between 2016 and 2020 in Google Earth Engine (GEE) and ESRI ArcGIS Pro 2.6.0, primarily using K-means clustering.

Suitable images had low cloud cover and limited snow cover.

We selected six bands representing the visual and infrared wavelengths from the images for classification.

We added three further bands to the image in the form of normalised difference snow index (NDSI), normalised difference vegetation index (NDVI), and normalised difference water index (NDWI).

We topographically corrected the images.

We conducted a principal component analysis of the images in GEE and the first three components.

We used a hierarchical approach to image classification. A first order land classification of "land", "ice", and "water" informed the subdivision of each of these classes in a second, more detailed, analysis of the dominant land cover types. To produce these broad land classes, we used a K-means clustering algorithm in GEE to split each image into 75 (K value = 75) discrete clusters.

The clusters are determined using the spectral information of each image, based on 500,000 randomly selected sampling points.

We assigned each of these sections a first order class in ArcGIS by visually inspecting the image they were derived from.

In some cases, we could not easily assign a cluster a first order class. This was usually because a cluster had conflated shadow with dark seawater. To address this, we split these clusters using a slope threshold of 3o, with pixels <3o being assigned as water.

Where this process resulted in obvious misclassification we used a random forest classifier to differentiate between water, land and ice. Some pixels were covered entirely by very dark shadows or clouds and, therefore, we could not classify them; these were assigned "No data".

We used this first order land classification to subset each image in GEE accordingly, and then to cluster these resulting images into 40 discrete groups (K = 40).

We interpreted these clusters to manually assign each of them a final land classification. Our first-order land class was subset into five classes "Bedrock", "Coarse/wet sediment", "Fine & dry sediment", "Vegetation". The water class subset into. "Water" and "Turbid water", while the ice class subset into "Ice" and "Wet ice".

In cases where clouds partially obscured land, we assigned pixels to the more general class of "Land (non-differentiated)".

We produced ten land classes that describe eight distinct surface types, plus partially obscured land (Land (non-differentiated)), and surfaces totally obscured by clouds or shadows (No data).

Data collection:

Data was produced and analysed in Google Earth Engine and ESRI ArcGIS Pro 2.6.0

Data quality:

Accuracy assessed using the methods of:

Olofsson, P., Foody, G. M., Stehman, S. V, and Woodcock, C. E.: Making better use of accuracy data in land change studies: Estimating accuracy and area and quantifying uncertainty using stratified estimation, Remote Sens. Environ., 129, 122-131, https://doi.org/10.1016/j.rse.2012.10.031, 2013.

Olofsson, P., Foody, G. M., Herold, M., Stehman, S. V, Woodcock, C. E., and Wulder, M. A.: Good practices for estimating area and assessing accuracy of land change, Remote Sens. Environ., 148, 42-57, https://doi.org/10.1016/j.rse.2014.02.015, 2014.

We found our classification to have an overall accuracy of 95.9 % overall. When just the proglacial land classes are considered, we have an accuracy of 77.0 %.