Greenland-wide accelerated retreat of peripheral glaciers in the twenty-first century

The long-term response of Greenland’s peripheral glaciers to climate change is widely undocumented. Here we use historical aerial photographs and satellite imagery to document length fluctuations of >1,000 land-terminating peripheral glaciers in Greenland over more than a century. We find that their rate of retreat over the last two decades is double that of the twentieth century, indicating a ubiquitous transition into a new, accelerated state of downwasting. .Observations of glacier response to climate changes prior to the satellite era are sparse. Here the authors use historical aerial photographs to document change in peripheral glaciers in Greenland since 1890, providing enhanced confidence that recent changes are unprecedented on a century timescale.

The long-term response of Greenland's peripheral glaciers to climate change is widely undocumented.Here we use historical aerial photographs and satellite imagery to document length fluctuations of >1,000 land-terminating peripheral glaciers in Greenland over more than a century.We find that their rate of retreat over the last two decades is double that of the twentieth century, indicating a ubiquitous transition into a new, accelerated state of downwasting.
Peripheral glaciers and ice caps (GICs) that are distinct from the Greenland Ice Sheet constitute just ∼4% of Greenland's total glaciated area but contribute a disproportionally large portion (∼14%) of the island's current ice loss 1 .In the first two decades of the twenty-first century, Greenland's GICs lost mass at a rate of -35.5 ± 5.8 Gt per year, a rate outpaced only by Alaskan GICs 2 .However, before the satellite era, measurements of fluctuations in Greenland's GICs remain sparse [3][4][5] .Thus, there is only very limited long-term, century-or longer-scale, context for recent loss [3][4][5] .
Beginning in the 1930s, Denmark's extensive mapping efforts in Greenland led to the collection of >200,000 aerial photographs covering the island's coast.The recent rediscovery of this collection, along with photos obtained by the US military during World War II and Cold War eras, has provided invaluable pre-satellite era GIC observations 3,4 .In this Brief Communication, we combine these twentieth-century historical photos with satellite imagery to extend the limited time frame of GIC records across several regions in Greenland, thus placing their contemporary retreat into a lengthened temporal context.Specifically, we map the ice front positions of 821 glaciers in South, North and West Greenland since the mid-twentieth century and extend to present (2021) 364 existing records of GIC length change that begin in the early to mid-twentieth century from the Southeast 4 , Northwest and Northeast 3 (Fig. 1, Supplementary Table 1 and Supplementary Figs.5-11).This work thereby more than doubles the number of GICs in Greenland with detailed air photo-based twentieth-century records of length fluctuations from ∼500 (refs.3-5) to ∼1,320.To further extend the records, we also document ice extents based on geomorphic evidence in South and West Greenland during the Little Ice Age (LIA) maximum (ca.∼1890-1900) [6][7][8] .Finally, we assess GIC response to changes in climate by comparing regional frontal change rates with long-term (1784-2020) meteorological station records 9 , as well as catchment-specific climate and surface mass balance (SMB) output from the Modèle Atmosphérique Régional (MAR) regional climate model 10 (Figs. 1 and 2, Supplementary Table 2 and Supplementary Fig. 15).Together, these records of frontal variations capture regional changes in climate and GIC mass balance over ∼130 years.We focus our study on GICs that terminate on land and that are distinct from the ice sheet (that is, physically detached or dynamically disconnected; Methods).These GICs are fed by local precipitation, are primarily sensitive to atmospheric climate variability, and respond quickly to climate shifts 3,11 .
Over the full ∼130 years of observation, Greenland's landterminating GICs have undergone substantial and widespread retreat (Supplementary Table 1).However, the acceleration in recession in the twenty-first century stands out as distinct and suggests that response to recent warming has been ubiquitous despite the range of climates and https://doi.org/10.1038/s41558-023-01855-6 Our data also show that recent retreat is largely unprecedented in over a century.Of the GICs studied here (with at least two observations, one in each century), the highest rate of retreat on record occurred within the twenty-first century for ∼87%.Further, of the six regions with records beginning around the turn of the twentieth century (that is, ∼1890-1910), the regional frontal change rate in the most recent observational period is unmatched in more than a century in all but two regions (the South and Northeast; Fig. 1 and Supplementary Table 1).Previous studies of Greenland's GICs that focused on the Southeast, Northeast and Northwest showed that early twenty-first-century changes were exceeded in magnitude by changes in the early twentieth century (in the Southeast; during a period of Early Twentieth Century Warming; ETCW), or following the end of the LIA (in the Northeast) 3,4 .GIC characteristics across Greenland.Across all eight regions, the mean frontal change rate (M) during the twenty-first century (∼2000-2022: M = −14.8± 5.4 m per year) is roughly double that of the twentieth century (∼1890-1999: M = −7.7 ± 5.1 m per year; Welch two-sample t-test, t(39) = −4.65,P < 0.001; Fig. 1).In addition to length change in absolute terms, we calculate twenty-first-century length change in relative terms, as the percentage of length lost relative to late-twentieth-century (1978-1987) length, using the common dataset between regions (Supplementary Fig. 3) 12 .Over the last two decades, we find that South Greenland's GICs lost ∼18.5% of their total late-twentieth-century lengths, while GICs in the North, Central-west, Southwest a , Southwest b , Southeast, Northwest and Northeast lost ∼4.9%, ∼9.5%, ∼7.7%, ∼10.7, ∼9.8%, ∼8.7% and ∼4.8%, respectively (Supplementary Fig. 12).https://doi.org/10.1038/s41558-023-01855-6 .Post-LIA retreat rates have also been nearly matched by recent loss in the South and have been far exceeded in the Centralwest and two Southwest regions (Fig. 1 and Supplementary Table 1).We note, however, that these periods are not directly comparable in terms of measurement uncertainty, climate, or glacier hypsometry.First, the uncertainty in the reported post-LIA rates is higher because constraint on the timing of LIA deglaciation is poor across Greenland, and GIC front position is inferred from undated geomorphic evidence 3 .Second, high rates of retreat post-LIA and during the ETCW can in part be explained due to differing glacier area-elevation configurations (for example, some glacier fronts would have been situated on flatter and lower elevation terrain than at present) 3,4 .This implies a higher sensitivity following the LIA maximum as comparatively more glacier surface area would have been exposed to ablation with a comparable rise in temperature.Lastly, unlike the Arctic's amplified and protracted contemporary warming, driven unequivocally by human influence 13 , ETCW in the Arctic was short lived (∼1919-1932) and ascribed to unusual atmospheric circulation and transport of warm air to the Arctic [14][15][16] .Thus, we consider the high rates of retreat in the twenty-first century even more striking, as Greenland's GICs are in a much more retracted and elevated position today than following the LIA.In addition, given continued rise in atmospheric greenhouse gases and further warming, the already exceptional twenty-first-century rates of retreat reported here are likely to be exceeded in coming decades (for example, ref. 17).
We also find that regional trends in climate and SMB from the regional climate model MAR agree well with our frontal change rate measurements, confirming that Greenland's land-terminating GICs are rapidly responding to climate changes via adjustment of their front positions (Figs. 1 and 2, Supplementary Table 2 and Supplementary Fig. 15).For example, we find that GICs in the South, North, West, Southeast and Northwest experienced a more negative SMB anomaly and higher retreat rate in the most recent observational period compared with that of the first decade of the twenty-first century, while GICs in the Northeast experienced a less negative SMB anomaly and lower retreat rate in the most recent observational period (as predicted by ref. 3  under a positive phase of the North Atlantic Oscillation; Figs. 1 and 2; Supplementary Fig. 14).In addition, in the most recent decade, there is a notable divergence in rate of retreat and SMB between GICs in Northeast and Southwest Greenland (Figs. 1 and 2, and Supplementary Figs. 14  and 15).While summer temperature has generally increased across all regions over the past several decades (with varying magnitudes of change; mean +0.04 °C per year from 1970), there is regional and interdecadal variability in precipitation (Figs. 1 and 2, and Supplementary Figs. 14 and 15).Thus, the East-West divergence probably reflects higher snowfall in East Greenland and drier conditions in West Greenland in the most recent decade (Northeast, Southwest a and Southwest b snowfall and summer temperature anomaly from MAR between 2010 and 2019: +49, −67 and −73 mm water equivalent (w.e.), and +1.5, +1.7 and +1.6 °C, respectively; Fig. 2 and Supplementary Figs. 14 and 15).Likewise, a slowdown in mass loss from Icelandic, Scandinavian and Eastern Greenland glaciers was observed over the last 5-10 years in a recent global-2 and regional-scale study 18 and was ascribed to wetter and/or cooler conditions in the North Atlantic region 2,18 .However, despite this, we find a high rate of frontal change in the Southeast between 2010 and 2021, suggesting the increase in accumulation in East Greenland was not sufficient to impede retreat of GICs there (which could, in part, be because of differing hypsometry between glaciers in the Northeast and Southeast; Supplementary Figs.13b and 15).
In summary, the ubiquitous recent acceleration in retreat of land-terminating GICs across Greenland's climate zones adds to a growing body of evidence for accelerated twenty-first-century glacier change across the Arctic (for example, refs.1,19-28)-a striking consequence of anthropogenic climate change 29,30 .The unique century-plus temporal context provided by this study suggests that the magnitude of GIC retreat over the last two decades is largely exceptional for Greenland, but also confirms that these glaciers respond to climate on short time scales, underscoring the urgency for immediate climate action to limit global temperature rise e.g. 17 .

