Area, volume and ELA changes of West Greenland local glaciers and ice caps from 1985–2020

Abstract

The behaviour of mountain glaciers on decadal time scales is a useful indicator for assessing climate change. Although less monitored and studied than the ice sheet, local glaciers and ice caps along the coast of Greenland are substantial contributors to meltwater runoff and sea level rise. This study analyses the cumulative area, ice mass and Equilibrium Line Altitude (ELA) change that occurred on 4100 glaciers and ice caps in West Greenland from 1985 to approximately 2020, using remotely sensed data and including glaciers smaller than 1 km² in the calculations. The glaciers involved in the study decreased in area by 1774 ± 229 km² which corresponds to almost −15%. Their surface elevation decreased on average by 20.6 ± 3.9 m, corresponding to a rate of −0.5 ± 0.1 m w.e. a⁻¹. The ELA shows a median regional rise of 150 m with marked local variability and higher median rise in the northern part of the study area. Strong regional gradients in ELA of individual glaciers are found, both towards the ice sheet and in areas where local orography affects precipitation. The observed high spatial variability of changes suggests that more monitoring on sub-regional level is needed.

1. Introduction

Mountain glaciers and their mass balance are useful indicators of regional changes in the climate system (Bojinski and others, 2014; IPCC, 2019). Glaciers have a shorter response time compared to ice sheets (Noël and others, 2017; Jakob and Gourmelen, 2023), making them an essential element in assessing the cryosphere's response to anthropogenic climate forcing (Haeberli and others, 2000; Rounce and others, 2023). In recent decades, there has been an increasing effort towards global estimation of glacier mass balance and their role within the water cycle. The measured global ice losses (Zemp and others, 2015) have highlighted the role of glaciers as one of the major contributors to global mean sea level rise (Zemp and others, 2019; Hugonnet and others, 2021). Various solutions for monitoring glacier mass changes from space have been put forward (Berthier and others, 2023), still, many remote areas lack high-resolution data over space and time.

The sensitivity of glaciers to climate conditions can be studied by analysing the variations in their Equilibrium Line Altitude (ELA) (Ohmura and others, 1992). ELA is defined as the altitude at which a glacier's net mass balance is zero (Cogley and others, 2011). Estimating the regional ELA is typically a demanding task, however, this information can be retrieved by using simple topographical parameters (Braithwaite and Raper, 2009) or by modelling it on a regional scale based on climate variables (Žebre and others, 2021; Del Gobbo and others, 2023).

Local Glaciers and Ice Caps outside the Greenland Ice Sheet (GrIS), at times also referred as Peripheral Glaciers (RGI region 05, Greenland Periphery), account for 13% of the global glacierized area and 10% of its total volume (excluding ice sheets, Farinotti and others, 2019). Rastner and others (2012) identified approximately 19 400 glaciers and ice caps in Greenland that account for 11 ± 2% (Khan and others, 2022) of current ice loss, despite making up only 4% of the Greenland ice covered area. Estimations of local glaciers and ice caps recent ice losses show substantial uncertainty but prove an acceleration of retreat rates trend, with values between 23.8 ± 2.0 Gt a⁻¹ (2010–2020, Jakob and Gourmelen, 2023) and 42.3 ± 6.2 Gt a⁻¹ (2018–2021, Khan and others, 2022). Widespread acceleration has also been observed in land terminating glaciers around Greenland, with a mean frontal change rate of −14.8 ± 5.4 m a⁻¹, roughly double that of the previous century (Larocca and others, 2023). Reconstructions of long-term average loss of glaciers and ice caps since the Little Ice Age (LIA) have confirmed that mass loss rates increased over the last two decades, despite the significant variability associated with local feedback (Carrivick and others, 2023). Local glaciers and ice caps are expected to lose between 2016 ± 129 and 3907 ± 108 Gt of ice by 2100 (Machguth and others, 2013). Remotely sensed variations in area or volume of glaciers and ice caps are available during different periods for South (Bjørk and others, 2012), Central East/West (Citterio and others, 2009; Jiskoot and others, 2012; Bjørk and others, 2018; Huber and others, 2020), North (Saito and others, 2016; Wang and others, 2021) and whole Greenland periphery (Bolch and others, 2013; Khan and others, 2022; Carrivick and others, 2023; Jakob and Gourmelen, 2023; Larocca and others, 2023).

Local glaciers play an important role for both society and local ecosystems. They provide water for reservoirs that generate hydropower as well as serving as recreational areas. However, their remote location often makes systematic monitoring challenging. To our best knowledge only six local glaciers and ice caps in the entire Greenland periphery are monitored annually or seasonally in the field. Of these, only Lyngmarksbræen in Qeqertarsuaq (Disko Island) and Qasinnguit, close to Nuuk, are located on the Central West coast. Machguth and others (2016) have compiled a comprehensive database with all available field measurements from past monitoring of local glaciers. The absence of adequate data coverage and field measurements is evident, which implies that numerous processes involving glaciers and ice caps in Greenland remain incompletely understood.

The aim of this study is to assess the recent changes and the current state of Central West Greenland glaciers. This is achieved by measuring variations of all local glaciers and ice caps from 1985 to approximately 2020 with a focus on area (i), surface elevation and mass balance (ii), and ELA (iii).

