The role of near-terminus conditions in the ice-flow speed of Upernavik Isstrøm in northwest Greenland

Abstract

Upernavik Isstrøm, the largest contributor to sea-level rise in northwest Greenland, has experienced complex and contrasting ice-flow-speed changes across its five outlets over the last two decades. In this study, we present a detailed remote-sensing analysis of the ice dynamics at Upernavik's outlets from 2000 to 2021 to evaluate the details of these changes. Previous research suggested that the presence or absence of floating ice tongues strongly influences Upernavik's ice dynamics. We use several lines of evidence to document the presence of floating ice tongues, and find that, while several outlets experienced ice-tongue formation and/or loss during the study period, these changes do not explain observed fluctuations in ice-flow velocity. Further exploration of ice-dynamic forcings using a flowline model suggests that changes in basal slipperiness near the terminus have a strong impact on upstream ice dynamics and can explain the velocity variations. Our results suggest that speed fluctuations at Upernavik's outlets may be seasonally and interannually controlled by bed conditions near the terminus, and highlight the need for further research on the influence of basal conditions on complex tidewater glacier dynamics.

1. Introduction

In the current warming climate, Greenland's marine-terminating glaciers have undergone significant changes in ice flow. Northwest Greenland has lost mass as a result of recent changes in glacier dynamics. This region has most of Greenland's deep-water outlet glaciers, and from 1972 to 2018 it contributed 4.4 mm to global sea-level rise, the highest contribution out of any region in Greenland (Mouginot and others, 2019; Wood and others, 2021). Between 1972 and 2018, ice dynamics in the northwest contributed to 86% of its total mass loss compared to surface mass balance, with the percentage increasing every decade from 1998 to 2018 (Mouginot and others, 2019). Controls on ice-flow acceleration have been linked to several processes, including thinning of floating ice tongues, retreat into deeper or wider troughs and changing subglacial drainage, and ocean warming (Nick and others, 2009; Moon and others, 2015). The changing conditions of a floating ice tongue can influence velocities at the terminus, which then propagate upstream (Moon and Joughin, 2008; Nick and others, 2009). Loss of floating tongues is thought to be a driver for the acceleration of other large glaciers in Greenland, including Jakobshavn Isbræ (also known as Illulissat Glacier and Sermeq Kujalleq) (Joughin and others, 2004; Thomas, 2004; Vieli and Nick, 2011). Retreat into varying trough geometry alters the glacier terminus width and thickness, which directly relates to how fast the terminus and upstream ice flows (Nick and others, 2009). Subglacial drainage due to warmer air temperatures and surface melting can lubricate the glacier bed and also cause short-term acceleration, thinning and retreat (e.g. Zwally and others, 2002; Fried and others, 2015; Straneo and others, 2015).

One of the largest and fastest-retreating ice streams in northwest Greenland is Upernavik Isstrøm. Model results by Haubner and others (2018) showed acceleration at Upernavik was responsible for 80% of its mass loss between 1995 and 2012. Upernavik Isstrøm has five marine-terminating glaciers, several of which have contributed significantly to Greenland's recent ice-mass loss. We will refer to the five distinct outlets as U0 to U4 (Fig. 1), from north to south, although previous literature has focused on the southernmost four. The two northern outlets have had the largest recent impact on sea level; as a result of acceleration and retreat, U1 underwent an increase in discharge of ~141% from 1996 to 2013, while U2 increased by 73% from 1993 to 2018 (Mouginot and others, 2019).

Figure 1. Upernavik terminus positions and flowlines. The change in terminus position for Upernavik's outlets from 2000 to 2019 using data from PROMICE (Andersen and others, 2019) and our own 2020 and 2021 polylines over a Landsat 8 image from 2 September 2021, courtesy of the US Geological Survey. Upernavik outlet flowlines obtained from averaged ITS_LIVE (Gardner and others, 2019) velocities from 2000 to 2018.

Upernavik Isstrøm is important not just due to its size, but also because it has experienced complex and contrasting behaviors in velocity between its individual outlets. There are observed periods of deceleration at the southern outlets while acceleration occurred at the northern outlets (Khan and others, 2013; Larsen and others, 2016). The northern outlets reside in a deeper region of the fjord than the southern outlets and are likely in contact with warm Atlantic Water at depths between 250 and 450 m (Andresen and others, 2014; Vermassen and others, 2019; Muilwijk and others, 2022). The dynamics of these deeper northern outlets may be sensitive to ocean forcing, as enhanced melt at times of warm-water intrusion may reduce terminus thicknesses and explain some terminus velocity fluctuations (Rignot and others, 2010; Chauché and others, 2014; Wood and others, 2021). However, acceleration at Upernavik also occurred during periods when Atlantic Water was not as warm, suggesting warm ocean water is not the only control on ice flow (Vermassen and others, 2019).

The loss of the floating ice tongue at U1 has also been proposed as a possible reason for its acceleration (Khan and others, 2013; Larsen and others, 2016). The relationship between ice-tongue presence and tidewater glacier ice-flow dynamics is complex. The acceleration of Jakobshavn Isbræ is likely due to the loss of buttressing from the disintegration of its long floating ice tongue (Joughin and others, 2004; Thomas, 2004; Holland and others, 2008; Vieli and Nick, 2011). Slower or stable velocities are observed when Greenland glaciers terminate with floating ice tongues, including at Jakobshavn Isbræ (Nick and others, 2010; Moon and others, 2012). Contrasting behaviors between Helheim and Kangerlussuaq Glaciers in southeast Greenland suggest velocity changes may be controlled by calving behavior, which is influenced in turn by glacier geometry and floating ice-tongue extent (Kehrl and others, 2017). U1 and U2, the two Upernavik outlets that have accelerated the most recently, had evidence of floating ice tongues while Upernavik's other outlets did not (Enderlin and others, 2013; Khan and others, 2013; Larsen and others, 2016).

Changes in resistive stresses along the glacier bed and fjord walls can also influence velocities (Howat and others, 2005, 2008; Pfeffer, 2007). For example, the acceleration of Jakobshavn Isbræ may have been due to weakened shear margins rather than solely from the loss of back-stress from the breakup of its floating ice tongue (Van Der Veen and others, 2011). Joughin and others (2012) suggested that basal water pressure may have driven acceleration along with the change in ice buttressing. At Nioghalvfjerdsfjorden (79N) glacier in northeastern Greenland, modeling showed that the loss of a 76 km ice tongue did not significantly affect the glacier's dynamics, though that work assumed some enhanced sliding along the ice tongue's margins (Rathmann and others, 2017). Furthermore, changes in basal slipperiness produced realistic seasonal velocities in that model. Recently, Downs and Johnson (2022) simulated Upernavik's response to perturbations in basal drag and submarine melt rates to investigate mass loss. They concluded that while subaqueous melt of the terminus was the primary driver of mass loss, the glacier remained sensitive to basal drag despite its rapid retreat, and changes in basal drag likely caused acceleration and thinning upstream (Downs and Johnson, 2022).

