Formation of mega-scale glacial lineations far inland beneath the onset of the Northeast Greenland Ice Stream

Rapidly-flowing ice streams drain the interior of the Greenland Ice Sheet, currently accounting for around half of its annual mass loss. The Northeast Greenland Ice Stream (NEGIS) is one of the largest, recognisable almost 600 km inland, and extends close to the central ice divide. Numerical ice sheet models are unable to accurately reproduce the configuration of the NEGIS, but understanding its bed properties and spatial and temporal evolution is critical to predicting its future contribution to sea-level change. Here, we use swath radar imaging to create a high-resolution Digital Elevation Model of the bed close to where the NEGIS initiates. Surprisingly, this reveals a landscape interpreted to include mega-scale glacial lineations (MSGLs) that are often assumed to be indicative of rapid ice stream flow (100s m yr⁻¹), under present-day flow velocities of only ∼ 60 m yr⁻¹. Given that MSGLs are thought to form under much higher flow velocities, their presence so far inland at an onset zone raises important questions about their formation and preservation under ice streams, as well as past configurations of the NEGIS. Elongate bedrock landforms outside the current shear margins also suggest that the NEGIS was wider than its present configuration at some point in the past.

1 Introduction

The Greenland Ice Sheet (GrIS) is currently the largest contributor to global mean sea level rise from the polar ice sheets (Otosaka et al., 2023). Its mass loss has increased sixfold over the last 40 years (Mouginot et al., 2019), around 50 % of which is accounted for by dynamic discharge via ice streams (Shepherd et al., 2020). One of the most prominent is the Northeast Greenland Ice Stream (NEGIS), a 600 km long feature draining ∼ 17 % of the GrIS (Krieger et al., 2019) via fast-flowing marine-terminating glaciers (Fig. 1a). The ice stream is characterised by several tributaries, the longest of which sharply narrows as it extends inland towards the central ice divide, and is also unusual in that there is no obvious topographic steering (Holschuh et al., 2019; Franke et al., 2020), meaning that the onset is not contained within a valley. The onsets of other ice streams in Greenland that are topographically unconstrained, such as Petermann Glacier (Chu et al., 2018), are often much wider and exhibit convergent flow into a main trunk, and lack such clearly defined shear margins, making the NEGIS unique in its geometry. Rather, the origins of the NEGIS have previously been attributed to subglacial water produced by exceptionally high geothermal heat fluxes in its onset zone (Fahnestock et al., 2001; Rysgaard et al., 2018; Smith-Johnsen et al., 2020), although this is disputed (Bons et al., 2021). Despite its importance, numerical ice sheet models are unable to fully replicate the NEGIS's unusual flow configuration (Aschwanden et al., 2016). This discrepancy between numerical ice sheet models and modern-day ice flow velocities is thought to be partly due to the missing parameter of ice crystallographic preferred orientation within the models (Gerber et al., 2023; Stoll et al., 2025), which can have an effect on the inferred basal drag or friction coefficient (Rathmann and Lilien, 2022), but may also be related to the representation of basal conditions.

Figure 1 Overview of the study site in northeast Greenland. (a) NEGIS ice flow velocity (Gardner et al., 2022) with survey flightlines centred over the EastGRIP ice core site (red dot). Inset map shows ice flow velocity across Greenland on a logarithmic colour scale. (b) Survey flightlines in swath mode at the onset of NEGIS with an outline (light grey) of the generated DEM, and ice flow velocity in the background. Black dots show the locations of seismic data points from Christianson et al. (2014).

Whilst the modern flow configuration and surface velocity of the ice stream is well constrained by satellite observations (Joughin et al., 2018) (Fig. 1a), much less is known about the bed properties and the spatial and temporal evolution of the NEGIS, the understanding of which is crucial to elucidating the driving factors which influence the ice stream. Basal properties of the onset of the NEGIS have been constrained using single seismic lines, which were used to infer the, at least local, presence of saturated, high-porosity, deformable subglacial sediment (Fig. 1b) (Christianson et al., 2014; Riverman et al., 2019). The presence of a streamlined bed, with potential elongate subglacial landforms in the upper reaches of the NEGIS, has also been hinted at indirectly via off-nadir reflections from widely spaced radar lines (Franke et al., 2020). In terms of the evolution of the ice stream, a recent study has revealed that the shear margins and flow configuration in the onset zone of NEGIS were only established from around 2000 years ago (Jansen et al., 2024), questioning the assumption that the ice stream configuration had been stable throughout the Holocene (Fahnestock et al., 2001). Here, the amplitudes of folded internal reflection horizons from radio-echo sounding data indicate a record of the deformation of the ice stream, where the shear strain rates localised in narrow shear margins in the upstream part of the NEGIS ∼ 2000 years ago (Jansen et al., 2024). Indeed, evidence from internal stratigraphy has also shown rapid reconfigurations in ice flow trajectories and drainage basin areas in northeast Greenland over short (centennial) timescales during the Holocene (Franke et al., 2022a).