2. Methods

2.1. Study area and dataset

This study focuses on local glaciers and ice caps located between the settlements of Upernavik (72.79° N), and Qeqertarsuatsiaat (63.10° N), west of the Greenland Ice Sheet. This area is the most populated zone of the country, containing more than 65% of the total inhabitants. According to the Randolph Glacier Inventory v7.0 (2023), 4093 glaciers and ice caps larger than 0.01 km² lie within the study area. These were subdivided into six geographical zones from north to south (Fig. 1), named as following: (1) Sigguup Nunaa (Nunavik Peninsula, also known and reported as Svartenhuk Halvø), (2) Nuussuaq, (3) Qeqertarsuaq (Disko Island), (4) Sisimiut, (5) Maniitsoq (Sukkertoppen), and (6) Nuuk. Although Greenland local glaciers smaller than 1 km² are almost negligible in terms of overall area (2%) and volume (<5%) (Carrivick and others, 2023), we chose to include them in the analysis, as they provide additional insights at a sub-regional level. Ice caps in the study area are all included in the general calculations, except for Kangaamiut Sermiat and Tasersiap Sermia, the two largest Maniitsoq (5) ice caps, which are separated in some of the overall statistics because of their larger size.

Figure 1. The study area in West Greenland (a); Northern Central (b) and Southern Central (c) coast of West Greenland with local glaciers and ice caps extents and zones used in this study (1–6). The number of glaciers per area and their position on the map are taken from the Randolph Glacier Inventory v7.0 (RGI Consortium, 2023). The base layer and GrIS area are taken from Styrelsen for Dataforsyning og Infrastruktur (SDFI) vector maps.

The full dataset used in this study is shown in Table 1. Present area of glaciers and ice caps was analysed using orthorectified Sentinel-2 and SPOT 6–7 imagery available from the Danish SDFI (Styrelsen for Dataforsyning og Infrastruktur). The glacier extent from 1985 was retrieved from PROMICE ice margin vector product (Citterio and Ahlstrøm, 2013), slightly updated in zone (1) using Korsgaard and others (2016) aerial imagery. Glaciers and ice caps surface elevation change was measured using two datasets: (i) 1985 Digital Elevation Model (DEM) from Korsgaard and others (2016) and (ii) 259 ArcticDEM strips (Porter and others, 2022). ArcticDEM strips were manually selected (mostly from August 2020 and 2021, see Fig. S1 for reference) based on three criteria: (i) minimal or absent cloud coverage, (ii) as close to the end of the ablation season as possible (with 20 August serving as reference date), (iii) as recent as possible – beginning in 2021 and going backwards if no suitable strip was found. ELA and Accumulation Area Ratio (AAR) of Greenland glaciers were determined using Sentinel 2-L2A satellite imagery for 2018–22, accessed from Sentinel Hub (Copernicus, 2023).

Table 1. Data sources used in this study for Area (A), Volume (V) and Equilibrium Line Altitude (ELA) analysis with references, resolution, main method used and year of the data collection

All the datasets are publicly available. Acronyms used throughout the manuscript are available in Appendix A.

^a A few ArcticDEM Strips are from 2019 and previous years, back to 2014 (see Fig. S1).

^b Info about the exact dates used for every zone is shown in Table S2.

2.2. Area change

The 1985 local glaciers and ice caps polygons (Citterio and Ahlstrøm, 2013) necessitated the addition of a few local glaciers in some areas of zone (1), Sigguup Nunaa. This was done manually using the Korsgaard and others (2016) 1985 orthophoto as a reference (Fig. S2). Glaciers and ice caps snow and ice areas from August 2019–2020 were retrieved from SDFI raster files (Table 1) using an Object Based Image Analysis (OBIA) approach, available in the open source Orfeo toolbox (Grizonnet and others, 2017). The 4-bands (RGB + NIR) raster dataset was segmented using the Large-Scale Mean Shift algorithm (Michel and others, 2015) with a Radius of the spatial neighbourhood of 10, a range radius of 15 and a tile size of 500 pixels both for X and Y axis. Zonal statistics were used on outputs of segmentation, retrieving information of mean and std dev. values for each band of the original raster.

Zonal statistics output was used as input for a semi-automatic libSVM classifier (Chang and Lin, 2011) that was manually trained using approximately 200 sampled points representing different classes (i.e. ice, snow, land, water). The classified vector dataset was then reprojected to WGS-84 to match with the 1985 data and merged into a single shapefile, using the Randolph Glacier Inventory (RGI Consortium, 2023) polygons to filter out any external area from further analysis. The 1985–2020 area change was calculated using vector difference, linking every area loss with the former 1985 glacier ID to recover relative area loss. As relative area loss is calculated using 1985 polygons that have no drainage basins, the calculations report only one value per group of glaciers sharing edges.

2.3. Surface elevation change and mass balance

