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, Reference Mouginot2019; Wood and others, Reference Wood2021). 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, Reference Mouginot2019). 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, Reference Nick, Vieli, Howat and Joughin2009; Moon and others, Reference Moon, Joughin and Smith2015). The changing conditions of a floating ice tongue can influence velocities at the terminus, which then propagate upstream (Moon and Joughin, Reference Moon and Joughin2008; Nick and others, Reference Nick, Vieli, Howat and Joughin2009). 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, Reference Joughin, Tulaczyk, MacAyeal and Engelhardt2004; Thomas, Reference Thomas2004; Vieli and Nick, Reference Vieli and Nick2011). 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, Reference Nick, Vieli, Howat and Joughin2009). 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, Reference Zwally2002; Fried and others, Reference Fried2015; Straneo and others, Reference Straneo and Cenedese2015).
One of the largest and fastest-retreating ice streams in northwest Greenland is Upernavik Isstrøm. Model results by Haubner and others (Reference Haubner2018) 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, Reference Mouginot2019).
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, Reference Andersen2019) 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, Reference Gardner, Fahnestock and Scambos2019) 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, Reference Khan2013; Larsen and others, Reference Larsen2016). 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, Reference Andresen, Kjeldsen, Harden, Nørgaard–Pedersen and Kjær2014; Vermassen and others, Reference Vermassen2019; Muilwijk and others, Reference Muilwijk2022). 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, Reference Rignot, Koppes and Velicogna2010; Chauché and others, Reference Chauché2014; Wood and others, Reference Wood2021). 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, Reference Vermassen2019).
The loss of the floating ice tongue at U1 has also been proposed as a possible reason for its acceleration (Khan and others, Reference Khan2013; Larsen and others, Reference Larsen2016). 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, Reference Joughin, Tulaczyk, MacAyeal and Engelhardt2004; Thomas, Reference Thomas2004; Holland and others, Reference Holland, Thomas, De Young, Ribergaard and Lyberth2008; Vieli and Nick, Reference Vieli and Nick2011). Slower or stable velocities are observed when Greenland glaciers terminate with floating ice tongues, including at Jakobshavn Isbræ (Nick and others, Reference Nick, Van der Veen, Vieli and Benn2010; Moon and others, Reference Moon, Joughin, Smith and Howat2012). 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, Reference Kehrl, Joughin, Shean, Floricioiu and Krieger2017). 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, Reference Enderlin, Howat and Vieli2013; Khan and others, Reference Khan2013; Larsen and others, Reference Larsen2016).
Changes in resistive stresses along the glacier bed and fjord walls can also influence velocities (Howat and others, Reference Howat, Joughin, Tulaczyk and Gogineni2005, Reference Howat, Joughin, Fahnestock, Smith and Scambos2008; Pfeffer, Reference Pfeffer2007). 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, Reference Van Der Veen, Plummer and Stearns2011). Joughin and others (Reference Joughin2012) 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, Reference Rathmann2017). Furthermore, changes in basal slipperiness produced realistic seasonal velocities in that model. Recently, Downs and Johnson (Reference Downs and Johnson2022) 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, Reference Downs and Johnson2022).
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, Reference Larsen2016). 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, Reference Downs and Johnson2022). 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, Reference Gardner, Fahnestock and Scambos2019). 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, Reference Andersen2019). 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, Reference Andersen2019).
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, Reference Melton2022). 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, Reference Morlighem2021) 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, Reference Gardner, Fahnestock and Scambos2019). We supplemented these mosaics with more recent data, including averaging data from the ITS_LIVE Global Glacier Velocity Point Data Access (Gardner and others, Reference Gardner, Fahnestock and Scambos2019) and GoLive image pairs (Scambos and others, Reference Scambos, Fahnestock, Moon, Gardner and Klinger2016), 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−1 for the mosaics and 18 to 99 m a−1 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, Reference Joughin, Smith, Howat and Scambos2015). This is available sporadically across the time period of interest with 200 and 500 m resolutions (Joughin and others, Reference Joughin, Smith, Howat and Scambos2015). 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, Reference Joughin, Smith, Howat and Scambos2015). 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−1.
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, Reference Porter2022) 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, Reference Joughin2008; Amundson and others, Reference Amundson2010, Kehrl and others, Reference Kehrl, Joughin, Shean, Floricioiu and Krieger2017; Melton and others, Reference Melton2022). Larsen and others (Reference Larsen2016) 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 (Reference Larsen2016) 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, Reference Wild2022; Bentley and others, Reference Bentley2023). When IceBridge thickness data were available, we conducted a hydrostatic elevation analysis to identify floating termini. For this, we used
where ρ i = 917 kg m−3 is the density of ice, ρ w = 1026 kg m−3 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, Reference Jenkins and Doake1991). 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, Reference Paden, Li, Leuschen, Rodriguez–Morales and Hale2010). 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, Reference Paden, Li, Leuschen, Rodriguez–Morales and Hale2014). 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, Reference Bindschadler2011; Enderlin and others, Reference Enderlin, Howat and Vieli2013; Li and others, Reference Li, Dawson, Chuter and Bamber2022).
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, Reference Blatter1995; Pattyn, Reference Pattyn2003), which allows both flow regimes while having considerably reduced complexity compared to the full Stokes equations (e.g. Shapero and others, Reference Shapero, Badgeley, Hoffman and Joughin2021). Following Shapero and others (Reference Shapero, Badgeley, Hoffman and Joughin2021), the model equations can be written in terms of an action functional, J, that must be minimized
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
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
where g is gravity and s is the surface elevation. The final term in the action is then
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, Reference Shapero, Badgeley, Hoffman and Joughin2021). 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, Reference MacAyeal1993). 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, Reference Morlighem2021). 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, Reference Joughin, Smith and Howat2018). The optimization was done using the Rapid Optimization Library (Ridzal and others, Reference Ridzal, Kouri and von Winckel2017). The inversion used Tikhonov regularization, with strength determined using an L-curve analysis (Calvetti and others, Reference Calvetti, Morigi, Reichel and Sgallari2000). 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, Reference Joughin, Smith and Howat2018) and the Greenland Ice Mapping Project surface used in Bedmachine (~2007; Howat and others, Reference Howat, Negrete and Smith2014; Morlighem and others, Reference Morlighem2021) 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, Reference Choi, Seroussi, Morlighem, Schlegel and Gardner2023), 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, Reference Downs and Johnson2022).
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, Reference Shapero, Badgeley, Hoffman and Joughin2021). 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, Reference Morlighem2021; Porter and others, Reference Porter2022). 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, Reference Langen2015). 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, Reference MacGregor2016) 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
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−5 K kPa−1 (e.g. Cuffey and Paterson, Reference Cuffey and Paterson2010).
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−3. 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, Reference Gagliardini, Durand, Zwinger, Hindmarsh and Le Meur2010), 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).
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, Reference Andersen2019). 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, Reference Gardner, Fahnestock and Scambos2019; Scambos and others, Reference Scambos, Fahnestock, Moon, Gardner and Klinger2016).
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, Reference Nick, Vieli, Howat and Joughin2009), 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.
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, Reference Joughin, Smith, Howat and Scambos2015; Gardner and others, Reference Gardner, Fahnestock and Scambos2019).
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.
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, Reference Joughin, Smith, Howat and Scambos2015; Gardner and others, Reference Gardner, Fahnestock and Scambos2019), 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 (Reference Khan2013) and Larsen and others (Reference Larsen2016) 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).
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, Reference Paden, Li, Leuschen, Rodriguez–Morales and Hale2010; Paden and others, Reference Paden, Li, Leuschen, Rodriguez–Morales and Hale2014).
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, Reference Morlighem2021). 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−1 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−1. Estimated ArcticDEM error based on nearby rock outcrop elevations was consistently under 3 m between 2011 and 2021.
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, Reference Morlighem2021) and ArcticDEM from 2011 to 2021 (Porter and others, Reference Porter2022). 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, Reference Wild2022; Bentley and others, Reference Bentley2023). 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−1 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.
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.
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−1 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 (Reference Melton2022) and Amundson and others (Reference Amundson2010) 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 (Reference Melton2022) 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, Reference Joughin, Tulaczyk, MacAyeal and Engelhardt2004; Thomas, Reference Thomas2004; Howat and others, Reference Howat, Joughin, Tulaczyk and Gogineni2005, Reference Howat, Joughin, Fahnestock, Smith and Scambos2008), 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, Reference Rathmann2017). Larsen and others (Reference Larsen2016) 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 (Reference Moon2014) 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 (Reference Vijay2019, Reference Vijay2021) and Solgaard and others (Reference Solgaard, Rapp, Noël and Hvidberg2022). 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, Reference Moon2014; Larsen and others, Reference Larsen2016; Vijay and others, Reference Vijay2019; Solgaard and others, Reference Solgaard, Rapp, Noël and Hvidberg2022). 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 (Reference Larsen2023) 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, Reference Cuffey and Paterson2010; Van Der Veen and others, Reference Van Der Veen, Plummer and Stearns2011), 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, Reference Van Der Veen, Plummer and Stearns2011; Lilien and others, Reference Lilien, Joughin, Smith and Gourmelen2019), 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 (Reference Rathmann2017) 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, Reference Rathmann2017). In addition, the results produced by our simplified flowline model are consistent with map-view shallow-shelf modeling of Upernavik. Downs and Johnson (Reference Downs and Johnson2022) 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, Reference Downs and Johnson2022). 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, Reference Moon2014), and many glaciers are influenced by both the subglacial hydrology and terminus-driven changes throughout the season. Solgaard and others (Reference Solgaard, Rapp, Noël and Hvidberg2022) found that classifications often varied based on distance to the terminus. On Upernavik specifically, Larsen and others (Reference Larsen2023) 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, Reference Solgaard, Rapp, Noël and Hvidberg2022), 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, Reference Downs and Johnson2022) 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, Reference Walter2012; Krug and others, Reference Krug, Durand, Gagliardini and Weiss2015). 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, Reference Walter2012; Krug and others, Reference Krug, Durand, Gagliardini and Weiss2015), 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, Reference Joughin, Smith and Schoof2019; Zoet and Iverson, Reference Zoet and Iverson2020). 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, Reference Straneo and Heimbach2013). 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 (Reference Morlighem, Wood, Seroussi, Choi and Rignot2019). 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, Reference Andresen, Kjeldsen, Harden, Nørgaard–Pedersen and Kjær2014; Muilwijk and others, Reference Muilwijk2022).
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, Reference Morlighem, Wood, Seroussi, Choi and Rignot2019). 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, Reference Holland, Thomas, De Young, Ribergaard and Lyberth2008; Rignot and others, Reference Rignot, Koppes and Velicogna2010, Khan and others, Reference Khan2013). 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, Reference Rathmann2017; Downs and Johnson, Reference Downs and Johnson2022). 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, Reference Moon2014).
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, Reference Morlighem, Wood, Seroussi, Choi and Rignot2019), 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.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/aog.2023.76.
Acknowledgements
This work was supported by funding from the Natural Sciences and Engineering Research Council of Canada (NSERC) grants: the Canada Excellence Research Chair CERC-2018-00002 and the Discovery Grant RGPIN-2021-02910. Additional financial support from the University of Manitoba Graduate Fellowship (UMGF) made this study possible. We are very grateful for the support from the community of the Centre for Earth Observation Science (CEOS) at the University of Manitoba.
Author contributions
Kelsey Marie Voss performed all data analysis and wrote most of the paper. Karen E. Alley oversaw the data analysis and contributed to writing the paper. David A. Lilien performed all model inversions and aided in the model set-up and simulations while also contributing to writing the paper. Dorthe Dahl-Jensen contributed to the subglacial hydrology discussion and overseeing the study logistics.
Open access