In order to further investigate the basal properties that lie beneath the onset of the NEGIS, this paper presents the first high resolution (25 m horizontal resolution, ∼ 10 m vertical resolution) Digital Elevation Model (DEM) of the subglacial topography, reconstructed from swath radar, which is capable of identifying individual subglacial bedforms. The DEM covers a 40 × 60 km² area surrounding the EastGRIP ice core site (Fig. 1b), where the ice thickness reaches 2.5 km and ice flow velocities are around 60 m yr⁻¹. Previously, mapping of in situ subglacial bedforms under other active ice streams at such a high resolution has been achieved with conventional nadir-focused radar techniques, which derive bed elevation from closely-spaced radar lines and use interpolation (King et al., 2016; Bingham et al., 2017; Schlegel et al., 2022). These methods, whilst providing excellent observations, have an inherent limitation in their cross-track resolution and spatial coverage, in that direct point measurements of the bed are only possible in the along-track direction. Instead, with swath radar processing, there is no requirement to interpolate between these point measurements, as the direction of the arrival of energy in both the along-track and across-track directions can be estimated from the sequential acquisitions. Therefore, the topography can be reconstructed in high resolution over kilometre-wide strips, and the imaging of subglacial bedforms is possible at a lower logistical cost (Holschuh et al., 2020) and over a much wider spatial area.

Importantly, swath radar data surveys enable a view into the small-scale (10s to 100s of metres) subglacial topography and landforms that reside beneath active ice sheets. This has a range of advantages, from identifying geomorphological features, to enabling the incorporation of small-scale variations in basal topography within ice flow models, which can produce modelled ice motion patterns that diverge substantially from those utilising smoothly varying topography (Law et al., 2023). Two other similar surveys of swath radar data have been generated in Antarctica (Holschuh et al., 2020; Hoffman et al., 2023), which investigate the subglacial geomorphology beneath Thwaites Glacier and Hercules Dome. These data, alongside our DEM, provide crucial insight into the formation of subglacial landforms, which previously have generally only been accessible and studied within the context of deglaciated landscapes. For example, when elongate subglacial landforms are identified on a deglaciated landscape, they are utilised as a proxy to reconstruct previous ice sheet extent and dynamics (Margold et al., 2015a; Stokes et al., 2015). Mega-scale glacial lineations (MSGLs), for example, are generally associated with high flow velocities under ice streams (King et al., 2009; Stokes et al., 2013; Spagnolo et al., 2014). In addition, the elongation of these landforms has often been said to follow a down-flow morphological continuum, where they evolve from flow-perpendicular landforms at the onset (for example, ribbed moraines) to elongated flow-parallel MSGLs, as a product of increasing ice flow velocity and changing basal conditions, such as cumulative strain at the bed (Ely et al., 2016; Zoet et al., 2021; Vérité et al., 2021; Vérité et al., 2023). Similar transitions in evolution are also observed laterally, such as when crossing from the centre of an ice stream to its shear margin and beyond (Stokes and Clark, 2002a; Vérité et al., 2021). However, the inferred link between MSGLs, ice streaming, and the duration of their deformation related to cumulative strain at the base of the ice is largely qualitative with very few quantitative modelling studies (Jamieson et al., 2016; Ely et al., 2022; Vérité et al., 2024).

In addition, different formational mechanisms for MSGLs have been proposed and there is, as of yet, little consensus in the literature. Early ideas tended to focus on their construction via processes of subglacial deformation under high ice velocities (a velocity-duration product) (Clark, 1993). Others have proposed an erosional mechanism through the catastrophic discharge of turbulent subglacial meltwater (Shaw et al., 2008) or via the erosion of pre-existing sediments to reveal a streamlined surface of megaridges (Eyles et al., 2016). Indeed, some have focussed on the grooves between the ridges and suggested that MSGLs may be the erosional remnants or intervening ridges from a “groove-ploughing” mechanism, where roughness in the ice base passes over a bed of soft, saturated sediments (Clark et al., 2003). A more recent hypothesis invokes a rilling instability (Fowler, 2010), where the water flowing between ice and deformable subglacial sediment is unstable, and will form linear streams/rills separated by intervening ridges. Related to this hypothesis is the instability theory of the coupled flow of ice, water and till, which has successfully modelled the formation of a continuum of subglacial bedforms, including MSGLs (Ely et al., 2022). Although most work on MSGLs has utilised observations from the palaeo-record, geophysical data of their presence under extant ice sheets (e.g. King et al., 2009) offers a powerful tool to test some of these various formational hypotheses.