The DEMs used in this study to calculate the surface elevation change required pre-processing. 1985 DEMs were re-reprojected to match WGS-84 Polar Stereographic North of modern data and their reliability masks (Korsgaard and others, 2016) were used for filtering nonreliable pixels (i.e. values below 40 on a scale from 1 to 100). This generated the presence of significant voids. The manually selected ArcticDEM Strips were merged keeping the data closer to August 2020 when overlapping and then upscaled to match 1985 data resolution.

Surface elevation change is quantified by comparing DEMs using the Multi Scale Model to Model Cloud Comparison (M3C2) algorithm in CloudCompare (Lague and others, 2013), which also allowed the elaboration of the distance uncertainty used for the error calculation (see Table S1). Voids within the 1985 DEM, often coinciding with glaciers and ice caps accumulation areas, were filled using a three-step interpolation approach based on the xDEM package for Python (xDEM contributors, 2021). A linear interpolation was used to fill smaller voids and jagged void margin with a distance-weighted average of the surrounding pixels values. Subsequently, local (ii) and regional (iii) hypsometric interpolations were applied. These are based on the idea that there is a relationship between the elevation and the elevation change in the M3C2 raster, which correlates strongly for a specific glacier, and weakly on regional scales (McNabb and others, 2019). Therefore, the local hypsometric interpolation, tailored for every individual glacier, was applied to fill the largest voids. Subsequently, remaining voids on glaciers with a limited coverage, were filled using the regional hypsometric interpolation.

The M3C2 xDEM-interpolated raster (e.g. Fig. S3) is used to compute zonal statistics with two different vector-based glacier datasets: (i) one coming from the RGI (2023) with drainage basins for mass balance (common area) and (ii) the 1985 polygons for the calculation of total mass losses. Mass balance conversions from geodetic measurements were undertaken using a density conversion factor of 850 kg m⁻³ as suggested in the literature (Huss, 2013) and already done in similar studies (Huber and others, 2020; Carrivick and others, 2023) involving local glaciers and ice caps. Glacier thickness data available from Millan and others (2022) were used to estimate the relative volume loss per glacier, assuming that their data are representative of the present.

2.4. Equilibrium line altitude and accumulation area ratio

The ELA marks the elevation on the glacier's surface where mass balance equals zero. AAR is the ratio of the area of the accumulation zone to the area of the glacier, bounded by values 0 and 1 (Cogley and others, 2011). In this study we refer to the effective ELA (effELA; Žebre and others, 2021) indicating the altitude of the equilibrium line of a glacier; simply referred as ELA, if not differently specified.

Here, we have assessed the AAR from 2018 to 2022, based on snow cover classification at the end of the ablation season, using thresholding and re-classification of Sentinel-2-L2A Band 08 (Near infrared). The snow-covered regions were converted into vectors and then clipped using 2020 glaciers area, as a reference for AAR. The 5-year average AAR was calculated for every glacier and ice cap, here referred as AAR_2018–22.

Modern-day ELAs were determined using a modified version of the AAR method (Pellitero and others, 2015) that involved a specific AAR_2018–22 for each glacier, instead of a fixed value. To prevent potential underestimation of ELAs, we have limited our calculations to glaciers and ice caps with AAR > 0 for all 5 years of analysis and an area greater than 1 km². This approach avoids errors that might result from the inclusion of glaciers with ELA above their topography (i.e. AAR = 0) or from the erroneous classification of snow in smaller glacial bodies bounded by steep terrain.

Retrieving 1985 ELAs from imagery classification was not possible. The precision of snow classification is reduced when using single-band black and white imageries, such as the Korsgaard and others (2016) orthophotos used in this study. Furthermore, the available data was limited to a single year with above-average temperatures (Fig. S4). Lastly, the photos captured in 1985 were mostly taken during the second part of July, not including the whole melting season. Based on the available measurements (Machguth and others, 2016; Cappelen and others, 2021) and studies (e.g. Kelly and Lowell, 2009; Abermann and others, 2017), we assume that glaciers in West Greenland were close to the steady state during the 1980s (Braithwaite and Muller, 1980). Consequently, the approach by Pellitero and others (2015) was used to calculate ELA_AAR with AAR of 0.67 (Braithwaite and Muller, 1980; Larocca and Axford, 2022) and ELA_AABR with AABR (Area-Altitude Balance Ratio) of 2.24 ± 0.85 (Rea, 2009; Carrivick and others, 2023) for all 1985 glaciers and ice caps. For these calculations, RGI v7.0 areas (2023) are used as the 1985 glacierized area lacks drainage basins. The xDEM interpolated raster with surface elevation change (section 2.3) is subtracted from the modern Digital Elevation Model to recover the 1985 surface elevation without voids.

2.5. Uncertainties assessment

