Ice acceleration and rotation in the Greenland Ice Sheet interior in recent decades

The study investigates whether presumedly stable interior regions of the Greenland Ice Sheet (GrIS) have experienced multi-decadal changes in flow speed and direction. Constraints on long-term interior dynamics are sparse due to limited high-quality pre-GPS observations. Many mass-balance assessments assume negligible ice-dynamic change in the high-elevation interior, and satellite-derived velocities struggle to detect summer accelerations above ~1250 m in central West Greenland due to poor signal-to-noise. Meanwhile, GrIS mass loss has accelerated via enhanced surface melt and dynamic discharge from marine-terminating glaciers like Jakobshavn Isbræ, which underwent floating tongue disintegration starting in 1997, followed by retreat, acceleration, and thinning. This raises the question of how far inland dynamic perturbations propagate and whether interior flow has changed. Leveraging historical EGIG (1959–1967) theodolite measurements and PARCA (1993–1995) GPS observations, the authors reanalyze and resurvey selected sites upstream of Jakobshavn Isbræ to test for inland trends in velocity and azimuth, assess potential drivers, and evaluate implications for ice-sheet dynamics and catchment geometry.

Prior work highlights: (1) Many GrIS mass-balance studies assume stable interior dynamics, with limited constraints from pre-GPS observations. (2) Satellite-derived velocities can have difficulty detecting inland changes at high elevations due to noise limitations. (3) Jakobshavn Isbræ experienced dramatic changes following disintegration of its floating tongue in the late 1990s, including retreat, acceleration, and thinning; similar responses are observed at other marine-terminating glaciers. (4) PARCA provided 1990s GPS velocities along the 2000 m contour, including sites overlapping this study. (5) Inland thinning can propagate as a diffusive kinematic wave; studies such as Felikson et al. assessed inland bounds for thinning propagation based on outlet geometry. (6) Previous inland accelerations have been attributed to different mechanisms: DS14 (site S10) inferred enhanced basal sliding transmitted via longitudinal coupling with distinct seasonal signals; WS21 found accelerations near Jakobshavn with no seasonality and attributed them to changes in local driving stress from geometry. (7) Studies on crevasse mechanics indicate local stress-field reorganization near crevasses and potential decoupling between surface and basal flow directions; borehole measurements show depth-varying azimuth related to ice properties and bed geometry. (8) Thermomechanical modeling and borehole data suggest temperate basal layers exist in parts of the Jakobshavn catchment, potentially influencing deformation and velocity.

Historical reanalysis: The authors reanalyzed 14 EGIG stakes (closest to Jakobshavn Isbræ) originally surveyed in 1959 and re-measured in 1967. To ensure comparability, 1959–1967 positions were converted from the Hayford ellipsoid to WGS84; projected coordinates used EPSG:3413. Azimuths were recalculated using spherical trigonometry (eq. 2), and velocities from position changes over time (eq. 3). Uncertainties were estimated with a Monte Carlo approach perturbing positions ±11 m (uniformly) 1e6 times, taking 95% CI half range as δv (eq. 4). All historical values were reproduced except T1 due to position inconsistencies; thus T1 historical values were excluded from inter-comparisons. Modern resurvey: From 2020–2022, 11 GPS stations were deployed at/near historical positions: EGIG T1–T5, T127a, T127–T129, and PARCA cd08, cd38, spanning ~1800 km². Position logs at 30 s sampling with Javad GPS receivers/antennas on elevated scaffolds; data processed with GIPSY-OASIS using JPL final products, VMF1 troposphere mapping, FES2014b ocean tide loading; coordinates in EPSG:3413. Observation durations varied from ~3.6 to >400 days; 6 sites had coverage >100 days (T3, T5, T127, T128, cd38, cd08). Velocity components were determined from linear fits to position–time series, with magnitude from component combination; azimuths via eq. (2). Velocity uncertainties were computed from residual-based fractional uncertainties (eq. 5). Present-day azimuth uncertainty was set to ±1°. Plate motion (13–16 mm/yr) was negligible relative to ice motion and not corrected. Spatial comparability: New observations were placed as close as possible to 1959 sites. To minimize biases from spatial gradients, comparisons were restricted where new vs historical points were within ~1 ice thickness (H) based on BedMachine ice thickness; some historical positions advected to 1967–1991 were farther than this threshold and were excluded from direct comparison. Altimetry and geometry change: Multi-mission altimetry datasets provided annual rates of surface elevation change for 1995–2011 and 2011–2020. Cumulative changes (1995–2020) were calculated to assess thinning and surface slope changes. Surface slopes were computed over three horizontal scales centered on each site: ±2H (~7 km), ±5H, and ±15H (~54 km) to bound local vs regional slope variability. The authors identified a knickpoint—transition from flat to steep surface—along a flowline through T1–T5 using a 2018/2019 DEM and altimetry time series to infer inland migration and estimate kinematic wave speeds. Dynamic attribution analyses: Seasonality was assessed by comparing normalized daily velocities for summer vs non-summer periods via unpaired two-tailed t-tests (1% significance) for four summer definitions; no significant seasonal differences were detected. Basal sliding influence through longitudinal stress-gradient coupling was evaluated using the Van der Veen formulation along a flowline near T4; coupling stresses were compared to gravitational driving stress, testing a hypothetical 10 m/yr sliding perturbation 1 km downstream of T4. The magnitude of observed strain-rate changes between T3 and T4 was compared to that required to explain the acceleration via coupling alone. Sensitivity of velocity to geometry: Using a simple-shear deformational model (η=3), surface velocity u = ub + (2A/(η+1))(ρ g H sin S)^(1/η) H, with ub=0, the authors explored theoretical percentage velocity changes for δH (±3%) and δS (±9%). Observed δH and δS at each site were derived from altimetry-updated 2018/2019 DEMs and compared against observed accelerations across the three slope scales. Creep instability case study: To explore whether warming/softening of near-bed ice could produce observed accelerations, a best-estimate temperature profile at T4 was extracted from Aschwanden et al.'s thermomechanical simulation (PISM-based). A temperate basal layer of 26 m (~1.4% of 1882 m thickness) was inferred initially. Assuming an enhancement factor E=2.5 to match 1959/1967 velocities, the temperate-layer height was manually increased by 18 m (to 43 m, 2.3% of thickness) and the resulting deformational velocity profile recalculated, showing a ~9.2% surface velocity increase comparable to observations. The authors note uncertainties due to smoothed bed topography in simulations and absence of in situ temperature profiles at the sites.

