Main

Satellite observations of current surface flow of the ice sheets in Greenland (GrIS) and Antarctica reach back to the 1970s1,2,3. Reconstructions of earlier ice-stream activity are based mostly on their geomorphological imprints in formerly glaciated regions but are limited to the period since the Last Glacial Maximum4,5. These reconstructions show that ice streams were activated and deactivated in different locations and at different times during the post-Last Glacial Maximum deglaciation (22–7 thousand years (ka) ago). Patterns of former ice flow in regions that are still covered by ice sheets may be revealed by folds in internal reflection horizons (IRHs) that are visible in high-quality radio-echo sounding (RES) data6,7,8,9. Most IRHs are considered to be buried former snow layers deposited on the surface of the ice sheets10. Englacial deformation can disturb the conformity and continuity of IRHs and lead to a non-monotonically increasing age distribution with depth, which is a challenge for ice-core analysis11.

Folds in the GrIS and Antarctic ice sheet range from centimetre-scale folding observed in ice cores12,13,14 to folds that affect nearly the entire ice column15,16. RES data have revealed the presence of overturned and sheath folds and plume-like structures overlying broad areas of disrupted basal stratigraphy8,9,17,18. A range of processes has been proposed to explain the diversity of folding, including basal freeze-on9,18,19,20, contrasts in ice rheology12, changes in basal resistance21,22 and convergent flow17.

In this Article, we combine recently acquired and older RES data over the northeast GrIS to reconstruct former flow patterns. Our findings reveal the Holocene activity of two now-extinct ice streams in an area close to the present-day ice divide.

Radiostratigraphy

We consider RES profiles in northeast Greenland upstream of the northern catchment of the Nioghalvfjerdsbrae (79° N Glacier (79NG); Fig. 1a–d) acquired by NASA’s Operation IceBridge (OIB) and by the Alfred Wegener Institute, Helmholtz Centre for Polar and Marine Research (AWI). The RES profiles in Fig. 1c (Set 1, black profiles in Fig. 1a,b) were acquired with AWI’s multi-channel ultra-wideband (UWB) radar system in 2018. Ice surface velocity is almost zero in the west of Set 1 close to the ice divide and increases eastwards to as much as 15 m a−1 (Extended Data Fig. 1). The survey was designed to closely investigate a set of folds already detected in earlier surveys23. Twelve RES profiles were flown, at 7.5 km spacing, and oriented perpendicular to the 100° true North trend of the fold axes (Extended Data Figs. 1 and 2 and Supplementary Figs. 1–12). We mapped cylindrical folds that reach their highest amplitudes at the centres of the RES profiles (Extended Data Fig. 3). The fold shapes and orientations do not mimic the underlying bed topography (Supplementary Figs. 13 and 14 and Extended Data Fig. 4). Similar folds have been observed at multiple locations in the GrIS8,9,18,21, including upright cylindrical folds in the active convergent flow region of Petermann ice stream in northwest Greenland17 (Fig. 1f).

Fig. 1: Overview of the study site and key RES observations.
figure 1

a,b, Northern Greenland’s bed topography (a)44 and ice surface velocity (b)27 with the locations of two sets of RES profile lines (Set 1 in black and Set 2 in red). The black rectangle indicates the map outline of Fig. 2d,e. The major outlets in northeast Greenland are indicated with black arrows (79NG, 79° North Glacier (Nioghalvfjerdsbrae); ZI, Zachariae Isstrom; SG, Storstrømmen Glacier). c, Selected RES profiles of Set 1. d, A subset of three radargrams of Set 2 with a repeating radar signature. e, A radargram of Set 2 with a characteristic radar signature in the upstream region f, A radargram from the onset of the Petermann ice stream. g, Radargrams at the onset of NEGIS (see b). h, A radargram showing the NEGIS shear margin fold patter (see g). The white arrows in f and g indicate today’s general ice-flow direction, looking downstream.

Full size image

The RES profiles of Set 2 (Fig. 1a,b,d and Supplementary Fig. 15) are composed of OIB and AWI profiles acquired with different RES systems. They extend from the central divide eastwards into the northern catchment of the 79NG. The radargrams reveal a recurring signature of vertical stripes of low reflectivity in the lower two-thirds of the ice column (Fig. 1d,e (yellow and white markers), Extended Data Fig. 5 and Supplementary Figs. 16–21), which results from the loss of reflected energy from steeply dipping internal layers whose steepness increases with depth24. The pattern of these features and their appearance in multiple radargrams (Fig. 1d and Extended Data Fig. 5) bears a strong resemblance to similar features within the shear margins of the Northeast Greenland Ice Stream (NEGIS; Fig. 1g,h and Supplementary Figs. 24–29; fig. 7a in ref. 25 and fig. 3 in ref. 15).