The accuracy of glacier area classification was evaluated by deriving a confusion matrix from the classifier results (see section 2.2), using the Orfeo toolbox (Grizonnet and others, 2017). Ground truth points were generated manually as no existing suitable dataset was available for validation. The metrics derived from the confusion matrix for ice, snow and other surfaces (i.e. land and water) are based on the following: True Positives (TP) – Negatives (TN) and False Positives (FP) – Negatives (FN) (Olofsson and others, 2014). These are used to calculate (i) Precision (P), that is the ratio between TP and the sum of TP and TN that measures the quality of positive predictions, and (ii) Recall (R), that is the ratio between TP and the sum of TP and FN that measures the classifier performance to find all positives (Marochov and others, 2021). The harmonic mean of P and R, F1 score, is typically used as one of the main performance metrics (Marochov and others, 2021). Among the outputs of the confusion matrix calculation there is also the overall accuracy index, which measures the quality of correct predictions including true positives and negatives (Olofsson and others, 2014). We used the overall accuracy value to calculate the absolute area uncertainties reported in results. Concerning the uncertainties of 1985 polygons, an error estimation is calculated using a 10 m buffer on 1985 polygons, as defined in Citterio and Ahlstrøm (2013) original work on these areas. The total error for area differences between the classifier and the 1985 data was calculated by combining the two errors in quadrature.

The accuracy of local glaciers and ice caps surface elevation change was estimated using the uncertainty analysis tools available in the xDEM package (xDEM contributors, 2021) and described in detail in Hugonnet and others (2022). We calculated the normalised median absolute deviation (NMAD) for each zone of the study area, using binning of DEM slopes. This approach is based on stable terrain data as glaciers are masked out and provides a robust measure of the error even in the presence of heteroscedasticity and outliers (Hugonnet and others, 2022). The weighted mean of the NMAD value for all slope bins in each zone (i.e. 1–6) was used as overall elevation error, in meters. Errors in the mass change estimations consider also an uncertainty range for the density conversion factor, corresponding to ± 60 kg m⁻³ (Huss, 2013). Uncertainties coming for the elevation change (V) and density (ρ) are combined for the total mass change (M) uncertainty, using the error propagation formula (1) (Taylor, 1997):

(1)$$\delta _{\Delta {\rm M}} = \vert {V\rho } \vert \sqrt {{\left({\displaystyle{{\delta V} \over V}} \right)}^2 + {\left({\displaystyle{{\delta \rho } \over \rho }} \right)}^2} $$

ELA accuracy estimation can be performed on the snow classification only, without being applied to the subsequent ELA calculation. This was done by using the same confusion matrix used for areas, despite this time evaluating the binary snow-no snow.

3. Results

3.1. Area and volume loss

From 1985 until the end of the 2020 ablation season, the 4094 glaciers and ice caps analysed in this study lost 1774 ± 229 km² of glacierized area, corresponding to 15 ± 2% (Fig. 2). Overall, all six sub-zones showed an area change of −15% or more, ranging between −16 ± 2% and −27 ± 2%. The sole exception is Maniitsoq (5), where the relative losses were lower, at 8 ± 2%. All largest area losses were observed in correspondence to the frontal areas, with Qeqertarsuaq (3), Nuussuaq (2) and Sigguup Nunaa having a few cases of net area gain related to surging glacier activity, observed on 12 glaciers in total (Fig. S5).

Figure 2. The 1985–2020 local glaciers and ice caps area loss (centroids) in Central (a) and Southern Central (b) West Greenland. Zones (1–6) are from north to south with pie charts' size proportional to the total glacierized area of each zone in 1985. Pie charts are divided into glacier area loss (red), bare ice area (darker light blue), and 2020 snow covered glacier area (light blue). Maniitsoq (5) major ice caps (i.e. Kangaamiut Sermiat and Tasersiap Sermia) are not included in the pie chart.

Relative area loss for different size classes of glaciers varies between an average value of 62% for small glaciers (<1 km²) and a much lower value of 13% for glaciers > 50 km² (Fig. S6). The variability of area loss across glaciers is not affected by latitudinal gradient or proximity to the GrIS, nor by glaciers median elevation (Fig. S7).

There are 279 very small glaciers (i.e. area < 0.5 km²) that disappeared between 1985 and 2020, with an average former median elevation of 864 m a.s.l. Most of them were small cirque glaciers. Surging glaciers that show a net area gain were excluded from total area calculation but are shown in Figure 2 (blue). Their contribution (+14.9 km²) is < 1% of the total area loss. Of these more than 10 km² are related to Kuannersuit and another surging glacier in Qeqertarsuaq (3).

The area classification accuracy, evaluated using a confusion matrix (section 2.5), showed a Precision value of 0.96 for both snow and ice and a Recall of 0.98 for ice and 0.97 for snow. The F1 score was 0.97 for ice and 0.96 for snow. We obtained an overall accuracy of 0.96, ±2% of the classified area (2020).

Surface elevation changes measured from 1985 to 2020 show a wide range of values for small local glaciers (<1 km²), and losses closer to the area-weighted average of 20.6 ± 3.9 m for the larger ones (Fig. 3). This average is computed over the entire glaciers' surface, resulting in a rate of −0.6 ± 0.1 m a⁻¹. Most of the surface elevation change comes from glacier fronts, where values up to −85 m are reached in all zones of the study area except for Sisimiut (4), where maximum losses are between 55 and 60 m. The average error of surface elevation change measurements is 3.9 m, with a std dev. of 1.38 m across different slope intervals (see Table S1).

Figure 3. The 1985–2020 glaciers and ice caps surface elevation changes in Central (a) and Southern Central (b) Western Greenland. Ice mass in Gt is calculated using 850 kg m-3 as density conversion factor (Huss, 2013). Total ice mass remaining and relative losses are estimated using Millan and others (2022) thickness data.

