Multi-decadal retreat of marine-terminating outlet glaciers in northwest and central-west Greenland

The retreat and acceleration of marine-terminating outlet glaciers in Greenland over the past 2 decades have been widely attributed to climate change. Here we present a comprehensive annual record of glacier terminus positions in northwest and central-west Greenland and compare it against local and regional climatology to assess the regional sensitivity of glacier termini to different climatic factors. This record is derived from optical and radar satellite imagery and spans 87 marine-terminating outlet glaciers from 1972 through 2021. We find that in this region, most glaciers have retreated over the observation period and widespread regional retreat accelerated from around 1996. The acceleration of glacier retreat coincides with the timing of sharp shifts in ocean surface temperatures, the duration of the sea-ice season, ice-sheet surface mass balance, and meltwater and runoff production. Regression analysis indicates that terminus retreat is most sensitive to increases in runoff and ocean temperatures, while the effect of offshore sea ice is weak. Because runoff and ocean temperatures can influence terminus positions through several mechanisms, our findings suggest that a variety of processes – such as ocean-interface melting, mélange presence and rigidity, and hydrofracture-induced calving – may contribute to, but do not conclusively dominate, the observed regional retreat.

1 Introduction

The Greenland Ice Sheet has lost significant mass over the last few decades (Enderlin et al., 2014; Shepherd et al., 2020) as many of its glaciers have retreated (Hill et al., 2017; Howat and Eddy, 2011; King et al., 2020; Moon and Joughin, 2008; Murray et al., 2015), and ice flow has accelerated (Joughin et al., 2010; Moon et al., 2012; Rignot and Kanagaratnam, 2006). Recent Greenland ice loss has contributed to rates at times approaching 1 mm a⁻¹ of global sea-level rise (Shepherd et al., 2020), with the contribution from northwest and central-west Greenland combined representing nearly half of the cumulative contribution from Greenland to sea-level rise since 1972 (Mouginot et al., 2019). Although surface mass balance has dominated Greenland's mass loss in the past 2 decades, over half of the mass loss in northwest and central-west Greenland is currently due to ice discharge (Mouginot et al., 2019), which has accelerated since 2000 in this region (King et al., 2020).

Changes in ice discharge are often related to changes to glacier terminus positions, with terminus retreat into deeper water driving acceleration and upstream thinning (Howat et al., 2008; Joughin et al., 2008b). Because of the relationship between the terminus position and calving rate, initial perturbations that increase calving can trigger further glacier retreat. While glacier retreat and acceleration are generally linked to changes at the terminus, it remains unclear which processes are most responsible for controlling perturbations to calving rates and subsequent terminus retreat (Straneo et al., 2013).

The recent acceleration and retreat of Greenland outlet glaciers have been attributed to warmer ocean temperatures (Fahrner et al., 2021; Holland et al., 2008; Howat et al., 2008; Morlighem et al., 2016; Rignot et al., 2012; Slater et al., 2019; Wood et al., 2021), to changes in the characteristics of sea ice and mélange – icebergs bound by a sea-ice matrix (Amundson et al., 2010; Carr et al., 2013; Cassotto et al., 2015; Foga et al., 2014; Joughin et al., 2008a; Moon et al., 2015), and to increased melting and crevassing associated with warming air temperatures (Benn et al., 2007; Nick et al., 2013; van der Veen, 1998). In turn, responses to these forcings are modulated by geometric factors associated with individual glaciers such as bed topography and fjord width (Carr et al., 2015; Catania et al., 2018; Felikson et al., 2021; Schild and Hamilton, 2013); for many glaciers, these modulating factors necessitate detailed records of terminus position changes in order to identify the importance of different forcing mechanisms on decadal-scale outlet glacier changes across a large area.

Most previous studies of Greenland outlet glacier terminus positions have been temporally or spatially limited. Some studies cover the entire ice sheet for over a decade but map termini only decadally or in non-consecutive years (Howat and Eddy, 2011; Moon and Joughin, 2008). Other studies map termini more frequently but only for a small sector of the ice sheet (Bjørk et al., 2012; McFadden et al., 2011; Moon et al., 2015) or for a few specific glaciers (Holland et al., 2016; Joughin et al., 2008b; Larsen et al., 2016; Motyka et al., 2017; Schild and Hamilton, 2013). Murray et al. (2015) mapped terminus positions at high spatial and temporal resolutions but only for a single decade. More recent studies have attempted to fill these observational gaps with ice-sheet-wide analyses of annual terminus positions spanning multiple decades (Fahrner et al., 2021; King et al., 2020; Wood et al., 2021), although they come to differing conclusions about the drivers of observed terminus retreat.

