Ocean‐Forcing and Glacier‐Specific Factors Drive Differing Glacier Response Across the 69°N Boundary, East Greenland

The Greenland ice sheet has become a significant contributor to global sea level rise over the last 40–50 years. Approximately half of Greenland's mass loss since 1992 was due to increased ice discharge from marine‐terminating outlet glaciers. Here, we present high temporal resolution (∼monthly) time series of ice frontal positions for 24 marine‐terminating outlet glaciers along Greenland's east coast between 2013 and 2020. The glaciers are located north and south of 69°N, which has previously been identified as a potential divide in glacier response to climate forcing. Frontal positions are compared to ice velocity, atmospheric and oceanic data, allowing investigation of change at both seasonal and interannual timescales. Our results reveal 19 of 24 study glaciers underwent net retreat. We find marked differences in interannual patterns of frontal position between glaciers located north and south of 69°N. South of 69°N, glaciers underwent multiyear retreat initiated in 2016, which we attribute to over‐winter calving, resulting from warmer ocean waters and repeated ice‐mélange break‐up. North of 69°N glaciers show either limited or gradual changes in frontal positions and retreat patterns are characterized by more year‐on‐year variability between glaciers. Although similar atmospheric conditions occur across both regions, glaciers north of 69°N experience minimal change in ocean conditions, and are strongly influenced by glacier‐specific factors. Our results show that 69°N continues to represent a boundary between different glacier responses and climate forcing, which is likely to persist under current conditions.

The external controls outlined above are modulated by glacier-specific factors, particularly the bed and fjord geometry, which can strongly enhance/suppress glacier response to forcing (e.g., Carr et al., 2015;Moon et al., 2012;Warren & Glasser, 1992). It is important to identify the factors driving enhanced ice discharge and how this varies between glaciers and regions to quantify near-future Greenland Ice Sheet losses and its contribution to sea level rise in response to climate change (e.g., Cowton et al., 2018;Fahrner et al., 2021;King et al., 2018;McFadden et al., 2011;Wood et al., 2021).
Much of the Greenland Ice Sheet underwent widespread and rapid outlet glacier retreat in the 2000 s (e.g., Goliber et al., 2022;King et al., 2018), but glaciers along Greenland's east coast showed markedly different behavior during this time (e.g., Carr et al., 2017;Seale et al., 2011;Walsh et al., 2012). Specifically, Seale et al. (2011) examined patterns of frontal retreat for 32 glaciers along the east Greenland coast between 2000 and 2009 and found a distinct contrast between glaciers located south of 69°N compared to those to the north. The southern glaciers underwent widespread, synchronous retreat between 2001 and 2005, accompanied by acceleration and thinning (e.g., Howat et al., 2008;Murray et al., 2015;Rignot & Kanagaratnam, 2006), followed by a period of slowdown and stabilization through 2009 (e.g., Murray et al., 2015;Seale et al., 2011). In contrast, the glaciers north of 69°N showed slow and limited terminus retreat or dynamic change during the same time period (e.g., Murray et al., 2015;Seale et al., 2011;Walsh et al., 2012). The transition between these two contrasting retreat regimes occurred at ∼69°N, with this latitude marking the northern extent of warm subtropical waters transported within the Irminger Current Murray et al., 2010;Seale et al., 2011). Therefore, it was suggested that: (a) variability in coastal heat transport was the dominant control on the difference in glacier dynamics observed in eastern Greenland between 2000 and 2009 (Seale et al., 2011); and (b) glaciers along Greenland's east coast showed a positive correlation between ocean temperatures and terminus retreat (Cowton et al., 2018;Murray et al., 2010). The period of slowdown and stabilization of glaciers in the southeast has persisted until at least 2016, likely resulting from regional frontal positions alternating between net annual retreat and advance since 2010 (e.g., Bunce et al., 2018;King et al., 2018).

10.1029/2022JF006857 3 of 24
A key exception to this pattern was Kangerlussuaq, one of east Greenland's largest glaciers, which entered a new phase of rapid retreat and acceleration from 2016 (Brough et al., 2019). This was triggered by the presence of warm ocean water proximal to Kangerlussuaq, which, along with warm air temperatures, caused multiple mélange break-ups, thus allowing calving to persist through two consecutive winters (Bevan et al., 2019). Given the observed substantive changes in mélange and ocean conditions at Kangerlussuaq and the new phase of accelerated retreat of the glacier, it is timely to reassess whether the dichotomy in glacier behavior around 69°N persists under the recent phase of warmer atmospheric and ocean conditions. In particular, it is vital to assess how far north the influence of these warmer conditions extends (i.e., whether the 69°N boundary has migrated northwards), and the resultant impact on glacier retreat, ice dynamics and ice surface thinning, with potential consequences for global sea level change.
Here, we extend and enhance previous work (i.e., Seale et al., 2011) by assessing, either side of the 69°N dichotomy, whether: (a) glacier response to climate forcing persists under warmer atmospheric and ocean conditions; (b) the dynamic response of east Greenland glaciers differs; (c) patterns of seasonal behavior, and the factors controlling it vary, with particular focus on cold season forcing factors; and (d) glacier-specific factors impact individual glacier behavior and/or contribute to the dichotomous response across the boundary.
To achieve this, we present high temporal resolution (∼monthly) time series of ice frontal positions (excluding winter) and velocity for 24 tidewater glaciers along Greenland's east coast between 2013 and 2020. The glaciers studied extend between 67 and 75°N ( Figure 1), but exclude the majority of tidewater glaciers along the Blosseville Kyst due to a preponderance of surge type glaciers (Jiskoot et al., 2003;Walsh et al., 2012), and further north (e.g., Zachariae Isstrøm and Nioghalvfjerdsfjorden), which have already been investigated (e.g., An et al., 2021). Specifically, we assess how spatial and temporal patterns of change in glacier dynamics correspond to atmospheric and oceanographic forcing. We subsequently assess whether significant differences in glacier-specific forcing factors exist across the 69°N boundary, which could potentially account for observed differences in glacier retreat rates and dynamic response to forcing.