The data presented here comprise the first high-resolution subglacial topography survey beneath the onset of an ice stream in the deep interior of the GrIS, located around 600 km from the grounding line. The classification of the now visible geomorphological features provides insight into both the basal conditions of the ice stream and its spatial evolution. The identification of MSGLs beneath the onset also raises important questions around their formation and preservation beneath ice streams, given the relatively low ice flow velocities.

Figure 2 (a) Individual waveforms from a multi-antenna radar array used to image both nadir and off-nadir return echoes from the ice-bed interface. (b) Time delays for the same phase at single antenna elements are used to geolocate the off-nadir return echoes.

2 Methods

2.1 Generation of swath data

The Alfred Wegener Institute's ultra-wideband radar system (AWI-UWB) (Hale et al., 2016; Franke et al., 2022b) was used to map the cross-track topography at a horizontal resolution of 25 m, in swaths ∼ 2 km wide (Fig. 2). The majority of the flightlines were oriented perpendicular to the direction of ice flow, in order to reduce the possibility of artefacts that could be misinterpreted as flow-parallel subglacial landforms (Fig. 1b). Historically, the ability to resolve the detail of subglacial topography from airborne radar surveys has been dictated by along-track spacing, survey line spacing and along-track focusing, as conventional bistatic antenna airborne radar is unable to localise energy in the cross-track direction. Therefore, to produce a high-resolution digital elevation model of the bed, we apply swath processing techniques to ultra-wideband radar data collected by the AWI-UWB radar system in 2018 and 2022, surrounding the EastGRIP ice core site. Using a multi-element antenna array (Fig. 2a) means that information from sequential acquisitions, combined with phase differences in arrivals between receiver elements (Jezek et al., 2011; Holschuh et al., 2020), can be used to estimate the direction of arrival for energy in both the along-track and across-track directions (Paden et al., 2010) (Fig. 2b). Application of the MUltiple Signal Classification (MUSIC) algorithm in the Open Polar Radar toolbox (Open Polar Radar, 2023) to SAR-processed radar data therefore allows mapping of the subglacial topography in three dimensions along a single flight line.

Figure 3 (a) Cross-track surface tracking at each nadir bed position produced from the MUSIC algorithm (Paden et al., 2010). (b) The tracked surface is range migrated to convert the two-way travel time values to depth.

Along-track and fast-time averaging was applied to the SAR-processed data to improve the signal-to-noise ratio when tracking the bed return and cross-track surface. In the along-track direction, the depth range where the bed return was expected was selected using bed elevation values from BedMachine v.5 (Morlighem et al., 2017, 2022). The nadir bed return is picked in order to produce seed points from which the cross-track reflector is digitised, to geolocate off-nadir reflection information (Al-Ibadi et al., 2017). In the cross-track surface picking (Fig. 3a), the bed return power was fitted to a gaussian distribution in order to overcome the broad distribution and multiple peak bed returns from which the maximum is picked. Tracking of the bed pick is limited using a guided window based on a theoretical hyperbola of a flat bed, as there are often englacial reflections close to the bed return. A −2.5° static angle correction was also applied to account for tilting in the tracked cross-track surface, in order to align and overlap each swath accurately. For each trace, after cross-track surface tracking, the hyperbola of the bed return is range migrated, including the refraction of air to ice, to convert the surface tracking to depth (Fig. 3b). Here, upwarping at the edges of each swath introduced some artefacts into the data, as spreading of the bed return energy reduced the accuracy of the surface tracking. Whilst this could have been excluded from the final DEM, reducing the swath width prevented overlap of each swath, meaning that landforms could not be mapped distinctly. The range migrated data produced a point cloud, which is projected from the reference frame to geographic coordinates (latitude and longitude) using the true heading and flight trajectory from the radar data. The final DEM raster is then interpolated from the georeferenced point cloud dataset using inverse distance weighting, at a pixel spacing resolution of 25 m and a vertical resolution of ∼ 10 m. The resulting product is a three-dimensional tomographic image of the ice base, with swath widths of ∼ 2 km.

2.2 Characterisation of the subglacial landscape