In this paper, we analyze glacier change at a high spatiotemporal resolution by constructing a multi-decadal (49 years) record of annual terminus positions for 87 marine-terminating outlet glaciers in northwest and central-west Greenland. This record allows us to identify the behavior of individual glaciers as well as regional trends in the magnitude and timing of glacier retreat. We compare this regional behavior with climate data – sea-surface and subsurface temperatures; sea-ice concentration; and ice-sheet surface mass balance, precipitation, melt, and runoff – to assess the relative influence of different forcing mechanisms on multi-decadal and annual terminus retreat in this sector of Greenland. It is important to note that we do not account for the effect of geometric factors such as bed topography on modulating the retreat of individual glaciers, as we instead focus on terminus retreat in connection with regional climate trends.

2 Data and methods

We used 455 synthetic aperture radar (SAR) mosaics and optical satellite images to trace annual terminus positions for 87 marine-terminating outlet glaciers in northwest and central-west Greenland (68.9 to 78.2^∘ N) (Fig. 1) from 1972 through 2021.

2.1 Satellite images

We digitized terminus positions from SAR images acquired by the European Copernicus Programme's Sentinel-1A and Sentinel-1B and the Canadian Space Agency's Radarsat-1 satellites for 13 of the 49 years in our record. These radar satellites are able to image the surface regardless of clouds or darkness. For the most recent 7 years (2014/15 through 2020/21), we traced terminus positions in mosaics of Copernicus Sentinel-1A and Sentinel-1B images (Joughin et al., 2016a). These mosaics are typically created from images acquired in the early February of their respective years and have a 25 or 50 m product-dependent resolution. Radarsat-1 mosaics were used to trace terminus positions for six winters (2000/01, 2005/06–2008/09, 2012/13) (Joughin et al., 2015; Moon and Joughin, 2008). These mosaics are formed from images collected from October through March and have nominal resolutions of 20 m.

For the remaining 36 years in our record, we used imagery from all Landsat missions to map terminus positions. We also used Landsat imagery to map individual glacier termini where they were missing from or indiscernible in the SAR imagery. We prescreened images in the USGS Global Visualization Viewer (GloVis) to confirm that glaciers of interest were not obscured by clouds and that the images were georeferenced well. We continued to use Landsat 7 images after the instrument's scan-line corrector failure in 2003 as Landsat 5 images were typically unavailable over our study area during this period, although we kept only images which retained a sufficient number of data to map the terminus. In those images, we digitized across the scan-line corrector gaps when they crossed a glacier terminus (Supplement Fig. S1); comparison of temporally close images with different data gaps indicates that this method only marginally increases errors. In all, the scan-line corrector gaps affected 348 terminus traces (9.7 % of the dataset). Due to winter darkness, we could not select images from the same time of year as the SAR mosaics; instead, we chose images as close to winter as possible, with a preference for spring over autumn to capture a more winter-like state to reduce the effects of seasonal variation (Fig. 2). As a result, the majority of the Landsat images we used were collected in spring (March–May). Because of the difficulty in finding sunlit, cloud-free imagery, however, we had to use data from other periods, so all months except January and December are represented. For some glaciers there are several years with missing data. The image resolution ranges from 15 m (Landsat 8 panchromatic band) to 60 m (Landsat 1–5 Multispectral Scanner). We digitized termini using the panchromatic band for Landsat 7 and Landsat 8 and using a single band that provided high image contrast (typically band 2) for Landsat 1–Landsat 5.

2.2 Terminus positions

Using ArcGIS, we manually traced annual terminus positions from the radar and Landsat imagery. This dataset builds on a preexisting dataset covering six winters between 2000 and 2013 (Joughin et al., 2015; Moon and Joughin, 2008). The study area ranges from Saqqarliup Sermia (68.9^∘ N, 50.3^∘ W; ∼ 35 km southwest of Jakobshavn Isbræ) to Bamse Gletsjer (78.2^∘ N, 72.7^∘ W) (Fig. 1). Individual glacier names and coordinates are detailed in Supplement Table S1. To focus our analysis on glaciers that produce substantial ice discharge, we limited our analysis to marine-terminating outlet glaciers that are at least ∼ 1.5 km wide at the terminus and are flowing at a rate of at least ∼ 1000 m a⁻¹. We traced each glacier's terminus position once per hydrological year (1 September through 31 August; Ettema et al., 2009) in every year for which suitable imagery was available. We used winter or near-winter imagery whenever possible as indicated by Fig. 2.

Errors in terminus position may arise from both the imagery used and the digitization process. The primary sources of errors introduced by the image data are the uncertainties in position after orthorectification and georeferencing. To reduce such errors, candidate images were compared with control images and discarded if they were noticeably offset or distorted. Manual digitization also introduces errors, which are exacerbated by poor image resolution and image artifacts (such as shadows or an indistinct transition from glacier to mélange) (Joughin et al., 2015). If a terminus position was ambiguous in one image, it was flagged during tracing and compared with close-in-time images from the same satellite platform or other satellite platforms when possible. Digitization errors are typically comparable to the image resolution; for example, a 25 m error for a Landsat 7 image with a 30 m resolution (Moon et al., 2015)