- Modern in situ GPS measurements (2020–2022) at 11 sites upstream of Jakobshavn Isbræ show increased velocities relative to historical EGIG (1959–1967) and PARCA (1993–1995) data. Increases are 6–19 m/yr, corresponding to 5.5–14.6% across all stations. The smallest percentage increases occur furthest inland (e.g., 5.4% at cd08; 8.3% at T5); the largest at T129 (14.6%). - EGIG sites exhibit consistent azimuth increases of 3–4.5°, indicating a northward deflection of flow beyond uncertainties. PARCA sites cd08 and cd38 show no statistically significant long-term azimuth change (<0.5° within uncertainties). - The majority of the regional velocity perturbation appears to have occurred post-1995 (PARCA period), as comparable magnitude increases over shorter intervals are observed at PARCA-overlap sites. - No seasonal signal in velocity was detected at any site between 2020 and 2022 (no significant summer vs non-summer differences). Nearby GC-NET PDDs (12–14 days) are far below levels associated with seasonal basal sliding elsewhere (40–120 PDD at Swiss Camp in late 1990s). - Longitudinal coupling stresses constitute ~6% of gravitational driving stress under no-sliding assumptions; a 10 m/yr basal sliding perturbation 1 km downstream of T4 would be required to raise total driving stress by 10%, inconsistent with observed along-flow strain-rate changes (Δε ≈ 8.2×10⁻¹¹ yr⁻¹), implying coupling alone cannot explain the acceleration. - Altimetry (1995–2020) indicates thinning of ~2.6–8.2 m at sites (0.1–0.4% of thickness) and regional steepening. However, geometry-induced changes in driving stress are generally insufficient to explain observed accelerations across most sites and slope scales; only at some scales for T2 and cd38 could geometry changes approach observed speeds. - A knickpoint (transition to steeper surface) migrated inland by ~82.2 km along the T1–T5 flowline between 2005 and 2019, implying a local kinematic wave speed of ~5.9 km/yr (range 4.5–9.1 km/yr accounting for ±5-year timing uncertainty), much slower than alpine down-glacier kinematic waves. - New, large transverse crevasses (~25 m wide) have appeared within the study area since the 1960s, bisecting historical traverse routes, suggesting a shift in local ice dynamics that may affect stress fields and potentially flow direction. - A simple case study at T4 shows that increasing the temperate basal layer height by 18 m (from 26 m to 43 m; 1.4% to 2.3% of thickness) can produce a ~9.2% surface velocity increase (from ~103.9 to ~113.2 m/yr), consistent with observed acceleration, supporting a creep-instability mechanism.