Image sources, rectification and digitization of glacier fronts
Glacier length data were acquired from historical aerial photographs (1943, 1951, 1953, 1958 and 1969), a 2-m-resolution late-twentieth-century orthophotograph (1978-1987; Supplementary Fig. 3) 12 , and Landsat, Sentinel and RapidEye satellite imagery (2000, 2001, 2002, 2008, 2010, 2013, 2014, 2017, 2021 and 2022).Data collection years were chosen on the basis of availability of quality late-summer, low-cloud and low-snow-cover imagery.The historical air photos were scanned and digitized from paper copies on a Canon photogrammetric precision scanner at a resolution of 600 dpi (Supplementary Fig. 1).The original photos are archived at the Danish National Archives in Copenhagen, Denmark.Vertical-only air photographs were selected for use.Once in digital form, the historical air photos were oriented, rectified to a standard shapefile using the image fiducials, and cropped in ArcGIS.In the South, Agisoft Professional was used to generate orthophotos from a set of overlapping historical images (Supplementary Fig. 2).In the West, a spline transformation was employed to rectify all historical air photos.All orthophotos were georeferenced in ArcGIS using tie points near the glacier fronts for co-registration to the 2-m orthomosaic, which we use as the reference dataset (Supplementary Fig. 3).
Glacier fronts in the twentieth-century historical imagery were manually digitized in ArcGIS as polylines at the transition from ablated ice to bare bedrock.In the South, twenty-first-century glacier fronts were digitized via the Google Earth Engine Digitization tool 31 in Landsat and Sentinel satellite imagery and exported as georeferenced shapefiles.All digitized glacier fronts (including relevant metadata such as sensor and acquisition date) were then combined in ArcGIS for each glacier.In the West, twenty-first-century glacier fronts were manually digitized in ArcGIS.In the North, glacier fronts in the historical (the twentieth-century orthomosaic) and Landsat imagery were manually digitized in ArcGIS as single point features at the glacier front, along the centreline.
When possible, in the South and West, the maximum extent of glaciers during the LIA was mapped using geomorphic evidence (that is, moraines and trimlines) as observed in the 2-m orthomosaic product.
We follow ref. 3 and do not digitize the fronts of GICs that show characteristics of surging (such as looped moraines, collapsed accumulation zones, crevasse squeeze ridges in forefields, and potholes).We also removed GICs from further study if surge-type characteristics were recognized during data processing or analysis (for example, substantial advances >50 m per year).We note, however, that evidence of surging may not always be easily identifiable, and thus, a few surge-type GICs may remain in our analysis.