Glacier Frontal Positions
Terminus positions were manually digitized from all available Level 1T pansharpened (15 m) Landsat 8 Operational Land Imager (OLI) satellite imagery between 2013 and 2020 using the Google Earth Digitization Tool (GEEDiT; Lea, 2018). Each Landsat 8 image was visualized within a web-browser and the glacier termini were manually digitized where an ice front was visible. Where multiple images were available for the same day as a result of overlap in imaging tracks, we measured the ice front using the first image acquired unless there was any discernible change, in which case both were included. The presence of cloud and/or year-round mélange also precluded the mapping of part or all of the glacier terminus in some images, with 40%-80% of all images available for each glacier being mapped. We only include in our analysis fully mapped termini. Specifically, we obtained an average of 239 (σ = 94) terminus traces for each of the 24 glaciers in our region of interest (ROI; Figure S1 in Supporting Information S1). On average, 30 (σ = 7) terminus traces were obtained per year, with an average sampling frequency of 13 days (σ = 6) for each glacier. Owing to poor solar illumination, no observations were made during late-November through mid-February. In general, glaciers in the northern part of our ROI (above 69°N) had between two and three times the number of observations than those in the southern part of our ROI (below 69°N), due in part to the convergence of the satellite flight path. Following the approach of Brough et al. (2019), we attribute uncertainty in ice front position mapping to error in both geolocational accuracy of imagery and precision in manual digitization of the ice fronts and was ±4.4 m, which is less than the pixel resolution (15 m) of our imagery.
Changes in frontal position were assessed using the curvilinear box method in the Margin change Quantification Tool (MaQiT; Lea, 2018). This method used a reference box of fixed width and upstream extent that intersects with contiguously mapped glacier termini, and mean retreat was calculated by dividing the change in reference box area by its width. This method therefore captured spatially asymmetric retreat and advance of a calving margin . Mapped glacier termini were subsequently exported from GEEDiT in vector format as GeoJSON files and converted to ESRI Shapefiles using the MaQiT (Lea, 2018). To factor out the absolute magnitude of individual glacier terminus change, changes in frontal position were also normalized using min-max feature scaling (Equation 1). All glaciers are therefore normalized between 0 (most retreated terminus position) and 1 (most advance terminus position), thus facilitating a more direct comparison in the pattern and timing of frontal position fluctuation between glaciers (e.g., Fahrner et al., 2021

Outlet Glacier Velocities
Ice surface velocities were compiled from ice-sheet-wide products from the Making Earth Science Data Records for Use in Research Environments (MEaSUREs) program and the Programme for the Monitoring of the Greenland Figure 1. Location of studied east Greenland marine-terminating glaciers (colored by latitude). Official glacier names are used where available (Bjørk et al., 2015). Bathymetry overlay indicating topography below sea level (Morlighem et al., 2017a(Morlighem et al., , 2017b Ice Sheet (PROMICE) ( Table S1 in Supporting Information S1). Monthly velocity mosaics were available from December 2014 onwards for the MEaSUREs product and provided a total of 60 velocity maps (Joughin, 2020;Joughin et al., 2010Joughin et al., , 2018. These ice velocity products were derived from a combination of Synthetic Aperture Radar data measured by TerraSAR-X (TSX), TanDEM-X (TDX), Sentinel-1A and Sentinel-1B satellites, and optical satellite imagery from Landsat 8, and have a nominal spatial resolution of 200 m. The PROMICE data set was produced from compositing all 6 and 12-day Sentinel-1A and Sentinel-1B image pairs over a 24-day period, and had a nominal spatial resolution of 500 m . A total of 164 velocity maps were available from September 2016 onwards, providing a sampling frequency between 6 and 12 days.
To increase the temporal coverage of the ice-sheet-wide velocity datasets, particularly at the start of our study period, we also used velocity fields provided by the Global Land Ice Velocity Extraction from Landsat 8 (GoLIVE) data set Scambos et al., 2016) and selected glacier site velocity maps from Interferometric SAR (InSAR) data set (Joughin et al., 2010. MEaSUREs and GoLIVE have been used successfully in previous work to create a consistent velocity time series for Greenland outlet glaciers (e.g., Catania et al., 2018). GoLIVE velocity maps were generated from pairs of Landsat 8 panchromatic images acquired from May 2013 onward, and had a nominal spatial resolution of 300 m. To maximize temporal coverage, we used all 33 Landsat-8 scenes which covered our study area and used all velocity fields generated from scenes separated by 16 days. To improve accuracy in the vector displacements, we used the masked velocity product and accounted for significant cloud cover and undetected offset errors in the data set  by removing scenes where: (a) the maximum velocity was ≥40 m d −1 (e.g., faster flow than observed at Sermeq Kujalleq [Jakobshavn Isbrae]); and/or (b) the average velocity of all valid pixels within the scene corresponded to a value that was greater than 10% of the maximum recorded pixel velocity in the scene. After this filtering process, we were left with 1,386 velocity maps.
The InSAR-derived data set was produced from 11-to 33-day image pairs measured by the TSX and TDX satellites, and had a nominal spatial resolution of 100 m . A total of 253 unique velocity maps were available from January 2013 onwards.
Following previous studies, velocities were sampled at a central location ∼5 km up the glacier from the most retreated frontal position of each glacier (e.g., McFadden et al., 2011;Walsh et al., 2012). Following data extraction, data time series were postprocessed to remove poor quality or unrealistically high/low velocity measurements where: (a) velocity errors > velocity magnitude; and/or (b) velocity magnitude was either larger than twice the mean of the time series or smaller than the mean of the time series divided by two. This postprocessing removed an average of 11 (σ = 3) measurements per glacier (3 ± 1%). Velocity time series were also normalized using the method outlined in Section 2.1 and Equation 1. We obtained an average of 361 (σ = 88) velocity measurements for each of the 24 glaciers in our ROI ( Figure S2 in Supporting Information S1), with an average of 45 (σ = 25) velocity measurements obtained per year. Where provided, relative velocity errors were calculated using the data set error values for each velocity field, and resulted in a mean error of ±5.0% (σ = 4.5%) for our study glaciers (Table S2 in Supporting Information S1).

Surface Air Temperature
Records of daily surface air temperature (SAT) were determined for Aputiteeq (67°47′N, 32°18′W), Mitt. Nerlerit Inaat (70°45′N, 22°39′W) and Daneborg (74°18′N, 20°13′W) meteorological stations ( Figure 2). Data were provided by the Danish Meteorological Institute at hourly resolution (Cappelen, 2021), and were only used to calculate daily averages where: (a) no more than six consecutive records were missing in a day; and/or (b) no more than nine records in total were missing in a day (e.g., Cappelen, 2011;Carr et al., 2013). Meteorological stations are located between 14 and 277 km from each glacier (Figure 2), length scales over which surface temperatures are generally consistent in Greenland (Box et al., 2009;Schild & Hamilton, 2013). For regional analysis, monthly SAT was calculated if there were more than 22 daily average measurements per month (Carr et al., 2013). Subsequently, anomalies in monthly mean temperatures were calculated relative to the 2002-2020 baseline period, 2002 being the earliest common date for all three meteorological stations.

Meltwater Runoff
To compare glacier behavior and meltwater availability, runoff for each study glacier catchment was determined using modeled output, with a 20-km horizontal resolution, from the Modèle Atmosphérique Régionale (MAR) 10.1029/2022JF006857 6 of 24 v3.11 (Tedesco & Fettweis, 2020) and individual glacier catchment basins   (Figure 2). For each glacier, total daily runoff was calculated by summing the values at each grid point across the glacier basin. For regional analysis, we also calculated cumulative annual runoff for each basin. Where glaciers are part of compound glacier catchments (Unnamed A and Unnamed B, Polaric B and Apuliliip Apusiia, Roligre Brae and Eielson), catchment wide runoff was assigned to each of these glaciers as no definitive split could be determined, but were only included once for regional analysis.

Ocean Temperature
Records of monthly ocean potential temperatures were extracted from the Level 4 TOPAZ4 Arctic Ocean Physics Reanalysis product (Xie et al., 2017), supplied by the Copernicus Marine Environment Monitoring Service (https://marine.copernicus.eu). The reanalysis product assimilates a combination of satellite and in situ observations using the Hybrid Coordinate Ocean Model (HYCOM) with a coupled sea ice model and has a spatial resolution of 12.5 km (Xie et al., 2017). As data do not extend into fjord waters, they were sampled from 50 by 50 km boxes that were located on the continental shelf ( Figure 2c). As the data provide information on ocean temperatures on the continental shelf and do not account for the complex processes within the glacier fjords or at the calving front, they are instead used to give a broad-scale indication of temperature change with depth in the water column on the continental shelf (Carr et al., 2017). Ocean potential temperatures were sampled for depths of 5 and 200 m in order to capture changes within the different water masses, specifically the 5 m depth data were chosen to represent the surface mixed layer, which is primarily composed of Polar Surface Waters (PSW) and is thus referred to as PSW henceforth, and 200 m depth data were chosen to include the upper layers of AW (e.g., Straneo et al., 2012;Sutherland et al., 2014). The PSW can impact glacier dynamics via its impact on the ice mélange (e.g., Bevan et al., 2019) and/or undercutting of the calving front, leading to calving at the water line (Benn et al., 2007), while warmer ocean temperatures at 200 m depth can cause ice loss via enhanced subaqueous melt across the calving face and/or via their interaction with subglacial plumes (e.g., Straneo et al., 2012;Straneo & Heimbach, 2013). Ocean potential temperature root-mean-squared-errors range between 0.32 and 1.04°C at depths of 0, 100, 300, 800, and 2,000 m (Xie & Bertino, 2020). For regional analysis, anomalies in monthly mean temperatures were subsequently calculated relative to the 1991-2019 baseline period, 1991 being the earliest available date in the TOPAZ record and 2019 being the last available date. To provide confidence in the use of modeled TOPAZ data, these are compared to observed temperature profile data obtained from Airborne eXpendable Conductivity Temperature Depth (AXCTD) casts from NASA's Oceans Melting Greenland mission (Fenty et al., 2016;OMG, 2019). Though the modeled and observed data represent different spatial and temporal scales, a significant correlation (r = 0.61, p = <0.01) between these two independent datasets is observed (Figures S3-S11 in Supporting Information S1). For the specific depths used in this study both show a significant correlation ( While the differences in spatial and temporal scales preclude a more direct comparison of the datasets, these results provide confidence that the use of modeled TOPAZ data is appropriate for East Greenland.

Mélange and Sea Ice Conditions
Given the difficulties in distinguishing between sea ice and ice mélange purely from remote sensing measurements, we followed previous studies (e.g., Fried et al., 2018;Moon et al., 2015) and referred to all ice immediately seaward of the terminus as mélange, even though in some scenarios this mélange may be more akin to seasonal sea ice. Following Moon et al. (2015), we characterized the timing of mélange coverage by tracking mélange conditions in each fjord throughout our observational period using our Landsat imagery. For each image, the near-terminus region of each glacier (typically within 5 km) was classified as: (a) (likely) rigid when there was (near-)complete mélange coverage and little relative motion of the mélange between near-time images; (b) potentially rigid (mixed) when there was evidence of fractures in the mélange and/or some mélange motion between near-time images; or (c) unlikely rigid (open) when there was open water adjacent to the calving front and/or extensive mélange motion between near-time images (e.g., Davison et al., 2020;Moon et al., 2015).
As well as these individual measurements, periods of consistent mélange conditions were classified where the observed conditions persisted for at least 14 days and there were three or more observations during the time period (Moon et al., 2015). It was assumed that if two observations of the same mélange condition were made, the condition was maintained between observations. This method allows for the tracking of mélange conditions and corresponding terminus change, irrespective of the time of year. From our classification of consistent mélange conditions, we compiled a total 412 mélange/terminus observation intervals, with an average of 17 (σ = 3) intervals per glacier. The data set comprises of 167 rigid periods, 72 mixed periods and 173 open periods. Following Moon et al. (2015), we remove observation windows that coincide with no noticeable terminus change, set to twice our uncertainty in ice front position (<10 m, i.e., >2σ), which removes 10%, 17%, and 4% of our observations for rigid, mixed and open conditions, respectively.
To better constrain broad regional sea ice patterns, this visual record was augmented with data from sea ice charts provided by the US National Ice Center (NIC) (https://usicecenter.gov/Products/ArcticData). The sea ice charts are compiled from a range of remotely sensed and directly measured data sources, have a spatial resolution of up to 50 m, and are provided at subweekly to biweekly temporal resolution (e.g., Carr et al., 2013). The estimated accuracy of sea ice concentrations is ±10% (Partington et al., 2003). Sea ice fraction values were sampled as close to the glacier terminus as possible from a polygon mirroring the width of the glacier terminus and, depending on fjord length, up to 5 km perpendicular to it (e.g., Carr et al., 2013). For regional analysis, data were resampled to a monthly resolution by taking the average of all measurements within the monthly time period (two to four measurements). These data were subsequently used to calculate the number of ice-free months, defined here as months with sea-ice concentrations of <10% (e.g., Carr et al., 2017).

Glacier-Specific Factors
We analyzed differences in glacier-specific factors across the 69°N boundary, specifically catchment area; width; ice thickness at the grounding line; mean ice velocities over the 8 years study period; mean ice flux over the 8 years study period; ice surface slope; bed slope; access of AW to the glacier terminus; and the change in fjord width between the least and most retreated terminus position. Catchment area was calculated using the basin outlines of Mouginot and Rignot (2019). Glacier width was calculated from the average of all terminus positions for each glacier, except for Kangerlussuaq, where we used the average width across the fjord, between the least and most retreated frontal positions during the study period, as terminus traces did not reach the fjord walls. We calculated the mean ice thickness, bed depth, ice surface slope and bed slope within the area occupied by the glacier termini during the study period (i.e., March 2013-October 2020). Ice thickness and bed depth were calculated from the thickness and basal topography layers of BedMachine v3 (Morlighem et al., 2017a(Morlighem et al., , 2017b, respectively, while the ice surface slope and bed slope were calculated from BedMachine's surface and bed layers, in degrees. Mean velocities were averaged for the 8-year study period using the monthly data described in Section 2.2. Mean ice flux was calculated by multiplying glacier width, grounding line thickness and mean ice velocity. Within the envelope of frontal positions occupied by each glacier during the study period, we identified whether a reverse (i.e., inland) sloping bed was present using BedMachine v3 data, and the inland change in fjord width, using the categorization of Carr et al. (2014) and Bunce et al. (2018). Finally, we categorized the study glaciers based on fjord geometry and water properties, following Wood et al. (2021): where DW indicates the presence of AW; CR indicates shallow ridges; SC indicates a shallow cold fjord with polar water; and NC indicates noncatergorized, due to lack of data .

Terminus Change
Between 2013 and 2020, 19 of 24 study glaciers underwent retreat (Figure 3a). The magnitude of retreat varied between glaciers, ranging from 65 m at Courtauld to more than 3,500 m at Unartit and Kangerlussuaq. Five glaciers underwent advance during the study period, with the largest advance of ∼100 m occurring at Vestfjord (Table 1).
We observed marked differences in interannual patterns of frontal position changes between glaciers located to the north and the south of 69°N (Figure 4a). Glaciers to the south of 69°N showed a coherent regional pattern  Table 1). In 2020, the 11 glaciers south of 69°N showed greater variability than in previous years, with 63% of glaciers advancing and 36% retreating, but regionally net advanced occurred ( Figure 5; Table 1). As a result of these interannual patterns of behavior between 2013 and 2020, the 11 glaciers south of 69°N retreated on average ∼1100 m, with an average retreat rate of ∼160 m a −1 (Table 1). Two southern glaciers, Apuliliip Apusiia and Christian IV, showed minor advanced (Table 1). In contrast, glaciers above 69°N exhibited much lower average retreat rates (30 m a −1 ) for the study period than those in the south (∼160 m a −1 ), but showed higher variability between individual glaciers, meaning there was no clear, coherent pattern of behavior for glaciers north of 69°N ( Figure 5; Table 1). Furthermore, in some years, such as 2018, the overall retreat signal was dominated by a single glacier (Daugaard-Jensen), which masked small advances of multiple glaciers across the region (Figures 4c and 5c; Table 1). Above 69°N, retreat rates were highest in 2019, where 11 out of 13 glaciers showed net retreat (Nunatakgletscher changed by less than 2 σ of observational uncertainty), with a regional average retreat rate of ∼70 m a −1 .
Intra-annual or seasonal variations were evident on all glaciers (Figure 4; Figures S12-S35 in Supporting Information S1), although there were differences in the timings of these variations between glaciers and years. For glaciers above 69°N, the onset of retreat tended to occur between May and July, and the onset of advance tended to occur between August and October (Figure 4 and Figure S28 in Supporting Information S1). For glaciers  Figures S38-S41 in Supporting Information S1 for further detail on regional variability and normalization).
below 69°N retreat onset occurred over a broader period of April to August, and was often coincided with the weakening or removal of the mélange (Figure 4 and Figure S27 in Supporting Information S1). The onset of advance occasionally started in August for a small number of glaciers, such as Unnamed A and Courtauld, but tended to occur between September and November for glaciers below 69°N, although we may miss the retreat/ advance transition if it occurred during the winter season where no margins were detected due to poor solar illumination. The magnitude of seasonal frontal position variations also varied by glacier, with seasonal cycles ranging from 80 m to over 3,000 m with a mean of 477 m for all glaciers (Figure 3c; Table 1). There was no significant difference at the 95% level in the means of the seasonal cycles between glaciers north and south of 69°N. Rather, the magnitude of seasonal advance and retreat cycles showed significant correlation at the 99% level to glacier velocity (Figure 6), indicating that seasonal cycles are broadly similar between glaciers and regions. Notwithstanding the broad seasonal patterns, there was a noticeable change in seasonal behavior between summer 2016 and spring 2017 for several glaciers below 69°N (Figures 4b and 4c). During this period, seven of the 11 glaciers showed over winter retreat and/or a reduced readvance through the early part of 2017. A similar pattern was also observed in 2017/2018 in four of the 11 glaciers below 69°N. These changes were particularly marked at Kangerlussuaq, which retreated by ∼7.5 km between June 2016 and June 2018 and accelerated by 15% during this period (Figure 4 and Figure S21 in Supporting Information S1). Glaciers above 69°N showed no noticeable change in their seasonal behavior during this period, where winter through spring advance occurred (Figures 4b and 4c).

Ice Velocity Change
The difference in patterns of terminus change above and below 69°N is reflected in the regional velocity records (Figures 4e and 4f). Ice velocities on glaciers south of 69°N generally decreased between 2013 and 2015, coincident with limited terminus position changes (Figures 4b and 4c). Ice velocities then increased during late 2016 and through to 2020, as glaciers began to retreat more rapidly (Figures 4b and 4c). Summer (July) velocities in 2020 were between ∼5 and ∼60% higher than in 2016 for seven out of the 11 glaciers. The remaining four glaciers located south of 69°N showed either negligible change (Apuliliip Apusiia, Styrtegletsjer) or a velocity decrease (Christian IV and Courtauld) for the same period. For glaciers above 69°N, there was a regional increase in velocity between 2013 and mid 2019, coinciding with the general retreat of these glaciers (Figures 4b and 4c). Starting in late 2019, surface velocities decreased north of 69°N and in 2020, velocities are comparable to those in 2014/2015 (Figures 4e and 4f).

Air Temperatures and Runoff
Here, we evaluate whether the interannual contrasts in glacier dynamics on either side of 69°N are influenced by trends in environmental forcing factors. Air temperatures and their anomalies are statistically correlated at the 99% level for the three meteorological stations across our study period (Figures 7b and 7c), indicating that air temperature changes are regional. A period of anomalously warm air temperature from March 2016 through to December 2016 occurred at all three stations, with anomalies relative to the 2002-2020 mean (2002 being the earliest common date for all three meteorological stations) peaking at +2.9°C in October at Aputiteeq (furthest south) and +7.1°C in November at Daneborg (furthest north). Air temperatures were also anomalously warm at Aputiteeq in February (+4. Relative cumulative runoff patterns are similar between glaciers north and south of 69°N, but vary in magnitude (Figures 7d and 7e). For glaciers below 69°N runoff was high in 2014 (∼16 Gt), 2016 (∼18 Gt), and 2019 (∼19 Gt) and reduced (∼10-13 Gt) in the remaining years. Runoff was lowest in 2018 (∼10 Gt). For glaciers above 69°N, cumulative annual runoff increased from ∼10-11 Gt in 2013 through 2015 to ∼16 Gt in 2016. Cumulative runoff reduced to ∼6-8 Gt in 2017 and 2018 before peaking at 18 Gt in 2019. Cumulative runoff was again reduced in 2020 at 10 Gt (Figures 7d and 7e).

Ocean Temperatures
Near-surface waters on the continental shelf were exceptionally warm in 2016 (Figures S47 in Supporting Information S1). Modeled ocean potential temperatures at a depth of 5 m for waters below 69°N were above average between May 2016 and February 2017: they were over +3.6°C warmer than the 1991-2019 mean between July and August 2016 (Figures 8b and 8c). Similar anomalously warm periods of ocean temperature at 5 m depth occurred between June 2018 and December 2018 and between April 2019 and November 2019 for glaciers below 69°N, with temperatures peaking at +2.7°C and +2.8°C warmer than average respectively. A less intense period of anomalously warm water at 5 m depth also occurred between July 2014 and December 2014 for glaciers below 69°N. For waters at 5 m depth above 69°N, there was also a period of above average temperature between June 2016 and December 2016, peaking at +3.9°C above average temperature in August 2016 (Figures 8b and 8c). Shorter and less intense periods of anomalously warm water at 5 m depth also occurred in spring through summer in all years apart from 2015 for glaciers below 69°N.   Figure S42 in Supporting Information S1 for individual glacier records for (d)).
Information S1). A particularly extensive period of anomalously warm water has occurred between May 2018 through to the end of our record in November 2019, peaking at +3.1°C above average in January 2019 (Figure 8c).

Terminus Correspondence to Mélange and Sea Ice
We first assessed whether our monthly sea ice data showed a difference in the seasonality of regional sea ice patterns above and below 69°N (Figure 8d). For glaciers above 69°N, sea ice usually reaches 100% coverage between November and January, before starting to decline between May and July. Sea ice coverage generally reaches a minimum in September and October when ice free (<10% ice coverage) periods occur in six out of the 8 years (2013, 2014, 2016, 2018, 2019, and 2020). For glaciers below 69°N sea ice generally reaches 100% coverage between February and April, before starting to decline between May and June. Sea ice coverage generally reaches a minimum in August and September, but regional ice-free conditions only occur for two (2016 and 2018) out of the 8 years (Figure 8e). Between 2016 and 2020 there is a noticeable reduction in regional sea ice below 69°N, where winter maximum ice coverage occurs for shorter periods and summer minima has less sea ice coverage (Figures 8d and 8e; Figure S45 in Supporting Information S1).
Following Moon et al. (2015), we also assessed the correspondence between terminus position change and the ice mélange by classifying intervals of consistent (≥14 days) rigid, mixed and open mélange between our   Figure S46 in Supporting Information S1 for break-down of frontal change data for individual glaciers). 95 days with mean retreat of 271 m (Figure 9). Regionally, glaciers showed the same type of terminus response to rigid and open conditions and only vary in their magnitude of change, with glaciers below 69°N advancing/ retreating further than those above 69 o N (Figure 9). However, some differences occurred in terminus behavior under mixed conditions for glaciers located north versus south of 69°N. For glaciers below 69°N during mixed conditions, 16% of observation intervals coincided with glacier retreat and 84% of observation intervals with glacier advance, with the terminus advancing 289 m on average. For glaciers above 69°N, 45% of mixed condition intervals coincided with periods of retreat and 55% of observation intervals with periods of advance, with the terminus retreating 9 m on average. Thus, glaciers located south of 69°N were far more likely to advance under mixed mélange conditions than those to the north and the magnitude of their terminus position change was far higher (Figure 9).

Glacier-Specific Factors
We observed no significant differences between our northern and southern study glaciers for catchment area, terminus width, ice thickness, bed depth, ice surface, or bed slope (Table 2). However, we did observe a significant difference in 8-year mean ice velocities, and hence ice flux, with glaciers south of 69°N flowing significantly faster than those to the north (Table 2). A greater proportion of glaciers south of 69°N were located on reverse bed slopes (5 out of 11) versus north of 69°N (3 out of 13). Overall, reverse bed slopes were associated with greater glacier retreat: of the five largest glacier retreats, four glaciers retreated down reverse slopes (Table 2). However, this was not universal, as glaciers on reverse slopes also underwent moderate to limited retreat and Vestfjord had the greatest advance, despite being on a reverse slope ( Table 2). The fjord bathymetry/water properties categorization was very similar either side of 69°N and the glaciers that underwent the greatest retreat were categorized as either "non-categorized" or terminating in "deep warm water", while those that advanced were "calving ridges," "shallow cold" or "noncategorized" (Table 2). Inland changes in fjord width were very similar on either side of 69°N and showed no clear relationship to glacier retreat magnitude, with almost all glaciers retreating between parallel fjord walls (Table 2).

Contrasting Glacier Dynamics North and South of 69°N
Our high temporal resolution data enable us to quantify the seasonal cycle in glacier frontal position for glaciers located north and south of 69°N. Our findings demonstrate that all east Greenland glaciers exhibit seasonal advance and retreat cycles proportional to glacier velocity ( Figure 6), but there is a distinct difference between the interannual trends in frontal position north and south of 69°N throughout our study period (Figure 4). We find that glaciers above 69°N showed either limited or gradual changes in frontal position and ice velocity over time (e.g., Figures 4b and 4c; Table 1). This contrasts with glaciers south of 69°N, where widespread retreat and glacier acceleration was observed between 2016 and 2019, following a period of relative stability in mean annual frontal position between 2013 and 2015 ( Figure 4; Table 1). The observed contrasting behavior in glacier dynamics between glaciers above and below 69°N demonstrates a persistent difference in regional glacier behavior and dynamics from at least 2000 along Greenland's east coast (e.g., Howat et al., 2008;Murray et al., 2010;Murray et al., 2015;Rignot & Kanagaratnam, 2006;Seale et al., 2011). Furthermore, our results indicate that the recent observed retreat at Kangerlussuaq (e.g., Bevan et al., 2019;Brough et al., 2019) was not isolated, and was part of a wider regional phase of retreat suggesting that these glaciers may have been responding to regional changes in environmental forcing(s). Similar observations have also been recently observed at peripheral glaciers along Greenland's south east coast (Liu et al. (2022)).

Climatic and Oceanic Controls on Glacier Changes
Previous work (Seale et al., 2011) identified ocean forcing as a primary control on outlet glacier retreat in east Greenland and the underlying driver of the boundary in glacier behavior either side of 69°N. Our high temporal resolution data enable us to assess the exact timings of the various forcing factors and the resultant changes in glacier dynamics, and hence to determine how observed interannual ocean work drives glacier change in east Greenland. Specifically, our work highlights the role of winter-time conditions in driving net glacier retreat and acceleration: the onset of retreat south of 69°N was coincident with a change in the seasonal pattern of glacier advance and retreat, with seven of the 11 glaciers showing over-winter retreat and/or a reduced terminus readvance through the early part of 2017 (Figure 4c; Figure S36 in Supporting Information S1). Furthermore, these glaciers showed over-winter increases in velocity (Figure 4f; Figure S36 in Supporting Information S1).
Changes in terminus position and velocity south of 69°N were coincident with a period of marked climatic and oceanic changes that persisted through winter, including anomalously warm atmospheric temperatures and modeled near-surface PSW (5 m depth) and AW (200 m depth) oceanic temperatures; high runoff; and a marked decline in sea ice/mélange (Figure 8). These conditions differed from previous periods when atmospheric and/or modeled near-surface (5 m depth) ocean temperatures were similarly high (e.g., 2013/2014) because warm conditions continued throughout the winter months (December/January/February; Figure 8). The continuously warm atmospheric and ocean temperatures appear to have delayed and/or disrupted the formation of mélange, which was typically associated with the seasonal onset of glacier advance (Figure 9 and Figure S36 in Supporting Information S1). Thus, delayed mélange formation may have allowed retreat to continue throughout winter, and into early spring (Figure 4c). This corresponds to the proposed mechanisms  . Glacier width was calculated from the average of all terminus positions for each glacier apart from Kangerlussuaq (KAN), where an average across-fjord width between the most advanced and retreated terminus position is used, as terminus traces do not extend the full fjord width. Mean ice thickness, bed depth, ice surface slope and bed slope were calculated within the area occupied by our terminus positions during the study period and ice and bed geometry data were obtained from BedMachine v3 (Morlighem et al., 2017a(Morlighem et al., , 2017b. Mean 8-year velocity was calculated using monthly interpolated velocity data ( Figure 4) and mean 8-year ice flux was calculated by multiplying glacier width, grounding line thickness and mean ice velocity. Following Carr et al. (2014), we categorized the bed slope and fjord width change experience by each glacier during the study period. Glaciers were also categorized by fjord geometry and water properties, following Wood et al. (2021): where DW indicates the presence of AW; CR indicates shallow ridges; SC indicates a shallow cold fjord with polar water and; NC indicates noncatergorized.The final row shows the P-Value for the Wilcoxon test, which was used to identify significant differences between glaciers above and below 69 o N. A P-value of ≤0.05 indicates a significant difference and are highlighted bold.
driving retreat at Kangerlussuaq over winters 2016/17 and 2017/18 (Bevan et al., 2019) and suggests that delayed and/or weakened mélange formation in response to warmer near-surface ocean temperatures (5 m) is a persistent control on ice loss for the glaciers located south of 69°N. Furthermore, warmer near-surface ocean temperatures and reduced ice mélange occurred after high air temperatures and runoff in the summer of 2016. These may have preconditioned the glaciers to retreat through surface melt driven thinning and/or increased terminus melt due to increased subglacial discharge and enhanced plume flow (e.g., Christoffersen et al., 2012;Jenkins, 2011;Motyka et al., 2003)-the latter of which may have also contributed to the weakening/melting of mélange. Therefore, a combination of forcings may have acted together to decrease terminus stability and promote above average seasonal retreat in winter 2016/17, in turn leading to net terminus retreat. Our seasonal data highlight the importance of winter conditions in driving glacier retreat, which has received much less scientific attention than summer forcing, and should be investigated elsewhere on the Greenland Ice Sheet.
In contrast to the reduced advance in winter 2016 through to spring 2017, the majority of glaciers south of 69°N experienced sustained spring-time readvances in 2018 and 2019 (Figure 4c). High annual retreat rates in 2018 and 2019 resulted from large summer through autumn retreat ( Figure 4c) that coincided with periods of anomalously warm ocean temperatures near the surface and at depth, and with reduced sea ice coverage ( Figure 8). Thus, we suggest that retreat in 2018 and 2019 was driven by warmer ocean temperatures. Surface runoff south of 69°N was also high during 2019, ranking second highest on record within the 1948-2019 period (Tedesco & Fettweis, 2020). However, it was at a minimum for our study period in 2018, which coincides with the period of maximum annual retreat observed in our study area (Figure 5a), suggesting that surface runoff is not the primary driver of the 2016-19 retreat.
One noticeable exception to this pattern of summer through autumn retreat was observed at Kangerlussuaq. As previously reported (Bevan et al., 2019;Brough et al., 2019), the glacier underwent a further phase of sustained calving during winter 2017 and summer 2018, with calving coinciding with mélange break-up and dispersal ( Figure S21 in Supporting Information S1; Bevan et al., 2019). Between 2018 and 2020, the annual terminus position remained roughly stable and the glacier returned to winter advance and summer retreat cycles. The glacier also stopped accelerating in 2018, but velocities remained elevated relative to their preretreat values ( Figure S21 in Supporting Information S1). The return to multimonth terminus advance over winter in 2018 and 2019 coincided with the mélange remaining in place over winter ( Figure S21 in Supporting Information S1) and has lead to the suggestion that the behavior of Kangerlussuaq is predominantly controlled by the presence, or lack thereof, of mélange (Bevan et al., 2019). Furthermore, the observed retreat between 2016 and 2018 left the ice front/grounding line at the base of a bedrock ridge (Bevan et al., 2019;Brough et al., 2019), which may have acted to, at least temporarily, promote stability (e.g., Gudmundsson et al., 2012;Schoof, 2007;Weertman, 1974) and dampen the glacier's response to subsequent environmental forcing.
In contrast to glaciers below 69°N, interannual retreat of glaciers above 69°N was more modest, with an average retreat of ∼200 m. Despite experiencing similar atmospheric conditions, modeled oceanic temperature increases and changes in sea ice/mélange above 69°N were not as marked or prolonged, and deeper ocean temperatures rarely exceeded 0°C (Figures 7 and 8; Figure S44 in Supporting Information S1). The more gradual and subdued glacier changes observed above 69°N relative to glaciers further south may therefore reflect the fact that these glaciers were not as exposed to the inflow of warmer oceanic water from the North Atlantic (e.g., Straneo & Heimbach, 2013). Furthermore, despite being subject to similar forcing, our high temporal resolution frontal position data show that retreat patterns were asynchronous and were characterized by more year-on-year variability between glaciers ( Figure 5c; Table 1). Rather than showing a regional response to environmental forcing, retreat patterns observed on glaciers above 69°N appear to be more closely governed by glacier-specific factors which act to enhance/ suppress the glacier response to external forcings (e.g., McFadden et al., 2011;Walsh et al., 2012). However, a clear exception was observed in 2019, where 12 out of the 13 glaciers underwent annual retreat (Figure 5c; Table 1). As noted previously, ice-sheet-wide runoff was the second highest on record within the 1948-2019 period in 2019 (Tedesco & Fettweis, 2020); thus, it is plausible that enhanced surface runoff in 2019 (e.g., Figure 7) may have been sufficient to overcome glacier-specific factors, resulting in the observed region-wide retreat. Nonetheless, our overall findings affirm the ongoing importance of oceanic forcing in driving the difference in behavior of east Greenland glaciers above and below 69°N and emphasize the greater role of glacier-specific factors in modulating the response of glaciers north of 69°N to external forcing. We expand on these glacier-specific influences below.

Influence of Glacier-Specific Factors
For the majority of glacier-specific factors, we see no significant differences between glaciers located either side of 69°N (Table 2). This suggests that glacier specific-factors are not a major driver of this divide in glacier behavior and that it is primarily due to differences in external oceanic forcing, particularly the northern extent of warm waters within the Irminger Current (e.g., Howat et al., 2008;Murray et al., 2010;Seale et al., 2011). However, we do see signficant differences in mean ice velocities, with the southern glaciers (mean: 5.3 m d −1 ) flowing three times faster than those in the north (mean: 1.8 m d −1 ). A greater proportion of glaciers south of 69°N was located on reverse bed slopes (Table 2), which is consistent with theory suggesting that such a scenario should facilitate more rapid retreat via a series of ice dynamic feedbacks when moving into deeper water. However, reverse bed slopes were associated with a range of changes in terminus position, including glaciers that underwent the most retreat and advance (Table 2). Thus, the impact of a reverse bed slope on glacier dynamics is variable within our study area, as has been observed elsewhere in Greenland (e.g., Bunce et al., 2018), and cannot alone account for observed differences in glacier dynamics across 69°N (e.g., Schoof et al., 2017;Sergienko, 2022).
Fjord bathymetry and water characteristics were comparable on either side of 69°N, but may explain the retreat patterns observed at specific glaciers. For example, three of the five glaciers (Apuliliip Apusiia, Christian IV and Vestfjord) that showed small advances between 2013 and 2020 terminate on stabilizing bedrock/basal ridges . In these configurations, the glaciers are protected from warmer waters at depth until they are forced from these ridges, and undercutting by the oceans only affects the short floating extension of ice beyond the ridge (e.g., Wood et al., 2021). As well as variations in basal topography, changes in fjord geometry may also enhance/suppress glacier response to forcing (e.g., Jamieson et al., 2012;O'Neel et al., 2005;Raymond, 1996). This may explain the limited terminus retreat observed at Unnamed A, Frederiksborg and Courtauld: despite following similar interannual patterns of advance and retreat to other glaciers below 69°N (Figure 4c), their magnitude of retreat was smaller and none of these glaciers retreated beyond their average seasonal fluctuations (Table 1). All three of these glaciers remained on lateral pinning points in their respective fjords throughout our study period. A similar pattern of reduced inter-annual retreat for glaciers confined to lateral pinning points has also been observed in northwest Greenland (e.g., Carr et al., 2013). Variations in local bed and fjord geometry may therefore help explain the observed variability of glacier terminus change, and corresponding changes in velocity, for glaciers in close proximity. These observations corroborate the hypothesis that glacier retreat, especially when viewed over short timescales, is at least in part controlled by local variation in glacier-specific factors (e.g., Bartholomaus et al., 2016;Bunce et al., 2018;Carr et al., 2015;Catania et al., 2018;Cowton et al., 2018;Fahrner et al., 2021;Moon et al., 2012;Warren & Glasser, 1992).

Impact of Changes in Glacier Dynamics on Ice-Sheet-Scale Mass Loss and Sea Level Rise Contribution
Variability in terminus position, ice velocity and therefore ice discharge have major implications for ice sheet mass loss and hence contribution to sea level rise. For example, during the phase of rapid retreat and acceleration observed on south and central eastern Greenland glaciers between 2001 and 2005, Greenland-wide ice discharge increased from ∼440 Gt a −1 to between ∼500 and 520 Gt a −1 in 2005 (King et al., 2018;. However, for the decade following this increase, ice discharge has been steady to declining in these regions (King et al., 2018;. Much of this increase in ice discharge between 2001 and 2005 resulted from the retreat, corresponding doubling in speed, and subsequent diffusive thinning of Helheim Gletsjer and Kangerlussuaq (Howat et al., 2007;Luckman et al., 2006;Rignot & Kanagaratnam, 2006;Stearns & Hamilton, 2007).
Although we see similar frontal position change at Kangerlussuaq between 2016 and 2018 as occurred between 2004 and 2005, increases in ice velocity (∼15-30%) and ice discharge (∼12%) were not as marked, with peak 2005 velocities and ice discharge remaining the highest observed at Kangerlussuaq (Bevan et al., 2019;. The changes observed at Kangerlussuaq also occur at the basin scale, where ice discharge from the central east basin of the ice sheet increased from ∼76 Gt a −1 in July 2016 to ∼82 Gt a −1 in July 2018 and was ∼81 Gt a −1 in July 2020 . Whilst recent changes in ice dynamics have resulted in a phase of increased sector ice discharge, they are below the peak ice discharge of ∼88 Gt a −1 recorded for the central east region during July 2005 . Nonetheless, these recent increases in ice discharge from the central east region, coupled with similar increases in ice discharge observed in the south east (e.g.,  appear to have offset declining ice discharge in the central west region-largely due to declining ice discharge from Sermeq Kujalleq (Jakobshavn Isbrae) likely resulting from ocean cooling (Khazendar et al., 2019)-and since 2018, the north west region . The annual rate of ice discharge sharply increased between 2000 and 2005, and has since then stayed at the elevated rate of discharge observed in 2005, which is ∼10% above the 2000 discharge rate . This sustained period of enhanced ice discharge has resulted in the ice sheet losing mass even in years where surface mass balance is high (e.g., King et al., 2020;Sasgen et al., 2020). However, rather than indicating stability in ice-sheet-wide dynamics, these changes highlight that the centers of dynamic mass loss show spatial and temporal variability in response to past and current oceanic and atmospheric forcing (e.g., King et al., 2020;Murray et al., 2015;Porter et al., 2018). Nonetheless, while centers of ice loss shift (e.g., the south east), some regional pattens persist, such as the divide observed at 69°N since at least 2000. Overlain on these regional signals, we have the response of individual glaciers, which can be modulated due to glacier specific factors (Section 4.3). Thus, it is important to assess patterns of glacier response to forcing at a range of scales from individual glaciers to the entire ice sheet.

Conclusions
Our intraannual time series of 24 east Greenland marine-terminating glaciers shows that while all glaciers exhibit seasonal advance and retreat cycles proportional to glacier velocity, there are distinct differences in behavior between those located north and south of 69°N. Glaciers south of 69°N underwent substantial multiyear retreat between 2016 and 2019, while those to the north showed minimal or gradual change in terminus position and velocity. The retreat to the south of 69°N coincides with warming of the ocean water column, while no substantial oceanic warming was observed to the north. Our data demonstrate that the contrasting glacier behavior observed either side of 69°N, first noted in the 2000s (e.g., Seale et al., 2011;Walsh et al., 2012), has persisted, most likely due to oceanic conditions, and impacts glacier dynamics at seasonal and interannual timescales.
South of 69°N, winter calving is a key control on multiyear retreat. Specifically, the delayed/disrupted formation of mélange that coincides with warm atmospheric and ocean temperatures in winter 2016 allowed termini to retreat throughout the 2016 winter instead of experiencing a seasonal readvance. By spring 2017, termini south of 69°N were therefore in a more retreated position as temperatures warmed, allowing continued retreat during summer/autumn. This forced the termini of multiple glaciers into regions of their fjords that were topographically conducive to retreat. This resulted in a multiyear terminus response, despite the majority of glaciers experiencing seasonal readvances in 2017/18 and 2018/19. North of 69°N, retreat patterns appeared to be more closely governed by glacier-specific factors, which acted to enhance/suppress the glacier response to external forcings. Overall, glaciers south of 69°N experienced multiyear retreat likely initiated by warm winter conditions, and those to the north of 69°N exhibited gradual terminus response, modulated by topography.

Data Availability Statement
Terminus data generated in this study are archived for download via Brough (2023), https://doi.org/10.5281/ zenodo.6904219. We also used a number of openly available datasets used in this study (also summarized in Table S1 in Supporting Information S1). Ice velocity maps of Scambos et al. (2016)