The findings demonstrate that significant dynamic changes—both acceleration and rotation—have occurred >100 km inland in the Jakobshavn Isbræ catchment, challenging the assumption of stable interior flow. The absence of seasonal velocity signals and limited PDDs, combined with the lack of supraglacial hydrologic activity and the distance from the terminus, argues against a dominant role for changes in local basal hydrology and sliding at the sites. Analytical estimates indicate longitudinal stress-gradient coupling from enhanced sliding downstream is insufficient to explain the observed accelerations, given the small coupling lengths (≈4–10H) relative to the spatial extent of observed changes and the measured strain-rate differences. Geometry-driven changes (thinning and slope steepening) contribute to driving stress but are generally too small—across multiple length scales—to account for the observed accelerations at most sites, in contrast to some previous interpretations. The observed inland migration of a surface knickpoint suggests a diffusive kinematic thinning wave moving up-glacier since the mid-2000s, which may reorganize stresses and facilitate dynamic responses inland. The consistent northward azimuth shifts at EGIG sites, coincident with acceleration, and the emergence of large transverse crevasses, point to a fundamental change in local ice rheology and stress regime. The authors posit a limited creep instability—warming and softening confined to the lower ~15–20% of the ice column—as a plausible mechanism capable of enhancing internal deformation and modifying depth-varying flow direction, thereby explaining both acceleration and azimuth rotation without large changes in local driving stress. A sensitivity experiment at T4 indicates that modest increases in temperate basal layer thickness can produce observed surface speedups. Such a mechanism is consistent with inland responses to terminus perturbations and suggests that the deep interior may be more dynamically sensitive than previously recognized. The observed azimuth shifts near the catchment boundary hint that Jakobshavn’s catchment may be evolving inland following its floating tongue collapse, with implications for flux-gate assumptions in mass-budget studies.

By resurveying historical EGIG and PARCA sites with in situ GPS, the study reveals a 5–15% acceleration and 3–4.5° northward rotation of ice flow at inland sites (>100 km from the terminus) upstream of Jakobshavn Isbræ. These changes likely initiated after ~2005 and have propagated inland as indicated by a knickpoint migration with a kinematic wave speed of ~5.9 km/yr (range 4.5–9.1 km/yr). Contrary to some previous interpretations, local geometry changes and longitudinal coupling alone do not fully explain the observed accelerations, and there is no evidence for seasonal basal hydrology-driven speedups at these elevations. The authors propose a limited creep-instability mechanism—warming and softening of near-bed ice layers—as a process capable of explaining both acceleration and azimuth rotation. This suggests that dynamic perturbations at Jakobshavn’s terminus can influence non-channelized interior flow far inland, and that the interior may be more dynamically responsive than assumed. The results underscore the importance of contemporary in situ GPS observations; reliance on 1990s PARCA velocities as proxies for present-day conditions can bias mass-budget and gravity analyses that assume constant inland discharge. Future work should obtain borehole temperature profiles, refine bed and thermal models to quantify temperate-layer evolution, and expand GPS networks to resolve spatial patterns of acceleration and rotation across catchment boundaries.

- Spatial comparability is limited where new and historical observation points are separated by >1 ice thickness, potentially introducing spatial gradients; some historical velocities (e.g., 1967–1991) could not be compared for this reason. - Observation durations varied widely (3.6 to >400 days) with non-continuous records due to logistics; 5 of 11 sites had <100 days coverage, which may affect precision of velocity and azimuth estimates. - Historical azimuth for PARCA sites could not be recalculated due to undocumented methods, limiting interpretation of azimuth changes at cd08 and cd38. - Catchment boundaries are uncertain, complicating attribution of dynamics to specific basins. - Geometry-change attribution depends on altimetry products, DEM choice, and slope length scales (2H, 5H, 15H); local variability at short scales affects inferred driving-stress changes. - No in situ ice temperature profiles exist at the surveyed sites; the creep-instability case study relies on a modeled temperature profile and simplified deformation assumptions (simple shear, E=2.5), and smoothed bed topography likely underestimates temperate-layer thickness. - The linkage between observed knickpoint migration and direct transmission of a front-originating thinning signal cannot be definitively established. - The impact radius of crevasse-related fracture mechanics on azimuth is poorly constrained; crevasse presence/visibility may also be influenced by firn property changes. - T1 historical position inconsistencies prevented inclusion of its historical values in inter-survey comparisons. - Plate motion was deemed negligible and not corrected, though this introduces a minor unmodeled component (~13–16 mm/yr).