Frontal change rate and late-twentieth-century GIC length measurements
In the South, frontal change rates were calculated using the multi-centreline method via the Margin change Quantification Tool (MaQiT) 31 .This method allows for the spatial variability in margin change to be quantified by defining multiple centrelines spanning the width of the margin (every 10 m across), which results in several one-dimensional measurements of change across the width of the entire margin 31 (Supplementary Fig. 4).MaQiT also requires a shapefile of the main centreline for each glacier.Main glacier centrelines were manually digitized using GEEDiT in Sentinel imagery.In the West, frontal change rates were calculated using the centreline method and the Glacier Termini Tracking tool for ArcGIS 32 .In the North, frontal change rates were calculated via the centreline method using the Euclidean distance between point features, divided by the number of years between observations.We divide the full observational periods, ∼1900-2021, ∼1890-2022 and 1978-2022 in the South, West and North regions into five (South), six (Central-west and Southwest b ), eight (Southwest a ) and four (North) subperiods, and calculate the regional frontal change rate (m per year) for each subperiod by averaging the front change of all measured individual glaciers (Fig. 1 and Supplementary Table 1).
The late-twentieth-century lengths of all GICs were measured using the only shared dataset between regions: the series of high-resolution (2 m) orthorectified mosaics derived from aerial photographs, which, depending on the region, captures GIC status between 1978 and 1987 (ref.12) (Supplementary Fig. 3).GIC lengths were manually digitized in ArcGIS along the flowline from the glacier headwall (or the centre point in the case of ice caps) to their late-twentieth-century fronts.Relative front change in the twenty-first century is reported as the percentage of the late-twentieth-century length lost (the total twenty-first-century front change (m) divided by the late-twentieth-century length (m), multiplied by 100) (Supplementary Fig. 12).
To extend the existing Southeast 4 , Northwest and Northeast 3 frontal change datasets to present, we digitize the GICs 2021 fronts, perpendicular to the centreline using Sentinel imagery.For each glacier, we measure the distance along the centreline between the late-twentieth-century front (via the orthomosaic product) and the 2021 front, then subtract the reported total front change (m) between the late twentieth century and the first decade of the twenty-first century.We then divide the remaining length by the number of years between the first twenty-first-century observational period and 2021 (Supplementary Fig. 4).

