1. Introduction
Approximately 17% of the present-day Greenland Ice Sheet is drained via the Northeast Greenland Ice Stream (NEGIS) (Krieger and others, Reference Krieger, Floricioiu and Neckel2019), which extends ∼600 km into the interior, reaching almost to the ice divide. The NEGIS is generally thought to be topographically unconstrained (Christianson and others, Reference Christianson2014; Holschuh and others, Reference Holschuh, Lilien and Christianson2019; Franke and others, Reference Franke2020), as the onset of fast-flowing ice is not confined to, or channelled by, a distinct valley or trough. Similarly unconstrained ice streams in Greenland, such as Petermann Glacier, however, exhibit convergent flow into a main trunk (Chu and others, Reference Chu, Schroeder, Seroussi, Creyts and Bell2018) and lack the clearly defined shear margins that make the NEGIS unique in its geometry.
The complex nature of the NEGIS has given rise to various theories about the onset of its fast flow so far into the interior of Greenland. For example, the NEGIS has been thought to originate as a result of high geothermal heat flux beneath the onset zone (Fahnestock and others, Reference Fahnestock, Abdalati, Joughin, Brozena and Gogineni2001; Rysgaard and others, Reference Rysgaard, Bendtsen, Mortensen and Sejr2018). However, modelling of this potential geothermal hotspot in order to initiate and sustain the NEGIS has provided unrealistically high values of geothermal heat flux (Smith-Johnsen and others, Reference Smith-Johnsen, de Fleurian, Schlegel, Seroussi and Nisanciolgu2020), leading to this theory being disputed (Bons and others, Reference Bons2021). The presence of water-saturated dilatant till has also been proposed to facilitate ice stream flow, suggesting a control through subglacial water routing (Christianson and others, Reference Christianson2014). Other work ascribes the characteristics of the ice stream to the interaction of ice crystal fabric anisotropy and dynamics (Gerber and others, Reference Gerber2023), where the crystal orientation fabric considerably affects the ice stiffness at the onset of the NEGIS. Knowledge of the stability and recent dynamics of the ice stream has also evolved. Originally thought to likely be a steady state feature of the ice sheet (Joughin and others, Reference Joughin, Fahnestock, Macayeal, Bamber and Gogineni2001), the current NEGIS has been shown to be a relatively young, transitory ice stream, where the shear margins have only consolidated 2000 years ago (Jansen and others, Reference Jansen2024), sometime after major deglaciation following the Last Glacial Maximum (Roberts and others, Reference Roberts2024) and subsequent drainage basin reconfiguration throughout the Holocene (Franke and others, Reference Franke2022a).
More widely, multiple factors can potentially influence the location of an ice stream and the onset of its fast flow, for instance, topographic focusing or sudden changes in bed gradient, basal roughness, subglacial geology, geothermal heat flux and subglacial water routing (Winsborrow and others, Reference Winsborrow, Clark and Stokes2010). In the absence of a strong topographic control, ice streams preferentially occur where basal stresses are lower, that is, where subglacial meltwater routing and geology are favourable, after which bed roughness, geothermal heat flux and topographic steps are thought to be the next most important controls (Winsborrow and others, Reference Winsborrow, Clark and Stokes2010). It has also been proposed that a relationship, although not directly causal, exists between faster ice flow and decreased basal roughness in Antarctic ice streams (Rippin and others, Reference Rippin, Vaughan and Corr2011, Reference Rippin2014); however, fast ice flow can also occur over a rough, hard bed (Falcini and others, Reference Falcini, Rippin, Krabbendam and Selby2018).
In the case of the NEGIS, basal roughness evolves downstream, changing from a smooth bed at the onset, to a rougher bed where the ice stream widens, correlating with the slowing of the ice stream and lateral variations in the ice surface velocity (Franke and others, Reference Franke2021). Subglacial water routing has also been established to play a significant role in the ice stream, through the reduction of basal shear stress and facilitation of basal slip (Joughin and others, Reference Joughin, Fahnestock, Macayeal, Bamber and Gogineni2001; Karlsson and Dahl-Jensen, Reference Karlsson and Dahl-Jensen2015; Franke and others, Reference Franke2021), whereas the subglacial geology is poorly constrained. The delineation of subglacial geological boundaries across Greenland has been determined from a synthesis of geophysical datasets by MacGregor and others (Reference MacGregor2024); however, the regions underlying the NEGIS have been highlighted as ones where the subglacial geology is highly uncertain.
Despite the acquisition of large amounts of bed elevation data across Greenland in recent years (MacGregor and others, Reference MacGregor2021), the subglacial topography of the NEGIS has yet to be fully characterised, aside from the more focused studies of Franke and others (Reference Franke2020; Reference Franke2021). In our study, we aim to investigate the subglacial topography along a length of approximately 400 km of the NEGIS from its onset to its divergence into outlet glaciers, in order to ascertain the characteristics of the subglacial conditions over which the NEGIS flows. We examine how the variations in topography beneath the NEGIS are related to ice flow, as well as whether there is an evident geological signal that would impact the subglacial topography. Therefore, we use a geomorphometric approach, utilising the metrics of hypsometry, spatial topographic roughness and valley morphometry, to improve our understanding of how or if topography is exerting a control on the NEGIS, as well as other processes such as geology or hydrology, which may leave topographic signals.
2. Data and methods
In this section, we first describe the airborne radio-echo sounding (RES) data utilised in this study, and the subsequent derivation of ice thickness and bed elevation to produce an updated version of the BedMachine digital elevation model (Morlighem and others, Reference Morlighem2017). We then use these datasets to elucidate the characteristics of the subglacial topography, through hypsometric analysis, spatial topographic roughness and valley morphometrics.
2.1. Airborne RES data
We use the airborne RES data of the EGRIP-NOR-2018 and NEGIS-FLOW-2022 surveys (Franke and others, Reference Franke2019; Carter and others, Reference Carter2023) (Fig. 1), which span from the onset area of the NEGIS to the divergence of the ice stream into its three outlet glaciers: Nioghalvfjerdsfjord Glacier (79NG), Zacharie Isstrøm (ZI) and Storstrømmen Glacier (SG). The EGRIP-NOR-2018 survey covers the EastGRIP ice core site with a line spacing of 5 km, with 10 km spacing further up- and downstream (Franke and others, Reference Franke2020) and flightlines oriented both perpendicular and sub-parallel to ice-flow direction. This survey was subsequently extended by NEGIS-FLOW-2022, which added downstream flightlines at spacings of 6–25 km, to cover a total area of ∼46,000 km2, reaching 200 km upstream and 260 km downstream of the core site.
Figure 1. (a) Surface velocity map (Gardner and others, Reference Gardner, Fahnestock and Scambos2022) of the northeastern section of the Greenland Ice Sheet, highlighting the NEGIS. The shear margins of the ice stream are highlighted, drawn where there is a sharp change in the velocity gradient (white dashes). The flightlines of the AWI surveys (Franke and others, Reference Franke2019; Carter and others, Reference Carter2023) and selected lines from Operation IceBridge surveys (CReSIS, 2024a) (black, red and blue lines) are centred around the EastGRIP ice core site, and span from the onset to where the ice stream diverges into multiple outlets. (b) Inset map shows the surface velocity of the Greenland Ice Sheet, with the outlet glaciers of NEGIS labelled (79NG, ZI, and SG).
These data were collected with the Alfred Wegener Institute (AWI) multi-channel ultra-wideband (UWB) airborne radar sounder, mounted on the Polar6 aircraft (Alfred-Wegener-Institut Helmholtz-Zentrum für Polar- und Meeresforschung, 2016), with a centre frequency of 195 MHz and bandwidth of 30 MHz. A detailed description of AWI’s UWB radar system is given by Rodriguez-Morales and others (Reference Rodriguez-Morales2014); Hale and others (Reference Hale2016); and Arnold and others (Reference Arnold2020), as well as a methodological description of the data acquisition and processing in Franke and others (Reference Franke2020) and Franke and others (Reference Franke2022b). Additional selected RES data, aligned with the flightline direction of the AWI surveys, were obtained from Operation IceBridge (OIB) (CReSIS, 2024a) (Fig. 1) in order to expand the area of analysis within the survey region.
2.2. Derivation of bed topography
The surface and bed reflection of the RES data were delineated using the seismic processing and interpretation software ParadigmTM by Aspen Technology Inc. To convert from two-way travel time to ice thickness, we used an ordinary relative dielectric permittivity of 3.15. Thicknesses were then subtracted from the GIMP DEM surface elevation model (Howat and others, Reference Howat, Negrete and Smith2014) to calculate bed elevation. Surface elevation measurements derived from the aircraft’s laser altimeter were not utilised in order to make the method reproducible and consistent with other datasets. A crossover analysis was performed in order to ascertain the range uncertainty across the bed picks, as described in Franke and others (Reference Franke2020), which takes the root mean square errors of the range resolution and dielectric constant. For our ice thickness values ranging between ∼1000–3300 m, the mean crossover deviation is 12.23 m. The bed elevation point data from NEGIS-FLOW-2022 were integrated into a new version of BedMachine (v6.2, Fig. 2) (Morlighem and others, Reference Morlighem2017), which showed clear differences in the representation of the trough edges as compared to BedMachine v5 (Appendix A, Fig. A1).
Figure 2. (a) Subglacial topography (Morlighem and others, Reference Morlighem2017) beneath the NEGIS, divided into three regions: upstream, middle and downstream.
From this new bed topography dataset, the study area was divided into three regions of similar morphology (Fig. 2). The upstream region encompasses the onset of the ice stream upstream of the EastGRIP ice core site and is characterised by lower elevation topography with flatter relief. Moving downstream, the middle region can be delineated from an upward step in the subglacial topography (∼300 m) situated immediately downstream of EastGRIP. It is characterised by large, broad and shallow subglacial troughs which trend parallel to modern ice flow. The transition from the middle to the downstream region occurs where there is a shift to a more alpine appearance of the topography, where the relief of the topography is higher.
2.3. Hypsometric analysis
In order to analyse the frequency distribution of elevations across the study area, hypsometric curves were calculated from BedMachine v6.2. The hypsometry was derived for each region (Fig. 2), as the percentage distribution of elevations. The distribution of elevations within an area can highlight the dominance of particular processes within a landscape (Jamieson and others, Reference Jamieson2014). For example, the presence of troughs cut to a particular elevation would result in a distinctive peak in the hypsometry, as the flat floor of a trough relates to a large percentage of the bed, which is situated in a particular elevation band.
2.4. Spatial topographic roughness
We use the bed elevation point data from the AWI 2018 and 2022 surveys, as well as selected survey lines from OIB (Fig. 1), for our basal roughness calculations. We evaluate spatial topographic roughness using the standard deviation of vertical change in bed elevation over a given horizontal distance (root mean square height [RMSh]) (Rippin and others, Reference Rippin, Bamber, Siegert, Vaughan and Corr2006). A rougher bed is represented by a higher RMSh, as it implies a greater spread between higher and lower elevations (Falcini and others, Reference Falcini, Rippin, Krabbendam and Selby2018).
For each data point on a linearly detrended profile, we calculated RMSh using a moving window size of 960 m (64 data points), in order to investigate the effects of small-scale (sub-kilometre) topography. The along-track resolution of all flightlines was ∼15 m, meaning that resampling was not required to achieve a constant point spacing. However, if data gaps along the flightlines were present due to missing bed picks, for example, RMSh values were omitted from the analysis as the point spacing was no longer consistent. As the interpretation of roughness parameters is highly directionally dependent (Rippin and others, Reference Rippin2014; Falcini and others, Reference Falcini, Rippin, Krabbendam and Selby2018, Reference Falcini, Krabbendam, Selby and Rippin2021; Eisen and others, Reference Eisen, Winter, Steinhage, Kleiner and Humbert2020), the analysis of the topographic roughness is divided into across- and along-flow profiles.
Anisotropy, or directionality, of topographic roughness quantifies the difference between roughness parallel and orthogonal to ice flow (Smith, Reference Smith2014). This can be applied to bed elevation roughness data at points where along- and across-flow profiles intersect, as described in Falcini and others (Reference Falcini, Krabbendam, Selby and Rippin2021), where the ratio of these roughness parameters was calculated following Smith and others (Reference Smith, Raymond and Scambos2006). The value of the anisotropy ratio is closer to −1 when across-flow roughness is higher than along-flow, is 0 when roughness is isotropic, and is closer to 1 when along-flow roughness is higher than across-flow. Values of anisotropy calculated where the bed has elongated bedforms, such as MSGLs or megagrooves, would therefore tend towards −1, as the bed is being streamlined and reducing the roughness parallel to ice flow. Areas with mixed landforms, for example, drumlins or cnoc-and-lochan terrain, would demonstrate more isotropic values (Falcini and others, Reference Falcini, Krabbendam, Selby and Rippin2021). In the context of ice streams, elongated landforms and smoother surfaces can infer either deformable sediment or a smooth bedrock surface which can be preferentially eroded. Rougher topography, in contrast, may infer a bedrock interface where the erosional patterns are more isotropic in terms of its fracture and joints, and glacial erosion is less able to be reinforced in the same place (Falcini and others, Reference Falcini, Krabbendam, Selby and Rippin2021). The RMSh values at each intersection of two flightlines were used to calculate this anisotropy ratio.
2.5. Valley morphometrics and classification
The morphology of valley cross-sectional profiles was analysed utilising the workflow of Paxman (Reference Paxman2023), which extracted valley profiles from the 2018 and 2022 AWI RES surveys, as well as a wider selection of all OIB survey lines within the study area (Greenland_P3 surveys from 1997, 1999, 2002, 2007, 2010, 2012, 2013, 2014 and 2015) (CReSIS, 2024b). As valley cross-sectional profiles in RES data may be oblique to its ideal orthogonal profile, this effect was corrected for by determining the angular difference between the azimuth of the flightline and the valley strike (perpendicular to valley orientation), as determined from BedMachine v5 (Paxman, Reference Paxman2023). A trigonometric correction was used to adjust the profile as if the survey had been flown orthogonal to the valley orientation, but if the difference between the flightline azimuth and valley strike was >60°, the valley was removed from the analysis. Morphometric indices relating to valley depth, top width, V-index and curvature ratio were calculated for each profile in order to produce quantitative characteristic metrics. Valley cross-profile form is a useful approach in the context of understanding a glaciated landscape, as the morphometry of a valley is characteristic of glacially eroded (U-shaped) or fluvially eroded (V-shaped) valleys, and can be utilised to distinguish between them. These indices were compared using a random forest classifier to a training dataset of fluvial and glacial valleys elsewhere on Earth, from which a classification score was assigned to each valley, ranging from 0 (most akin to fluvial) to 1 (most akin to glacial). Our analysis focused only on valleys with scores of <0.25 or >0.75, as these reflect ‘higher confidence’ classifications. The valleys were also compared in terms of their width/depth (W/D) ratios, in order to differentiate between broad, shallow valleys and relatively narrower and deeper valleys, which may have different mechanisms of formation or different geological conditions (Fig. 3).
Figure 3. Morphometric classification of valley profiles identified from AWI and Operation IceBridge RES surveys. Colour scale relates to the classification of more fluvial (0, blue) to more glacial (1, red). Circle size represents the W/D ratio of each valley profile (see radius given in legend).
2.6. Isostatic rebounding of topography
To ascertain the elevations at which pre-glacial topography would reside, the bed topography was isostatically rebounded relative to a hypothetical ice-free world sea level, calculated as caused by the removal of the modern ice sheets without adjustment for thermosteric or ocean dynamical effects (Paxman and others, Reference Paxman, Austermann and Hollyday2022a). The topography change resulting from ice-sheet unloading was added to BedMachine v6.2 (Paxman and others, Reference Paxman, Austermann and Hollyday2022b), and a second set of (rebounded) hypsometric curves was calculated as the percentage distribution of elevations (Fig. 5), as detailed above. This analysis was carried out to investigate the potential formational processes of the landscape, as the relative elevations of the pre-glacial topography to sea level can give indications of marine influence, or lack thereof.
3. Results
The subglacial topography of the NEGIS can be divided into three geomorphologically distinct regions (upstream, middle and downstream, outlined in Fig. 2), which are characterised by distinct patterns of spatial topographic roughness, hypsometry, and valley morphometry.
3.1. Upstream region
The topographic region upstream of the EastGRIP ice core site (Fig. 2) underlies the onset zone of the NEGIS, where the ice flow velocity begins to accelerate to ∼50 m a−1 within the shear margins of the ice stream. The topography here is low elevation (ranging from −300 to 200 m, Fig. 4a), with the majority lying below sea level, sloping slightly inland. A few isolated topographic high points are present, spanning approximately 10–30 km across, however the overall relief is relatively low and smooth. The narrower range of elevations within the hypsometric distribution (432 m, from 5th to 95th percentile) (Fig. 4a) indicates that there is little topographic variation within this region compared with that further downstream. Smaller peaks within the hypsometric curve can be attributed to a shallow topographic low point and low-relief hills. Few valleys (29 profiles) were identified within this region (Fig. 3), all with shallow depths ranging between 102–187 m, thus suggesting that the bed is not strongly incised near the NEGIS onset.
Figure 4. Hypsometric curves of both modern and isostatically rebounded elevations for the (a) upstream, (b) middle and (c) downstream regions. The isostatically rebounded subglacial topography beneath the NEGIS is calculated as if the full ice-sheet load has been removed (Paxman and others, Reference Paxman, Austermann and Hollyday2022b), the elevations of which are relative to a hypothetical ice-free world sea level.
In the isostatically rebounded topography, which represents how the elevations would change with the ice sheet completely removed, the entire study area rebounds above sea level (Fig. 4). The range of elevations for the upstream region shifts to 400–900 m above sea level, but the hypsometric curve retains a similar shape (i.e. elevations are shifted near-uniformly across the region) (Fig. 4a). After isostatic rebounding, these elevations align more closely with the middle and downstream regions (Fig. 4a–c), with the main peak of the upstream region now coincident with the main peak of the downstream region, centred at ∼600 m above the ice-free world sea level.
In terms of kilometre-scale topographic roughness, in both along- and across-flow directions, the majority (59%) of RMSh values are low, between 0–3 m, with no obvious transition in roughness across the shear margins of the ice stream (Fig. 5). Some local higher roughness values are evident over distinct high points in the landscape. In the across-flow direction, a downstream transition to higher roughness values occurs ∼25 km upstream of a ∼300 m step in the topography at the boundary of the upstream and middle regions. In the along-flow direction, on the western side of the ice stream, the transitional zone can be seen to be more spatially variable. Higher roughness is evident upstream of the topographic step, which decreases moving upstream. The lowest roughness values are contained within the shear margins of the NEGIS (Fig. 5b).
Figure 5. Spatial topographic roughness values (RMSh) calculated using bed elevation point data from AWI and Operation IceBridge RES surveys, separated into (a) across-flow and (b) along-flow flightline orientation. (c) Spatial topographic roughness anisotropy, calculated at the crossing points of perpendicularly oriented flightlines. Values are closer to −1 when across-flow roughness is higher than along-flow, 0 when roughness is isotropic, and closer to 1 when along-flow roughness is higher than across-flow. (d, e) Histograms illustrating the distribution of across- and along-flow RMSh values for the upstream and middle regions.
The anisotropy of the spatial roughness is strongly negative surrounding EastGRIP (Fig. 5c), meaning that the roughness across flow is higher than parallel to flow, as the features at the bed are likely longer than they are wide. Outside of the shear margins, more so on the western side, the anisotropy is more isotropic to positive, with limited anisotropy evident very far upstream.
3.2. Middle region
The transition between the upstream and middle region can be delineated almost immediately downstream of EastGRIP, where there is a step of ∼300 m in the subglacial topography. The elevations within this middle region range from −400 to 400 m, and the hypsometric curve is negatively skewed with three major peaks (Fig. 4b). Once isostatically rebounded, the hypsometry still demonstrates a strong negative skew, however, the main peak resides at a higher elevation (∼750 m) than the main peak of the upstream and downstream regions (Fig. 4b). The region is characterised by two large flow-parallel subglacial troughs, which are evident as the lowest peak in the hypsometric curves of this region (Fig. 4c). The largest of the troughs spans up to 35 km in width and is up to 600 m deep, whilst the smaller trough is ∼400 m in depth, and 30 km wide (Fig. 6). Here, the ice stream is widening and increasing in ice flow velocity up to approx. 130 m a−1, and higher ice velocities are focused on the eastern side of the ice stream.
Figure 6. Across-flow radargrams and velocity profiles progressing downstream in the middle region, illustrating the change in surface ice flow velocity (Gardner and others, Reference Gardner, Fahnestock and Scambos2022) towards the east, coinciding with the large subglacial trough. Contour lines of the ice flow velocity are spaced at 50 m a−1 (black) and 10 m a−1 (grey). Numbers on top of graphs A, B, and C refer to the RES profiles.
The morphology of the valley profiles in this region would indicate a mix of both ‘fluvial’ and ‘glacial’ valleys (Fig. 3). Sixty-nine valley profiles can be identified, however, some of these are consecutive profiles of the same trough, resulting from flightlines crossing multiple times. ‘Glacial’ valleys here mostly have high W/D ratios (from 8 to 80), meaning that they are relatively wide and shallow, with the highest ratios evident within the largest subglacial trough. ‘Fluvial’ valleys tend to have lower W/D ratios and are more randomly distributed across the region.
Spatial roughness within the middle region is more variable, with standard deviations of RMSh values of 6.22 m (across-flow) and 4.65 m (along-flow) as compared to the upstream region (3.74 m along-flow and 2.31 m across-flow) (Fig. 5d, e). A distinct, sharp transition to higher roughness values (Fig. 5) is evident at the topographic step that separates the two regions. However, there is no notable lateral difference in roughness across the shear margins in the middle region, as was also the case for the upstream region. In the along-flow direction, areas of low roughness occur most commonly in the topographic lows, or subglacial troughs, with the higher ridges corresponding with higher roughness (Fig. 5b). Across-flow, the largest trough contains a large range of RMSh values from 0.1 to 58.2 m, however, the adjacent trough has a strong signal of high roughness. Further downstream, towards the transition to the downstream region, there is some differentiation in roughness values across the shear margins that can be seen (Fig. 5a). In terms of the roughness anisotropy, stronger negative anisotropy is evident on the topographic highs in the centre of the ice stream, with low to no anisotropy in the deepest parts of the troughs (Fig. 5c).
3.3. Downstream region
Downstream of the large subglacial troughs, the subglacial geomorphology evolves into higher relief topography (Fig. 4c), with a more alpine appearance evidenced by multiple smaller valleys with complex planform geometries (Fig. 2). The valleys have a range of orientations, in contrast to the middle region where the troughs are aligned parallel to ice flow. At this point, the ice stream diverges into its three outlet glaciers (79NG, ZI, and SG), with a noticeable increase in ice flow velocity immediately downstream of this region.
The hypsometric curve for this region demonstrates a unimodal, near-normal histogram, with elevations ranging between −300 to 1000 m (Fig. 4c). Compared to the middle region, the ranges of elevations in terms of major peaks are similar, and the singular peak coincides with the highest elevation peak in the middle region (Fig. 4b, c). The shape of the hypsometric curve remains the same after adjusting for isostatic rebound, but the main peak now sits 400 m higher at 600 m, and the entire region resides above the ice-free sea level (Fig. 4c).
Many ‘glacial’ valleys can be identified within this region (Fig. 3); however, these have distinctly lower W/D ratios than those in the middle region (ranging between 5–29), meaning that they are relatively narrower and deeper. Some valleys also appear to be channelling the diverging ice stream flow, at the furthest point downstream (Fig. 2). Similarities can also be observed between the valleys identified here and those in the subglacial highlands to the east of the study area, in terms of their highly ‘glacial’ morphometry and comparable W/D ratios.
Analysis of the spatial roughness within the downstream region is limited, due to the orientation of OIB flightlines not being aligned with ice-flow direction. However, the sparse flightlines that are oriented along the along-flow direction show similar roughness values to the middle region, with some small areas of very low roughness evident in the topographic lows (Fig. 5b).
4. Discussion
The characteristics of the three distinct geomorphological regimes beneath the NEGIS can be clearly differentiated through the changes in hypsometry, spatial roughness and valley morphometry that occur moving progressively downstream from the onset zone (Table 1). The potential origins of the differences between these regions, as well as their effects on ice dynamics, are discussed below.
Table 1. Summary of geomorphological characteristics of each region.
4.1. Origins of the regional geomorphological variability
4.1.1. Upstream region
The region surrounding the onset of the NEGIS is characterised by low roughness values (Fig. 5), an absence of valley incision (Fig. 3), and low relief, which would indicate a likely area of soft sediments (Rippin and others, Reference Rippin2014; Siegert and others, Reference Siegert2016). Parts of the central area of Greenland have been identified as having low roughness, with sediments potentially being responsible for the observed smooth topography in places (Rippin, Reference Rippin2013; Cooper and others, Reference Cooper, Jordan, Schroeder, Siegert, Williams and Bamber2019). In the context of Antarctic ice streams, low roughness has been attributed to the presence of unconsolidated sediments, often marine in origin (Bingham and Siegert, Reference Bingham and Siegert2007; Rippin and others, Reference Rippin, Vaughan and Corr2011, Reference Rippin2014). However, it can also indicate a streamlined hard bed, or exposed crystalline bedrock (Jeofry and others, Reference Jeofry2018; Munevar Garcia and others, Reference Munevar Garcia, Miller, Falcini and Stearns2023). Geophysical surveys around the EastGRIP ice core site indicate the presence of saturated sediments (Christianson and others, Reference Christianson2014; Riverman and others, Reference Riverman2019), and sedimentary basins containing MSGLs have been identified beneath the onset zone (Carter and others, Reference Carter2025). The strongly negative anisotropy of the roughness in this region (Fig. 5c) also corroborates the presence of elongated subglacial features beneath the ice stream (Falcini and others, Reference Falcini, Krabbendam, Selby and Rippin2021). Therefore, there are strong indications that the upstream region is underlain by an area of soft sediments, the unconsolidated nature of which suggests a relatively recent origin, rather than a much older Precambrian sedimentary basin.
A possible origin for subglacial sediments is a marine source (Rippin, Reference Rippin2013). Areas situated below sea level could have formerly hosted seaways and been subjected to sediment deposition, as is interpreted to be the case for the region underlying the Siple Coast ice streams in West Antarctica (Studinger and others, Reference Studinger2001). However, when the ice-sheet load is removed from the bed topography of Greenland through isostatic rebounding, the elevations of the interior (including the upstream region described in this study) are significantly above the palaeo-sea level, with the lowest elevations situated at 400 m (Fig. 4a). Even in the case of partial isostatic rebounding, for example, with the presence of an ice cap on the eastern subglacial highlands (Paxman and others, Reference Paxman, Jamieson, Dolan and Bentley2024a), there is insufficient subsidence for the elevations of the sediment-covered region to be situated close to or below sea level. Therefore, we rule out a marine origin for these sediments.
A second potential source of sediments within the interior of Greenland could be glacio-lacustrine, resulting from deposition within large proglacial meltwater lakes. The genesis of the Petermann mega-canyon, which spans the interior of Greenland (Bamber and others, Reference Bamber, Siegert, Griggs, Marshall and Spada2013), has been posited to have been excavated via repeated catastrophic outburst floods, the mechanism of which also necessitates the presence of large proglacial lakes which occupy the interior (Keisling and others, Reference Keisling, Nielsen, Hvidberg, Nuterman and Deconto2020). The modelled evolution of the ice sheet throughout the Pliocene and Pleistocene produces ice-dammed proglacial lakes filling a large overdeepened subaerial basin in central Greenland during phases of deglaciation, with subglacial outlets in both the northwest and northeast (Keisling and others, Reference Keisling, Nielsen, Hvidberg, Nuterman and Deconto2020). Therefore, it could be plausible that sediments from these large proglacial lakes remain in this region, as such lakes interrupt the delivery of meltwater and sediment to proglacial zones, and act as highly efficient sediment traps (Carrivick and Tweed, Reference Carrivick and Tweed2013). Lake sediments have been shown to be preserved beneath the Laurentide Ice Sheet as well as beneath the ice sheet in northwest Greenland (Briner and others, Reference Briner, Axford, Forman, Miller and Wolfe2007; Paxman and others, Reference Paxman, Austermann and Tinto2021).
A third sediment source is glacio-fluvial outwash, which may have been deposited in the interior at a time when erosive mountain valley glaciers were present as part of an ice mass restricted to the southern and eastern highlands, as is inferred to have existed during the late Miocene and/or Pliocene (Bierman and others, Reference Bierman, Shakun, Corbett, Zimmerman and Rood2016; Paxman and others, Reference Paxman, Jamieson, Dolan and Bentley2024a). Some outlets of this modelled ice cap, which is geomorphologically and climatologically constrained (Paxman and others, Reference Paxman, Jamieson, Dolan and Bentley2024b), coincide directly with the upstream region, which generally sits lower than the surrounding topography, with the only possible outlet to the northwest (Fig. 7). Therefore, this region would have been a natural location for sediment being shed from the interior side of the eastern highlands to be accommodated, in front of the margin of an early highland ice cap. A comparable analogue for this mode of sediment deposition would be the Northern Patagonia Icefield, where outlet glaciers drain a large temperate icefield, terminating in proglacial lakes and forming glaciofluvial outwash plains (Glasser and others, Reference Glasser, Harrison and Jansson2009; Loriaux and Casassa, Reference Loriaux and Casassa2013).
Figure 7. (a) Modelled ice cap present on the eastern subglacial highlands (constrained from the configuration of subglacial valley networks) during the past warmer climates (e.g. the late Miocene/Pliocene) (Paxman and others, Reference Paxman, Jamieson, Dolan and Bentley2024b). (b) Potential glaciofluvial sediment outwash pathways from the ice cap into the upstream region of the NEGIS.
4.1.2. Middle region
A distinct transition from the low roughness values (upstream region), to higher and more variable spatial roughness values (middle region), occurs at the topographic step which lies just downstream of the EastGRIP ice core site. This change marks the hypothesised limit of the sediment infill, as the rougher topography implies a discontinuity of the smooth sediments that lie upstream. This resembles findings at Pine Island Glacier and Thwaites Glacier (West Antarctica), where higher roughness relates to outcropping bedrock without significant sediment cover (Rippin and others, Reference Rippin, Vaughan and Corr2011; Schroeder and others, Reference Schroeder, Blankenship, Young, Witus and Anderson2014). At this transitional zone between the upstream and middle regions of the NEGIS, the subglacial landscape comprises a mixed bed assemblage (Carter and others, Reference Carter2025), where low roughness values coincide with sediment, and higher roughness relates to areas of bedrock outcrops, identified as crag and tails (Fig. 8). Areas of low roughness do occur again much further downstream, within the topographic low points of the large subglacial troughs (Fig. 5), which could suggest that sediments have also collected in these basins.
Figure 8. (a) Across- and (b) along-flow RMSh values overlain on high-resolution bed topography data generated from swath radar (Carter and others, Reference Carter2025). Polygons delineate areas of likely sediments and exposed bedrock outcrops.
As well as rougher terrain, the middle region is characterised by the presence of two large subglacial troughs oriented in the direction of ice stream flow. These troughs are likely to be glacial in origin (Fig. 3), with classical U-shaped cross-sectional profiles and high width-to-depth ratios. Their linear structure suggests that these troughs potentially exploit inherited tectonic structures, such as a fault, as a line of geological weakness would act as a conduit for focused ice flow and selective erosion (Patton and others, Reference Patton, Swift, Clark, Livingstone and Cook2016; Paxman and others, Reference Paxman2017). The effect of these troughs on the dynamics of the NEGIS is explored in further detail below.
4.1.3. Downstream region
There are distinct morphological similarities between the valleys in the downstream region of the NEGIS and those identified to the south in the eastern subglacial highlands (Paxman, Reference Paxman2023) (Fig. 3), such as their classification as highly ‘glacial’ in morphology, and comparable W/D ratios. The normally-distributed hypsometric curve is also characteristically alpine (Brocklehurst and Whipple, Reference Brocklehurst and Whipple2004; Jamieson and others, Reference Jamieson2014), which can be clearly differentiated from the other two regions. This suggests that these landscapes had a common origin, meaning that the downstream region could comprise the outer, northernmost continuation of the eastern subglacial highlands, which were incised beneath erosive mountain glaciers during times of restricted ice extent (Paxman and others, Reference Paxman, Jamieson, Dolan and Bentley2024a). Some valleys underlying the NEGIS demonstrate slightly greater W/D ratios than the highlands; this could be attributed to subsequent continued erosion from the ice stream. A geological transition could also demarcate the middle and downstream regions, as the subglacial geology underlying the highlands has been identified as the Caledonian Orogen (MacGregor and others, Reference MacGregor2024), and it would follow that the downstream region beneath the NEGIS represents the limit of this province.
4.2. Subglacial geological provinces
Greenland’s subglacial geology, particularly beneath the interior of the ice sheet, is poorly constrained, and has been mostly derived from interpolation of geological mapping of its ice-free margins. MacGregor and others (Reference MacGregor2024) utilised geophysical data to improve upon previous mapping (Dawes, Reference Dawes2009), and refine the boundaries of its subglacial geological provinces. However, three subglacial regions were left unreconciled with marginal geology, one of which also underlies the NEGIS (Fig. 9a). The upstream and part of the middle region coincide with the unknown ‘Central Region A’, whereas the downstream region and other half of the middle region are proposed to be underlain by the Independence Fjord Basin (Fig. 9a). However, the exact positions of the province boundaries in this area are designated as ‘low confidence’ (MacGregor and others, Reference MacGregor2024). Based on our new findings, we propose adjustments to these geological province boundaries, as the first-order differences in regional geomorphology beneath the NEGIS suggest varying geologies.
Figure 9. (a) Subglacial geological provinces that underlie the NEGIS (MacGregor and others, Reference MacGregor2024), plotted over the subglacial topography (Morlighem and others, Reference Morlighem2017). (b) The red area with a white dashed outline illustrates the proposed adjustment to the subglacial geological province boundaries, which would expand the area designated as the Caledonian Orogen to also underlie the downstream region. The white hashed area shows the region of sedimentary infill, which provides a low-friction subglacial environment.
As highlighted above, the strong similarities of the valley morphology between the downstream region and eastern subglacial highlands may imply a common lithology. The abrupt transition from wider, shallower valleys in the middle region to relatively narrower and deeper valleys in the downstream region likely demarcates the geological boundary, as valley formation and morphology can be controlled by the underlying lithology. For example, in fluvial settings, valley width is two to three times wider in erosion-prone lithologies, such as sedimentary bedrock as compared to relatively erosion-resistant basalt (Schanz and Montgomery, Reference Schanz and Montgomery2016; Langston and Temme, Reference Langston and Temme2019). Therefore, the clear difference in W/D ratio between the middle and downstream region is likely to be attributable to geological differences. Since the Caledonian Orogen underlies the highlands, we suggest that the subglacial province boundary should be modified so that the northern end of the Caledonian Orogen extends inland to incorporate the alpine-like topography of the downstream region of the NEGIS (Fig. 9b). Therefore, we suggest that the boundary between the Caledonian Orogen and Independence Fjord Basin coincides with the transition from the middle to downstream regions. We would also note here that the Independence Fjord sedimentary basin also differs from the unconsolidated sediments present in the upstream region, whilst also being distinct from the Caledonian Orogen, and might explain the presence of much larger subglacial troughs in the middle region due to its potentially greater erodibility.
4.3. Relationship between bed geomorphology and ice dynamics
At the onset of the NEGIS, bed morphology does not appear to have a controlling effect on the location of the shear margins, as the low roughness values do not vary significantly from within to outside of the ice stream in both the along- and across-flow directions. However, spatial topographic roughness beneath the NEGIS evolves downstream, with a distinct increase evident at the transition from the upstream to the middle region (Fig. 5), which could be interpreted to reflect the downstream limit of the sediment cover. This is opposite to what might be expected with increasing ice flow velocity, as roughness often decreases with increasing flow speed (Rippin and others, Reference Rippin, Vaughan and Corr2011, Reference Rippin2014), and low roughness has been attributed to the facilitation of ice streaming and the presence of deformable sediment (Siegert and others, Reference Siegert, Taylor and Payne2005; Rippin, Reference Rippin2013). However, increasing roughness values with increasing ice flow velocity have been observed to occur over a rough, hard bed, such as at the Institute and Möller Ice Streams, and Minch Palaeo-Ice Stream (Falcini and others, Reference Falcini, Rippin, Krabbendam and Selby2018). Fast ice flow, in the case of high roughness, could occur where there is a thick basal temperate ice layer, which would reduce drag (Krabbendam, Reference Krabbendam2016). An alternative explanation could also be an enhanced availability of water at the ice-bed interface, if the sediments in the upstream region act as an aquifer (Li and others, Reference Li, Aitken, Lindsay and Kulessa2022), which terminates or discharges at the boundary between the upstream and middle regions.
Another significant controlling effect on the flow dynamics of the ice stream is the macroscale (10s of km) subglacial topography. To date, the NEGIS has been thought of as topographically unconstrained (Christianson and others, Reference Christianson2014; Holschuh and others, Reference Holschuh, Lilien and Christianson2019; Franke and others, Reference Franke2020), as the onset of the fast flow is not contained within a valley or trough. However, the large subglacial troughs in the middle region exert a clear influence on the flow dynamics of the ice stream, by channelling the fast flow towards the eastern side (Fig. 6). These valleys topographically steer the ice, but this focussing of faster flow occurs within the confines of the shear margins, which are themselves not contained by one all-encompassing valley or topographic ‘low’. Even so, as ice can preferentially flow through the linear valleys in the middle region, the speed of that flow may propagate longitudinally upstream into the basin area, enabled by the low friction characteristics of the sedimentary basin which lies upstream (Fig. 9b). This could potentially explain the narrow band of fast flow experienced all the way inland, and the long, linear nature of the NEGIS, the shear margins of which are reinforced by directional hardening caused by the formation of anisotropic ice crystal orientation fabric (Gerber and others, Reference Gerber2023).
Another notable aspect of the NEGIS is its currently changing dynamic behaviour, as highlighted by its recent (last 35 a) acceleration, which has occurred both in the onset zone and its outlet glaciers (Grinsted and others, Reference Grinsted2022; Khan and others, Reference Khan2022). However, between the upper reaches of the NEGIS and its outlets near the ice margin, there is a distinct zone where surface velocities have not accelerated (Grinsted and others, Reference Grinsted2022). This zone directly overlies the ‘downstream region’ mapped in this study, characterised by alpine-like topography (Fig. 10). Presently, it is unclear whether this is a coincidence or there are unknown causal mechanism(s) linking the spatial variation in bed topography and the dynamical changes of the ice. Nevertheless, the coincidence of the disconnection between the two areas where velocity has been observed to increase in recent decades (i.e. near the outlet and near the onset zone) and the spatial change in the characteristics of bed topography is clearly evident and requires further investigation.
Figure 10. Recent dynamical acceleration of the NEGIS (Grinsted and others, Reference Grinsted2022) overlying the subglacial topography. The disconnect between changes at the outlet glaciers and within the interior is evident as an area of zero acceleration in the downstream region. Values of zero acceleration are transparent in the colour scale.
5. Conclusions
In this study, we have investigated the detailed subglacial topography of the NEGIS, from its onset to where it splits into the three distinct outlet glaciers 79NG, ZI, and SG, to describe and demarcate three geomorphologically distinct regions beneath the ice stream. To ascertain the characteristics of the ice stream bed, we used the topographic analysis techniques of hypsometry, spatial roughness, and valley morphometry. In the upstream region at the onset of the ice stream, the low roughness, low relief, and lack of valleys indicate the presence of sediments. The origin of these sediments is unlikely to be marine incursion, but instead glacio-fluvial outwash from an early ice cap on the eastern highlands (Paxman and others, Reference Paxman, Jamieson, Dolan and Bentley2024a) and/or proglacial lake infill (Keisling and others, Reference Keisling, Nielsen, Hvidberg, Nuterman and Deconto2020). A distinct transition to the middle region, with an abrupt increase in basal roughness, marks the downstream limit of this sediment infill. Two large, flow-parallel subglacial troughs characterise this middle region, which affect the ice stream geometry through topographic steering. Downstream, as the ice stream starts to diverge towards its current termini, the topography evolves again into smaller, alpine-like valleys, which are akin to the eastern subglacial highlands of Greenland. Here, we propose that the changes in valley morphometry are likely to be attributable to mechanical differences in the underlying geological provinces, as the morphological similarities between the downstream region and the adjacent eastern highlands would suggest that this downstream region is part of the Caledonian Orogenic Belt that comprises the highlands. This enables us to place an important constraint on the sub-ice geological structure of a region that was previously poorly understood.
Whilst the upstream regime appears to have little effect on the location of the onset of the ice stream and its shear margins, the downstream subglacial topography has a distinct impact on ice stream geometry, by preferentially steering ice flow through one of the wide troughs. The low friction of the upstream regime may enable the propagation of fast ice flow longitudinally upstream, which could be one part of the explanation for the long, linear nature of the ice stream. On the basis of our data, we argue that the NEGIS is more influenced by basal topography than previously suggested.
Data availability statement
The OPR Toolbox is available at https://gitlab.com/openpolarradar/opr/. UWB radar data of the AWI UWB 2018 radar campaign is available on PANGAEA at https://doi.org/10.1594/PANGAEA.928569, and of the AWI UWB 2022 radar campaign at https://doi.org/10.1594/PANGAEA.986099. Bed picks from both AWI UWB radar campaigns are available on PANGAEA at https://doi.org/10.1594/PANGAEA.907918 (2018) and https://doi.pangaea.de/10.1594/PANGAEA.963555 (2022). Spatial roughness data (https://doi.org/10.1594/PANGAEA.984362) and valley classification results (https://doi.org/10.1594/PANGAEA.984371) are also available on PANGAEA. The upcoming BedMachine v.6 will include the bed elevation datapoints presented here.
Acknowledgements
The authors gratefully thank Mathieu Morlighem for incorporating the new AWI RES data from 2022 into a new version of BedMachine and for providing this dataset for analysis in this study. The authors thank the AWI polar aircraft technicians for their support in the field during the 2022 UWB radar flights as well as the Kenn Borek flight crew. For the NEGIS-FLOW airborne campaign, the authors acknowledge support via the AWI funding grant AWI_PA_02128. The authors thank the entire EastGRIP team during the 2022 field season for scientific and logistical support. EastGRIP is directed and organised by the Centre for Ice and Climate at the Niels Bohr Institute, University of Copenhagen. It is supported by funding agencies and institutions in Denmark (A. P. Møller Foundation, University of Copenhagen), USA (US National Science Foundation, Office of Polar Programs), Germany (Alfred Wegener Institute, Helmholtz Centre 475 for Polar and Marine Research), Japan (National Institute of Polar Research and Arctic Challenge for Sustainability), Norway (University of Bergen and Trond Mohn Foundation), Switzerland (Swiss National Science Foundation), France (French Polar Institute Paul-Emile Victor, Institute for Geosciences and Environmental research), Canada (University of Manitoba) and China (Chinese Academy of Sciences and Beijing Normal University). Furthermore, the authors acknowledge the use of software from Open Polar Radar generated with support from the University of Kansas, NASA grants 80NSSC20K1242 and 80NSSC21K0753, and NSF grants OPP-2027615, OPP-2019719, OPP-1739003, IIS-1838230, RISE-2126503, RISE-2127606, and RISE-2126468. The authors would like to thank Aspen Technology, Inc. for providing software licenses and support.
Author contributions
Charlotte M. Carter wrote the manuscript and conducted the main part of the data processing and analysis. Olaf Eisen and Daniela Jansen were PI and Co-I of the airborne campaigns. Daniela Jansen coordinated and conducted the fieldwork alongside John Paden. Charlotte M. Carter, Steven Franke, Guy J. G. Paxman, Stewart S. R. Jamieson, and Michael J. Bentley interpreted the radar and bed topography data. All co-authors discussed and commented on the manuscript.
Funding statement
Charlotte M. Carter was funded by the Alfred Wegener Institute’s INSPIRES III programme. Steven Franke was funded by the Walter Benjamin Programme of the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation; project number 506043073). Guy J. G. Paxman was supported by a Leverhulme Trust Early Career Fellowship (award number ECF-2021-549).
Appendix A
Figure A1. Difference (v5 minus v6.2) in elevation between BedMachine v5 (Morlighem and others, Reference Morlighem2022) and the updated BedMachine v6.2, which includes the new bed elevation data obtained from the AWI UWB 2022 flightlines (black lines).
Open access