In order to characterise the subglacial landscape, digital geomorphological mapping of the landforms was carried out in QGIS v3.26.3 (QGIS, 2025) (Fig. 4). Multiple hillshade directions and contour lines were utilised to define the edges of landforms (Chandler et al., 2018), which were manually delineated across the DEM. The landforms were initially categorised, without interpretational bias, as geomorphic ridges, lineations, channels, and basins. Lineations (Fig. 4a) were identified as linear, low relief, and flow-parallel features. Geomorphic ridges, in contrast to lineations, are landforms that are less consistent in their width and less elongate, meaning that they were more circular or irregularly shaped. Channels are incised valleys that cut across the landscape as linear to sinuous depressions (Fig. 4e). Basins were mapped as flatter areas in topographic low points (Fig. 4d), often between larger ridges, with a much smoother appearance compared to the areas of geomorphic ridges.

Figure 4 Examples of subglacial landforms identified beneath the onset of the NEGIS. The initial non-interpretational mapping of lineations (red lines), geomorphic ridges (black outlines), channels (dark blue thalwegs, light blue rims) and basins (purple outlines) is shown in each panel. These landforms were then interpreted as described below, examples of which are shown in each panel: (a) Mega-scale glacial lineations. (b) Crag and tails. (c) Barchan/erosional drumlins. (d) Sedimentary basins. (e) Meltwater channels. (f) Shear margin moraine.

Elongation ratios were calculated in QGIS as the length of a minimum bounding rectangle divided by the width. As the swath processing introduced some artefacts into the DEM, due to upwarping at the edge of each swath creating flow-perpendicular rises in elevation, only landforms which are clearly identifiable as geomorphological features, and could definitively be observed crossing multiple swath widths, were mapped.

2.3 Comparison of landforms to velocity data and shear margins

The NASA MEaSURES ITS_LIVE velocity mosaic (Gardner et al., 2022) was used to derive the streamlines of the current velocity field, as well as the overlying surface velocity related to each landform. The mean velocity value within each polygon was used as a representative value for each individual landform. The shear margins of the ice stream were delineated where the shear strain rate was at its maximum value from multiple cross-profiles. The shear strain rate was derived from TerraSAR-X ice velocity data (Hvidberg et al., 2020). Relative locations of subglacial landforms within or outside of the shear margins were determined using the outermost shear margin where the shear margins overlap, and the polygon delineating the landform had to be contained entirely within the shear margin boundaries.

3 Results

High resolution subglacial geomorphology at the onset of the NEGIS

A detailed image of the subglacial landscape is shown in Fig. 5, with our mapping of the landform assemblage shown in Fig. 6a (see also Fig. 4). Low relief, highly elongated lineations are interpreted as MSGLs (Fig. 4a) and are almost entirely constrained within the modern shear margins of the NEGIS (Fig. 5). The elongation ratios of these lineations range between 7.7 : 1 and 59.2 : 1, with the more elongate features generally located beneath the central trunk of the NEGIS (Fig. 6b). MGSLs can be identified as linear, low relief, and flow-parallel features, which generally have elongation ratios greater than 10 : 1, and are generally consistent in terms of width (Stokes and Clark, 2002b; Stokes et al., 2013; Ely et al., 2016). The majority are oriented in the direction of the overlying surface flow, although some close to the shear margins are oriented at an angle up to 13° away from the trajectory of the current flow field. Heights of the MSGLs have a mean value of 66.5 m, and range between 21 m (minimum) and 148.4 m (maximum). The MSGLs range from 968 to 7397 m in length, with a mean value of 3025 m (Fig. 6c), but some are truncated by ridges that have a rougher appearance in the landscape. In addition, there are instances where the ridges appear to act as seed-points for the MSGLs, which then continue downstream into the sedimentary basins and have a much smoother appearance (see identified MSGLs in Fig. 5).

Figure 5 Subglacial topography produced from swath radar measurements at the onset of NEGIS (Carter et al., 2025). The hill-shaded DEM (referenced to WGS84) allows interpretation of an assemblage of subglacial landforms, characteristic of ice streaming, such as mega-scale glacial lineations, crag and tails, and barchan drumlins. Location of the survey shown in Fig. 1 and the red box shows the location of Fig. 7. Arrows on top indicate True North and direction of ice flow. The line delineating the location of the shear margin (maximum shear strain rate) is dashed, as this represents the central point of the broader transition between slow and fast flow (the shear zone).

Figure 6 Subglacial landform mapping related to the overlying ice dynamics. (a) Geomorphological mapping of the subglacial topography, with surface ice flow velocity (Gardner et al., 2022) and streamlines (Hvidberg et al., 2020), which are flow-parallel lines, illustrating the direction of flow of the current steady-state velocity field. Inset polar histogram shows the orientation distribution of the MSGLs. (b) Scatter plot showing the elongation ratios of the MSGLs and geomorphic ridges compared to surface flow velocity. (c) Histogram of MSGL length. (d) Histogram of the elongation ratios of the geomorphic ridges, both within and outside of the modern shear margins of the NEGIS.