The mass loss in the entire study area is 194.5 ± 36.8 Gt or an average loss rate of 5.6 ± 1.1 Gt a⁻¹. This corresponds to a loss of 19 ± 3% of the estimated ice volume of 1985. The location with the largest relative mass loss is Sisimiut (4), where 39 ± 5% in volume loss was measured. On the contrary, the relative mass loss is smaller where the largest ice caps are, in Maniitsoq (5), with 12 ± 2%. Table 2 shows the mass balance rate for the six West Greenland's glaciers and ice caps zones, divided into different size classes. The cumulative specific mass balance between 1985 and 2020 is −17.5 ± 3.6 m w.e., with an area-weighted average rate of −0.5 ± 0.1 m w.e. a⁻¹.

Table 2. The 1985–2020 area-weighted average mass balance rate for different local glaciers and ice caps size classes

Data are divided by zones of the study area (1–6). The average rate is computed according to the year that happened most frequently among ArcticDEM Strips, which is 2020 (49% of the 259 Strips used). Of the remaining Strips, 33% are from 2021, 6% from 2019, 12% from 2014–18 (see Fig. S1 for months and years used).

3.2. AAR and ELA variability

Smaller local glaciers (<1 km²) have a wide range of AAR values. Larger glaciers generally show higher AARs, with glaciers and ice caps > 50 km² keeping a substantial snow cover at the end of each ablation season. For three years out of five AAR of larger glaciers was > 0.5. AARs show a statistically significant (p < 0.001 every year) positive linear correlation with median elevation. The coefficient of correlation r ranges from 0.28 (in 2018) to 0.59 (in 2020), indicating varying degrees of strength. Highest correlation values are from glaciers < 1 km² and lowest from glaciers and ice caps > 10 km². Weighted area averages of AAR show that the snow cover compared to the glacierized area is 47%, on average (from 2018 to22), considering overall glacierized area at the end of the ablation season (Fig. 4).

Figure 4. Local glaciers and ice caps AAR with median elevation, divided into five size classes in 2018–2022. Average AAR for all glaciers and ice caps is weighted considering their area. Exact dates of the imagery, taken at the end of the ablation season, can be found in Table S2.

The snow classification on which AAR calculations are based was evaluated using a confusion matrix, showing Precision, Recall and F1 score in the range between 0.98 and 1.00, depending on the year and the zone of the study area involved. The overall accuracy of snow classification was also between 0.98 and 1.00.

2018–22 ELA_AAR calculations were performed on 933 local glaciers and ice caps out of the total 4094 included in the study area. Approximately 2900 glaciers and ice caps were below the 1 km² area threshold, while another portion had at least one year during the 2018–2022 period where no AAR data was available, that is AAR = 0 or the presence of cloud cover invalidated the snow cover classification (see section 2.4). Comparison with 1980s ELA to estimate ELA rise was undertaken on the same set of glaciers and ice caps.

Modern-day ELA_AAR shows a wide range of values, spreading from below 500 m in Maniitsoq (5) to above 2000 m in Sigguup Nunaa (1). In terms of geographical distribution, the lowest ELAs are found in the coastal regions, with a clear coast to ice sheet increase evident in all regions, except Nuussuaq (2) (Fig. 5), where the highest ELAs are closer to the GrIS. No north-south gradient is discernible from single glaciers and ice caps ELAs values across different zones, despite the median ELA per area is decreasing further towards the south (Table 3 and Figure 5).

Fig. 5. Equilibrium Line Altitude (ELA) variation between 1980's and 2018–22: Kernel density estimation plot of ELA calculated from Accumulation Area Ratio (AAR) in 2018–22 (grey) and 1980s (dashed line) and their median value per area. Map of 2018–22 ELAAAR in Central (a) and South Central (b) West Greenland. Nuuk (6) has no data available for 2022. Qeqertarsuaq (3) has only partial coverage in the 2021 dataset. Glaciers with no residual snow (AAR = 0) or fully covered (AAR = 1) between 2018–22 are not included in ELA calculations.

Table 3. Median ELA_AAR calculated for 2018–22 and 1980s on the same set of glaciers and ice caps (* glaciers and ice caps > 1 km² with 5-years of AAR > 0 with available satellite data with low cloud cover) and median ELA_AAR and ELA_AABR (All glaciers and ice caps in the 1980s) for each zone

Differences between ELA_AAR and ELA_AABR calculated for all glaciers and ice caps for the 1980s are minimal. Results from the 0.67 AAR method always fall within the ELAs interval resulting from the 2.24 ± 0.85 AABR method. The 1980s ELA_AAR calculated for the same set of glaciers and ice caps used for modern-day ELA_AAR show a median rise of 154 m, ranging between + 98 m for Maniitsoq (5) and + 256 m for Sigguup Nunaa (1). Sisimiut (4) and Nuuk (6) exhibit a shift in distribution while maintaining the same pattern, whereas the other areas exhibit both a shift and a change in distribution shape. ELAs > 1500 m increased from 8% in 1985 to 21% in 2018–22.

