Mass Loss of Glaciers and Ice Caps Across Greenland Since the Little Ice Age

Glaciers and ice caps (GICs) are important contributors of meltwater runoff and to global sea level rise. However, knowledge of GIC mass changes is largely restricted to the last few decades. Here we show the extent of 5327 Greenland GICs during Little Ice Age (LIA) termination (1900) and reveal that they have fragmented into 5467 glaciers in 2001, losing at least 587 km3 from their ablation areas, equating to 499 Gt at a rate of 4.34 Gt yr−1. We estimate that the long‐term mean mass balance in glacier ablation areas has been at least −0.18 to −0.22 m w.e. yr−1 and note the rate between 2000 and 2019 has been three times that. Glaciers with ice‐marginal lakes formed since the LIA termination have had the fastest changing mass balance. Considerable spatial variability in glacier changes suggest compounding regional and local factors present challenges for understanding glacier evolution.

• Total volume loss of at least 587 km 3 since the Little Ice Age (LIA) termination, equating to 499 Gt and to 1.38 mm sea level equivalent • Glacier mass balance from 2000 to 2019 is three times more negative than since the LIA but five times more negative in the North region • Lake-terminating glaciers have experienced the greatest change in rate of mass loss

Supporting Information:
Supporting Information may be found in the online version of this article.
10.1029/2023GL103950 2 of 9 precipitation (Machguth et al., 2013). However, satellite-derived data sets have revealed that there has been considerable variability in the absolute amount and rate of change of mass loss between Greenland regions (Hugonnet et al., 2021;Khan et al., 2022).
These recent decadal-scale changes to Greenland GICs lack a longer-term (pre-satellite era, centennial) context. A longer-term context is important because projections of Greenland GICs into the future require base line data sets, such as past glacier extents and past glacier ice surfaces, to hindcast and to spin-up model simulations from. Such calibration can increase confidence in glacier evolution models over a timescale that is representative to that of their projections.
The aims of this study are therefore to (a) quantify the past glacier extent and past ice surface of Greenland GICs during the Little Ice Age (LIA; 1900: Kjaer et al., 2022) termination, and (b) make a critical analysis of the variability in mass loss rates by region and by glacier type. We achieve these aims by employing a novel large-scale methodology to calculate glacier mass changes over a longer-term (centennial) timescale since the LIA.

Little Ice Age Moraine Mapping
We mapped LIA extents for 5467 GICs larger than 1 km 2 in the RGI v6.0 inventory (mean year 2001; GLIMS, 2021) primarily using a hillshaded image of the ArcticDEM (Porter et al., 2018) with 2 m horizontal resolution. Reference was also made to color optical-wavelength satellite images from PlanetScope (Planet, 2021) and, where that was not available, from 15 m resolution ESRI ArcGIS World Imagery. Moreover, we aided our mapping using visual interpretations of landforms in the 2 m AeroDEM orthomosaic (Korsgaard et al., 2016), and using landcover classifications of recently disturbed sediment made from 10 m resolution Sentinel 2 images from July and August 2021. Whilst the number of GICs larger than 1 km 2 is roughly one third of the ∼17,000 GICs within the RGI v6.0 inventory (GLIMS, 2021) with no connection to the GrIS (Rastner et al., 2012), they account for 99% of total GIC area.
We mapped LIA glacier extents by extending the RGI v6.0 outlines (GLIMS, 2021) down-valley to the crests of prominent moraine ridges interpreted to represent the LIA termination. That interpretation was based on position, typically several kilometers down-valley from contemporary GICs, and morphology and surface character. Moraine ridges were often sharp-crested, steep sided, and typically devoid of vegetation ( Figure 1a, Figures S2 and S3 in Supporting Information S1). If the LIA moraine comprised a set of multiple ridges, we deliberately mapped the innermost crest, which typically appeared more subdued in relief, more vegetated, and likely formed during older late-Holocene advances . Outlines from LIA lateral moraine crests were extended along trimlines using the aforementioned data sets. Our choice of the innermost LIA moraine ridges maintains consistency, permitting replicability, and provides a conservative estimate of glacier extent during the LIA termination. We were unable to map LIA extents where geomorphological evidence was absent such as for some water-terminating GICs and for some outlet glaciers in the extreme north. Overall, this mapping protocol together with our observations of some empty cirques above the regional LIA ELA (cf. Carrivick et al., 2019) means that our analysis provides a minimum estimate of LIA area, and of volume loss and mass loss since the LIA. Our calculations of glacier areal extent changes are to 2001 as that is the mean date of images used to construct the RGI v6.0 glacier outlines.