In this study, we begin by investigating the hypothesis that the changing conditions of floating ice tongues are the main driver of significant terminus acceleration at Upernavik over time. We use remote-sensing data to update the records of all outlets since 2013, then assess major changes over the period 2000–2021. We find significant changes in velocities since 2013, indicating that changes to the system have continued since the last detailed analysis of remote-sensing observations (Larsen and others, 2016). We use tabular icebergs, plume–polynyas and slope-break patterns to determine when a floating ice tongue existed at each of Upernavik's outlets. We estimate and analyze retreat, calving type, subglacial drainage, elevation profiles and thinning. We then compare these data with velocity data to analyze whether velocity fluctuations can be explained by changing conditions of floating ice tongues.

The results of this analysis suggest that floating ice-tongue changes cannot explain the variations in velocity observed at Upernavik's outlets. We therefore use a simple flowline model to explore the stress balance of Upernavik's two largest outlets, U1 and U2. Our modeling analysis suggests that changes in basal conditions at the terminus may plausibly explain the observed ice-flow-speed variations, and this result is consistent with 3-D modeling studies investigating similar questions (Downs and Johnson, 2022). The combined results of our remote-sensing analysis and simple modeling exploration suggest that Upernavik's large velocity fluctuations may be explained by changes in bed conditions near the terminus, highlighting the need for further observational and modeling studies focused on basal conditions at tidewater glaciers.

2. Data and methods

The majority of our data and results were extracted along outlet flowlines (Fig. 1) to estimate bed topography and near-terminus bathymetry, elevation and ice flow. The flowlines were generated using stacked ITS_LIVE mosaics from 2000 to 2018, which required an input coordinate at the upper-glacier starting point for each outlet and used along-flow velocities to generate the flowline (Gardner and others, 2019). For the remote-sensing analysis, we used ~30 km long flowlines to capture near-terminus ice dynamics for the three central and long trunks, U1, U2 and U3, and ~15 km for the smaller U0 and U4. For our ice-flow model, we used ~75 km long flowlines for U1 and U2 to capture ice-flow sensitivity from the Greenland Ice Sheet.

2.1 Terminus and bed observations

Landsat 7 and 8, courtesy of the US Geological Survey, provide a detailed record of visual changes up to 2021. Terminus positions were available from 1999 to 2019 from PROMICE, derived from Landsat, Aster and Sentinel-2 imagery for all outlets (Andersen and others, 2019). We supplemented the existing PROMICE record by digitizing the terminus positions for all outlets once per year from 2019 to 2021, using Landsat imagery and the same method as PROMICE. The end of the melt season was determined by PROMICE as the most retreated position of the glacier from July through November (Andersen and others, 2019).

While it is difficult to gather evidence about specific subglacial hydrologic patterns using remote sensing, some information on channelization and drainage events can be gleaned by tracking plume–polynyas. Plume–polynyas are open-water areas that form when subglacial channels release buoyant meltwater that rises to the glacier front (Melton and others, 2022). We manually inspected all available and cloudless Landsat 7 and 8 imagery from 2000 to 2021 and developed a record of all the occasions where plume–polynyas were visible to relate to both floating ice-tongue evidence and seasonal velocity changes.

We used NASA's Operation Icebridge BedMachine Greenland, Version 4 (Morlighem and others, 2021) as an estimate to evaluate the depth of the fjord near each outlet and the bed topography of each outlet so we could assess the water bodies in contact with the ice front and how future retreat may affect the glacier geometry. For U0, U1 and U2, BedMachine elevations in the fjord are inferred from gravimetry, while more accurate swath sonar measurements are used to constrain bathymetry at U3 and U4.

2.2 Velocity

We used a variety of ice-flow data available at Upernavik to investigate the timing and magnitude of acceleration. Most of these data came from annual grids from the NASA MEaSUREs ITS_LIVE project, which are updated through 2018 (Gardner and others, 2019). We supplemented these mosaics with more recent data, including averaging data from the ITS_LIVE Global Glacier Velocity Point Data Access (Gardner and others, 2019) and GoLive image pairs (Scambos and others, 2016), which extended the record to 2021. GoLIVE image pairs from four Landsat path/rows covering the Upernavik region were interpolated onto a common grid and stacked. A mosaiced product for each summer season was created by taking the mean of each pixel stack. Annual ITS_LIVE velocities are available at a 240 m resolution while GoLIVE velocities are available at a 300 m resolution. As they are both derived using panchromatic optical imagery, which is only available during months with sufficient sunlight at Upernavik's latitude, we used these data as an estimate of summer velocities from approximately May to September. Error estimates were provided by ITS_LIVE and range from 1 to 326 m a⁻¹ for the mosaics and 18 to 99 m a⁻¹ for the point data.

NASA's MEaSUREs project also offers a Greenland Ice Sheet Velocity Map for winter velocities using InSAR speckle tracking data (Joughin and others, 2015). This is available sporadically across the time period of interest with 200 and 500 m resolutions (Joughin and others, 2015). Annual winter velocities are derived from varying sources from approximately September to May, including RADARSAT-1, ALOS, TerraSAR-X/TanDEM-X and Sentinel-1A and Sentinel-1B (Joughin and others, 2015). The winter velocities are derived from the fall of one year to the early spring of the next. Error estimates for the MEaSUREs dataset range from 0.5 to 115 m a⁻¹.

2.3 Evidence for floating ice tongues

In order to test the hypothesis that ice tongues are the primary control on ice flow, we needed a variety of evidence to assess the presence of floating ice. We used a multi-dataset approach that relied on identifying tabular icebergs throughout the entire time period, extensive comparison of ice elevations and thicknesses in a hydrostatic analysis for floating ice and ArcticDEM elevations (Porter and others, 2022) to identify slope-break and horizontally sloped floating ice at the terminus. Our analysis focused on the largest outlets, U1 and U2, where floating ice tongues have been observed during two periods in the 2000s and the 2010s.

The first method we used to determine the presence of floating ice tongues was observations of tabular icebergs throughout the study time period from 2000 to 2021. Marine-terminating glaciers in deep fjords can produce tabular icebergs, which are considered evidence of (near) floatation because full-thickness tabular bergs can only separate and remain upright from a terminus if the ice is sufficiently thin (Joughin and others, 2008; Amundson and others, 2010, Kehrl and others, 2017; Melton and others, 2022). Larsen and others (2016) manually identified tabular icebergs being produced from U3 not necessarily from floatation but due to a long, lightly grounded horizontal ice-surface slope and shallow bed topography near the toe. Though tabular icebergs were tracked by Larsen and others (2016) from 2009 to 2010 for U2 for its floating ice tongue, we looked for tabular icebergs and evident changes in calving behavior and type throughout the entire time period and for all outlets by manually identifying upright icebergs using Landsat 7 and 8 imagery in order to provide supporting evidence indicating that floating ice tongues were present.

The second method we used to determine the presence of floating ice tongues was the hydrostatic analysis for floating ice, which has been shown to be a reliable method in many settings (e.g. Wild and others, 2022; Bentley and others, 2023). When IceBridge thickness data were available, we conducted a hydrostatic elevation analysis to identify floating termini. For this, we used

(1)$$Z_{\rm s} = \left({1-\displaystyle{{\rho_{\rm i}} \over {\rho_{\rm w}}}} \right)H, \;$$