A palaeo ice stream analogous to the Petermann ice stream

We traced three IRHs with estimated ages of 45, 52 and 60 ka (ref. 26) in the lower part of the ice column that are continuously visible in all radargrams (Fig. 2a–c). We reconstructed the shapes of the three stratigraphic horizons by interpolating surfaces between the IRHs in our RES profiles (Methods). This reveals two parallel cylindrical anticline–syncline pairs (Fig. 2b,c) whose geometry was validated using an OIB RES section oriented obliquely to the fold axes (Supplementary Fig. 13). The internal layers are folded far up into the base of the Holocene ice column (Fig. 2a). The steepest upper parts of the axial planes strike east–west with a shallowing of the dip at the base that suggests top-to-the-north bedrock-parallel shearing (Fig. 2b). The reconstructed folds show a strong resemblance to the folds at the onset region of the active Petermann ice stream, which are suggested to be products of convergent flow and lateral shortening17. In contrast to the Petermann folds, they are systematically overturned towards the north. They have smaller amplitudes than their Petermann counterparts, and their fold axes lie oblique (by ~25°), rather than parallel to present-day flow direction17 (Fig. 2d and Extended Data Fig. 4b).

Fig. 2: Folded radiostratigraphy in northeast Greenland.
figure 2

a, Radargram of Set 1 with three traced internal reflections (green, pale blue and orange; 45 ka, 52 ka and 60 ka, respectively) and axial planes of two synclines (red and yellow) and one anticline (dark blue). The base of the Holocene ice column (~11.5 ka) is highlighted with a white line. b, The reconstructed 3D internal horizons and axial surfaces. c, The individual horizons. The location of the radargram in a is indicated with yellow arrows. d,e, Overview of the locations of these Petermann-type folds (two cylindrical anticlines of Set 1) and the NEGIS-type shear margin folds (tight fold sequences) superimposed on the ice surface velocity27 (d) and superimposed on the gradient of that velocity on a logarithmic scale (log10) (e). The background map in c shows the bed topography (Supplementary Fig. 17b).

Full size image

Folds similar to the shear margins folds of NEGIS

We find characteristic upright, short-wavelength (~100–500 m) folds with amplitudes up to ~100 m extending almost from the ice divide to the northern catchment of the 79NG in the RES stratigraphy of Set 2 (Figs. 1b,d and 2d,e, Extended Data Fig. 5 and Supplementary Figs. 15–21 for details). Especially in the downstream portion of Set 2, these bear a remarkable resemblance to folds in the shear margins of the active NEGIS (Fig. 1d,e,g,h, Extended Data Fig. 6 and Supplementary Figs. 24–29). Set 2 reveals further similarities between the folds in our survey area and the NEGIS-type folds, including downwarping internal layers at the outer margins of the fold sequences, kink folds in the region between them (Extended Data Fig. 6) and characteristic reflections in RES profiles oriented obliquely to the fold axes (Supplementary Fig. 27). The folds at the active shear margins of NEGIS extend almost up to the ice surface, and the folds in Set 2 clearly disrupt the stratigraphy of the base of the Holocene (~11.5 ka; Extended Data Fig. 6). Furthermore, the folds in Set 2 coincide with a southwest–northeast-trending zone of slightly elevated surface velocity gradient that extends far upstream (Fig. 2e and Supplementary Fig. 22). A second trail of folds to the southeast is less pronounced, and far upstream, we find only solitary large folds, including one that was previously interpreted as bedrock (Supplementary Fig. 23).

Folding sequences and ice flow