Geomorphic ridges were subsequently interpreted as three different categories: crag-and-tails (Fig. 4b), barchan drumlins (Fig. 4c), or a shear margin moraine (Fig. 4f). They are mostly interpreted as bedrock features, resembling crag and tails (Fig. 4b, metrics in Table A1). These are elongate features with a steep stoss side on a bedrock crag and lee side tails oriented in the direction of overlying ice flow. They are differentiated from MSGLs by their higher relief and with stoss-side crags having a rougher appearance, with a smoother tapering lee-side tail that is usually composed of unconsolidated sediment (Dowdeswell et al., 2016; Nitsche et al., 2016). There are some potential barchan, or erosional, drumlins (sensu Eyles et al., 2016), which are identified based on the grooves carved into the bedrock ridges, bisecting a larger “parent” drumlin into smaller outcrops that are grouped together in a crescentic shape (Fig. 4c). The bedrock ridges are clustered in groups, particularly in the southwest of the survey area, and are present both within and outside of modern shear margins, more evidently on the western side. The distribution of elongation ratios of these ridges demonstrates that they show some degree of elongation outside of the shear margins, as well as inside (Fig. 6d). However, in the absence of seismic data directly overlying these areas of hypothesised bedrock, this interpretation may have a higher level of uncertainty. Some geomorphic ridges (e.g. in the centre of Fig. 4e), with a long axis 45° oblique to the MSGLs and shear margins, may represent deformed ribbed bedforms, which are a transitional bedform resulting from the deformation of flow-transverse ridges (Vérité et al., 2023, 2024).

A long, low relief, narrow feature aligned closely with the south-eastern shear margin, but clearly misaligned with ice flow direction (see yellow arrows in Fig. 5), is interpreted as a lateral shear margin moraine (Stokes and Clark, 2002a). This is based on its proximity and alignment with the shear margin, its continuity across multiple swaths, and its relief of between 30 to 40 m that is prominent in the surrounding landscape (Fig. 4f). The ridge appears to be dissected in places, but this is likely due to artefacts at the edge of each overlapping swath (see Methods and Fig. 3a). This dissection may be a function of the ridge's orientation, which is oblique to the flightlines, as mapped MSGLs are less subject to this effect and appear continuous across multiple swaths. Riverman et al. (2019) also identified a shear margin moraine with similar dimensions of ∼ 500 m width and ∼ 50 m height, and orientation nearly parallel to ice flow, although it is not clear if we image exactly the same bedform.

Other features include depressions within which are likely to be softer/unconsolidated sediments, rather than hard bedrock, as might be suggested by their smoother appearance (Roberts et al., 2019; Dowdeswell et al., 2016). Previous seismic surveys crossing the topographic low point beneath the north-western shear margin demonstrate acoustic impedance values that correspond to soft, deformable till, indicating the presence of saturated subglacial sediment in this region, including where the MSGLs are identified (Christianson et al., 2014; Riverman et al., 2019). These sediments contained within the basins (Fig. 4d) appear as flatter areas in topographic lows, which are often downstream of bedrock features, the largest of which has an area of 49 km². Some large (between 0.5–3 km wide) incised channels with undulating thalwegs are also identified, often deviating around topographic highs and cutting across the landscape both inside and outside of the ice stream (Fig. 4e). These are interpreted as large meltwater channels, and potentially tunnel valleys which indicates incision by subglacial meltwater (Livingstone and Clark, 2016; Kehew et al., 2012), although the possibility remains that they are part of a preserved fluvial network.

4 Discussion

4.1 Implications for the basal conditions and spatial evolution of the NEGIS