where ρ _i = 917 kg m⁻³ is the density of ice, ρ _w = 1026 kg m⁻³ is the density of seawater, H is the ice thickness (derived from radar data as described in the next paragraph) and Z _s is the hydrostatic ice-surface elevation (e.g. Jenkins and Doake, 1991). As the terminus of Upernavik's outlets lies in the ablation zone, we assumed the entire thickness is solid glacial ice. If Z _s matched or exceeded the measured surface elevation, this was taken as evidence of floatation, while Z _s below the measured surface elevation indicated a grounded terminus.

We constructed a record of ice thicknesses from 2010 to 2017 for each outlet using a mixture of data products from NASA's Operation IceBridge Multichannel Coherent Radar Depth Sounder (MCoRDS). Where available, we used ice thicknesses directly from the pre-traced MCoRDS L2 Ice Thickness Version 2 (Paden and others, 2010). Only 2013 had coverage for all outlets. The MCoRDS L1B Geolocated Radar Echo Strength Profiles Version 1 dataset was used to fill in gaps, as this allowed us to take our own echogram picks when the L2 product could not detect the ice thickness (Paden and others, 2014). Our echogram picks were selected by manually picking the strongest reflector, which we assume indicated the ice bottom. We obtained thickness measurements for U0 in 2011 and 2013. U1 had five IceBridge flightlines with good coverage available in 2010, 2011, 2012, 2013 and 2014. Two IceBridge flightlines with good coverage were available for U2, in 2013 and 2017. We were able to obtain picks for six flightlines available for U3 in 2010, 2011, 2012, 2013, 2015 and 2017. U4 only had one IceBridge flightline available, in 2013 (Figs S6–S16).

The third method we used to determine the presence of floating ice tongues was a detailed elevation analysis using ArcticDEM. Though there are multiple sources of ice-surface elevations for Upernavik, the most recent edition of ArcticDEM provided detailed elevation data for each outlet over recent years. ArcticDEM includes 2 m resolution strips from 2010 to 2021. We estimated vertical error by using rock outcrop elevation points near the glacier termini from the ArcticDEM high-resolution mosaic and calculating the mean squared difference from all the overlying elevation strips. Slope break near the terminus of glaciers has been used as evidence of the grounding line: ice with relatively steep slopes is considered to be grounded, and an abrupt shift to nearly horizontal slopes marks the transition to floating ice (Bindschadler and others, 2011; Enderlin and others, 2013; Li and others, 2022).

In addition to analyzing slope breaks, we used ArcticDEM to analyze ice-surface slope change over time and used the ice-surface elevations to investigate thinning rates for each outlet as a measure of glacier change in relation to acceleration and reaching floatation. Thinning rates from ArcticDEM were calculated from year to year and from 2011 to 2021 to compare between glaciers as U1 and U0 have very sparse and incomplete 2010 data.

2.4 Numerical modeling

To supplement remote-sensing observations, we used a 2-D, flowline ice-flow model to help us understand how the presence of floating ice tongues, bed properties and stresses relate to changes in ice flow. We chose to examine the two largest outlets, which both experienced changes in ice-tongue presence and had variable velocities during the study period. We first established reasonable parameters for basal drag using an inverse 2-D model on the earliest elevation and ice thickness datasets available. We then ran diagnostic, flowline simulations of U1 and U2 in multiple years to test whether changes in floating ice tongues and changes in driving stress (caused by thickness changes) can explain the observed changes in velocity. Since these simulations were unable to reproduce the observed velocities, we introduced sidewall drag (i.e. the friction between the glacier tongues and their margins) and ran further simulations varying sidewall and basal drag.

2.4.1 Model physics

Given the fast but confined flow experienced by the outlets of Upernavik, a model of these outlets should be able to capture both internal deformation and basal sliding. We thus modeled the system using Blatter–Pattyn higher-order equations to describe ice flow along a central flowline (Blatter, 1995; Pattyn, 2003), which allows both flow regimes while having considerably reduced complexity compared to the full Stokes equations (e.g. Shapero and others, 2021). Following Shapero and others (2021), the model equations can be written in terms of an action functional, J, that must be minimized

(2)$$\begin{align}J = & \displaystyle{n \over {n + 1}}\mathop \smallint \limits_{\rm \Omega } \mathop \smallint \limits_0^1 HA^{{-}1/n}\sqrt {\vert {\dot{\varepsilon }( {\boldsymbol u} ) \vert^2 + H^{{-}2}} \vert \partial {\boldsymbol u}/\partial \zeta \vert ^2} ^{1/n + 1}{\rm d}\zeta \;{\rm d}x \cr & + {\rm friction}\hbox{-}{\rm gravity}\hbox{-}{\rm terminus, \;\;}\end{align}$$

where n is the flow exponent (assumed to be 3), A(T) is the temperature-dependent fluidity rate factor from Glen's law, u is the velocity, $\dot{\varepsilon }$ is the horizontal strain rate, ζ is the normalized vertical coordinate (z/H) and x is the model domain. The first term in Eqn (2) describes viscous effects, and we are left to describe other forces. We take

(3)$$\eqalign{{\rm friction\;} = \mathop \smallint \limits_\Omega C_{{\rm bed}}\vert {{\boldsymbol u}( {\zeta = 0} ) } \vert ^2{\rm d}x + \displaystyle{m \over {m + 1}}\mathop \smallint \limits_{\rm \Omega } \mathop \smallint \limits_0^1 HC_{{\rm side\;}}\vert {\boldsymbol u} \vert ^{{1 \over m} + 1}{\rm d}\zeta \;{\rm d}x}$$

where m is a friction exponent (assumed to be 3) and C _bed and C _side are friction coefficients. The first term in Eqn (3) describes friction at the ice–bed interface, while the second describes friction from the contact between a glacier and its side walls. Here we have assumed that the domain is a flowline, since otherwise the outer integral in the second term must only be evaluated on the sidewalls of the domain. The resistive forces of Eqns (2) and (3) are balanced by the driving effects of gravity and the pressure difference at the terminus, given by

(4)$${\rm gravity\;} = {-}\mathop \smallint \limits_{\rm \Omega }^\;\mathop \smallint \limits_0^1 \rho _{\rm i}gH\nabla s\cdot {\boldsymbol u}{\rm d}\zeta \;{\rm d}x$$

where g is gravity and s is the surface elevation. The final term in the action is then

(5)$${\rm terminus\;} = \mathop \smallint \limits_{\rm \Gamma }^\;\mathop \smallint \limits_0^1 ( {\rho_{\rm i}-\rho_{\rm w}} ) {\boldsymbol u}\cdot {\boldsymbol n}{\rm d}\zeta \;{\rm d}\gamma $$

where n is the outward-pointing normal vector, and the integral is evaluated over the lateral boundary of the domain.

2.4.2 Inference of basal drag