Three different mechanisms for the origin of folds as in Set 1 have been proposed: basal freeze-on18, variations in basal resistance21 and response to convergent ice flow17. Folds resulting from basal freeze-on or basal-resistance effects are expected to trend perpendicular to the direction of ice flow17, thus, north–south in our study region. Similarly, present-day ice flow in the region of the folds is not convergent27. The current local flow regime can therefore not explain the occurrence of the observed folds. All three folding mechanisms17,18,21 require ice to flow in a given direction for a long period, but the inferred direction of ice flow is different for basal freeze-on and variations in basal resistance than for folding by convergent flow. According to the direction of the fold axes (which is determined by the orientation of the fold hinges) and the direction of overturning, ice flow should have been northwards if the folds were created by freeze-on or variations in basal resistance (Extended Data Fig. 4a,b). Fold formation by convergent flow would imply a flow direction parallel to the fold axes, which would thus be eastwards. A constant northward flow over a longer period is unlikely considering that outlet glaciers in the north are much smaller than in the east (Extended Data Fig. 4a,b). Furthermore, folding due to freeze-on and variations in basal resistance would require highly complex basal properties (for example, recurrent alternations of melting and freezing to create pairs of anticlines and synclines) for which there is no evidence (Extended Data Fig. 4c,d and Supplementary Section 1.3).

Therefore, on the basis of the three-dimensional (3D) geometry, glaciological setting and orientation of the folds, we consider the following two-stage scenario as the most likely one to explain initial fold formation and subsequent deformation (Fig. 3a). In the first stage, a convergent flow regime similar to the one in the onset zone of today’s Petermann ice stream created large upright folds. The axes of folds formed in such settings align with the positions of today’s major outlets in northeast Greenland (Figs. 1 and 4). The greater number of folds downstream (Extended Data Fig. 3) indicates an increase in the intensity of folding due to an increase in horizontal shortening. In the second stage, the flow regime had changed to its present configuration. Ice flow is oblique to the previous flow direction but no longer convergent, resulting in wholesale bedrock-parallel shearing of the existing folds (Fig. 3a). Applying the current flow field to restore the shearing and reconstruct a set of upright folds requires a period of at least 1.0–2.5 ka (Supplementary Fig. 31 and Supplementary Table 1). This is a minimum estimate as we do not know the exact geometry of the axial planes before shearing (Supplementary Fig. 31). Ongoing accumulation without further fold growth would have led to compression of the folds, reducing their amplitudes. Along-flow stretching and shearing during this second phase of deformation appears to have affected the downstream portion of the mapped features more strongly due to the higher ice-flow velocities there.

Fig. 3: Conceptual model for the ice-dynamic reconstruction from fold patterns.
figure 3

a, The formation of overturned cylindrical folds of Set 1. The upper row illustrates the formation of upright folds due to a convergent palaeo flow regime. The row below indicates a change from the palaeo flow regime to the present flow regime, which passively shears the upright folds. The schematic ice-flow regime is indicated in the left column, the effects of the flow mechanisms on the folds are illustrated with a 3D horizon in the middle column and the schematic layer deformation in a 2D radargram is in the right column. b, The reconstructed sequence of flow regimes producing the observed RES stratigraphy patterns in Set 2. An active NEGIS-type ice stream with localized ice flow (orange) created the tight folds at its shear margins (red). Later, ice flow ceased and subsequent accumulation buried the folded layers beneath younger flat layers (blue), compressing them. Outline of Greenland from QGreenland, NSIDC45.

Full size image
Fig. 4: Illustration of the palaeo ice-stream locations in northeast Greenland.
figure 4

The exact positions of the two palaeo ice streams during their active lifetimes depend on the temporal evolution of the ice-flow field since their deactivations (Extended Data Fig. 7). Magnitudes of present ice-flow velocities in Greenland27 are in the colour-coded background, and superimposed are the former convergent flow regime (blue arrows) and the NEGIS-type ice stream (orange arrow).

Full size image

The similarities of the RES signatures in Set 2 to radargrams observed at NEGIS suggest that the observed features are remnants of the shear margins of a now-extinct ice stream (Fig. 3b, Extended Data Figs. 5 and 6 and Supplementary Section 2). As folds act as passive tracers of ice deformation, remnants of past ice-stream activity are still visible even if the shear zone is no longer active. These observations suggest that the folding in Set 2 is the product of a ~20-km-wide palaeo ice stream with distinct shear margins similar to the currently active NEGIS (Fig. 3b).

