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 two decades has 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 10 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 around 1996. The acceleration of glacier retreat coincides with the timing of sharp shifts in ocean surface temperatures, duration of 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 15 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.

Recent Greenland ice loss has contributed to rates at times approaching 1 mm a -1 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 25 dominated Greenland's mass loss in the past two 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 up-stream thinning (Howat et al., 2008;Joughin et al., 2008b). Because of the relationship between terminus position and calving rate, initial perturbations that increase calving can trigger further glacier retreat. While glacier 30 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).
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 40 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 45 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. 50 In this paper, we analyze glacier change at 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 datasea surface and subsurface temperatures, sea-ice concentration, and ice-sheet surface mass balance, precipitation, melt, and runoffto assess the relative influence of 55 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 correspondence with regional climate trends.

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

Satellite Images
We digitized terminus positions from SAR images acquired by the European Copernicus Program's Sentinel-1A/B 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 65 the surface regardless of clouds or darkness. For the most recent seven years (2014-2015 through 2020-2021), we traced terminus positions in mosaics of Copernicus Sentinel-1A/B images (Joughin et al., 2016a). These mosaics are typically created from images acquired in early February of their respective years and have 25 or 50-m product-dependent resolution.
Radarsat-1 mosaics were used to trace terminus positions for six winters (2000-2001, 2005-2009, 2012-2013) Moon and Joughin, 2008). These mosaics are formed from images collected from October through March and have 70 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 75 instrument's scan-line corrector failure in 2003 as Landsat-5 images were typically unavailable over our study area during this period, although we only kept images which retained a sufficient amount of data to map the terminus. In those images, we digitized across the scan-line corrector gaps when they crossed a glacier terminus ( Supplementary Fig. 1); comparison of temporally close images with different data gaps indicates that this method only marginally increases errors. In all, the scanline corrector gaps affected 348 terminus traces (9.7% of the dataset). Due to winter darkness, we could not select images 80 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 fall 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 85 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-5.

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 Moon and Joughin, 2008 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 least ~1000 m a -1 . We traced each glacier's terminus position once per hydrological year (September 1 through August 31 (Ettema et al., 2009)) in every year for which suitable imagery was 95 available. We used winter or near-winter imagery whenever possible as indicated by Figure 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 100 shadows or an indistinct transition from glacier to mélange) . If a terminus position was ambiguous in one image, it was flagged during tracing, and compared with close-in-time images from the same or other satellite platforms 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 115 to get 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 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 120 periods. Since we mostly focus on decadal scale trends, issues of seasonal sampling should not greatly affect our results.

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 130 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 135 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 get an approximate length change. This measurement should 140 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 in the beginning of the record, we did not interpolate prior to the first observation in the record. The largest temporal gap interpolated is nine 145 years (1975-1976 through 1983-1984), at Sermeq Avannarleq (#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. 150 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 155 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. Although this method always finds one break point in the time series, even if the number of break points is greater or there is no break point at all, we adopted it in order to identify the single-most important year for each glacier when a change in behavior may have occurred.

Climate Data 160
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) 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 165 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) 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 6km x 6km grid over the period 1979-2020. 175 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 1/3° x 1/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° x 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

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