In order for the model to produce velocities that matched observations, it was necessary to apply a realistic, spatially variable basal drag. We inferred the basal drag coefficient, C _bed, using a 3-D model implemented in Icepack (Shapero and others, 2021). The model physics were described by the Blatter–Pattyn equations (in the form of Eqns (2–5) above), and the unknown C _bed was inferred using standard glaciological inverse methods (e.g. MacAyeal and others, 1993). This inversion covered the entire Upernavik catchment with 300 m horizontal resolution. The model was solved on a single vertical layer with second-order Gauss–Legendre trial and test functions. The surface and bed were extracted from Bedmachine (Morlighem and others, 2021). The necessary boundary conditions were constraints on the velocity on the sides of the domain; we took the standard approach of fixing the velocity to observations at all edges except the terminus, and applying the hydrostatic pull at the terminus (Eqn (5)). In effect, this meant that the sidewall drag (second term in Eqn (3)) was subsumed into the boundary condition and no explicit sidewall drag was calculated.

The inversion sought to minimize the root mean square misfit between the modeled velocity and the MEAsUREs multiyear velocity product (Joughin and others, 2018). The optimization was done using the Rapid Optimization Library (Ridzal and others, 2017). The inversion used Tikhonov regularization, with strength determined using an L-curve analysis (Calvetti and others, 2000). The inversion results were used to determine C _bed for the basal boundary condition of the flowline models described below.

The mismatched timestamps between the MEAsUREs multiyear product (1995–2015; Joughin and others, 2018) and the Greenland Ice Mapping Project surface used in Bedmachine (~2007; Howat and others, 2014; Morlighem and others, 2021) potentially cause some change in surface elevation to be erroneously interpreted as effects of changes in basal shear stress. However, this problem is likely relatively minor, as the timestamp of the DEM lies within the dates spanned by the data used in the velocity product. Moreover, the issue of such mis-matched timestamps for inversions is largely unavoidable except by newly developed time-dependent inverse methods (e.g. Choi and others, 2023), which are beyond the scope of this work, and instead most studies take a similar approach to that used here (e.g. Downs and Johnson, 2022).

2.4.3 Flow-line modeling

We simulated changes along the flowlines of U1 and U2 using the open-source, finite-element model Icepack (Shapero and others, 2021). We used Icepack's hybrid flow solver, i.e. the Blatter–Pattyn higher-order equations (Eqns (2–5)). The model uses terrain-following coordinates with one vertical layer; despite this low resolution, the model is still able to resolve realistic vertical variations in velocity by using fourth-order Gauss–Legendre trial and test functions in the vertical.

We ran diagnostic simulations (i.e. snapshots) of flowlines for U1 and U2, derived from ITS_LIVE mosaics averaged from 2000 to 2018. We extracted ArcticDEM surface elevations for the years of interest for our model runs and BedMachine Greenland bed topography along the flowline (Morlighem and others, 2021; Porter and others, 2022). For years with floatation, we determined where the likely ice bottom was compared to the bedrock by using a reverse hydrostatic calculation. We used the hydrostatic ice bottom rather than the bed elevation for the thickness calculation whenever the ice bottom was shallower than the bed. This was most important for U2 in 2013, when the ice tongue was significantly thinner near its collapse.

We used a depth-varying temperature to determine the fluidity rate factor A. Annual mean surface temperature was obtained from HIRHAM from 1980 to 2014 (Langen and others, 2015). We assumed the reanalysis value at the surface varied linearly to the pressure melting temperature at the bed. While this assumption likely has a slight warm bias, the bed is likely thawed (MacGregor and others, 2016) and this is the simplest assumption that follows the expected trend of the temperature profile without requiring a separate model. We solved for the pressure melting temperature, T _m, through

(6)$$T_{\rm m} = T_{{\rm tp}} + \gamma _{\rm p}( {\rho_{\rm i}gH-p_{{\rm tp}}} ) , \;$$

where the triple-point temperature and pressure, T _tp and p _tp, were 273.16 K and 611.73 kPa and the Clausius–Clapeyron constant for pure ice and air-free water, γ _p, was 7.42⦁10⁻⁵ K kPa⁻¹ (e.g. Cuffey and Paterson, 2010).

We required boundary conditions for the bed, inflow and outflow. At the inflow, we fixed the horizontal velocity to match the surface velocity for all depths. The outflow boundary is handled by Eqn (5), where we took ρ _w = 1024 kg m⁻³. At the bed, we used a spatially variable C _bed, with values evaluated directly from the 3-D inversion described above. It was necessary to use this nested model approach (boundary condition for the 2-D model extracted from a 3-D one) due to limitations in the available data; imperfect coverage of the elevation and velocity products prevented simulating all outlets in all years. Though not technically a boundary condition, for initial simulations, we used C _side = 0 in Eqn (3). Initial flowline model runs (simulating 2011 for U1 and 2013 for U2) produced velocities much faster than observations, despite using a basal drag from the inversion procedure described above. These outlets flow through confined fjords, where sidewall drag may be a significant component of the force balance (Gagliardini and others, 2010), so we incorporated a non-zero C _side in Eqn (3). Because the 3-D model used for inversions spanned the width of the outlet glaciers and continuing into the areas of thin ice surrounding them, the effects of the sidewalls were naturally incorporated, but that resistance must be explicitly included in the flowline model. To determine realistic values for the sidewall drag, we performed inversions along each flowline to infer the sidewall drag needed to reproduce the observed velocities. The optimization procedure used the same optimization and regularization techniques as described above for the 3-D modeling. We conducted the sidewall drag inversion only along our flowlines for the initial model years: 2011 for U1 and 2013 for U2. The inversion procedure sought to optimize C _side to minimize the misfit between the modeled and observed velocities. Using the inferred sidewall drag, modeled velocities closely matched the observations. This inferred C _side was then used as input for subsequent simulations.

Based on the available data and the most interesting changes, we ran five simulations:

• • U1 2011: with no evidence of floatation

• • U1 2014: with evidence of re-floatation

• • U2 2013: with evidence of floating ice tongue

• • U2 2015: with no evidence of floatation; retreat

• • U2 2018: with evidence of re-floatation; thinning

The surface and bed geometry and the position of the calving front were updated for each simulation based on the remote-sensing observations, allowing us to test how velocities responded to varying geometry. We explored the sensitivity of the surface velocity in each model simulation to varying basal and sidewall drag coefficients as well as observed retreat, thinning and bed topography.

While the Blatter–Pattyn equations provide a good approximation to the full Stokes equations, the flowline modeling used here makes significant simplifications to the complex reality of the physical system. Because of these simplifications, results of the modeling should be interpreted cautiously; more detailed analysis of the limitations of the modeling can be found in section 4.3 below.

3. Results

3.1 Velocities and terminus positions

The terminus positions of Figure 1 and the Hovmöller velocity diagrams of Figure 2 show some consistent behaviors. As is typical of marine-terminating glaciers, velocities for all outlets are faster near the glacier terminus than upstream. We also note that the termini of the outlets are either retreating or remaining relatively steady throughout the record, with no distinct periods of sustained readvance (with the possible exception of a very slight readvance at U1 during the second half of the record).

Figure 2. Upernavik flowline summer velocities. (a–e) Hovmöller summer velocity diagrams from 2000 to 2021 for the flowlines of U0–U4, respectively, ending at the terminus position in each year, determined using data from PROMICE (Andersen and others, 2019). Velocity scales vary between outlets. All velocity data from 2000 to 2018 are from ITS_LIVE velocity mosaics and from 2019 to 2021 the data are from stacked GoLIVE velocities (Gardner and others, 2019; Scambos and others, 2016).