By reconstructing folds from structural information archived in the RES stratigraphy, we have revealed the extents and characteristics of two palaeo ice streams (one Petermann type and one NEGIS type) in the northern 79NG catchment, both of which originally extended far into the ice-sheet interior (Fig. 4). It is difficult to ascertain exactly when the ice streams were active. Assuming it to be likely that the deformation processes are younger than the deformed radar stratigraphy, both ice streams can be shown to have been active into the Holocene, whose basal reflector (~11.5 ka) is deformed by both sets of folds (Fig. 2a and Extended Data Fig. 6b). The lateral shearing of the Petermann-type folds is on the order of 10–20 km (Fig. 3a). At the current surface velocity of ~3–5 m a–1, the lateral shearing must thus have operated for at least ~2.5 ka since formation of the folds. This is consistent with the youngest ice layers not showing substantial folding. Another possibility could be that the two palaeo ice streams do not represent synchronous separate events but instead evolutionary stages of a larger, evolving ice-stream system.

Ice streams come and go

Numerical modelling studies suggest that ice streams with highly localized shear margins can spontaneously appear by self-reinforcing processes within ice sheets28,29,30,31,32. These studies all require an initial instability to activate the feedback loop that enforces flow localization. This instability can be generated by coupling of the temperature field and ice flow using simplifications such as the shallow ice approximation29, but yields more consistent results by implementing longitudinal stresses (membrane stress approximation30). Schoof and Mantelli31 show similar patterns emerging due to a feedback of temperature and sliding at the bed. The study by Syag et al.32 presents temporally variable ice-stream systems, which switch between ‘catchment’ areas feeding the stream, depending on perturbations in the driving stresses. All of these studies use simplified, artificial geometry set-ups for their simulations as they all focus on investigating mechanisms and do not aim to simulate any particular real ice streams. Despite this, their results support the view that spontaneous localization of self-reinforcing flow by any of a number of mechanisms is likely to be a process that leaves observable traces in real ice sheets. These processes require no further evidence to be consistent with our observations of ice-stream products deep in the interior of the GrIS.

Our study provides extensive observational evidence that ice streams are temporally and spatially variable components of the northeast GrIS. The mechanisms and boundary conditions required for flow regimes to reorganize and generate ice streams in new locations are not yet fully understood. The initiation and termination of fast flow could be considered a particular case of surging behaviour33, whose specific properties vary according to temporally changing combinations of climate, ice-sheet geometry and bed properties, for example, the rearrangement of water routing networks34. However, in situ observations of some of these properties and detailed dating of our palaeo ice streams are still lacking, making it difficult to fully test and quantify this possibility for the northeast GrIS.

Flow regime shifts on similar spatio-temporal scales to our observations have previously been demonstrated around the outlets of the Siple Coast ice streams in Antarctica35,36,37. One, the active Bindschadler Ice Stream (formerly Ice Stream D), presents a single large fold that has been argued38 to have experienced shearing after originally forming as an upright flow-parallel fold, very similar to the multiple folds documented here. The switching behaviour of the Siple Coast ice streams has been related to the competition for ice discharge between neighbouring basins during a period of regional ice-sheet thinning37. A similar line of argument is followed by Stokes et al.5 in a reconstruction-based study of the Laurentide Ice Sheet. They suggest that streaming was triggered at times of large ice volume and that it ceased abruptly during volume reductions. Hence, ice-stream initiation is particularly likely to occur at the beginning of deglaciations. Geological evidence39 and modelling studies40 that indicate substantial changes in the geometry and elevation of the GrIS over the Holocene are consistent with this hypothesis. On the basis of the coarse dating that is possible, we suggest that a Holocene north-to-south reconfiguration of the ice-sheet geometry in northeast Greenland either (1) deactivated both northern palaeo ice streams, where the NEGIS and the palaeo NEGIS-type ice stream were simultaneously active or (2) relocated the main far-inland reaching ice drainage system southwards by deactivating the northern palaeo ice streams and initiating localized ice flow at the onset region of the NEGIS (Fig. 4). The latter possibility suggests that the upstream part of the present-day NEGIS is relatively young.

The evidence for the activity of a former ice stream with well-defined shear zones is also of major importance for ongoing research into comparable present-day ice streams. Today’s NEGIS is a prominent feature whose origin and activity have been controversially related to exceptionally high geothermal heat flux (GHF) at its onset41,42. However, the existence of such a high GHF was challenged43. Our discovery of an extinct ice stream that would have had very similar properties, but was located well to the north of today’s NEGIS, adds to the range of observations that suggest localized high GHF is not necessarily a prerequisite for NEGIS-type flow. Our observation that the palaeo NEGIS-type and Petermann-type ice-flow systems extended much farther inland than the present flow field at this location (Fig. 4) indicates that effective ice drainage expands or retreats depending on the ice discharge conditions forced by the overall glaciological setting. As a consequence, the divide between the eastern and northwestern catchments may have shifted to the west when the palaeo ice stream was active and possibly back to the east again after its deactivation. This has consequences for our interpretation of deep drilling projects (NEEM, NGRP, GRIP and so on) that are now on divides or domes, but may not always have been so during the Holocene.