Figure 1 Map of study area, showing individual glacier locations (white) and points where ocean data were acquired (orange, numbered). Glacier labels are (a) every 10th glacier, numbered, and (b) specific glaciers named in the paper, abbreviated as follows: Jakobshavn Isbræ (JI), Alianaatsup Sermia (AS), Sermeq Avannarleq (SA), Store Gletsjer (ST), Sermilik Isbræ (SI), Kangilleq (KA), Kangerlussuup Sermia (KS), Rink Isbræ (RI), Upernavik Isstrøm (UI), Naajarsuit Sermiat (NS), unnamed glacier no. 26 (26), Alison Glacier (AL), unnamed glacier no. 36 (36), Kjer Gletsjer (KJ), Kong Oscar Gletsjer (KO), unnamed glacier no. 55 (55), unnamed glacier no. 66 (66), Tracy Gletsjer (TR), Verhoeff Gletsjer (VH), and Morris Jesup Gletsjer (MJ). Basemaps are (a) bed topography from BedMachine Greenland V3 (Morlighem et al., 2017b, a) and (b) ice-sheet velocity from MEaSUREs GIMP (Joughin et al., 2016b, 2018a).

In addition to the errors associated with digitization of the images, there is additional uncertainty introduced by sampling seasonal variability at different times of year. For the terminus positions mapped from Landsat imagery, it was not possible to obtain sunlit and cloud-free images over each glacier at the same time every year (Fig. 2). This timing variation complicates the year-to-year comparison at a glacier, which might include seasonal variability in the terminus position. Such seasonal variability could produce some short-term deviations in our data; for example, Moon et al. (2015) found a mean annual terminus range of 610 m for a subset of our study glaciers. However, the data collection season is largely consistent across each of our glaciers (Fig. 2), and because these seasonal errors are not independent, they tend to cancel out over longer periods. Since we mostly focus on decadal-scale trends, issues of seasonal sampling should not greatly affect our results.

Figure 2 Years in which a terminus position was observed (filled circles) for each glacier, colored by the season of the observation. If a terminus position was not observed, an estimated position was interpolated (open circles) from prior and subsequent observations to use in annual analyses. Glacier nos. 76–78 are not shown because they are not ice-sheet outlet glaciers and therefore are not included in this study.

2.3 Glacier change measurements

We calculated glacier area change over time using the box method (Moon and Joughin, 2008). Each glacier has a static, open-ended reference box (polygon) that approximately delineates the main region of ice flow. The box sides are roughly parallel to ice flow, and the “back” of the box is perpendicular to ice flow at an arbitrary location up-glacier of the extent of maximum observed terminus retreat (see example in Fig. 3). Where a terminus trace intersects the open end of the box, the polygon is closed, and the area of that polygon represents the glacier's reference area at that point in time. Repeating this process for each terminus trace for a glacier forms a time series of reference areas, from which we determine the annual area change. While the areas are arbitrarily determined by the box size, the area differences between successive terminus traces represent the annual gain or loss of area. By focusing on the area change between measurements rather than the absolute area of each measurement, we can ignore the arbitrariness of how the boxes are drawn, as well as any small regions of stagnant ice that may be included within the box boundaries. There is a small error associated with these area-change measurements because the boxes do not completely conform to the glacier sides.

As length change can provide a more intuitive measure of retreat than area change, we determine the nominal length change by scaling the area change by the average width of the box to obtain an approximate length change. This measurement should be interpreted as a proxy for length change rather than an exact measurement. Compared to the centerline method of measuring length change, this method is less sensitive to uncertainties in the centerline position at the terminus.

In years when no observations were made, we linearly interpolated between the prior and subsequent observations to estimate glacier length and area during the missing years (Fig. 2). For glaciers with missing observations at the beginning of the record, we did not interpolate prior to the first observation in the record. The largest temporal gap interpolated is 9 years (1975/76 through 1983/84), at Sermeq Avannarleq (no. 8). Most temporal gaps are in the 1970s and early 1980s, with additional gaps into the 1990s for some high-latitude glaciers, and there are no temporal gaps after 2002.

Because it is difficult to compare changes in area between glaciers of different sizes, we also determined the percent area change over time for each glacier. For each individual glacier area time series, we normalized the minimum observed glacier area to 0 and the maximum observed glacier area to 1 and linearly scaled the other measured areas between those set points. This method normalizes every glacier's area change to the same 0–1 scale, allowing straightforward comparison of size changes between different glaciers. Because the equivalent length is simply a scaled area, the results are identical for area and length.

To pinpoint the timing of changes in glacier area and length, we performed a break-point analysis on the time series for each glacier. We fit each time series with a piecewise linear function with two segments (Jekel and Venter, 2019); the break point between the two segments corresponds to a year in the time series that is the best fit for approximating the time series as a two-segment piecewise linear function. We performed an F test to determine whether the piecewise linear function was statistically a better fit than a single linear regression and discarded any break points that were not statistically significant (p<0.05).