However, beyond these broad similarities, the outlets show contrasting behavior in both terminus position and velocity. The most rapid retreat events are seen on U0 and U1, which were connected as a single outlet until 2006. U1 retreated abruptly in 2007 (totaling close to 5.5 km between 2000 and 2008) before stabilizing. U0 experienced a steady retreat up until this time, and then experienced a smaller abrupt retreat event around 2011. This stands in contrast to the steadier retreat of U2 and U4 over the same time period, with smaller stepped retreat events in 2010 (~1 km) and 2015 (~2 km) at U2 and 2008 (~0.7 km) at U4. Meanwhile, U3 maintained a relatively steady ice-front position throughout the record.

Velocity changes are similarly complex, and difficult to represent on a single diagram. In the Hövmoller diagrams in Figure 2, only U3 maintains an approximately steady velocity pattern along its central flowline throughout the record, which is consistent with the approximately steady terminus position. The four other outlets retreat throughout the record, sometimes smoothly and sometimes episodically, and also show fluctuating velocities on interannual timescales. However, the velocity fluctuations generally do not match up with time periods that terminus retreat took place. For example, a distinct acceleration began at U1 a few years prior to its large retreat in 2007. While many glaciers accelerate as they retreat (e.g. Nick and others, 2009), the terminus of U4 decelerates after its retreat in 2008. While velocity generally increases at U2 as the terminus retreats, stepwise retreat events do not seem to correspond closely to distinct changes in velocity. Figure 3 shows terminus velocities throughout the record ~400 m from the glacier front of each outlet. Note in particular the two distinct peaks in speed observed at the termini of U1 and U2. While these outlets are next to each other in the fjord and both exhibit variable terminus velocity, their velocity peaks occur at different times.

Figure 3. Terminus velocity comparison at Upernavik. A comparative seasonal terminus velocity plot of U0–U4 showing ITS_LIVE summer (S) and MEaSUREs winter (W) ice velocities from 2000 to 2021 (Joughin and others, 2015; Gardner and others, 2019).

Figure 4 shows velocities at four points along the flow lines used in the Hovmöller diagrams for each outlet, with data included for both summer and winter. U2 and U4 show fairly consistent patterns of summer velocities being slightly below the winter velocities. U0, U1 and U3 also show evidence of seasonal variability, but the patterns are inconsistent from year to year. In general, seasonal variations are of a much smaller magnitude than interannual variations for all outlets. Figure 4 also shows that all outlets experience interannual variability of the greatest magnitude close to the terminus.

Figure 4. Upernavik seasonal velocities. (a–e) Seasonal velocity from 2000 to 2021 for U0–U4, with evidence of floatation provided for U1, U2 and U3. The four points along each glacier were selected by identifying the terminus and end flowline point (upper-glacier) and selecting the lower-middle (LM) and upper-middle (UM) points evenly spaced in between. These points change over time as the terminus retreats. S and W indicate summer velocities derived from ITS_LIVE data and winter velocities derived from MEaSUREs data (Joughin and others, 2015; Gardner and others, 2019), respectively. Evidence for floatation is determined by available hydrostatic data, calving behavior and previous studies. Likely floatation is evidenced by elevation profiles. No floatation is evidenced by elevation profiles, calving behavior and/or when hydrostatic data show no floatation.

3.2 Calving behavior

For floating ice-tongue evidence, including calving behavior, hydrostatic elevation and elevation, we have summarized the results for each outlet in Figure S1. Significant transitions in calving behavior from tabular to non-tabular calving were seen for U1 and U2. U1 produced large tabular icebergs, typically from 0.5 to 2.5 km long, until the end of summer in 2007 when Khan and others (2013) and Larsen and others (2016) estimated its floating extension fully broke up. U1 then transitioned to primarily non-tabular calving. U2 also produced large tabular icebergs, around the same length as U1, from 2000 to the end of summer in 2014 (Fig. S2) and then transitioned to infrequent non-tabular calving. U3 produced tabular icebergs throughout the entire time period from 2000 to 2021, suggesting that it has been near or at floatation for the last two decades (Fig. S2). U0 and U4 did not produce any tabular icebergs and experienced minimal calving.

3.3 Hydrostatic analysis

Figure 5 displays the five flightlines (two for U1, two for U2 and one for U3) that showed evidence of floatation along some portion near the terminus, where the hydrostatic surface height calculated based on ice thickness matches the observed surface height. Though hydrostatic elevation lines in 2010, 2011 and 2012 did not show floatation (Figs S8–S10), U1 reached floatation according to the MCoRDS ice thickness profiles in both 2013 and 2014 (Figs 5a, b). There were higher terminus elevations past the slope-break toward the ocean in 2013 compared to the more horizontal slope in 2014, suggesting that the terminus may have been partially floating and partially pinned to the bed. Our hydrostatic results for U2 showed floatation in both 2013 and 2017 (Figs 5c, d). In 2013, there was a long floating extension of the terminus below 100 m in elevation that was mostly at hydrostatic elevation or very close. This floating extension was lost during a period of terminus retreat, and a small floating tongue reformed by 2017 when the terminus region thinned to hydrostatic elevation. U3 was not at the hydrostatic elevation in any profiles (Figs S12–S15) until the last observation in 2017, when the terminus was near or at floatation (Fig. 5e). U0 and U4 did not have hydrostatic evidence of floatation in any available year (Figs S6–S7, S16).

Figure 5. Upernavik hydrostatic profiles. (a) U1 hydrostatic elevation compared to the actual elevation on 18 April 2013. (b) U1 hydrostatic elevation compared to the actual elevation on 26 April 2014. (c) U2 hydrostatic elevation compared to the actual elevation on 18 April 2013. (d) U2 hydrostatic elevation compared to the actual elevation on 10 April 2017. (e) U3 hydrostatic elevation compared to the actual elevation on 10 April 2017. The horizontal axis for U1 and U2 covers the same glacier range over multiple years. The hydrostatic elevation, calculated using the hydrostatic equilibrium equation (1) with MCoRDS data (Paden and others, 2010; Paden and others, 2014).

3.4 Elevation and bed topography

Figure 6 shows available surface elevations along the central flowlines of each outlet from ArcticDEM during the study period, along with bed elevations from BedMachine (Morlighem and others, 2021). These profiles show that most of the outlets thinned throughout the record, with most of the thinning concentrated in the lower reaches of the outlets, closer to the terminus. The fastest thinning rates were observed at the northernmost three outlets, with average thinning rates of ~5.5, 3 and 6 m a⁻¹ at U0, U1 and U2, respectively. U3 showed almost no elevation change throughout the record, consistent with its unchanging velocity and terminus position patterns. U4 thinned at an average rate of 2.5 m a⁻¹. Estimated ArcticDEM error based on nearby rock outcrop elevations was consistently under 3 m between 2011 and 2021.

Figure 6. Upernavik bed and ice elevation profiles. (a)–(f) Bed and ice elevation profiles from 2010 or 2011 to 2021 for U0–U4. Data are obtained from BedMachine (Morlighem and others, 2021) and ArcticDEM from 2011 to 2021 (Porter and others, 2022). Error was calculated for ArcticDEM and upper and lower bounds are included in this plot.