4. Discussion

4.1. Uncertainties and potential errors

As reported in the PROMICE map of Greenland peripheral ice masses (Citterio and Ahlstrøm, 2013), the 1985 polygons used in this work may be less precise when referring to Sigguup Nunaa (1), mapped by the KMS (National Survey and Cadastre, Denmark). Missing glaciers and ice caps areas were added manually (see Fig. S2) using Korsgaard and others (2016) orthophotos for mapping and RGI as reference. Misinterpretation of seasonal or occasional snow may represent a limiting factor in the 1985 data, while in the SDFI imagery snow cover is not relevant for glacier margin. Other limitations on area evaluation are related to the classifier performance, and the limits of this type of analysis when dealing with debris-covered surfaces (Pfeffer and others, 2014). Nevertheless, this problem occurs in only a minor part of the dataset. If the 10 m buffer approach used for 1985 areas (Citterio and Ahlstrøm, 2013) is extended also for the 2020 data, the error estimation gives an overall uncertainty of 1.7% on the total area loss, against the current 1.9%. The total number of extinct glaciers (279) should be regarded as an estimate as it is subject to inaccuracy due to the small size of these glacial bodies, which implies the possibility of misclassification in the 1985 dataset.

Sources of uncertainty for the volume change arise from DEM dates, their alignment with the 1985 data and the unavailability of suitable density measurements for mass conversion. The timestamp of the ArcticDEM Strips used is variable. However, the percentage of data from 2020 ± 1 a is high (88%) and most of the Strips used are from August, a period when snow cover is generally at its minimum (see Fig. S1 for reference). When dealing with long-term geodetic surface elevation change to mass balance conversions, a fixed ice density conversion factor of 850 kg m⁻³ was shown to be a good average value (Huss, 2013). This was used in similar works in West Greenland (Huber and others, 2020; Carrivick and others, 2023). While the presence of firn is less relevant for smaller local glaciers and ice caps (Gillespie and others, 2023), this could be more important when dealing with the larger ice caps of Maniitsoq (5) or with the 1985 dataset.

ELA values reported in this study are representative of the effective ELA only, and the dates of the imagery used are available from Table S2. The effective ELA of a singular glacier has less climatic significance than the environmental ELA of an area (Ohmura and others, 1992; Žebre and others, 2021), although the two can be much closer in value when the median effective ELA of a region is considered. The presence of superimposed ice was not considered in this study and remains an unsolved problem. The 1980s ELA was considered as steady state ELA as shown from the mass balance measurements available for three local glaciers (Glacier 33, Tasersiaq and Qaparfiup sermia) within the study areas in the 1980s (Machguth and others, 2016). This monitoring was undertaken as part of regional assessment of hydroelectric potential of West Greenland (Braithwaite and Olesen, 1982; Braithwaite, 1986; Olesen, 1986) and confirms the choice of excluding the 1985 imagery as representative of the 1980s. The available reports mention that 1985 winter was warmer than previous years in the 1980s, but more importantly summer melting rates higher and mass balance measurements lowest in the set of observations A warming period started the 80's in West Greenland as reported by Box (2002) and Abermann and others (2017).

4.2. Changing glaciers

Relative area loss of West Greenland local glaciers and ice caps since 1985 decreases for increasing glacier size classes, as previously shown in Citterio and others (2009). The range of relative area losses is also wider in values for smaller local glaciers and narrows significantly for larger ones (>50 km²). Out of the group of analysed glaciers and ice caps, area increase was observed in only 12 surging glaciers, 5 in Qeqertarsuaq (3), 5 in Nuussuaq (2), and 2 in Sigguup Nunaa (1). Central West Greenland is already known as a cluster for surge-type glaciers (Yde and Knudsen, 2007; Citterio and others, 2009). All the surging glaciers we have detected were already mentioned in the comprehensive overview by Lovell and others (2023), with the most studied and largest surge episode being in Kuannersuit glacier between 1995–1998 (Yde and others, 2019). Regarding surface elevation change, all surge-type glaciers in the study area show a frontal area with evident mass gain, although two of the sites lack enough data to reconstruct the variations, caused by the 1985 DEMs voids. None of the surge-type glaciers exhibits a positive surface elevation change from 1985 if the entire glacier is considered.

We found a mass loss rate of 0.5 m w.e. a⁻¹ in Central West Greenland from 1985 to 2020, which can be compared with the reference glaciers (RG) of the World Glacier Monitoring Service (WGMS, 2021) in the Arctic region. According to 1985–2020 averages, our study area shows rates of mass loss similar to Svalbard and Jan Mayen (0.52 m w.e. a⁻¹, 2 RG), lower than Alaska (0.79 m w.e. a⁻¹, 3 RG) and Western North America (0.84 m w.e. a⁻¹, 7 RG) and higher than North Arctic Canada (0.31 m w.e. a⁻¹,4 RG) and Norway (0.37 m w.e. a⁻¹, 8 RG). Our observed mass loss rate is also higher than the 2010–2020 measures for the entire Greenland periphery (i.e. 0.27 ± 0.02 m w.e. a⁻¹; Jakob and Gourmelen, 2023).