Classification of the geomorphology beneath the onset of the NEGIS allows us to infer some aspects of the basal properties, as the differing regions of subglacial landforms would suggest heterogeneity of the bed. The geomorphology bears close resemblance to a “mixed bed” landform assemblage (Fig. 7), composed of a combination of hard bedrock outcrops interspersed with softer sediments and subglacial bedforms, such as those identified on palaeo-ice stream beds in west Greenland (Dowdeswell et al., 2014), the Amundsen Sea Embayment (Graham et al., 2009), western Northern America as part of the Late Wisconsin Cordilleran Ice Sheet (Eyles et al., 2018), and the North Sea (Roberts et al., 2019) (Fig. 7). These types of subglacial landscapes are characterised by the combined presence of crag and tails, drumlins, and highly elongate lineations, and thought to be formed under ice streams (Graham et al., 2009; Dowdeswell et al., 2014; Roberts et al., 2019; Eyles et al., 2018). Previous seismic data from around the EastGRIP ice core site (Christianson et al., 2014) showed saturated, high-porosity, presumably deforming sediment within the main trunk of the ice stream, and non-deforming and more compacted sediments outside of the main trunk. However, the presence of crag-and-tails would imply additional bedrock outcrops in the broader region of the onset now visualised by our new DEM. Christianson et al. (2014) infer from their observations that an underlying layer of dilatant till might explain the presence of the ice stream in this location and its lack of a major subglacial trough, but our observations of a mixed bed landform assemblage would suggest that the characteristics of the bed are perhaps more heterogenous than previously recognised.

Figure 7 Comparison of mixed-bed landform assemblages. (a) North Sea bathymetry data showing streamlined bedforms and bedrock highs (Roberts et al., 2019). (b) Bedrock features beneath the onset of NEGIS (inset in Fig. 5).

Within this mixed-bed assemblage, the distribution of the identified subglacial landforms can be used to infer changes in the location of the ice stream shear margins, using the crag-and-tails, which are commonly used as indicators of fast palaeo-ice flow (Dowdeswell et al., 2010; Nitsche et al., 2016; Dowdeswell et al., 2016). The degree of elongation of crag-and-tails outside the present shear margin is similar to those inside of the shear margins (Fig. 6c, Table A1), indicating a potential spatial evolution of the shear margins of the NEGIS. It suggests that they could have been subject to similar flow conditions as inside the shear margins of the current ice stream at some point in the past, as the current overlying ice velocity conditions outside of the shear margins are slower, and not streaming (10 to 25 m yr⁻¹, Fig. 6b). The crag-and-tails are more prevalent on the north-western side of the ice stream compared to the southeast, which correlates with internal folding of ice horizons that demonstrate a “jump” between an active and inactive section of the north-western shear margin, where the margin has shifted inwards (Jansen et al., 2024). The south-eastern shear margin, however, has a much more abrupt transition between the slow-flowing and fast-flowing ice inside of the ice stream trunk. The transitory nature of the north-western shear margin could, therefore, be related to a potential narrowing of the ice stream associated with the localisation of the shear margins, as streaming conditions would have overlain the ridges currently outside of the shear margins. Therefore, we propose that the ice stream may have been up to 10 km wider at some point in the past, or experienced convergent or a slightly different direction of flow prior to the localisation of the current shear margins (Bons et al., 2016; Jansen et al., 2024). The relative lack of crag-and-tails, as well as the preservation of a potential shear margin moraine evident along the south-eastern shear margin, suggest a more stable configuration on this side.

4.2 Implications for MSGL formation

Subglacial landforms used to reconstruct palaeo-ice streams, such as drumlins and MSGLs, have often been argued to represent a subglacial morphological continuum that evolves with cumulative basal slip speed (Stokes et al., 2013; Ely et al., 2016; Barchyn et al., 2016; Zoet et al., 2021; Ely et al., 2022; Vérité et al., 2024). Along this continuum, ice flowing at velocities of 10s to < 100 m yr⁻¹, i.e. similar to the velocity at the onset of NEGIS, has been associated with landforms oriented perpendicular to ice flow direction known as ribbed moraines, but merging into classic drumlins as velocity increases (Stokes et al., 2013). Hence, ice stream onset zones are generally thought to be associated with landforms that are transverse to ice flow, or oriented in the direction of ice flow but with low elongation ratios that are more characteristic of drumlins, rather than MSGLs (King et al., 2007; Anderson and Fretwell, 2008; De Angelis and Kleman, 2008). When modelled, as ice velocities increase to 100s m yr⁻¹, drumlins have been hypothesised to evolve and elongate into MSGLs (Stokes et al., 2013; Barchyn et al., 2016; Ely et al., 2022; Vérité et al., 2024), but there is little work that quantitatively links bedform elongation to ice flow/slip speed and the assumptions are mostly qualitative.