Along with revealing surface elevation change, the ArcticDEM profiles show spatial patterns in ice-surface slope. A distinctive break in surface slope from steeper to nearly flat near the terminus often represents the transition from grounded to floating ice (e.g. Wild and others, 2022; Bentley and others, 2023). This pattern is clearly visible in most of the elevation profiles in Figure 6, but specifically at U1 throughout the record with the exception of the period between 2008 and 2012. This pattern is also seen until the major terminus retreat and floating tongue disintegration in 2014, and it reforms as a new floating tongue in 2017.

The BedMachine bed topography shows that all five outlets are grounded below sea level near their termini. U0–U3 all have termini at or near bedrock highs with inland reverse bed slopes, while the bed of U4 has a consistently prograde slope. U1 and U2 are the most deeply grounded outlets, with termini 400–500 m below sea level. U3 is currently grounded ~300 m below sea level, but a reverse bed slope just upstream of the terminus leads to a bed grounded closer to the depths of U1 and U2.

3.5 Evidence of plume–polynyas

Distinct plumes were visible at the terminus of U0 through U3 in many years between June and August. We were limited by the lack of mélange at the terminus of U4 to identify plume–polynyas. U0 through U3 all appear to have sporadic subglacial hydrologic outflow patterns, with no clear patterns in plume–polynya appearance other than clustering around the mid to late 2000s and the mid 2010s. Plume–polynyas were observed at U1 and U2 both when floating ice tongues were present and when they were not.

U0 and U2 both experienced large subglacial drainage events, evidenced by larger-than-average plume–polynyas, often spanning across the terminus. U0 experienced several large outbursts over our record whereas U2 only experienced one shared drainage event in the same 2-week period in late June and early July of 2010 (Figs S3, S4) as observed in available imagery. Three large supraglacial lakes upstream of U2, which has the largest upstream lakes (~3 km in diameter) for this ice stream, drained during the 2010 event (Fig. S5).

3.6 Modeling

The inverted basal and sidewall drag produced modeled ice-flow velocities that closely matched observations at the time of model initialization, observable in Figures 7 and 8 (2011 for U1 and 2013 for U2). We investigated whether changes in driving stress alone (based on a new surface elevation) could explain the slowdown observed by 2014 at U1. However, modeled velocity remained much faster than 2014 observations when forced with the 2014 surface geometry. A simple calculation supports this conclusion. If we approximate the driving stress along the flowline as ρ _igHsinθ where θ is the ice-surface slope, then an observed average thinning rate of 3 m a⁻¹ would lead to a reduction in driving stress of ~1% over the 3-year time span between model runs. This cannot explain a 50% observed decrease in ice-flow velocity near the terminus during this time period. It is therefore likely that an increase in sidewall drag, a decrease in basal slipperiness or both are responsible for the observed velocity decrease.

Figure 7. U1 model inputs and ice-flow speeds. (a) U1 measured and modeled velocity using the same drag inputs (C _base, C _side) over time. (b) U1 measured and modeled velocity in 2011 and 2014 using modified basal drag coefficient, along with ice-surface elevation and bed depth. (c) U1 basal temperature in 2011 and 2014. (d) U1 basal drag inversion output (C _base) in 2011 and manipulated drag in 2014. (e) U1 sidewall drag inversion output (C _side) used in both model runs. The flowline starts at 0 m at the upper glacier and approaches the terminus point toward 75 km.

Figure 8. U2 model inputs and ice-flow speeds. (a) U2 measured and modeled velocity using the same drag coefficient inputs (C _base, C _side) over time. (b) U2 measured and modeled velocity in 2013, 2015 and 2018 using modified drag coefficient inputs, along with ice-surface elevation and bed depth. (c) U2 basal temperature in 2013, 2015 and 2018. (d) U2 basal drag inversion output (C _base) in 2013 and the manipulated drag coefficient in 2015 and 2018. (e) U2 sidewall drag inversion output (C _side) in 2013 which was also used in 2015, and the manipulated drag coefficient in 2018. The flowline starts at 0 m at the upper glacier and approaches the terminus point toward 75 km.

Figure S17 shows a scenario with increased sidewall drag that reproduces the observed velocities at U1 in 2014. The sidewall drag would have to have increased by more than 50% between 2011 and 2014, with little thinning and retreat, to decrease the modeled velocity enough to match the observations with spatial velocity inconsistencies (Fig. S17). Much smaller changes in basal slipperiness (Fig. 7d) were required for modeled velocities to match observed 2014 velocities at U1. Changes in the basal drag coefficient alone beneath small areas affected the ice flow upstream, and so small changes in slipperiness were sufficient for the model to reproduce the observed velocity in 2014.

U2 experienced more complex changes over time than U1, with terminus retreat and the loss of its floating ice tongue between 2013 and 2015, and then high thinning rates between 2015 and 2018. Following model initialization, we examined if the inferred 2013 basal drag coefficient and sidewall drag would yield velocities that matched observations in 2015 and 2018. Modeled terminus velocities were almost 4000 m a⁻¹ too fast in 2015 when only accounting for changes in glacier geometry (Fig. 8a). In 2018, the model more closely matched observations, with small deviations at the ice front. Reduced driving stress as a result of thinning likely caused some slowdown by 2018, but did not completely explain the slowdown at the terminus. For both 2015 and 2018, we investigated whether adjusting the sidewall drag coefficient would give a better match to observed velocities (Figs S18, S19). For 2015, similar to U1, the sidewall drag coefficient would have to have increased significantly (>30%) to cause the observed slowdown and the velocity result was still spatially inconsistent near the terminus (Fig. S18). Although slightly increased sidewall drag (10%) could help explain increased velocities near the terminus in 2018, this resulted in inconsistent upper-glacier velocities (Fig. S19). We then manipulated the basal drag coefficient using the same method as U1 and found that realistic changes in the basal drag coefficient could explain the slowdown in 2015 (Fig. 8d). To fully explain the 2018 U1 velocities, reasonable decreases in both sidewall drag and basal drag were required (Figs 8d, e).

4. Discussion

4.1 Floating ice tongue evidence and methods

Our methods for identifying floating ice-tongue presence yielded mostly consistent results, and we had most confidence in floating ice-tongue presence when results from multiple methods agreed. We found U1 to be floating from 2000 to 2007 and from 2013 onwards while U2 was floating from 2000 to 2014 and reached floatation again in 2017. Tabular icebergs appeared when both U1 and U2 were floating, which is consistent with iceberg records by Melton and others (2022) and Amundson and others (2010) when Helheim Glacier and Jakobshavn Isbræ were floating. Both hydrostatic elevation and ice-surface slope supported the tabular iceberg observations, indicating floating ice-tongue presence. When floatation occurred again for U1 and U2, tabular icebergs were not produced, yet hydrostatic elevation, a well-defined slope break and thinning to a horizontal slope strongly indicated the presence of a floating ice tongue. U3 produced tabular icebergs throughout the record while not floating, including during 2017 when a short-lived floating ice tongue was evidenced by hydrostatic elevation and surface slope. This may indicate that tabular icebergs alone are not a reliable proxy for floatation.