2.4 Climate data

Earlier work has shown that climate-related processes including terminus ablation and undercutting (driven by ocean warming) (Holland et al., 2008; Motyka et al., 2011; Slater et al., 2019; Rignot et al., 2012; Wood et al., 2021), mélange rigidity (driven by changes in ocean temperature, sea-ice concentration, and/or runoff) (Carr et al., 2017; Joughin et al., 2020; Moon et al., 2015; Sohn et al., 1998; Todd and Christoffersen, 2014), and enhanced hydrofracture (driven by changes in runoff and surface mass balance) (Benn et al., 2007; Nick et al., 2010) may affect terminus position. Hence, we acquired several ice-sheet and oceanographic datasets in order to compare our glacier terminus position changes with climatic factors. For each variable and dataset, we considered the annual and decadal mean values. For the ice-sheet variables, we considered the annual mean at a fixed location near the front of each glacier, as well as the population annual and decadal means. For the oceanographic data, due to lack of reliable long-term data in narrow fjords, we used offshore observations as a proxy for fjord conditions. We used eight points offshore (Fig. 1), selected to be representative of clusters of glaciers: Disko Bay, Uummannaq Fjord, Upernavik Icefjord, south Melville Bay–Wilcox Head, central Melville Bay, north Melville Bay–Cape York, Wolstenholme Bay–Thule, and Inglefield Fjord (coordinates for each point are detailed in Supplement Table S2).

Ice-sheet surface mass balance, snowfall, rainfall, meltwater production, and runoff were extracted from the Modèle Atmosphérique Régional (MAR), Version 3.11, on a 6 km × 6 km grid over the period 1979–2020 (Fettweis et al., 2017).

Sea-surface temperatures came from the NASA JPL Estimating the Circulation and Climate of the Ocean (ECCO) consortium ocean circulation model, Version 5 (Forget et al., 2015; Zhang et al., 2018), on a $\mathrm{1}/\mathrm{3}$^∘ × $\mathrm{1}/\mathrm{3}$^∘ grid over the period 1992–2017, and the merged Hadley-OI sea-surface temperature and sea-ice concentration dataset (Hurrell et al., 2008; Shea et al., 2020), on a 1^∘ × 1^∘ grid over the period 1972–2020.

We also obtained subsurface temperatures from ECCO for most of our study area. In the Disko Bay area we used field observations from the ICES Dataset on Ocean Hydrography (ICES, 2014), collected over the period 1977–2016 by the Greenland Institute of Natural Resources. For both datasets, subsurface temperatures were averaged over the lower 60 % of the water column in order to capture warm water at depth while excluding the near-surface mixed layer. For the ECCO data, we applied a temperature correction to bring the data in line with Wood et al. (2021) (see Supplement).

Sea-ice concentration came from the merged Hadley-OI dataset, as well as the NOAA/NSIDC Climate Data Record of Passive Microwave Sea Ice Concentration, Version 3 (Meier et al., 2017; Peng et al., 2013), on a 25 km × 25 km grid over the period 1978–2019. For the sea-ice concentration, we calculated the annual duration of the period when sea-ice concentration exceeds 15 %, as well as the seasonal mean sea-ice concentration.

3 Results

We produced a comprehensive multi-decadal dataset of terminus positions for glaciers in our study area and collected climate data both near the termini of these glaciers and in the ocean offshore of clusters of these glaciers.

3.1 Terminus positions

We created a dataset of 3606 annual terminus positions for our 87 selected glaciers from 1972/73 through 2020/21 (see example in Fig. 3). The median number of annual observations per glacier is 41, and nearly all glaciers were observed in 38 to 46 of the 49 years examined. Only three glaciers have fewer (33–34) observations; these glaciers are located just south of Thule Air Base and have limited imagery available in the 1980s and early 1990s. After interpolating area changes between glacier observations, the first year with either an observation or an interpolation available for each glacier is 1977/78. Therefore, our glacier analyses start in this year except where noted.

Figure 3 Illustration of terminus traces and box method for Kjer Gletsjer (no. 42). Each trace intersects the glacier's box (green) to form a closed polygon. The areas of sequential polygons are differenced to create a time series of glacier area change. The basemap is a Sentinel-1 SAR mosaic from 4–9 February 2020.

3.1.1 Terminus behavior