Ice Surface Reconstruction
To mitigate artifacts and high uncertainty in parts of the DEMs, due to snow cover and cloud cover in the images used to make them, we limited our analysis of surface elevation changes to glacier ablation areas only. Glacier ablation areas (for LIA and RGI v6.0 glacier outlines) were delineated using the Area-Altitude Balance Ratio (AABR) method, which was automated using code developed by Pellitero et al. (2015) and as adapted and used across several other world regions (Carrivick et al., 2019;Carrivick, Andreassen, et al., 2022;Carrivick, Tweed, et al., 2020;Lee et al., 2021). We assumed an AABR of 2.24 ± 0.85 as suggested by Rea (2009) for arctic glaciers but the suitability of this value varies on a glacier-by-glacier basis. We therefore aggregate our ELAs spatially, and we do not infer climate change from our ELAs. A LIA glacier surface was interpolated between point elevations extracted from the ablation area outline ( Figure S2 in Supporting Information S1).
Calculating the difference between our LIA surface and the ArcticDEM indicated the surface lowering that has occurred ( Figure 1b). Some positive elevation changes occurred, either due to glacier surges, convex glacier surfaces, or artifacts in the DEM. These causes could not be objectively discriminated and so we replaced all these positive values using a void filling spatial interpolation technique to negate them. Surface lowering maps were converted to a volume change by summing the elevation changes per area of interest (glacier, region) and multiplying by grid cell size. The difference in volume change computed with/without positives amounted to <−10% by region (Table S1 in Supporting Information S1). Importantly, our glacier mass balance estimates derive from these calculations of volume change across ablation areas, not from the ELA estimates. The parts of the ArcticDEM mosaic covering Greenland GICs were constructed from images dating between 2007 and 2017 and with a mean date of 2015, which we use for our volume change rate calculations. We take the LIA termination as year 1900 following Kjeldsen et al. (2015), and whilst acknowledging regional differences in the timing of the LIA termination (1850-1965: see Section 5.3 in Kjaer et al., 2022). Our rates of mass loss are most sensitive to the year chosen for the LIA termination.
More details of our data sets and mapping are given in Supporting Information S1 where we also show: (a) that the effect on our calculations of excluding RGI v6 glaciers <1 km 2 is near-negligible, typically 2% for areal extent and <5% for volume; (b) our LIA surface interpolation routine, (c) the sensitivity of our results to the choice of an AABR value for ablation area estimation, and to a date for the LIA termination, (d) our conversions of volume to mass and mass balance, and (e) our categorization of glacier types.