While Melton and others (2022) suggested observations of plume–polynyas (which appeared when the glacier was grounded and disappeared when the glacier was floating due to interactions between the ocean and the subglacial hydrologic system) as a proxy for the floating state of Helheim Glacier, we found contrasting observations at Upernavik. Out of available images, we observed plume–polynyas for both floating ice tongues and grounded termini. These contrasting observations may be related to differing degrees of floatation at Helheim and Upernavik, contrasts in ice-tongue length, variability in floatation degree across the glacier front and/or differences in the strength of the channelization of subglacial water flow. The variability in our results and contrasts with related literature highlight the importance of using a variety of evidence when investigating the floatation of marine-terminating glaciers, with the hydrostatic elevation and surface slope being the most reliable methods used in this study. When considering the hydrostatic elevation and using ice thickness to measure the floatation of glaciers, we suggest that assessing bed slope and bed depth is equally important to distinguish lightly grounded glaciers. Tabular iceberg calving and plume–polynya presence may serve as secondary evidence, as the geometric characteristics of the grounding zone may dictate different types of calving and meltwater channelization.

While our observational results provided new insights into floatation and thinning for Upernavik's outlets, we were limited by available ice thickness data. The rapidly changing thicknesses and floatation conditions at the termini of U1 and U2 are best constrained by ice-penetrating radar, highlighting the importance of future ice-thickness data collection. However, the multiple methods used and consistency with past studies give us confidence in our assessment of the presence of floating ice tongues at Upernavik.

4.2 Causes of complex ice-flow changes at Upernavik's outlets

4.2.1 Impacts of floating ice tongue changes on ice flow

Changes in velocity did not correlate with changes in floating ice tongues at U1 and U2 (Fig. 4), contradicting the hypothesis that the changing floating conditions were responsible for changing ice flow. The loss of resistive buttressing as a result of the disintegration of floating ice tongues or large calving events has been shown to be a driver of acceleration for many marine-terminating glaciers across Greenland (Joughin and others, 2004; Thomas, 2004; Howat and others, 2005, 2008), although it is likely that some ice tongues – even large ones – may have very little impact on a glacier's stress balance (Rathmann and others, 2017). Larsen and others (2016) suggested the acceleration of U1 between 2008 and 2009 was the result of the break-up of its floating ice tongue. Instead, we found the disintegration of the U1 and U2 floating ice tongues occurred after a sustained period of acceleration, including leading up to the partial retreat of the U2 floating ice tongue from 2009 to 2010. We also found deceleration of U1 following the disintegration of its floating ice tongue. This may indicate that the thinning and disintegration of both floating ice tongues was a response to acceleration, rather than its cause.

4.2.2 Impacts of subglacial hydrology on ice flow

Moon and others (2014) categorized the seasonal velocity fluctuations of Greenland's outlet glaciers into three types depending on their sensitivity to ice-front position or meltwater from 2009 to 2013, followed by further classification by Vijay and others (2019, 2021) and Solgaard and others (2022). Our seasonal velocity observations, extending published records through 2021, suggest that Upernavik's outlets are type 2 or 3, which is consistent with previous classifications (Moon and others, 2014; Larsen and others, 2016; Vijay and others, 2019; Solgaard and others, 2022). Glaciers of type 2 and 3 show seasonal velocity fluctuations related to changes in meltwater at the bed and the state of the subglacial hydrologic system. Larsen and others (2023) also linked velocity variations at Upernavik to changes in surface meltwater, and noted that this relationship shows complex variations related to distance from the terminus and response to terminus change.

These studies show that Upernavik's outlets are sensitive to subglacial hydrologic changes. Although seasonal changes are of a much smaller magnitude than observed interannual changes, this may also suggest sensitivity to changes at the bed that operate on interannual scales. This hypothesis is supported by exclusion of other controls through our numerical modeling results. Our experiments indicate that changes in glacier geometry leading to changes in driving stress could not explain the observed velocity changes. Both U1 and U2 experienced thinning over the time period of the model experiments, which should reduce driving stress. However, updating the model geometry did not produce observed velocity slow-downs. While shear-margin strengthening and weakening can also strongly impact ice-flow velocities (Cuffey and Paterson, 2010; Van Der Veen and others, 2011), the pace and direction of observed velocity changes between our model years at Upernavik are not consistent with this mechanism. Crevassing or rifting in the margins could cause rapid loss of sidewall drag and, consequently, velocity increases (e.g. Van der Veen and others, 2011; Lilien and others, 2019), but we know of no mechanism that could cause rapid strengthening in the shear margins to explain the magnitude of the velocity decreases observed at U1 and U2. Even the one modeling experiment that required a velocity increase, the 2018 run for U2, required both reasonable weakening in sidewall drag and a decrease in basal drag in order to match observations. All of this suggests that interannual changes in basal drag are likely to play an important role in changing ice-flow velocities at Upernavik.

We note that there is room for considerable future modeling and observational work to constrain the details of the changes in subglacial hydrology and basal slipperiness. Although we asserted changes in basal slipperiness (Figs 7d, 8d) that fall within realistic ranges, these are just experiments to assess sensitivity, rather than exact simulations.

Our results are also supported by evidence that near-terminus bed conditions can control velocity fluctuations at some similar marine-terminating glaciers. Numerical ice-flow modeling by Rathmann and others (2017) found that the dynamics at 79N Glacier were controlled most strongly by changes in basal slipperiness. This is despite 79N having an ice tongue 76 km long, which was shown to have very little impact on the stress budget at the terminus (Rathmann and others, 2017). In addition, the results produced by our simplified flowline model are consistent with map-view shallow-shelf modeling of Upernavik. Downs and Johnson (2022) showed that U1 and U2 were highly sensitive to changes in basal drag from the terminus to 45 km upstream, and acceleration extending into the upper glacier was observed as a result of these reductions in the basal drag coefficient near the terminus (Downs and Johnson, 2022). This relationship held for both rapidly retreating glaciers and those with stable termini locations. Similarly, we found that the surface velocities throughout U2 were sensitive to reductions in the basal drag coefficient near the terminus, even during rapid retreat. U1, which had a relatively stable terminus location during the modeling study period, displayed similar sensitivity to changes in basal drag.

While our results point toward changes in basal drag as the primary control on the large velocity fluctuations observed at Upernavik, it is worth noting that ice dynamics are often influenced by a combination of forcings. The classification of seasonal velocity fluctuations discussed above represents behavioral endmembers (Moon and others, 2014), and many glaciers are influenced by both the subglacial hydrology and terminus-driven changes throughout the season. Solgaard and others (2022) found that classifications often varied based on distance to the terminus. On Upernavik specifically, Larsen and others (2023) showed that seasonal changes caused by basal hydrologic evolution could be overshadowed near the terminus by impacts caused by terminus retreat. Furthermore, classifications are closely tied to the volume of surface meltwater available (Solgaard and others, 2022), which is controlled by regional climate patterns that may evolve through time. The large variations observed in velocity at Upernavik may therefore be related to variability in interannual weather patterns that are tied to basal conditions through surface meltwater availability. A thorough analysis of atmospheric conditions during the study period is required to assess these relationships in future studies.