Longer term mass loss rates were reported by Carrivick and others (2023) for Central West (0.37 m w.e. a⁻¹) and Southwest (0.30 m w.e. a⁻¹) Greenland glaciers and ice caps (ablation areas only), from LIA up until 2015. Their analysis showed much higher rates of loss for 2000–2019 period, with 0.78 m w.e. a⁻¹ for Central West and 1.06 m w.e. a⁻¹ for South West. Our analysis from 1985 to 2020 in Sigguup Nunaa (1), Nuussuaq (2) and Qeqertarsuaq (3) ablation areas only (i.e. Central West), show a mass loss rate of 0.59 m w.e. a⁻¹.

Despite being significantly lower in absolute terms compared to GrIS mass losses (Mankoff and others, 2021), the overall contribution brought by glaciers and ice caps of the study area is significant, with −5.56 ± 0.66 Gt a⁻¹ on average since 1985. The effect of these losses is particularly evident on small local glaciers (<1 km²), which have experienced relative thickness losses of up to 62%, according to Millan and others (2022) dataset. Moreover, several glaciers experienced frontal thickness losses exceeding 80 m and reaching 100 m in a few sites. Surface elevation change of West Greenland glaciers were retrieved from 1985 to 2012 in Huber and others (2020), showing a rate of −0.5 ± 0.2 m a⁻¹ on average, matching our results. Our mean surface elevation change rate of −0.6 ± 0.1 m a⁻¹ is close to that found by Hugonnet and others (2021) for the whole of Greenland periphery between 2000 and 2020 (i.e. −0.5 ± 0.04 m a⁻¹).

The relative volume loss data presented in this study is based on thickness calculations by Millan and others (2022), performed at global scale. Therefore, these values cannot be used to assess individual local glacier losses, especially if small (Fig. S8). For instance, Marcer and others (2017) measured a thinning of 0.60 ± 0.11 m a⁻¹ from 1985 to 2014 in the Aqqutikitsoq glacier (Sisimiut). In this work we observed a thinning rate of 0.57 m a⁻¹, which matches the previous findings. The same work included geophysical surveys that showed a total loss of a quarter of the glacier's 1985 volume in 2014. This was different from our estimation of a loss of 38% from 1985 to 2020.

In West Greenland, inland regions generally exhibit the highest ELAs compared to the coast. The geographical distribution of the ELA_AAR 2018–22 is relatable to the spatial distribution of precipitation observed in West Greenland from meteorological records (Weidick and others, 1992) and climatological reanalysis data (Box and others, 2023). Smaller amount of precipitation and higher ELAs proceeding inland primarily depend to orographic factors as well as the proximity to the GrIS, resulting in a dry climate with continental traits (Weidick and others, 1992). The meteorological and climatological causes that generate this situation are well described already in Ohmura and Reeh (1991), and are related to the combined role of topography and the southwesterly flows that direct most of the storms towards western Greenland. The moist flows directed towards western Greenland release most of their water vapour content due to the forced ascent of the flow induced by the presence of the first reliefs, leading to condensation and precipitation. Moving eastward, the air mass becomes progressively less humid. An additional factor in the drying of the air mass is related to the katabatic winds from the ice sheet, which make the climate near it more continental.

Glaciation level (i.e. critical altitude at which any mountain exceeding this limit and with space for glaciers will be glacierized; Humlum, 1985) trends in West Greenland are shown in earlier works (Humlum, 1985; Weidick and others, 1992), reporting year-to-year variability that can reach 100 m and an altitude increase towards inland areas following local climatology. In-depth exploration of deviations from the general trend owing to topo-climatic effects on precipitation is lacking, despite findings from Humlum (1985) research on Glaciation level in West Greenland demonstrate similar patterns to our own modern-day ELA observations, albeit our results being of higher elevation. The coast-to-interior temperature gradient in West Greenland is noticeable and its effect on glaciers and ice caps ELA is pronounced (Fig. S9). The topography shadow effect of near-coastal mountains and the increasing distance from the oceanic moisture source leads to lower precipitation areas towards the inland. This pattern is reflected by the higher ELAs observed further inland as shown by the 2018–22 ELA distribution map (Fig. 5) and was already evidenced in previous studies (Abermann and others, 2019).

According to Carrivick and others (2023) glaciers and ice caps in Central-West and South-West regions of Greenland have a median ELA of 1075 m and 896 m a.s.l., respectively, with a net rise of 99 and 73 m from the LIA until the present. In our study, a smaller sample of GICs > 1 km² yielded slightly higher values for the present (2018–22) compared to Carrivick and others (2023). This may be related to the different method used, as Carrivick refers to current glaciers median elevation. LIA ELAs from Carrivick and others (2023) are comparable on a regional level with the ELA from the 1980s in our findings, supporting the hypothesis of glaciers at steady-state in 1985. LIA-present ELA difference was also estimated using snowlines and geomorphological evidence in southern Greenland (Brooks and others, 2022), finding a rise of 122 ± 64 m.