Identifying the location, type and, if possible, time of activity of palaeo ice streams in the GrIS is of paramount importance for understanding ice-stream longevity and their past and future activity. Our observations provide a benchmark for the interpretation of other proxies of past ice-dynamic activity, for example, sedimentation rates and characteristics obtained from marine sediment records off the Greenland east coast. The reproducibility of our findings in models is particularly important to increase our confidence about whether, how and how suddenly the GrIS might reconfigure in response to ongoing mass changes and loss. Our work highlights that ice-stream activity is subject to temporal changes also on the larger ice-sheet scale and that the processes governing these changes need to be better understood and implemented in numerical models to accurately predict Greenland’s and Antarctica’s evolving contributions to sea-level rise in the future.

Methods

AWI UWB RES data profiles

In April and May 2018, we collected RES data upstream of the northern catchment of 79NG and at the onset of NEGIS with the UWB airborne RES system of AWI (black lines in Fig. 1a–c,g). The system consists of an antenna array with eight channels (co-polarized with the polarization plane oriented in flight direction), and the data were recorded at a frequency range of 180–210 MHz. For a comprehensive description of AWI’s UWB RES system and the survey design of the NEGIS survey, see refs. 6,46,47. We used the CReSIS Toolbox48 for RES data processing, which includes synthetic aperture radar and array processing16,46. For the geometrical analysis of the IRHs, we depth convert the two-way travel time RES data with a 2D dielectric constant (ε) model for air, with a relative permittivity εair = 1, and ice with εice = 3.15.

Construction of 3D horizons

For a 3D analysis of the IRHs, we build on the approach of Bons et al.17 and use the 3D model-building software MOVE (MOVE Core Application, Version 2019). The radargrams in depth domain were converted to the SEGY format with the ObsPy Python toolbox49 and imported into MOVE. Internal layers, which could be clearly identified in all RES sections, were traced manually. The traced layers were split into segments that show similar geometries (for example, anticlines and synclines) as segments of the same IRH in adjacent radargrams. For this purpose, the optimal orientation of the radar profiles is perpendicular (90°) to the fold axes, which approximately applies to the reconstruction of the 3D surfaces at Petermann17 as well as to the reconstruction of the Set 1 folds. Supplementary Fig. 30a shows how the line segments are interpolated between a pair of radargrams.

We express the degree of deformation of the 3D horizons with depth by comparing the area of the horizons for a specific region. The area with zero deformation (analogous to the horizons at initial deposition) is represented by a polygon with the extent of the referenced region (Supplementary Fig. 30c). The area of a distorted horizon is larger than the area of the polygon. The degree of deformation is then expressed by the ratio in percentage terms. The method generally underestimates the degree of deformation unless the layers between the folded horizons do not thicken.

Fold axial surfaces

We construct the axial surfaces of one anticline and two synclines that make up the cylindrical folds of Set 1 (Figs. 1c and 2). The method is analogous to the construction of the 3D horizons (Supplementary Fig. 30b). We trace the points of maximum curvature in the radargrams for two synclines and the anticline of the cylindrical fold to create a line segment and interpolate between the lines of each radargram.

Ice-flow reconstruction from axial traces

Due to the curvature of the fold axial traces, we define a minimum and a maximum distance for a measure of the displacement of the ice column (dSmin and dSmax) for five profiles of synclines 1 and 2, respectively (Fig. 2a). At each position, we use the absolute ice-flow velocity (vS) to calculate the velocity component parallel to the direction of shearing with the offset angle α: vy = vS cos(α) (Supplementary Fig. 31 and Supplementary Table 1). We assume that the ice-flow velocity has not changed over time and varies only slightly in space. Then we calculate the time needed for the displacement of the upper ice column tS, assuming that the deepest end of the axial trace is immobile: \(t_{\mathrm{S}} = \frac{{{\mathrm{dS}}}}{{v_{\mathrm{y}}}}\). The vector components of the velocity field are shown in Supplementary Supplementary Fig. 31a.