Frontal change rate uncertainties
Uncertainty in frontal change rates arise from several sources: the pixel size of the imagery, image rectification, and the manual digitization of the glacier fronts.We calculate the uncertainty in frontal change rate for each subperiod following ref.3: using the maximum digitizing uncertainty and the maximum mean error from the orthomosaic (that is, 8.6 m).The retreat rate error (1) is given as the root of the sum of squares of the digitizing and rectification uncertainties of the two observations used in the observation period: where pa and pb are the digitizing and rectification errors for the given observations, and ∆t is the time interval in years between the two observations 3 .We use a maximum digitizing error for each image source (Supplementary Table 1).This was measured by Bjørk et al. 4 , by digitizing objects with known positions with the orthomosaic as reference.The digitizing error for the orthomosaic is given as one pixel size, or 4 m (see Bjørk et al. 4 ; Supplementary Table 2).The rectification error is included in all frontal change rate error calculations that rely on the orthomosaic product (which is used as the reference for all other rectification).Bjørk et al. 4 calculated a mean error of 8.6 m for the orthomosaics (see Bjørk et al. 4 ; Supplementary Information Section 2.3).Thus, the frontal change rate uncertainties should be interpreted as the combined error in the measurements.We note that longer observational periods result in a relatively smaller frontal change rate error.Following ref. 3, we add ±5 years to the post-LIA observational periods to account for poor constraint on the timing of regional LIA deglaciation.
We emphasize that precise constraints on the regional timing of LIA deglaciation are sparse in the South and three western regions included in our study (Central-west and Southwest a,b ).However, we assign the onset to ∼1900 in the South and ∼1890 in the West, based on available historical accounts and lichenometric data [6][7][8]33,34 . In sothern Greenland, refs.7,8 give historical accounts of the fluctuations of several glaciers in the Julianehåb-Godthåb districts and conclude that most attained their maximum extensions between 1890 and 1900, after which a period of recession set in, irrespective of the type of glacier.For instance, Kiagtût sermia (official name Kiattuut Sermiat), located to the northeast of Narsarsuaq, retained its maximum extension as late as 1876, and still in 1890 showed only slight shrinkage 7 .

