Methodology:

Glacier Boundaries

Landsat-2 scenes were used to manually delineate glacier outlines for the 1970s. All other glacier outlines were produced from full-archive Landsat-5 (1984 - 2012) and Landsat-8 (2013 - 2020) imagery using a Normalised Difference Snow / Water Index (NDSI, NDWI) in Google Earth Engine. A cloud-free composite was produced for each year by extracting the median pixel value from the driest months (August to December) to minimise snow and cloud cover. Due to heavy cloud in these regions, there are some years within the time series with no usable imagery.

Elevation Change

All L1A ASTER tiles collected over these regions in the 20-year timespan from 1st January 2000 to 31st December 2020 were used to derive DEMs and co-registered to off-ice areas using the method of Nuth and Kääb (2011). DEMs were clipped to the range of elevation of the cordilleras, ± 100m of the median to remove erroneous values. A single value per pixel per year was extracted from all DEM tiles by collecting the median value. Using these single values per year, a linear fit through the surface elevation was then calculated over three time periods: 2000 - 2010; 2010 - 2020; and 2000 - 2020 to represent elevation change per year during these intervals.

Data collection:

ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer) L1A (AST_L1A)

Landsat-2 MSS L1TP

Landsat-5 TM Collection 1 Tier 1 TOA Reflectance (Google Earth Engine)

Landsat-8 Level 2 Collection 2 Tier 1 (Google Earth Engine)

Data quality:

All 1984 - 2020 glacier boundaries have an accuracy of 5%, with the 1975 Landsat-2 scenes at 11% owing to larger pixel sizes in the original dataset Volume change error is calculated as a median of off-ice areas within a 2 km buffer of glaciers and this error carried forward to calculate mass (following Brun et al. 2017) and mass balance.