The geomorphological signature beneath the onset of NEGIS has a clear ice streaming imprint with the presence of MSGLs, which is perhaps surprising given that ice velocities within the study area currently peak at 58 m yr⁻¹. All previous observations of MSGLs beneath active ice streams are associated with ice velocities > 100 m yr⁻¹, such as the landforms detected beneath the main trunks of Rutford Ice Stream (Schlegel et al., 2022) (380 m yr⁻¹) and Thwaites Glacier (Holschuh et al., 2020) (between 100 and 400 m yr⁻¹), which have formed under continuous and relatively stable high-velocity ice streaming conditions (Woodward and King, 2009; Alley et al., 2021; Gudmundsson and Jenkins, 2017). The setting of these ice streams contrast with the temporally and spatially variable NEGIS (Franke et al., 2022a), with its slower (< 60 m yr⁻¹) velocities at the onset. To our knowledge, the only other high-resolution DEM of a modern ice stream onset zone (King et al., 2007) revealed classic drumlin forms (elongation ratios 1 : 1.5 to 1 : 4) and a potential ribbed moraine under Rutford Ice Stream, West Antarctica, where velocities accelerate from 72 to > 200 m yr⁻¹. The landforms beneath the onset zone of the Rutford Ice Stream are much less elongate, but are forming beneath higher ice flow velocities than our study area of the NEGIS. Furthermore, the location of the MSGLs in our study, 600 km from the grounding line and only ∼ 200 km from the main ice divide, is also the furthest inland that MSGLs have been identified beneath a modern ice sheet. Even so, ice stream onset zones clearly can extend very close to ice divides in both palaeo (Margold et al., 2015a, b) and modern settings (Bamber et al., 2000). Indeed, there are numerous examples from the palaeo-record that show ice stream imprints extending very close to the inferred position of ice divides (Hodder et al., 2016; Rice et al., 2020; Dubé-Loubert et al., 2021; Rice et al., 2024). However, it is very rare for MSGLs to be identified/mapped so close to an ice divide.

There are a number of hypotheses that might explain the presence of MSGLs beneath the modern-day slow ice velocities observed at the onset of the NEGIS. Firstly, one possibility is that the MSGLs are akin to bedrock mega-grooves and, therefore, their length is not directly related to ice velocity (Newton et al., 2018). However, we reject this hypothesis, as the morphometric characteristics of the MSGLs here are different to those relating to bedrock mega-grooves (Newton et al., 2023). The heights of the MSGLs average 66.5 m, which are much greater than the depths seen in bedrock mega-grooves (5–15 m; Newton et al., 2023). Whilst some studies have shown the majority of amplitudes of MSGLs to have low heights (between 0.4–38.1 m; Vérité et al., 2024), MSGLs of the scale identified here have also been observed beneath Thwaites Glacier, where elongated bedforms composed of soft till locally exceed 100 m in height (Alley et al., 2021), as well as in past ice streams off the Antarctic Peninsula (Canals et al., 2000). The contrast of the rougher bedrock outcrops to the smoother appearance of the MSGLs would also indicate that this is a mixed bed assemblage (Roberts et al., 2019), as the landscape is akin to a transitional “mixed bed” zone, where streamlined bedrock features and MSGL occur adjacently as the subglacial bed changes between hard and soft material (Eyles and Doughty, 2016; Mulligan et al., 2019).

In addition, the role of subglacial meltwater and its potential contribution to MSGL formation in this context is difficult to constrain. Subglacial hydropotential modelling has indicated that water is concentrated into and along the shear margins (Christianson et al., 2014; Franke et al., 2021), and moves through this region in a distributed system with only limited regions of continuous channelisation (Riverman et al., 2019), suggesting that it has a limited role in bedform development. Observations from phase-sensitive radar measurements at the EastGRIP ice core site quantify basal melt rates at 0.19 ± 0.04 m a⁻¹ (Zeising and Humbert, 2021), which is suggested to arise from the subglacial meltwater system, but information on subglacial water volume is very limited. Subglacial drainage cascades have also been observed, albeit much further downstream, by differential satellite synthetic aperture radar interferometry (DInSAR) (Andersen et al., 2023). As the amount of subglacial water cumulatively increases downstream, smaller amounts of subglacial water are likely to exist at the onset further upstream, but below the detection limit of DInSAR.