Brief Communication
https://doi.org/10.1038/s41558-023-01855-6 Ref. 8 describes two glaciers at the head of Tasermiut fjord, northeast of Nanortalik.Accounts of the southernmost glacier, Sermitsiaq, suggest that it began to advance sometime before 1833, as at that time it already had an extent close to the historical maximum and held there until it began to thin ∼1900 (ref.8).The eastern-most glacier Sermeq reached its maximum extent by ∼1888, and rapidly receded after 1900 (ref.8).In addition, two ice lobes of the Julianehåb ice cap (Kujatdleq and Jespersens Brae) are thought to have attained their maximum extensions around the same time, between 1890 and 1900 (ref.7).In the West, ref. 7 reported fluctuations of ten peripheral glaciers on Nuussuaq and reported that their maximum extent occurred between 1880 and 1890.Similarly, the majority of the GICs investigated in the Southwest a region reached a position near their LIA maximum extent by the mid-1800s to 1900.GICs located south of Nuuk (and of the region Southwest b ) are reported to have reached their maximum LIA extents ∼1870-1880 (refs.6-8).These historical observations are further supported by lichenometric dating that indicates some outlet glaciers of the Sukkertoppen Iskappe reached their LIA maximum extents between 1870 and 1890 (refs.33,34).We also acknowledge that LIA moraines are probably different in age between individual glacier systems.As such, assigning a single date to the LIA maximum, or the onset of deglaciation for a large region is unrealistic.This imposes an additional and difficult-to-quantify uncertainty on the retreat rates that use the LIA maximum as one end of the observational period.We note that some studies, notably those that use cosmogenic exposure dating of glacial moraines, generally suggest an earlier and more variable LIA maximum timing.For example, in southern Greenland, ref. 35 provide evidence for a much earlier LIA glacial advance at ∼0.71 ka.In this case, an assigned age of ∼1900 for the timing of LIA deglaciation in the South would result in artificially high retreat rates for the first observational period.If the onset of LIA deglaciation occurred earlier than reported here, thereby making our estimated rate of retreat in the earliest observational periods artificially high, this would ultimately bolster the conclusion that twenty-first-century retreat is unprecedented on a century time scale.

Climatological data and glacier-specific characteristics
We compare the regional GIC frontal change rates with climatological data and modelled SMB.Specifically, we use summer temperature and accumulated precipitation (winter half year: October, November, December, January, February, March) from seven Danish Meteorological Institute (DMI) stations with long-term observations located closest to the regions of interest 9 .In addition, we extract air temperature (summer, °C), snowfall and SMB annual sum (mm w.e.) between 1950 and 2019 from the regional climate model, MAR 10 v.3.12.MAR v3.12 is forced at the boundaries with ERA5 atmospheric re-analysis.We use the monthly mean values downscaled to 1 km and take an average of the variables over each glacier catchment.
We also collect information on glacier-specific characteristics including elevation (minimum, median and maximum), area, aspect and slope from the Randolph Glacier Inventory (RGI) version 6.0 (ref.36) and calculate regional hypsometry using the ArcticDEM 32-m mosaic product 37 and GIC polygons from the RGI version 6.0 (ref.36) (Supplementary Fig. 13a,b).Although we focus our study on GICs that are physically detached from the main ice sheet and that terminate on land, we include some GICs that have a Connect value of 1 or 2 in the RGI, indicating a weak or strong connectivity level to the ice sheet.We also include a few GICs that have a TermType value of 1 in the RGI, indicating a marine terminating environment.These represent cases where either (1) the glacier or ice cap was initially marine-terminating and has since become land-terminating, and/or (2) the glacier or ice cap has several lobes, of which some are marine-terminating, and some are land-terminating.In the former case, it is possible that a change in the terminal environment (that is, marine-to land-terminating) has affected the rate of retreat over time; however, we are unable to account for this due to limited observational periods.In the latter case, change in one or more land-terminating fronts was measured, but the glacier or ice cap is labelled as dominantly marine-terminating in the RGI.Of the glaciers and ice caps included in the analysis with unique RGI IDs, 86.9% have a Connect value of 0; 10.4% have a Connect value of 1; 2.7% have a Connect value of 2; 95.5% have a TermType of 0; and 4.5% have a TermType value of 1 (Supplementary Fig. 13a).

Fig. 1 |Fig. 2 |
Fig. 1 | Comparison of regional GIC frontal change rates and climate over the past ∼130 years.a, Locations of GICs and nearby DMI stations.b, Frontal change rate (m per year) for GICs (n = 1,185) in eight regions divided into observational periods: Northwest (NW, n = 139) 3 , Central-west (CW, n = 59), Southwest a (SW a , n = 53), Southwest b (SW b , n = 61), South (S, n = 61), Southeast (SE, n = 28) 4 , Northeast (NE, n = 197) 3 and North (N, n = 587).Error bars represent the combined measurement error given the maximum digitizing uncertainty and the maximum mean error from the orthomosaic (Methods).The increased transparency in the six reported post-LIA rates signify their higher uncertainty due to poor constraint on the timing of LIA deglaciation across Greenland.Dashed black lines show the mean frontal change rate across all observations in An exception is in the Northeast, where retreat rates in the twenty-first century have not yet exceeded the rate of retreat in the early twentieth century (1910-1932: −22.1 ± 6.4 m per year; 2000-2013: −12.2 ± 1.8 m per year; and 2013-2021: −11.8 ± 2.4 m per year)