The majority of the glaciers in our study area retreated between 1977 and 2021. The cumulative area loss of all of these glaciers is 1067 km², equivalent to a cumulative retreat of 287 km. The individual area and length changes are plotted in Fig. 4 with glaciers with the greatest change broken out separately. We identify these dominant glaciers as those with a net change falling more than 2 standard deviations beyond the population mean, which yields four glaciers that dominate the area change: Jakobshavn Isbræ (no. 3; −130.9 km²), Alison Glacier (no. 35; −59.4 km²), Kjer Gletsjer (no. 42; −81.5 km²), and Tracy Gletsjer (no. 81; −58.7 km²). These glaciers together are responsible for 31.0 % of the total area loss. Five glaciers dominate the length change: Jakobshavn Isbræ (−16.9 km), Alison Glacier (−14.7 km), Kjer Gletsjer (−14.7 km), Tracy Gletsjer (−14.0 km), and one unnamed glacier (no. 36; −11.1 km). These glaciers are cumulatively responsible for 24.9 % of the total retreat. For the remaining glaciers, the mean area change is −8.9 km² and the mean length change is −2.6 km. Net area and length changes for individual glaciers are detailed in Supplement Table S3.

A total of 15 glaciers were stable in that the net changes in area were within 2 standard deviations of their respective observed interannual variability, and no glaciers advanced significantly over the observation period. The 15 stable glaciers are as follows: Alianaatsup Sermia (no. 7); Sermeq Avannarleq (no. 8); Store Gletsjer (no. 9); Sermilik (no. 11); Kangilleq (no. 12); Kangerlussuup Sermia (no. 16); Rink Isbræ (no. 17); Upernavik Isstrøm (no. 21); Naajarsuit Sermiat (no. 25); Kong Oscar Gletsjer (no. 51); Verhoeff Gletsjer (no. 86); Morris Jesup Gletsjer (no. 87); and unnamed glacier nos. 26, 55, and 66.

Figure 4 Cumulative (a) area change and (b) length change for large glaciers that dominate the total observed change. Glacier changes in each plot are labeled using the abbreviation scheme in Fig. 1. The color map follows the order of the glaciers (from south to north), so glaciers with similar colors are spatially closer together. For all other glaciers, cumulative (c) area change and (d) length change are provided but not distinguished by glacier. The final net (e) area change and (f) length change are reported as histograms.

3.1.2 Timing of change

Because a small number of glaciers tend to dominate the overall trends in glacier changes, we scaled each glacier's length and area to a 0–1 scale as described above in order to consistently compare the relative timing of each glacier's behavior (Fig. 5a). The resulting data are noisy, so to better identify overall trends, we computed the mean and standard deviation of the scaled glacier changes (blue overlay in Fig. 5a). Of the observed mean size change (area or equivalent length), ∼ 4 % occurs each decade before 1996 and ∼ 31 % occurs each decade after 1996. The population of glaciers notably readvanced in 2017 and 2018 but began to retreat again in 2019.

We performed a break-point analysis (see Sect. 2.3) for each glacier's area time series (Fig. 5b). The step change in retreat rates in around 1996 for the normalized time series is consistent with the time series break points for individual glaciers. The most common break-point years were 1996/97 and 1997/98, with additional smaller peaks in break points in the mid-1990s and mid-2000s. The break-point years for each glacier are detailed in Supplement Table S3.

Figure 5 Break-point analysis to identify the onset of increased terminus retreat rates. (a) Glacier cumulative area change is normalized such that for each glacier the greatest observed extent is 1 and the smallest observed extent is 0 so that glaciers of all sizes are placed on a common scale. Each glacier's cumulative change is then approximated as a two-segment piecewise linear function, and for cases where the piecewise linear function is statistically better than a linear regression, we compute (b) a histogram of the break-point year between the two segments.

3.2 Climate data

All of the climate variables that we surveyed showed long-term trends consistent with regional climate warming. Most of these variables, in most of the locations that we sampled, also showed a sharp change in their long-term trends between the 1990s and 2000s.

Figure 6 summarizes the data from the MAR climate model. While there is considerable interannual variability, the data indicate a net decrease in surface mass balance (Fig. 6a) and a net increase in both meltwater production (Fig. 6b) and runoff (Fig. 6c), with corresponding increases in the annual number of meltwater production days (Fig. 6d) and runoff days (Fig. 6e), that have occurred since the late 1970s. The mean surface mass balance was steady in the 1980s and 1990s but dropped by 0.48 m a⁻¹ between the 1990s and the 2000s, primarily due to substantial increases in meltwater production (+0.42 m a⁻¹) and runoff (+0.48 m a⁻¹). By contrast, there were only small changes to the mean snowfall (−0.01 m a⁻¹; Fig. 6f) and the mean rainfall (+0.02 m a⁻¹; Fig. 6g) from the 1980s to the 2010s. For the 2010s, while the annual surface mass balance anomaly was positive in 2013, 2017, and 2018, it reached its lowest values for the full record in 2019 and 2020. These extremes are reflected in coincident meltwater production and runoff anomalies. For this decade as a whole, the decadal mean surface mass balance anomaly is slightly more negative than that of the 2000s and the decadal mean meltwater production and runoff continued to increase.