We found a median ELA rise for local glaciers and ice caps in West Greenland of 154 m. This can be compared to a few Arctic reference glaciers with long term ELA data available (WGMS, 2021), from which 1980–85 to 2018–22 ELA rise can be calculated. Rise values comparable to our findings can be observed in Svalbard (+165 m, Austre Brøggerbreen; + 136 m, Midtre Lovénbreen), North Arctic Canada (+142 m, Devon Ice Cap NW) and Alaska (+112 m, Gulkana; + 174 m, Wolverine; + 353 m, Lemon Creek). As mentioned above, effective ELA measured on a small number of glaciers can only serve as an indicator for the range of values expected, but it's not as robust as an environmental ELA assessment.

The average ELA rise from this study ranges between 2.4 m a⁻¹ (Sisimiut, Maniitsoq) and 6.4 m a⁻¹ (Sigguup Nunaa) and can be compared to the common rates of ELA shift observed around the world (i.e. 2–5 m a⁻¹, Ohmura and Boettcher, 2022), except for the values of Sigguup Nunaa, which are slightly higher. ELA time series from singular arctic glaciers presented by Ohmura and Boettcher (2022) showed ELA rise rates of 7.6 m a⁻¹ (Mittivakkat, 1996–2018) in southeastern Greenland; 6.7 and 8.2 m a⁻¹ (White Glacier and Devon Ice Cap, 1979–2018) in the Canadian arctic; 4.2, 3.5, 7.6, 9.8 m a⁻¹ (Gulkana Glacier, Wolverine Glacier, Taku, Lemon Creek Glacier, 1979–2018) in Alaska. Data from Svalbard and Iceland showed lower ELA rising rates (between 1.1 and 2.2 m a⁻¹) and data from Scandinavia a wider range, from 0.4 to 8.4 m a⁻¹. ELA rise was measured also on larger numbers of glaciers with remote sensing ESS approach in the West European Alps by Rabatel and others (2013), with an average rate of + 6.5 m a⁻¹ (1984–2010).

The reduction in glacial coverage and the ELA rise observed in West Greenland provides clear evidence of unfavourable climatic conditions for glaciers in recent decades. At present, they are out of balance with the climate. Various studies in recent years have focused on the analysis of the evolution of certain meteorological-climatic parameters. In particular, Box and others (2023), using the Copernicus Climate Change Service Arctic Regional Reanalysis (CARRA), highlight that on peripheral glacial bodies, as well as on the GrIS, snowfall does not indicate any trend of increase or decrease in the 1991–2021 period. In the absence of further changes to other climatic parameters, this factor alone should in theory result in glaciers tending towards a stable state. However, the same authors point out that rainfall increased significantly (1 − p = 0.90) by 37% in the same period. It is demonstrated that rain directly enhances melting processes by up to 15–20% (Box and others, 2022), while indirect effects include snowpack heating through percolation and refreezing, and subsequent melt-albedo feedback.

Besides, Hanna and others (2021) found that Greenland coastal and inland weather station records have shown strong and significant warming from 1991 to 2019 in summer (~1.7°C), spring (~2.7°C), and winter (~4.4°C). In autumn, a significant warming trend is observed only in the eastern and northeastern areas. The most pronounced warming is observed in winter on the western coast and in northwest Greenland, up to ~6.5°C. Warming during the summer seems to be linked to variations in the Greenland Blocking Index (GBI). GBI is a measure related to atmospheric conditions over Greenland and represents the mean 500 hPa geopotential height over the region bounded by coordinates 60–80°N, 20–80°W. It quantifies the occurrence and strength of atmospheric high-pressure systems that tend to remain stationary when they occur. GBI increased in strength in spring and summer significantly (r ~ 0.8) in recent decades (Hanna and others, 2021) enhancing warm-air advections over the GrIS, and decreasing cloud-cover (Hofer and others, 2017). An increased prevalence of clear skies in summer enhanced solar radiation boosting daytime temperature maxima. On the contrary, winter temperatures, both since 1981 and 1991, exhibit a tendency to be higher during the night. This implies an increased influence of cloud-cover in hampering nocturnal cooling along coastal areas and might explain the higher observed warming trend (Hofer and others, 2017).

5. Conclusions

Our study on West Greenland Local Glaciers and Ice Caps underscores their diverse nature in terms of glacier type, size, aspect, and hypsometry. Measurements show high variability in the net changes from 1985, including area and volume shrinkage, alongside an elevation rise of ELAs. We observe a nearly 15% reduction in overall glacierized area, accompanied by an average mass loss rate of −0.5 ± 0.1 m w.e. a⁻¹ (5.6 ± 1.1 Gt a⁻¹), aligning with observed rates in the Arctic. Small local glaciers are the most affected by ice loss. The median ELA rise of 154 m corroborates existing understanding of regional climatic shifts. The distribution of ELAs currently displays a clear increase towards the ice sheet, with local variability being influenced by orographic effects on precipitation. These findings highlight the significance of ongoing monitoring and sub-regional analysis to reduce uncertainties regarding the future of these important components of Greenland's cryosphere. Further developments in this field will undoubtedly benefit from the increasing accessibility of these regions, advancements in remotely sensed imagery offering higher spatial and temporal resolution, growing scientific interest, and international collaboration.