Terminus Positions
We created a dataset of 3606 annual terminus positions for our 87 selected glaciers from 1972-1973 through 2020-2021 (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 195 Thule Air Force 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-1978. Therefore, our glacier analyses start in this year except where noted.

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 2 , equivalent to a cumulative retreat of 287 km. The individual area and length changes are plotted in Fig.  205 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 Isbrae (#3; -130.9 km 2 ), Alison Glacier (#35; -59.4 km 2 ), Kjer Gletsjer (#42; -81.5 km 2 ), and Tracy Gletsjer (#81; -58.7 km 2 ). These glaciers together are responsible for 31.0% of the total area loss. Five glaciers dominate the length change: Jakobshavn Isbrae (-16.9 km), Alison Glacier (-14.7 km), Kjer Gletsjer (-14.7 km), Tracy Gletsjer (-14.0 km), 210 and one unnamed glacier (#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 2 and the mean length change is -2.6 km. Net area and length changes for individual glaciers are detailed in Supplementary Table 3.

Timing of Change 225
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 re-advanced in 230 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 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-1997 and 1997-1998, with additional smaller peaks in break points in the mid-1990s and mid-2000s. The break-point years for each glacier are detailed in Supplementary Table 3. 235

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 245 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), since the late 1970s. The mean surface mass balance was steady in the 1980s and 1990s, but dropped by 0.48 m a -1 between the 1990s and the 2000s, primarily due to substantial increases in meltwater production (+0.42 m a -1 ) and runoff (+0.48 m a -1 ). By contrast, there were only small changes to the mean snowfall (-0.01 m a -1 ; Fig. 6f) and the mean rainfall 250 (+0.02 m a -1 ; 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. 255 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 days and 12 days, 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 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 one to two months between the 1990s and the 2000s, and remained relatively steady at the new shorter duration in the 2000s and 2010s. This 285 pattern is borne out in the seasonal sea-ice concentration as well at many locations . Figure 1.

Regression Analysis 290
We conducted a suite of multiple linear regressions of key climate variables against terminus retreat (Table 1), using data from 1992-2017, when 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, 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σ 295 errors, as well as the R 2 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 2 value is small, likely due to large interannual variability in the climate variables. When we drop runoff from the regression, the R 2 value drops by 0.098, suggesting that runoff accounts 300 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.

Discussion
Our measurements (Fig. 5) show a multi-decadal trend of regional glacier retreat with a step-change in terminus retreat rate 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 310 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 315 contributor to terminus retreat, with ocean temperatures accounting for much of the remaining variance. Offshore sea ice has a small and likely negligible effect.

Observations of Regional Retreat
The timing of our observed step-change in terminus retreat rate is consistent with previous studies of northwest and/or central-west Greenland that identified accelerated glacier retreat Catania et al., 2018;Fahrner et al., 2021;320 Howat and Eddy, 2011;Wood et al., 2021) and ice speedup and discharge (Joughin et al., 2018b;King et al., 2020) beginning in the mid-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 21 st 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 325 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. 69% of glaciers retreated in both decades and 9% were stable in both decades.
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 Isbrae in 2017 and 2018 (Khazendar et al., 330 2019); however, while Jakobshavn Isbrae 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 20 th 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 335 a relationship between those factors and regional glacier behavior, our regression analysis indicates that runoff anomalies have the strongest effect on terminus position.

Terminus Melting
Ocean subsurface 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 for individual glaciers (Holland et al., 2008;340 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 Slater et al., 2015Slater et al., , 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 345 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 stepchange 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 350 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. 355 Similar 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 seasurface 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). 360

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 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 365 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 of sea-ice concentrations , coincident with a step-change in regional glacier retreat rate in 1996; however, our regression analysis indicated that offshore sea-ice conditions have only a weak correspondence with terminus retreat. This suggests that 370 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, 375 causing iceberg melt to peak in late summer and early fall (Moon et al., 2018). At least for Jakobshavn Isbrae, 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 380 to more mobile mélange near the terminus. Thus, the combination of increased runoff, sea-surface temperature, and subsurface temperature together suggest a reduction of mélange presence and rigidity that might have increased calving and retreat.

Enhanced Hydrofracture
In addition to reducing mélange rigidity, increased meltwater production and runoff could also increase hydrofracture-385 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 days (Fig.  390 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.

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 395 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. 400 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.
Further research is needed to improve our understanding of the dominant processes contributing to terminus retreat and the resulting increases in ice discharge. 405

Data Availability
Some of the terminus positions are available on NSIDC as part of the MEaSUREs Annual Greenland Outlet Glacier Terminus Positions from SAR Mosaics dataset (https://nsidc.org/data/NSIDC-0642), and the remainder are being prepared 410 for delivery as part of the same dataset. Bed topography is from BedMachine Greenland V3 (https://nsidc.org/data/idbmg4).

Author Contribution
T.B. and I.J. conceptualized the project. T.B. carried out the terminus data collection, analysis, and visualization. I.J. prepared the SAR data products. T.B. prepared the manuscript, with contributions from I.J.

Competing Interests 425
The authors declare that they have no conflict of interest.