4.3 Study limitations

While our results and other studies (Downs and Johnson, 2022) suggest that Upernavik's two largest outlets are highly sensitive to changes at the bed and that changes in basal drag could plausibly explain the observed interannual variability in ice-flow speed, there are other mechanisms that have not been explored in detail. In particular, we did not examine the role of terminus back-stress imparted by land-fast ice and iceberg mélange. As rapidly flowing and actively calving glaciers, both U1 and U2 produce enough icebergs to often have mélange at or near their termini. This material can provide back-stress to glacier termini that can influence both calving rate and ice-flow speed at marine-terminating glaciers (e.g. Walter and others, 2012; Krug and others, 2015). However, to explain the observed large peaks and minima in terminus velocities (Fig. 3), changes in ice mélange would have to be sustained over several seasons, and we saw mostly consistent mélange conditions at the U1 and U2 terminus throughout the record. In addition, the changes typically observed at marine-terminating glaciers in response to ice mélange are of the same magnitude as the seasonal velocity changes we observed at Upernavik (e.g. Walter and others, 2012; Krug and others, 2015), and could not explain the large interannual changes. Finally, if changes in back-stress from the loss of the floating ice tongues did not explain the velocity changes, it is unlikely that mélange, which is typically thinner and less coherent than floating ice tongues, could have had a bigger effect. Nonetheless, a thorough analysis of this hypothesis is required to definitively rule out this mechanism.

As discussed previously, numerical flowline modeling also has many inherent limitations, and we addressed some of these by using a 3-D inversion for basal drag. While using a separate inversion for the basal drag in each year would have better allowed us to determine whether basal slipperiness changed through time, incomplete surface elevations precluded this approach. However, using the 3-D inversion to infer the basal drag prevents some of the non-uniqueness that is introduced by the double-inversion for basal and sidewall drag. The 3-D model resolves across-flow variations in velocity, effectively incorporating the process described by the sidewall drag. Thus, the remaining resistance that is inferred by the sidewall-drag inversion is likely to be assigned to the correct physical process, whereas we would have had no ability to separate the different sources of drag in a purely flowline model.

While this approach helps isolate effects of basal and sidewall drag, it is nevertheless susceptible to errors in the time-varying basal shear stress, which is controlled by the form of the sliding law. Though the drag coefficient likely remains constant, the basal shear stress varies as a function of velocity, ice geometry and other basal properties such as effective pressure. We used a linear Weertman sliding law that linearly relates the basal shear stress to the sliding velocity using a constant coefficient. This differs from the Coulomb-type laws that recent work suggests most accurately describes temporal variations in drag (Joughin and others, 2019; Zoet and Iverson, 2020). However, Coulomb sliding depends on the effective pressure, and thus incorporates an additional unknown parameter. We use the linear Weertman law because it avoids additional uncertain inputs, though it is an imperfect representation of changes in basal shear stress in response to the changing geometry and velocity. To the extent that the temporal variations in basal shear drag are not accurately captured by the simple sliding law used here, there is still some risk that we misinterpret changes in basal resistance as changes at the sidewalls.

In addition, our analysis provides no evidence constraining the nature of the changes at the bed that could lead to the interannual velocity changes. Our remote-sensing data allowed us to rule out changes in floating ice tongues as the cause of the interannual changes. Our flowline model for U1 and U2 enabled us to show that glacier geometry evolution and changes in shear-margin strength were also unlikely explanations. While a flowline model makes significant simplifications to the complex stress state of a glacier like Upernavik, 3-D modeling was prevented by gaps in the input and validation data and was thus out of the scope for this study. Despite its limitations, the flowline model is able to capture some of the key controls on outlet-glacier flow, particularly terminus retreat, thinning and basal and sidewall drag. Based on our elimination of other controls, we suggest that future work constraining the nature of the bed and investigating possible large-scale controls on basal slipperiness is needed at Upernavik.

4.4 Likely future changes at Upernavik's outlets

Acceleration of U1, the fastest-flowing outlet, is possible in the future if the glacier retreats past its shallow bed-ridge into the deeper part of its trough. Velocities at U2 have recently reached the same speed as U1, along with the beginning of new retreat, which indicates the potential for further acceleration of this outlet should it also retreat into the deeper part of its trough. Our results indicate U1 and U2 still have floating ice tongues, making them more sensitive to ocean temperature as a result of the larger area in contact with the ocean (Straneo and Heimbach, 2013). Ocean-driven retreat of both glaciers could cause the ice front to retreat inland by up to 11 km for U1 and 23 km for U2 by 2100 based on model results by Morlighem and others (2019). We recommend further focus on these two outlets, as the depth of their termini could allow continued contact with the warmest water mass in the fjord (Andresen and others, 2014; Muilwijk and others, 2022).

Though U0, U3 and U4 have been decelerating, U3 has the potential to accelerate if it retreats into its long and deep bed. Increased subglacial discharge or thermal forcing could cause this glacier to retreat up to 29 km (Morlighem and others, 2019). There is the potential for continued thinning and retreat of U0 and U4 with increasing atmospheric temperatures as increased subglacial drainage interacting with the terminus could affect their stability (Holland and others, 2008; Rignot and others, 2010, Khan and others, 2013). Further retreat will lead U4 to become land-terminating.

5. Conclusions

Our results show that Upernavik's five outlets have experienced large fluctuations in velocity over the past two decades, accompanied by changes in ice thickness, terminus position and floating ice-tongue presence. As tidewater glaciers are significant contributors to sea-level rise, it is important to understand which forcings control Upernavik's ice dynamics in order to improve predictions of future change.

The extended record of ice-tongue changes developed in this study show that ice-tongue changes are not coincident in time with large fluctuations observed in ice-flow speeds, suggesting that other physical mechanisms must be primary controls on Upernavik's ice dynamics. This conclusion is supported by flowline modeling of the two largest outlets, U1 and U2, which investigates the outlets’ sensitivities to changes in basal drag, sidewall drag, terminus retreat and thinning. Adjusting basal slipperiness yielded modeled velocities that were the most spatially consistent with observed velocities at both outlets. This sensitivity to bed conditions near the terminus is consistent with recent ice-flow modeling of Upernavik as well as other marine-terminating glaciers in Greenland (Rathmann and others, 2017; Downs and Johnson, 2022). Additionally, observations of velocity fluctuations on both seasonal and interannual timescales at Upernavik suggest that Upernavik's outlets are sensitive to meltwater availability at the bed (Moon and others, 2014).

The results highlight the importance of better understanding the basal hydrology and bed conditions at Upernavik Isstrøm. Upernavik's largest outlets have the potential to enter deeper sections of their bed if terminus retreat is sustained (Morlighem and others, 2019), which could further exacerbate ice-mass loss. If Greenland's air temperatures warm as climate change continues, it is likely to cause changes in the meltwater supply and thus ice flow. The contrasting behaviors observed at Upernavik's five outlets suggest that basal conditions can play a large role in tidewater glacier response to climate change, and that predictions of glacier retreat therefore depend strongly on an understanding of ice-base conditions.