Like the mean meltwater production and runoff, the mean number of melt days (Fig. 6d) and runoff days (Fig. 6e) increased between the 1990s and the 2000s by 7 and 12 d, respectively. However, instead of continuing to increase in the 2010s, the number of melt and runoff days decreased, suggesting that the net increase in meltwater production and runoff was due to more intense melt events rather than a greater number of melt events.

Figure 6 Annual and decadal mean anomaly of (a) surface mass balance, (b) meltwater production, (c) number of melt days, (d) runoff, (e) number of runoff days, (f) snowfall, and (g) rainfall, from the MAR climate model, averaged over all glacier fronts in the study area. The mean for each anomaly is reported in each panel.

Figure 7a shows that Disko Bay warmed by between 0.21 ^∘C (ICES) and 0.69 ^∘C (ECCO) between the 1990s and 2000s and cooled in 2016–2017 (the final years of the two subsurface temperature datasets). North of Disko Bay (Fig. 7b–h), the decadal subsurface ocean temperatures followed a continuously increasing trend at every location we sampled. The temperature warmed substantially (0.79–1.07 ^∘C) between the 1990s and 2000s at each of these locations, with a lesser magnitude of warming (0.01–0.45 ^∘C) between the 2000s and the 2010s. The temperature trends, especially between the 1990s and the 2000s, are largely imposed by the corrections from Wood et al. (2021) (see Supplement).

Figure 7 Annual and decadal mean deepwater (bottom 60 %) ocean temperature anomaly from ECCO (blue) and ICES (green). Each panel corresponds with an ocean point in Fig. 1. The mean for each anomaly is reported in each panel. The “decadal” means are not taken over the entire decade in the 1990s (for ECCO) and 2010s due to lack of data.

Figure 8 shows that sea-surface temperatures (SST) were relatively steady in the 1970s and 1980s. At some locations the SSTs dipped slightly in the 1990s, followed by a sharp rise (0.57–0.78 ^∘C from ECCO, 0.33–0.81 ^∘C from Hadley-OI) in SSTs in the 2000s at all sites. SSTs in the 2010s were slightly warmer or cooler than in the 2000s, depending on the location and data source, but remained consistently above the pre-2000s temperatures.

Figure 8 Annual and decadal mean sea-surface temperature anomaly from ECCO (blue) and Hadley-OI (orange). Each panel corresponds with an ocean point in Fig. 1. The mean for each anomaly is reported in each panel.

Figure 9 shows that the annual duration of the sea-ice season (defined as when sea-ice concentration exceeds 15 %) has shortened since the 1970s. In particular, from Upernavik Icefjord north to Wolstenholme Bay (Fig. 9c–g), while the decadal mean sea-ice season length was relatively steady from the 1970s through the 1990s, it decreased by 1 to 2 months between the 1990s and the 2000s and remained relatively steady at the new shorter duration in the 2000s and 2010s. This pattern is borne out in the seasonal sea-ice concentration as well at many locations (Supplement Figs. S2–S5).

Figure 9 Annual duration of sea-ice season (when sea-ice concentration is greater than 15 %) from Hadley-OI (orange) and NOAA (purple). Each panel corresponds with an ocean point in Fig. 1.

3.3 Regression analysis

We conducted a suite of multiple linear regressions of key climate variables against terminus retreat (Table 1), using data from 1992–2017, where we have data available for all variables. We ran two cases with all variables included: one with the corrected ECCO deepwater temperatures and one with the original uncorrected ECCO deepwater temperatures. The remaining cases systematically drop one of the variables at a time and then include only one type of variable (ocean temperatures, sea ice, runoff). For each case, we report the sensitivity of terminus retreat to each climate variable with 1σ errors, as well as the R² value of the regression.

While it is difficult to interpret the results because of multiple correlations between the different parameters, it is clear that runoff has the dominant effect. Overall, terminus retreat appears to be most sensitive to runoff, moderately sensitive to ocean temperatures, and least sensitive to sea ice. For each case, the R² value is small, likely due to large interannual variability in the climate variables. When we drop runoff from the regression, the R² value drops by 0.098, suggesting that runoff accounts for ∼ 10 % of the variance in terminus sensitivity. Aside from runoff, dropping any individual variable seems to have little effect. Although dropping ocean temperatures (surface and deepwater) from the regressions has a small effect, keeping only ocean temperatures accounts for about half as much variance (0.094) as keeping only runoff (0.186). The sea-ice parameter sensitivities are not significantly different from zero, and the effect of dropping them is weak.

Table 1 Sensitivities of terminus retreat to different climate variables for a suite of multiple linear regressions, with 1σ error.

4 Discussion