Results
There has been considerable fragmentation of GICs since the LIA termination in the form of tributaries separating from trunk glaciers. The 5327 GICs that we mapped LIA outlines for now number 5467; a 2.6% increase in number. Regionally, the total glacier extent during the LIA termination was between 76,342 and 84,378 km 2 , and so there has been at least 4,120 km 2 (5.1%) to 4,554 km 2 (5.7%) areal extent loss between 1900 and 2001. Areal extent losses between the LIA termination and 2001 were at least 17%-19% in the central west (CW) region, and at least 11%, 10%, and 9% in the south east (SE), south west (SW) and north west (NW) regions, respectively, but only ∼1% in the north (NO) region (Table S1 in Supporting Information S1).
GICs across Greenland have lost at least 528-646 km 3 volume since the LIA termination to year 2015. The greatest proportion of volume loss (25%-32%) has come from the north east (NE) region, but that region has 36%-40% of the total LIA glacier area. In contrast the SE, SW, and NW regions have undergone a disproportionate volume loss since the LIA with 8%, 11%, and 11%, respectively, comprising just 5%, 7%, and 8% of the total glacier area, respectively.
The median elevation rise of glacier termini per region has decreased with higher latitude on the west coast from 228 m a.s.l. (±12.5 m a.s.l.) in the SW to 169 m a.s.l. (±19.4 m a.s.l.) in the NW (minimum elevation; Figure 2). On the east coast the SE region is an anomaly to this latitudinal pattern, having the lowest terminus elevation  Figure 2). These regional medians reveal a higher rise in ELAs on the west of Greenland compared with those on the east coast and a N-S trend of decreasing ELA change with increasing latitude. However, there is considerable intra-region variability in ELA changes, particularly on the east coast ( Figure 2).
Using an ice density of 850 kg m −3 (Huss, 2013), the mass balance of all Greenland GICs in ablation areas since the LIA termination to 2015 has been at least −0.18 to −0.22 m w.e. yr −1 . That rate is at least −4.0 m w.e. yr −1 less than the −0.60 m w.e. yr −1 for 2000 to 2019 that we derived for Greenland GIC ablation areas only from the data set of Hugonnet et al. (2021). Thus, we suggest a conservative estimate of a 3 times faster rate of mass loss between 2000 and 2019 than the mean rate since the LIA termination. If we assume no mass loss from LIA glacier accumulation areas and apply our mass changes glacier-wide, then we estimate that the mass balance rate since the LIA for all GICs was −0.054 to 0.066 m w.e. yr −1 , which shows that the rate of −0.37 to 0.45 m w.e. yr −1 for 2000 to 2019 is an order of magnitude faster than the mean since LIA (Table S1 in Supporting Information S1).
Similarly, by examining regional glacier changes, we find that the Greenland-wide mean GIC mass balance rate over the past 20 years is at least twice as more negative as the LIA to 2015 mean. In the SW and NW regions the rate has become approximately three times more negative, from −0.3 to −0.76 and from −0.26 to −0.79 m w.e. yr −1 , respectively, and in the NO region the rate has become five times more negative from −0.11 to −0.57 m w.e. yr −1 (Figure 3). We note that the NO region was highlighted by Khan et al. (2022) for its rapid acceleration of glacier mass loss from 2003 to 2009. These inter-regional patterns are the same for absolute volume ( Figure  S6 in Supporting Information S1), that is, without consideration of the contributing glacier area and time duration. The GIC mass loss equates to 1.24-1.52 mm SLE since the LIA termination, that is, between 0.011 and 10.1029/2023GL103950 6 of 9 0.013 mm yr −1 . We compute that the SLE of GIC mass loss between 2000 and 2019 has been 0.08 mm yr −1 using the data set of Hugonnet et al. (2021). For comparison, the mean annual mass loss of the GrIS over the last 20 years has been an order of magnitude more at 2.7 mm yr −1 (Muntjewerf et al., 2020).
Overall, the spatial pattern of GIC mass balance has transitioned to become more spatially heterogeneous between 2000 and 2019 compared with between the LIA termination and 2015 ( Figure 3). Specifically, the range of mean mass balance in the three west regions has gone from −0.11 m w.e. yr −1 since the LIA termination to −0.28 m w.e. yr −1 during 2000-2019, and the range in the three east regions has gone from −0.18 to −0.52 m w.e. yr −1 (Figure 3). The NO region has experienced a −0.46 m w.e. yr −1 change in mass balance rate (Figure 3).
Marine-terminating GICs had the most negative mass balances since the LIA termination compared with lake-or land-terminating ( Figure 3). However, during the period 2000 to 2019 lake-terminating GICs had the most negative mass balances (Figure 3). We find 24 of these lakes formed between the LIA moraine and the contemporary ice front, thus changing the terminus environment from land-to lake terminating. As a result, lake-terminating GICs have experienced the greatest change in mass balance regime since the LIA termination (Figure 3).