Secondly, it could be argued that a plausible hypothesis to explain the presence of the MSGLs under this part of the NEGIS, is that the study area experienced higher ice velocities in the past, after which the MSGLs have simply been preserved. This might therefore suggest that ice velocities under this part of the NEGIS were much higher (> 100 m yr⁻¹) at some point in the past. However, it is unlikely that the NEGIS experienced higher ice velocities prior to the localisation of the shear margins 2000 years ago (Jansen et al., 2024). The folds in the isochrones observed at the NEGIS onset by Jansen et al. (2024) are consistent with folding due to convergent ice flow (e.g. similar to Petermann Glacier), which are then sheared where they are intersected by the shear margins, causing them to rotate and tighten. The timing of this intersection of the folds by the shear margin is constrained to 2 ka by both the offset of ∼ 55–75 km of the fold hinges (as they are advected with ice flow over ∼ 2000 years) as well as the cessation of fold amplification at 2 ka. When the shear margins localised in their current location, the ice stream was decoupled from the interior of the ice sheet, as the rotation of the orientation of the ice crystal basal planes lead to mechanical weakening in the shear margins (Gerber et al., 2023), which in turn facilitates localised deformation and the faster flow observed at present. Higher velocities would have produced higher shear strain within the margins, for which the analysis of the fold amplitudes by Jansen et al. (2024) did not show any evidence. Modelling of the ice stream dynamics of the northeastern sector of the Greenland Ice Sheet since the Last Glacial Maximum has also shown no evidence of higher velocities within the interior of the ice sheet, prior to the formation of the NEGIS (Tabone et al., 2024). Thus, we would consider the plausibility of these MSGLs being formed from a faster configuration of the NEGIS to be unlikely.

Thirdly, it could be argued that the MSGLs under the onset of the NEGIS were formed by the current slow ice velocities over much longer (multi-millennial) timescales (Clark, 1993), where the present NEGIS has been very stable in this location at this flow speed. However, it has been shown that the NEGIS is a relatively young ice stream (< 2000 years) and that this sector of Greenland has undergone drainage basin reconfiguration over the Holocene (Franke et al., 2022a; Jansen et al., 2024). Earlier configurations of now-extinct ice streaming in the northeastern part of the ice sheet during the Holocene did not reach so far into the interior (Franke et al., 2022a; Tabone et al., 2024), instead draining a more northerly part of the ice sheet, indicating that the NEGIS as seen today is a relatively recent feature established within the last 2000 years (Jansen et al., 2024). Thus, although we are unable to rule out that the MSGLs may have formed by slow flow over a much longer duration (> 2000 years), we view this as unlikely because they correspond to the margins of the present enhanced flow, rather than everywhere across the study area. Indeed, if that were the case more generally, then MSGLs would be far more common on both palaeo and modern ice sheet beds. As well, there is evidence from the palaeo-record that ice streams can switch on and off over just a few hundreds of years. During deglaciation of the Laurentide Ice Sheet, it has been shown that around 40 % of ice streams operated for < 2000 years (Stokes et al., 2016; Margold et al., 2018), although very few of these ice streams record MSGLs so close to the inferred ice divide.

Whilst we cannot rule out an episode of enhanced flow at this location in a previous glaciation, the sedimentary basins outside of the northwestern shear margin (Fig. 6) show little evidence of MSGL formation, similar to our observations within the shear margins. This would mean that the ice stream would have formed in the same configuration as observed today, potentially with a higher velocity, to produce the observed MSGLs. Even so, if this had occurred, this would then suggest that MSGLs can be preserved for 100s to 1000s of years under relatively slow ice velocities. Rather, we propose that these MSGLs within the current shear margins, are likely to have formed in situ under modern-day slow ice velocities (< 60 m yr⁻¹), over the 2000 years in which the NEGIS onset has existed at this location, questioning the paradigm that they are exclusively associated with fast ice flow. This potentially could be attributed to the presence of highly deformable sediments, which have high porosity and high dilatancy (Christianson et al., 2014; Riverman et al., 2019), where these MSGLs have been identified.

5 Conclusions

In conclusion, high-resolution bed topography data from swath radar enables key observations into the basal conditions of an ice stream, as well as its potential spatial evolution. The identification of the previously unseen subglacial geomorphology beneath the onset of the NEGIS allows insight into the processes under which these landforms form. The presence of MSGLs so far inland under the NEGIS is a surprising discovery, given that ice velocities are < 60 m yr⁻¹ and the study area is only ∼ 200 km from the ice divide. We deem it unlikely that the MSGLs were formed by slow ice flow over a much longer duration than 2000 years, but are unable to rule out the possibility that a previous episode of fast flow created the MSGLs, which would have since been preserved under ice flow conditions similar to present. However, we view this as unlikely given that the folding of the internal stratigraphy and shear strain rates of the margins indicate a slower, confluent flow regime preceding the formation of the shear margins. Rather, our favoured interpretation is that the MSGLs are the product of the relatively low ice velocities (10s m yr⁻¹) observed today, within the 2000-year timescale of the current configuration of the NEGIS. If correct, these observations may prompt a re-evaluation of some models of MSGL formation and their use as an indicator of high (> 100s m yr⁻¹) ice velocities, but further work, such as a multitemporal analysis, is required to examine whether the features are relicts from a previous flow regime or are actively forming under current conditions.

Appendix A

Table A1 Quantitative metrics of geomorphic ridges mapped.