Our measurements (Fig. 5) show a multi-decadal trend of regional glacier retreat with a step change in the terminus retreat rate from around 1996. Several mechanisms have been proposed as drivers of outlet glacier retreat, including ocean-warming-induced terminus ablation and undercutting, mélange rigidity, and enhanced hydrofracture, all of which can affect calving (Straneo et al., 2013). All of these mechanisms are tied to climatic warming in some sense, whether due to rising air or ocean temperatures, and associated changes in meltwater production, runoff, and sea-ice concentration. Our observed step change in terminus retreat rate was approximately coincident with sharp increases in meltwater production (22 %; Fig. 6b), runoff (26 %; Fig. 6c), deepwater temperature (0.21–1.07 ^∘C; Fig. 7), and sea-surface temperature (0.33–0.81 ^∘C; Fig. 8) and a sharp decrease in the duration of the sea-ice season (1–2 months; Fig. 9). Thus, any or all processes related to these anomalies could have contributed to the terminus retreat rates. Multiple linear regressions suggest that runoff is the dominant contributor to terminus retreat, with ocean temperatures accounting for much of the remaining variance. Offshore sea ice has a small and likely negligible effect.

4.1 Observations of regional retreat

The timing of our observed step change in the terminus retreat rate is consistent with previous studies of northwest and/or central-west Greenland that identified accelerated glacier retreat (Carr et al., 2017; Catania et al., 2018; Fahrner et al., 2021; Howat and Eddy, 2011; Wood et al., 2021) and ice speedup and discharge (Joughin et al., 2018b; King et al., 2020) beginning in the middle to late 1990s. Nearly all glaciers in our study region retreated between the onset of this step change and the end of our observations in 2021, in line with other observations of sustained retreat in the 21st century (Bunce et al., 2018; Fahrner et al., 2021; Howat and Eddy, 2011; Murray et al., 2015). Howat and Eddy (2011) found that between 2000 and 2010, nearly 100 % of glaciers in northwest Greenland retreated; although we considered a larger set of glaciers, we found that 85 % of our glaciers retreated and 15 % were stable over the same time period. Between 2010 and 2020, 74 % of glaciers retreated and 22 % were stable. In both decades, 69 % of glaciers retreated and 9 % were stable. Overall, between 2000 and 2020, 86 % of glaciers retreated and 14 % maintained stable terminus positions.

We observed a brief period of regional advance in 2017 and 2018, which was negated by regional retreat in 2019 (Fig. 5a, b). This broad pattern of readvance coincides with the readvance of Jakobshavn Isbræ in 2017 and 2018 (Khazendar et al., 2019); however, while Jakobshavn Isbræ remained in a relatively advanced position in 2019 and 2020 and retreated in 2021 (Fig. 4a, b), the region as a whole began to retreat again in 2019. The regional behavior is coincident with a sustained negative runoff and meltwater production anomaly (and positive surface mass balance anomaly) in 2017 and 2018 with values near the 20th-century mean anomaly. This negative anomaly was followed by the strong positive anomalies in the record in 2019 and 2020 (Fig. 5a). Although we lack sufficient ocean temperature or sea-ice data during this period to assess a relationship between those factors and regional glacier behavior, our regression analysis indicates that runoff anomalies have the strongest effect on terminus position.

4.2 Terminus melting

Ocean subsurface warming and surface warming in Baffin Bay and adjacent fjords have been cited as contributors to regional glacier retreat (Rignot et al., 2012; Slater et al., 2019; Wood et al., 2021) as well as to retreat for individual glaciers (Holland et al., 2008; Khazendar et al., 2019; Motyka et al., 2011; Rignot et al., 2010). Warmer ocean water increases melt at the calving face, and, in conjunction with subglacial discharge plumes, subsurface warming undercuts the terminus (Motyka et al., 2013; Slater et al., 2015, 2017), which could enhance calving (How et al., 2019; Luckman et al., 2015; Morlighem et al., 2019; Rignot et al., 2015).

Wood et al. (2021) found that ocean thermal forcing switched from a stable period to one of rapid warming between 1997 and 1998, consistent with the timing of the sharp step change in glacier retreat rate we observed. Based on the same ocean dataset, we found that ocean subsurface temperatures increased an average of 0.90 ^∘C coincident with our observed step change in glacier retreat rate (Fig. 7) and together with sea-surface temperatures were moderately associated with terminus retreat. Wood et al. (2021) identified deep glaciers sitting in warmer water as those that are retreating the most. Terminus response, however, tends to scale non-linearly with depth, and deeper termini should be more responsive to any change that induces retreat (Schoof, 2007). As evidence of this effect, we note the strong response to seasonal perturbations on three of Greenland's deepest outlet glaciers (Jakobshavn, Kangerlussuaq, and Helheim), which is large compared to the degree of melting at the terminus (Joughin et al., 2020; Kehrl et al., 2017). Thus, irrespective of the forcing that caused the initial perturbation, a greater response would be expected for the deeper glaciers (> 100 m depth), whether the forcing is due to ocean melting or some other process.