Mass Balance
Studies concerned with Greenland GIC glacier evolution, or of isostatic effects of glacier mass loss for example, will find it useful to consider that a total volume of between 528 and 646 km 3 has been lost since the LIA termination, equating to between 449 Gt and 549 Gt at a mean rate of 4.34 Gt yr −1 . That rate is one order of magnitude less than the 27 to 42 Gt yr −1 reported by Khan et al. (2022) for 2003 to 2021 for GICs in north Greenland. It is also an order of magnitude less than the 75 to 100 Gt yr −1 reported by Kjeldsen et al. (2015) for the entirety of the GrIS since year 1900. However, the areal extent and volume of the GrIS is 1.72 × 10 6 km 2 and 2.9 × 10 6 km 3 the GrIS (Box et al., 2022;their Table 1), which are two and four orders of magnitude greater than 90 × 10 3 km 2 and 11.8 ± 3.7 × 10 3 km 3 , respectively, for all of the GICs combined (Millan et al., 2022). Therefore, although the absolute magnitude of mass loss (and hence global mean SLE contributions) from GICs has been approximately one twentieth of that of the GrIS since the year 1900, the mass loss from GICs is disproportionately large, given their combined size, when compared with mass loss from the GrIS. The greatest uncertainty in our rate(s) of change is the date used for the LIA termination. An earlier date (e.g., 1880) would make a slower rate of change (i.e., 3.7 Gt yr −1 ; Table S3 in Supporting Information S1) and that would suggest a greater change in the rate of mass loss when compared with 2000-2019.
Numerical models of glacier evolution should become mindful of changing terminus environments. Given that we find 242 fewer marine terminating GICs in 2001 than during the LIA termination and that the number of lake-terminating GICs have increased slightly by 24 (Figure 3), we interpret that (a) many GICs that were marine-based during the LIA are now terminating on land, and (b) some GICs have developed ice-marginal lakes at their termini as those lakes have formed in the last few decades and inside LIA moraine ridges, for example, as shown across Greenland by How et al. (2021), Carrivick, How, et al. (2022). The influence of lake effects on glacier mass balance (e.g., Carrivick, Tweed, et al., 2020), in addition to climate forcing, could explain why the greatest changes of mass loss have occurred for lake-terminating glaciers, which is an association that we find exists to a greater or lesser extent in all Greenland regions (Figure 3).

ELAs
Absolute LIA ELAs tend to gain elevation with increasing latitude on the east coast but have quite similar median elevations per region on the west coast ( Figure 2). This latitudinal pattern of LIA ELAs is counter to what would be expected with a sole control of colder air temperatures northwards and demonstrates that effective (solid) precipitation is very important (cf. McGrath et al., 2017 for Alaska), perhaps as important as air temperature, for controlling mass balance sensitivity to climate change across Greenland. Specifically, precipi tation across Greenland is higher in the south due to warmer ocean waters, atmospheric temperatures and more evaporation (Box, 2002), thereby sustaining GIC accumulation zones at lower elevations in southern Greenland than in the north. There is also a continentality effect whereby LIA ELAs were higher inland than on the coast (Figure 2 inset), presumably reflecting similar precipitation gradients to those today (e.g., Taurisano et al., 2004).
Furthermore, topography and hypsometry (rather than solely glacier area) affects glacier response times (cf. Raper & Braithwaite, 2009;McGrath et al., 2017), especially considering that ice caps, which become more prevalent northwards (regions NW, NO, NE), have a top-heavy (skewed to higher normalized elevations) hypsometry that contrasts with that of valley glaciers (cf. Pfeffer et al., 2014).
Spatial heterogeneity in GIC changes ( Figure 2) has recently been noted by Khan et al. (2022) in north and north-east Greenland, by Cooper et al. (2022) in central-east Greenland, and also by Brooks et al. (2022) for their sample of 42 glaciers in south Greenland, even though the latter excluded glaciers with confounding properties such as debris cover, and marine-or lake-termini. Brooks et al. (2022) attributed this spatial variability in glacier responses to local topographically-modified microclimates. Some regions have experienced increased precipitation that can have a stabilizing effect on glacier mass loss (Bjørk et al., 2018). Such spatial diversity in glacier mass balances and ELA shifts must cast considerable doubt on using a small sub-sample of glaciers to reconstruct the spatial variability of climate changes (cf. Carrivick & Brewer, 2004).
Some glacier-specific rises in ELA since the LIA termination to 2015 have been >150 m but region medians are quite modest; most are just a few tens of meters, suggesting GICs are not changing geometry in pace with climatic shifts; that is, climate change has progressed faster than glacier response time. Therefore, we suggest that GICs have committed ice losses, just as the GrIS has (Box et al., 2022). In the central and northern regions this could be due to glacier thermal regime; polythermal and fully-cold-based glaciers recede and thin much more slowly than temperate glaciers, becoming in disequilibrium with climate (see paragraph 11 of Irvine-Fynn et al., 2011). The NO region has the least dispersed histogram of ELAs and therefore has greatest proportion of glacier areal extent distributed at a similar elevation meaning that it is the most sensitive region to a rise in ELA.