Similarly to ocean subsurface temperatures, sea-surface temperatures averaged over all of our subregions also increased (0.58 ^∘C) in concert with accelerated glacier retreat (Fig. 8). Fahrner et al. (2021) found a significant relationship between sea-surface temperatures and terminus change in northwest Greenland, whereas Murray et al. (2015) found no relationship between retreat and sea-surface temperature in this region. Elsewhere in Greenland, warming sea-surface temperatures have been linked with rapid terminus retreat (Howat et al., 2008).

4.3 Mélange rigidity

Models indicate that the absence of a rigid mélange may be more important than ocean-driven melting of the terminus in enhancing glacier retreat (Todd and Christoffersen, 2014). Observations show a correlation between terminus change and sea-ice and mélange conditions (Carr et al., 2017; Moon et al., 2015; Sohn et al., 1998) and that reduced sea-ice and mélange formation could have triggered retreat at several Greenland glaciers (Amundson et al., 2010; Howat et al., 2010; Joughin et al., 2008a; Sohn et al., 1998). The presence of a rigid mélange can exert sufficient force to inhibit calving, facilitating seasonal glacier advance (Amundson et al., 2010; Cassotto et al., 2015; Cook et al., 2021; Reeh et al., 2001; Robel, 2017). We observed a sharp reduction in both the duration of the sea-ice season (Fig. 9) and sea-ice concentrations (Supplement Figs. S2–S5), coincident with a step change in the regional glacier retreat rate in 1996; however, our regression analysis indicated that offshore sea-ice conditions have only a weak connection with terminus retreat. This suggests that the effect of offshore sea ice on terminus retreat is weak, although we note that the data we sampled are an imperfect proxy for near-terminus sea-ice conditions.

Runoff increases and ocean warming are both coincident with terminus retreat. While these factors also contribute to terminus melt, the majority of the ocean heat in the fjords goes into melting icebergs (thus weakening mélange) rather than the terminus face (Moon et al., 2018; Mortensen et al., 2020). Moreover, the runoff-enhanced iceberg melting lags discharge, causing iceberg melt to peak in late summer and early autumn (Moon et al., 2018). At least for Jakobshavn Isbræ, reductions in mélange rigidity appear to correspond to periods of warmer temperatures at depth in the fjord (Joughin et al., 2020). Warmer sea-surface temperatures, although unlikely to contribute to submarine undercutting of the terminus, may also inhibit sea-ice formation and reduce mélange rigidity. Additionally, the observed increase of 26 % in runoff and 11 % in the number of runoff days (Fig. 6c, e), coincident with the observed step change in terminus retreat rates, may have seasonally contributed to more mobile mélange near the terminus. Thus, the combination of increased runoff, sea-surface temperature, and subsurface temperature together suggests a reduction in mélange presence and rigidity that might have increased calving and retreat.

4.4 Enhanced hydrofracture

In addition to reducing mélange rigidity, increased meltwater production and runoff could also increase hydrofracture-induced calving. Crevasses filled with surface melt penetrate deeper than dry crevasses (van der Veen, 1998; Weertman, 1973), and if a crevasse near the terminus penetrates the full thickness of the ice, it can cause calving (Nick et al., 2010; Sohn et al., 1998). We observed a sustained 26 % increase in runoff (Fig. 6c) coincident with accelerated glacier retreat in 1996, which could contribute to increased filling of surface crevasses and subsequent hydrofracture, which would facilitate greater calving and retreat. Although runoff remained high in the 2010s, the duration of the runoff season decreased by 6 d (Fig. 6e) while sustained glacier retreat continued. These results suggest that to the extent that hydrofracture may have contributed to retreat, it is through increased runoff volume rather than seasonal duration.

5 Conclusions

We have built a comprehensive record of annual terminus positions for 87 marine-terminating outlet glaciers in northwest and central-west Greenland from 1972 through 2021. The majority of these glaciers retreated and lost area over the observation period, with retreat accelerating after 1996. We observed a brief regional readvance in 2017 and 2018, which was offset by losses in 2019. Ice-sheet climate data indicate that surface mass balance, meltwater production, and runoff increased substantially between the 1990s and 2000s, coincident with accelerated glacier retreat. Similarly, in most regions sea-surface and deepwater temperatures increased and sea-ice season duration decreased between the 1990s and 2000s. Our results indicate that increased runoff and ocean temperatures correspond with periods of increased terminus retreat. Runoff and ocean warming are expected to increase terminus retreat through a combination of terminus undercutting and reducing mélange rigidity; runoff could also increase calving rates through enhanced hydrofracture. Thus, some combination of these processes is likely responsible for the widespread observed retreat across northwest and central-west Greenland. In order to more conclusively reveal the primary drivers of retreat, future studies would benefit from additional information such as more spatiotemporally dense oceanographic measurements and estimates in fjords (e.g., Wood et al., 2021) and more widespread analysis of fjord mélange conditions (e.g., Joughin et al., 2020). Further research is needed to improve our understanding of the dominant processes contributing to terminus retreat and the resulting increases in ice discharge.