Global Comparison
Centennial-scale rates of GIC volume and mass losses since the LIA have been determined for other World regions, for example, a doubling in rate of mass loss for Vatnajökull, Iceland (Hannesdóttir et al., 2015), Patagonia (Glasser et al., 2011) and the Southern Alps, New Zealand (Carrivick, James, et al., 2020), and a ten-fold increase across the Himalaya (Lee et al., 2021). These studies have conducted their analysis via an inventory-style mapping approach, as herein, which mitigates problems of sampling bias and glacier-specific conditions that weaken the link between glacier geometry changes and climate fluctuations, including; (a) glacier size; larger glaciers lose the most absolute volume ( Figure S3 in Supporting Information S1) but typically lose the least volume proportional to their area, (b) terminus environment; marine-and lake-terminating lakes typically have higher mass loss rates than land-terminating glaciers (e.g., Lee et al., 2021;Mallalieu et al., 2021), (c) glacier surface characteristics; predominantly debris-cover effects, and (d) glacier surges, which are a rapid translation of ice mass usually into an ablation area often in the form of a bulge or kinematic wave. The Greenland-wide trebling (at least) in rate of mass loss identified in this study since the LIA is therefore in accordance with these other studies, but the NO region of Greenland stands out for a five times increase in rate, which comparing our results with the results of Khan et al. (2022) has apparently occurred entirely in the last two decades.

Conclusions
Greenland GICs are under-studied in comparison to the GrIS, yet are very sensitive to climate change and are important contributors to global mean sea level. In this study we have determined glacier areal extent changes, reconstructed a LIA ice surface, and calculated ice volume changes since then within the ablation areas of Greenland GICs. We compute volume loss rates that are proportionally higher than areal extent loss rates, which demonstrates the importance of glacier thinning to overall mass loss, especially in northern regions where ice velocities are especially low and glacier termini retreat slowly. Comparison of mass loss rates between the LIA to 2015 and 2000-2019 reveals that glacier mass loss has accelerated over the past 20 years by at least two to three times for most regions and by five times for the NO region. We find that marine influences on GIC mass loss have diminished since the LIA and that proglacial ice-marginal lake effects on glacier mass balances have apparently increased. There has been considerable spatial variability in glacier-specific mass changes and ELA changes, which leads us to caution against inferring regional trends from analyses on only a small number of glaciers. Given the sparse and disparate nature of palaeo-climate proxy data across Greenland, numerical modeling of glacier evolution could make use of our data set to help unravel the complex interplay between regional climate change and local topographic controls on GIC mass loss since the LIA.