Holocene thinning in central Greenland controlled by the Northeast Greenland Ice Stream

Domain: earth_science
Published: 2024-07-31
Editorial outcome: accept

Abstract
Ice-core records from the interior of the Greenland ice sheet suggest widespread thinning during the Holocene. However, the recurring underestimation of this thinning in numerical models raises concerns about both the veracity of such reconstructions and the reliability of glaciological models. Recent work suggests the 8000-year-old Northeast Greenland Ice Stream (NEGIS), including a now-extinct northern tributary, may have been an early influence on Greenland ice-sheet dynamics. Yet, the inaccurate reproduction of NEGIS-like dynamics in most models hampers investigation of whether this feature played a role in Holocene ice-sheet thinning. Here we show that grounding-line retreat in northeast Greenland triggers elevation changes at the northern summit via ice-dynamic effects modulated by the paleo NEGIS system. In our simulations, fast ice-stream flow caused by transiently imposed reduced basal shear stress following the northeast retreat explains 55% (±18%) of the estimated ice thinning, showing that ice-stream dynamics is one of the main drivers of the NGRIP Holocene surface elevation drop. Our findings show that the ice-flow in northeast Greenland plays a larger role in ice-surface elevation changes in central Greenland.

Introduction
Understanding the history of the Greenland ice sheet (GrIS) during past warm periods is crucial to constrain its response to current and future climate change1. Local information retrieved from proxy records, such as the abrupt lowering of the ice surface elevation by several hundreds of metres (~200 to 1000 m) over the past 11,000 years estimated at ice-core sites2,3, provides valuable insights into the local climate and/or ice-sheet evolution, as well as constraints for evaluating paleo ice-sheet models. Yet, the capability of correctly capturing this thinning represents a long-standing challenge for three-dimensional ice-sheet models2–5. Improved climate and isostatic rebound model representations have helped to reduce the early-Holocene (~11 ka ago) data-model mismatch at High Arctic, peripheral locations from 900 m to only 200 m3. Nevertheless, these features cannot explain the reconstructed thinning at ice-core sites at lower latitudes, as High Arctic ice-elevation changes are mainly driven by the polar amplification, glaciostatic effects resulting from the demise of the nearby North American ice sheets3, and changes in precipitation2. Thus, the reasons behind the prominent thinning in central Greenland, such as the ~200 m elevation drop estimated at the NGRIP site2, remain unclear.

The proximity of the ~600 km-long Northeast Greenland Ice Stream (NEGIS, defined here as all ice stream branches and all outlet glaciers) to the ice divide at the present6, as well as the presence of another well-defined NEGIS-like stream stretching north from the 79°N glacier far to the interior during the early Holocene7, suggest a possible far-reaching influence of ice-stream dynamics during the Holocene (Fig. S1). Yet, large uncertainties in basal processes driving the NEGIS-like flow hamper its proper model representation8. The inclusion of subglacial hydrology schemes9, as well as extremely high basal geothermal heat flow (GHF)10, has improved the ability of models to simulate the present NEGIS fast flow. Surface-velocity assimilation for basal shear stress estimation yields the most accurate basal friction approximation, although only where the bed is currently temperate, assuming a specific temperature profile. Yet, assumptions of steady-state basal dynamics are only valid for decades-long simulations. The dynamics of the ice sheet can change considerably as the timescale increases due to non-negligible variations in the drainage system and changes in its geometry7. Therefore, the simulation of the paleo NEGIS flow pattern on multi-millennial time scales is a challenging task. Accordingly, its potential role in the GrIS evolution, in particular in the dramatic early Holocene thinning, remains unexplored.

Here, we investigate the impact of the NEGIS dynamics on the Holocene surface elevation drop in north-central Greenland. We perform an ensemble of 3000 simulations of the GrIS evolution through the last deglaciation using a state-of-the-art, three-dimensional, thermomechanically coupled ice-sheet model11 specifically adapted to represent various configurations of the paleo NEGIS (see the “Methods” section for details). The accurate timing of the deglaciation is ensured by evaluating the model results against deglaciation timings inferred from marine sediment records from the Northeast Greenland continental shelf. We characterise the past northeast fast flow regime as a system of a northern, now-extinct paleo ice stream and a southern present-day-like ice stream that activated during the Holocene by imposing a time-dependent basal friction change beneath the NEGIS (Fig. S1). Such a stream configuration is in line with recent improvements in our understanding of this system7. Basal lubrication is promoted by reducing the bed-dependent basal friction coefficient and by imposing sufficiently high GHF12, which results in high sliding rates for the ice-stream system during most of the Holocene. We find that NGRIP thinning is initiated by the NEGIS onset, enhanced by upstream propagation following early Holocene ice-sheet retreat in northeast Greenland, and maintained by increased meltwater production and reduced basal friction. The far-inland development of the present-like NEGIS during the last 8 kyr, here fostered by imposing a time-dependent reduced basal friction coefficient along the ice stream, drives the dynamic thinning, explaining at least ~100 m of the surface elevation drop in north-central Greenland. Such a far-inland reach of dynamic thinning in northeast Greenland during a major retreat phase suggests that such processes may also amplify the response of this sector to future warming and ice loss.

Results
The Northeast Greenland Holocene retreat
Our transient simulations show a progressive inland retreat of the northeastern sector throughout most of the Holocene (Figs. 1 and S2). At the end of the Last Glacial Maximum (LGM), the grounded ice sheet extends over the continental shelf, in agreement with reconstructions13. The ice-stream configuration at the onset of the deglaciation shows three well-defined fast-flowing regions draining ice mass from the present-day location of the Nioghalvfjerdsbrae (79°N) and Zachariae Isstrom (ZI) outlet glaciers to the shelf break. Two topographically constrained ice streams flow over the northern and southern troughs, whilst a third stream flows eastward over the shallow bank. The location of this ice-stream system aligns with glacial lineations mapped in the area from bathymetry data13. An ice stream flowing eastward and draining the shallow bank between the troughs would support the hypothesis of the GrIS extending toward the shelf edge during the LGM14.

Deglaciation began ~15 kyr ago when a rapid annual air temperature rise drove ice-sheet retreat to today’s ice front by about 10 kyr ago (Figs. 1 and 2). In our simulations, the NEGIS ice front retreats from the outer to the inner coast between 11.7 and 10.2 kyr ago, followed by a further recession of several tens of kilometres (~110 ± 40 km) 7–5 kyr ago (Fig. 2d). Such a retreat history compares well with cosmogenic-based Holocene ice-front reconstructions14 (see Fig. 2d) and independent 14C-exposure dating of marine organic material15. Although the retreat at the present 79°N ice-shelf front is simulated ~1 kyr earlier than in data (Fig. 2d), the timing of the outer-to-inner shelf recession aligns closely with deglaciation timings of both moraines at today’s NEGIS ice front14,15 and marine locations in the northeast shelf derived from gravity core analysis (blue dots and triangles of Fig. 3, respectively).

The fast inland retreat simulated at the early Holocene is accompanied by the formation of a paleo-NEGIS stretching from the northern flank of present-day NEGIS outlet glaciers into the ice-sheet interior (Fig. 1). The catchment shows three main well-defined north-to-south outlets and replaces the system of fast glaciers simulated over the northeastern shelf during the LGM. Increased sliding velocities propagate upstream in response to thermo-mechanical and geometrical changes occurring during the retreat of the marine-based ice sheet, through an imposed basal friction coefficient reduction through time. At the onset of the deglaciation, ice-mass loss at the shelf edge due to rising air and oceanic temperatures (Fig. S3) weakens buttressing, accelerating inland ice flow. This results in glacier thinning, ice flow speed-up, and retreat of the grounding line, which further accelerates the ice flow, promoting ice-front recession. Although ice thinning at the front reduces basal friction (as effective pressure decreases), sliding significantly increases when a local increase in subglacial melt occurs due to increased frictional heating (Fig. S4). As the retreat continues, regions showing decreased basal friction conditions allow for flow acceleration farther upstream toward the ice divide (Fig. 2).

Our simulations first show a surface elevation increase in north-central Greenland from 12 to 9.4 kyr ago (Fig. 2a), mostly due to increased precipitation2 (Fig. 2f). Then, the velocity propagation associated with the imposed inland expansion of the paleo-NEGIS induces dynamic thinning of ice at the NGRIP site (Fig. S5). Northern summit ice is advected downstream, causing the surface elevation to decrease abruptly (Fig. 2a) until the end of the Holocene Thermal Maximum (HTM) (~5 kyr ago). For the last 5 kyr the northeastern ice front slowly regrows towards its present position, reducing the NGRIP surface ice velocity (Fig. 2b) and the rate of thinning (Fig. 2a). Yet, the decrease in basal sliding imposed at the present NEGIS (Fig. S1) ensures continuous drainage of the GrIS northern sector until the present (Fig. 1). The hypothesis of the presence of a large, active ice stream resembling that of today since the mid-Holocene is further supported by a recent study of past ice-flow trajectories connecting the EGRIP ice back to its source location16.

Influence of ice stream dynamics on the NGRIP thinning
The surface-elevation increase induced by the rise in precipitation at the early Holocene is rapidly followed by a significant drop persisting until the late Holocene (Fig. 2a). We simulate a score-weighted mean (see the “Methods” section) maximum elevation change of 105 ± 34 m, corresponding to 55 ± 18% of the reconstructed 190 m surface-elevation fall2 (Fig. 4). Yet, the model scenario showing the highest thinning within the score-weighted 1-σ uncertainty shows a maximum elevation drop of 139 m, indicating that up to 73% of the total estimated thinning is reproduced.

We propose that the mechanism behind the early-Holocene GrIS north-central thinning is its dynamic adjustment to the NEGIS-induced inland flow acceleration during margin retreat (i.e. the retreat of the ice front and associated drainage basin). This hypothesis is in agreement with studies from extinct ice streams in Canada and Iceland, suggesting that the dynamic switch of paleo-ice streams has crucial consequences on ice-sheet stability17,18. To test this hypothesis, we investigate the relationship between the NEGIS dynamics, as defined by its basal sliding state, and the nearby NGRIP elevation change. We find that the NGRIP Holocene thinning rate is well correlated with the NEGIS basal shear stress (R = 0.59; Fig. 4b, c), calculated as the average between all the NEGIS catchments (northern, central and southern branches), with the strongest dependence for the early-mid Holocene interval with respect to the late Holocene (Fig. S6). This implies that faster sliding in the proximal NEGIS onset zone enhances thinning at NGRIP, which intensifies as the ice-stream onset propagates inland. No similar correlation is found for any other GrIS basins in our experiments (Fig. S7), indicating that there is no other catchment where ice flow plays a significant role in the NGRIP thinning.

The bed-elevation-dependent basal friction coefficient strongly controls the ice-stream basal conditions, hence the NGRIP thinning (Fig. S8). The imposed changes in the basal friction coefficient (Fig. S1) represent any transient processes at the base of the ice sheet not explicitly represented in the model. These could include changes to bed roughness due to sediment transport and deposition, as well as changes in subglacial drainage and routing. Yet, basal lubrication from melting enhances the flow. In our simulations, the initial increase in velocity is associated with dynamic adjustment due to ice-margin retreat (Fig. S5), which causes a decline in the basal effective pressure and friction, and an increase in meltwater production stemming from increased frictional heating. Sustained basal melt rates at the early Holocene saturate the basal till within a few thousand years, promoting a progressive upstream decrease in basal traction and faster flow (Fig. S4). Here, a mean GHF of ~63 mW/m2 12 is required along the NEGIS flow line to reproduce the rapid NEGIS inland expansion upon ice-flow acceleration needed to drive a significant NGRIP elevation drop (Fig. S9). Despite even higher GHF values suggested in northeast Greenland as attributed to the passage of the Icelandic hotspot19, a GHF of just 63 mW/m2 allows our model to simulate a temperate bed at the NEGIS, in line with estimates20. The high basal melting rates modelled for the same region are also in agreement with recent work21.

Model simulations with a very fast ice stream result in a larger Holocene NGRIP thinning (Fig. 4b, grey points). Yet, these tend to present an overall low skill score (Fig. 4c), mainly due to mismatches in the simulated present NGRIP elevation, northeast geometry, and ice-ocean front position (Fig. S10). Here the NEGIS bed roughness parameter is set to a constant value for each experiment. The internal switch applied at the mid-Holocene (Fig. S1) and the time-dependent variables in the friction law (effective pressure, see Eq. (4) and scaling parameter λ depending on bedrock elevation, see Eq. (3)) introduce some time variability in the definition of the friction coefficient. However, these approximations, the coarse model resolution and the simple hydrology model employed here are otherwise insufficient to reproduce the complex, unsteady dynamic state of the ice stream due to continuous fluctuations in the basal water system22. The inclusion of physical processes important for describing present-day ice-stream behaviour, such as non-local basal water routing9 could allow for a better investigation of the complex NEGIS dynamics in the future. Such complex behaviour of the NEGIS observed at present might still reflect the response to large geometry adjustments that occurred during the last deglaciation22,23, likely associated with a substantial change in the streaming activity in northeast Greenland during the Holocene7.

Early-Holocene temperature in the Northeast
We investigate the role of forcing climatology in the Holocene retreat by running an additional 3000-member ensemble experiment for a different deglaciation temperature scenario. This suggests a high sensitivity of the northeast margin retreat to the climate forcing (Fig. 3). The Holocene temperature anomaly based on ice-core-adjusted climate model simulations24 is here replaced by independent paleoclimatic reanalysis reconstructions25. Both scenarios are based on the transient climate Trace-21ka experiment26, yet the diverse ice-core data and techniques used to improve the original climate outputs may result in temperature differences, enhanced in areas with poor proxy records for model evaluation. This is clear at the NEGIS ice front, where significant discrepancies are present at the early Holocene (Fig. S17, right panel, yellow and blue lines). We find that temperature reconstructions that peak earlier (~9 kyr ago) and higher (+3.5 °C with respect to preindustrial levels) in the northeast sector are crucial for a correct retreat (Fig. 3).

A comparison of our two ensembles points to an early-Holocene temperature transition that should be closer to that estimated for the Agassiz ice cap3 (Fig. S17, right panel). Yet, the latter shows an even sharper (+7 °C compared to pre-industrial levels) and earlier (10 kyr ago) glacial-interglacial temperature rise, which could be related to high-latitude amplification effects24. An early HTM (11–7 kyr ago) would also be supported by marine proxy records in the Northern Hemisphere27. Summer air temperatures at EGRIP in the early Holocene could have been 3 °C higher than today28, while northwest and central Greenland could have been even warmer (3–5 °C)29. Interestingly, our early northeast retreat is in contrast with the ice-front recessional pattern found in southwest Greenland associated with a late HTM30. This heterogeneity in time and magnitude of warming trends may be interpreted as the climatic response to regional effects that are reflected in the asynchronous ice-front retreat.

Drivers of the northeast retreat
Simultaneous air and ocean temperature increases took place around ~15 kyr ago (Fig. S3), but submarine melting only occurred 13.3 kyr ago when the oceanic temperature anomalies increased by 2 K from their glacial value (see the “Methods” section). Simulations with extremely low oceanic forcing (κ ~ 0) still simulate a retreat from the continental shelf, potentially causing rapid thinning at NGRIP (Fig. S8). This implies that in our model, initial deglaciation in the northeast is triggered by atmospheric warming. Elevated air temperatures intensify surface ablation, thin the floating tongues, accelerate ice flow, and lead to ice front retreat. Sensitivity tests using different climatologies highlight a primary control of atmospheric temperatures on inner shelf retreat timing (Fig. 3, e.g. Sondre Mellemland, Lambert Land).

Still, introducing a significant ocean temperature anomaly at ~13 kyr ago (see the “Methods” section and Fig. S3) is crucial for accurately modelling the Holocene retreat rate. Initial tests with water temperature data from Trace-21ka exhibited a consistent delay in the mid-inner shelf retreat (not shown). Abundant evidence from sediment core analysis supports the presence of Atlantic Waters (AW) on the outer continental shelf as early as ~13.4 kyr ago31,32 and on the inner shelf during the early Holocene (~11 kyr ago33,34). This hints at the potentially significant role of warmer ocean temperatures in the disintegration of the northeast ice sheet’s floating front and grounding-line migration during the last deglaciation.

The central role of the oceanic forcing in the northeast ice sheet’s retreat could be similar to that observed at the present. Ice depletion at the grounding line of 79°N and ZI due to increased water temperature is the principal cause of their current retreat35. In the case of ZI, air temperature warming decades ago reduced sea ice, promoting calving and floating tongue collapse, decreasing buttressing on the glacier36. However, the present ZI grounding-line recession is predominantly attributed to the intrusion of warmer North AW in the fjord35, an effect seen persisting until a topographic ridge stabilises the retreat37. Topography restrained AW inflows at the 79°N glacier front38, yet recent observations reveal warmer waters beneath the ice shelf and near the grounding line39, causing the ice tongue to lose most of its mass through submarine melting. The potential collapse of the floating tongue would have severe implications for the near-future stability of the glacier, as reduced buttressing prompts acceleration, thinning, and grounding-line retreat.

We propose that the northeast margin retreat occurred due to a combination of warmer air and oceanic temperatures. The abrupt air temperature increase (>15 °C) around ~15 kyr ago initiated the retreat through enhanced surface ablation, thinning the ice front and favouring grounding line retreat. Warming AW on the northeast shelf since ~13 kyr ago likely amplified the retreat through increased submarine melting at the grounding line, fostering ice margin instability and accelerating inland retreat. Sea-ice buttressing is not simulated here, therefore its potential impact cannot be evaluated.

Extended dynamic effect of the retreat on the interior of the ice sheet
The NGRIP lowering can be interpreted as the indirect response to margin retreat via geometry and subglacial water system alterations. In our simulations, the northeast ice sheet’s retreat during the last deglaciation has an extended dynamic effect on the ice-sheet interior (Fig. S5). Margin thinning and retreat, driven by the early Holocene temperature rise, induce velocity propagation at higher elevations, causing ice thinning, reducing effective pressure and increasing basal sliding (Fig. S4). Still, basal sliding increases considerably only when a significant basal water increase from enhanced basal frictional heating saturates the till. This reduces the effective pressure and basal friction, increasing velocity along the ice stream and allowing upstream paleo-NEGIS propagation (Fig. S11). Notably, the onset of the NGRIP dynamic thinning (and surface elevation drop) occurs independently of whether basal friction beneath the NEGIS is artificially reduced or not, since it is triggered by the dynamic adjustment and velocity propagation in the northeast caused by margin retreat (Fig. S5). In fact, the NGRIP thickness decrease occurs following the decrease in basal friction at EGRIP, which in turn follows the decrease in basal drag averaged at the NEGIS basin due to margin retreat (Fig. S12).

Still, the imposed reduction in the basal friction enhances the dynamic thinning and velocity increase at higher elevation throughout the rest of the Holocene. This is primarily due to an earlier saturation in the hydrological system and a structural enhancement of inland velocity once saturation occurs. The extensive and sustained dynamic effects of the retreat, evident in the best simulation, contrast starkly with the limited impact observed in an unperturbed test simulation where a constant cf = 1 is imposed (Figs. S4, S11, and S13). Indeed, in this unperturbed simulation, the constrained velocity increase toward the ice divide limits the NGRIP dynamic thinning for the duration of the Holocene. Therefore, although the NGRIP thinning is triggered by the northeast retreat and the dynamic adjustment associated with it, independently from the applied NEGIS basal friction coefficient, the sustained thinning during the Holocene is controlled by the high sliding environment at the NEGIS, here favoured by the imposed reduced basal friction coefficient.

It has been shown that, over time spans ranging from centuries to millennia, variations in surface elevation trigger a corresponding reaction in the basal water drainage system (and vice versa)22. As a result, adjustments in thickness and velocity can extend several hundred kilometres inland through wave propagation40. Long-term ice-flow acceleration can be propagated at almost the highest elevations in the ice sheet as a result of mass loss at the terminus41. This has been observed at ZI recently, where the ice front retreat of ~10 km has induced dynamic thinning along the NEGIS more than 200 km upstream42. This has also been observed in some glaciers in west Greenland43. Moreover, localised perturbations at the Antarctic ice-shelf front have been shown to propagate even more than 900 km upstream along the ice stream44. These findings align with our hypothesis that ice acceleration and ice front retreat in northeastern Greenland could have induced dynamic thinning more than 600 km from the ice front on centennial to millennial timescales.

Secondary effects on NGRIP thinning from other parameters
Ice anisotropy, controlled by ice temperature and stress variations, exerts a major control on the NGRIP drop down during the Holocene (Fig. S8). Stiffer ice (represented via a lower value of E in our model) in shear-controlled areas promotes slower ice flow, hampering upstream dynamic mass loss and facilitating the surface elevation increase at the early deglaciation. Although the enhancement factor varies spatially across the ice sheet (see the “Methods” section), we do not investigate how different streaming flow factors might affect the dynamics of fast flow. Still, our simulations show higher ice viscosity for ice stream regions (Estream = 1, Eshear > 1, see the “Methods” section and Table S1). This agrees with flow factors recently estimated at the NEGIS from radar and ice-flow modelling, suggesting stiffer ice along the stream flow and softer at the shear margins45. Although simulations with slower ice flow show greater accuracy in predicting ice thickness within the northeast catchment, they simulate too much ice in the ice sheet interior at present (Fig. S10). This highlights the need for spatial variability in most parameters (e.g. cf from Eq. (2) and fp from Eq. (10) in the “Methods” section) to fully capture the complexity of the ice-sheet dynamics.

Second-order effects on NGRIP thinning arise from parameterised physical processes like submarine melting, GIA, and surface ablation. The submarine melt rate primarily controls the timing and magnitude of the northeast ice front retreat triggered by the intrusion of warmer AW at the early Holocene along the NEGIS coasts31–34,46. A stronger ice-ocean heat flux exchange (see Eq. (7) in the “Methods” section) ensures an early onset of the northeast retreat, matching the deglaciation timings of sediment core locations from the outer and inner coast, generally inducing the ice front to retreat farther inland. However, a higher ice-ocean heat flux tends to reduce the NGRIP drop down overall. Cooler atmospheric temperatures during the Neoglacial associated with a reduced advection of Atlantic-sourced waters promote ice regrowth towards the present ice front, causing a reduction in dynamic thinning. If such a regrowth is not outpaced by a moderate/strong active ice stream sustaining ice mass flow from the interior, such as in the case of our best simulation, surface elevation tends to increase after the HTM, reducing the total Holocene NGRIP elevation drop.

Higher lithosphere relaxation times approximating higher mantle viscosities cause a delayed bedrock response to the ice unloading during the deglaciation. A delayed bedrock uplift in the northeast, as the ice front retreats inland and ice mass is lost, keeps the bedrock elevation low enough to favour earlier sliding in streaming areas (see Eq. (3) in the “Methods” section) and leads to a prompt inland propagation of velocities. This initiates early dynamic thinning in central Greenland, causing the NGRIP surface elevation to decrease prior to the model having time to respond to the increase in snow accumulation at the Holocene onset. A northeastern higher mantle viscosity would induce an earlier thinning in north-central Greenland, extending the drop down duration as suggested by data2. This would agree with the presence of a higher mantle viscosity in central Greenland47. Yet, a faster deformation (lower relaxation time) compares better with elevation data overall, as the model has a delayed response in the surface elevation increase around 11–10 kyr ago. A more complex GIA model accounting for lateral mantle viscosity variations would likely reduce the reconstruction-model mismatch related to such an uncertainty48.

Enhanced surface ablation at higher latitudes (see Eq. (6) in the “Methods” section) results in higher northeast ice mass loss and grounding line retreat. The geometric reorganisation upon ice-front recession propagates velocities farther upstream, leading to a higher NGRIP elevation drop.

Mismatch between modelled-reconstructed NGRIP surface elevation change
The NGRIP thinning history simulated in our full-deglaciation simulations compares well with that estimated from ice-core data2,4 (Figs. 2, 4a), in particular from ~10 kyr ago onwards. Although the model simulates the surface elevation change throughout the Holocene qualitatively well, it fails to capture the full NGRIP lowering both in terms of magnitude and timing.

The residual data-model mismatch in the early Holocene could stem from the underestimated precipitation due to atmospheric warming. A comparison with Holocene NGRIP accumulation data49,50 indicates an underestimation by ~0.08 m/yr (Fig. S14a), largely attributable to the Regional Climate Model’s climatological mean51. A supplementary test, by adjusting the reference precipitation for the 0.08 m/yr deviation in north-central Greenland, with an unvaried precipitation rate of change, shows minor improvement (~20 m) (Fig. S14). Still the response time of the NGRIP surface elevation to changes in accumulation at the early Holocene might have been faster than modelled. In fact, the rapid NGRIP elevation increase suggested by data occurs only ~1 kyr after the abrupt increase in accumulation (from 11.7 to 10.7 kyr ago, Figs. 2a, f and S14). Climatic reconstructions50 would support a faster early Holocene accumulation rise than the reanalysis precipitation signal used here25. A higher rate of change in accumulation could yield a faster elevation increase, possibly reducing the data-model temporal offset before 10 kyr ago.

Elevated basal ice temperatures due to higher GHF (Fig. S9) might hinder elevation gain at the early deglaciation by fostering sliding. Increased GHF aligns with borehole data indicating a temperate bed beneath NGRIP52. However, consensus is lacking on north-central Greenland’s basal thermal state20. The early-Holocene elevation rise could also be affected by uncertainties in post-glacial isostatic rebound. Our approach attributes bedrock uplift solely to GrIS mass loss, neglecting the GIA signal from Innuitian ice sheet deglaciation3, which could impact the isostatic rebound in central Greenland. Our simulations underestimate bedrock uplift in north-central Greenland by half (see Fig. S15 and Fig. 11 from ref. 5). As the North American ice sheet system disappeared almost completely 8 kyr ago5, including such a GIA signal might initiate an earlier NGRIP bedrock uplift, increasing the modelled surface elevation.

We attribute the temporal offset in the NGRIP lowering mainly to a simulated delay in dynamic thinning following margin retreat. The elevation decrease results from dynamic thinning outpacing the accumulation increase. In our simulations this occurs around 9.8 kyr ago, ~800 years after the onset of the reconstructed surface elevation decrease2. We attribute this discrepancy to a simulated delay in the early Holocene retreat along the inner northeast coast (Fig. 3, right panel). Although the retreat from the outer-middle shelf (before 11 kyr ago) is well simulated, retreat near the present-day ice front lags by centuries (millennia for Lambert Land and Sondre Mellemland). This implies that the impact of the northeast ice stream in north-central Greenland may have occurred earlier than our simulations indicate, thereby outpacing accumulation-driven elevation rise. Furthermore, the effect of the dynamic propagation from margin retreat strongly depends on the location where the paleo-ice stream initiates. It is possible that the streaming area was situated closer to the ice divide already 10 kyr ago than shown in Fig. 1. In our model, the ice stream onset is governed by the inland extent of regions of low effective pressure (Eq. (4)). The predefined NEGIS basal friction coefficient, altered temporally (Fig. S1, left panel), dictates the potential presence of an ice stream. However, actual stream activation occurs when effective pressure dips low enough to reduce basal friction and initiate fast flow, i.e., when basal water increases due to enhanced frictional heating (Fig. S4). In our simulations, this occurred around 10 kyr ago, as the northeast ice front still retreats from the inner shelf (Fig. 1), and the increase in basal shear stress produces basal water to lubricate the bed (Fig. S4). However, the presence of an active hydrological drainage system at the NEGIS already at the early deglaciation (~11 kyr ago), and its reorganisation associated with dynamic thinning, could have induced an earlier and far-reaching paleo ice stream onset. This cannot be accurately simulated by our model as the latter takes into account only local changes in the water system53. In support of an earlier paleo-stream activation, Franke et al. (2022)7 suggested that ice-stream initiation occurs likely at the beginning of the deglaciation, possibly driven by large geometry shifts, subglacial water routing variations22, and catchment piracy17. An early Holocene paleo ice stream extending well inland might have induced early ice advection from north-central Greenland, decreasing the NGRIP surface elevation earlier. Finally, we cannot rule out that the retreat of the northwest coast had a non-negligible impact on the NGRIP elevation drop down. In fact, although basal friction in northwest Greenland does not correlate well with the NGRIP thinning rate in our model (Fig. S7), it also inadequately represents northwest dynamics and retreat. We likely simulate a delayed retreat in northwest Greenland, as the applied temperature anomaly exhibits a damped early Holocene increase compared to reconstructions, which, on the contrary, indicate a rapid temperature rise before 10.5 kyr ago3. The early regional retreat might have initiated a stronger early Holocene dynamic thinning in the northwest, possibly affecting north-central Greenland’s elevation history. Moreover, we overlook the buttressing effect of the demise of the Innuitian Ice sheet, potentially amplifying retreat from the American continent. An earlier temperature increase at the early Holocene, might reduce this time discrepancy.

Uncertainties in the reconstructed Holocene thinning curves2 might be underestimated, as (1) they have been estimated by assuming a homogeneous climate across Greenland4 and (2) the uplift-corrected δ180 signal used to infer the site’s elevation changes does not consider the non-negligible influence of ice-load changes in the Innuitian ice sheet at the Holocene onset3. Such uncertainties might partly compensate for the missing elevation drop simulated at the early-Holocene.

For the period 8–0 kyr, our modelled evolution fits well within the reconstruction uncertainties. Yet, ice-core data show a greater elevation decrease throughout the last phase of the Holocene (e.g. the last ~4 kyr). Although we impose a decreased basal friction coefficient at the present NEGIS for the last simulated 8 kyr, we fail to reproduce the rate of the late-Holocene lowering. We attribute this discrepancy to the northeast ice front advance after the HTM, which slows down the ice flow in the northeast, reducing ice advection and the NGRIP elevation decrease. Still, our late-Holocene NGRIP elevation change is consistent with previous modelling efforts evaluated against RSL estimates4,5.

Here, we attribute at least half of the estimated NGRIP elevation drop to the dynamic thinning associated with the Northeast retreat as induced by the Holocene warming. Such dynamic thinning is the result of ice-flow acceleration in the northeast triggered by the formation and inland development of the paleo-NEGIS during the last deglaciation. Our results imply that a future northeast ice front retreat stemming from increasing temperatures might similarly provoke tens of metres of dynamic thinning in central Greenland. Given that the present-day NEGIS velocity pattern is not captured by most ice-sheet models54, our results suggest that inland dynamic thinning in north-central Greenland is underestimated by future projections and Greenland’s contribution to future sea-level rise may be greater than currently anticipated.

Methods
Ice-sheet-shelf model
We use the open-source state-of-the-art Yelmo ice-sheet model11,55 version 1.801 (https://github.com/palma-ice/yelmo/). Yelmo is a three-dimensional, thermodynamically coupled, higher-order ice-sheet-shelf model that is here run over the Greenland domain at 8 km resolution. The model has been recently improved to adopt the depth-integrated-viscosity approximation (DIVA) velocity scheme56,57, a higher-order stress approximation proven to be efficient for multi-millennial simulations. This solver is found to replicate well the flow found in the ice-sheet interior but also that of fast ice streams and of floating ice shelves57.

The basal shear stress (τb) is defined through the regularised Coulomb friction law58 as:

τb = βub = Cb ub |ub|q / (|ub|q + u0q)

where ub is the basal velocity, q is the exponent of the friction law and u0 is a reference velocity parameter (100 m/yr). The bed-dependent term Cb is defined as

Cb = λ cf Neff

where cf is a unitless coefficient representing the local bed characteristics. The bed-dependent term is decreased through a factor λ depending exponentially on the bedrock elevation (zb)55, as

λ = max(exp((zb - z1) / z0), λmin) if zb < z1
λ = 1 if zb > z1

with z1 = 400 m and z0 as a tunable parameter. This relationship ensures that ice flows faster in topographic depressions55. The parameter cf is set to a constant value everywhere, except in the northeastern basin. Here, cf is reduced to better capture the system of paleo and present ice-streams (Fig. S1). In this way, we attempt to replicate the recent findings from radar stratigraphy analysis suggesting that at least part of the Holocene was characterised by the presence of a NEGIS-like stream located to the north of the present NEGIS and now deactivated7. We outline this paleo NEGIS following the location of the internal deformation features (folds) found from radargrams. This allows us to define a stream that extends along the northern part of the northeast catchment from the 79°N glacier towards the ice divide (northern branch of Fig. S1). The present NEGIS (southern and central branches) has been set as a region characterised by present surface velocities >30 m/yr for latitudes <77.5°N and by velocities >50 m/yr for higher latitudes. This approach defines a long penetrating stream in correspondence with the present-day main branch and a more constrained streaming region in the upper tail. The central branch of the ice stream is kept active for the whole simulation with a fixed value (fmid). There is great uncertainty around the timing of activation/deactivation of this ice-stream system and the processes that might have caused it. Yet, significant restructuring of the ice stream’s shape and strength is more probable during the initial phase of deglaciation due to the large ice-sheet geometrical reorganisations upon margin retreat. Therefore we assume an abrupt relocation of the ice-stream flow from its past to its present position around 8 kyr ago. This is also supported by the recent hypothesis of the NEGIS showing a similar configuration to today since 8000 years ago16. To this end, cf is linearly reduced in the southern branch to fmin to favour the build-up of faster flow along the present-day NEGIS, even as the ice sheet re-advances during the Neoglacial. This ensures a further reduction of the NGRIP surface elevation (Fig. S13). Therefore, in our simulations, cf can be set to either 1, fmid or fmin depending on the considered NEGIS branch and time in the Holocene (Fig. S1). The parameters fmid and fmin are perturbed in our simulations (see Table S1). Although our approach approximates the past activation/deactivation of the northeast ice stream through an imposed relocation of the paleo ice stream to its present position, the root causes of such an internal switch are not investigated. The timing and the form set for the internal switch between the ice-stream modes have negligible effects on the NGRIP surface elevation lowering (Fig. S16). We investigate a potential effect also for different ice-stream configurations and we conclude that considering both streams active at the same time instead of an abrupt switch during the Holocene has a marginal impact on our results (Fig. S16). Also, it should be noted that in our simulations, the presence of an active northern branch for part of the Holocene is not crucial for triggering the NGRIP drop down. On the contrary, the activation of the elongated ice stream in the southern branch at 8 kyr ago, propagating several hundred km inland (Fig. 1), ensures a sustained dynamic adjustment in north-central Greenland driving the NGRIP drop down (Fig. S5, lower panel). Conversely, an unperturbed simulation restricts the inland velocity propagation to lower elevations (Fig. S11), limiting the effect of the stream activity at the NGRIP site (Fig. S13). The cf term reduced at the NEGIS (fmid) is the dominant term within the basal friction parameters to directly influence the NGRIP thinning rate at the early-mid Holocene, as suggested by the high correlation between the two (Fig. S8).

The effective pressure Neff is computed through the relationship for till hydrology53, such that

Neff = N0 δ (P0 - s) 10e C0 c (1 - s)

where P0 = ρi g H and s is the fraction of basal water in till (Wtil) over its maximum value held by the saturated till (Wmax til, here set to 2 m). Values of e0 and Cc are set as in ref. 53, whilst the fraction (δ) of the overburden pressure P0 is considered here as a tunable parameter. Wtil evolves in Yelmo as

∂Wtil / ∂t = (ρw / ρb) bg - Cd

where the basal mass balance, bg, is computed for grounded regions. Such a subglacial hydrological scheme only accounts for local changes in the water system and the till rate of drainage Cd is fixed and set to 0.001 m/yr.

The surface ablation is calculated by the insolation-temperature melt (ITM59) method, a simplified scheme that, contrary to PDD schemes, directly accounts for insolation changes and melt changes due to surface albedo. The potential surface melt is calculated as

Ms = Δt / (ρw Lm) * (τa (1 - αs) S + c + λT)

where c = c1 + b(lat - lat0), Lm is the latent heat of melting, αs is the surface ice or snow albedo, S is the incoming insolation, Δt is the day length and λ is the long-wave radiation coefficient (10 W/m2K). The terms b and c1 correct the surface ablation for high-latitude regions and the ice sheet uniformly, respectively. τa is the atmosphere transmissivity, calculated as τa = 0.46 + 0.00006 * (zs)0.5, where zs is the surface elevation60. The reference latitude lat0 is set to 65°N and both b and c1 are considered as tunable parameters.

A flow enhancement factor is introduced in Glen’s flow law to approximate the effect of ice anisotropy on the ice flow11. Specific flow factors depending on the flow regime can be assigned in the model. The enhancement factor for floating areas (Eshelf) is set equal to 0.7, that for streaming zones (Estream) is set to 1, whilst that for shearing-driven regions (Eshear) is considered as a tunable parameter. Subsequently, the shear and stream enhancement factors are averaged for grounded areas, with the weighting determined by the vertical shearing fraction computed at every grid cell11.

Bedrock elevation changes are computed via the Elastic Lithosphere Relaxing Asthenosphere (ELRA) scheme61. This simple but effective glacial isostatic adjustment (GIA) scheme assumes viscoelastic deformation of the crust from the elastic behaviour of the lithosphere and the ability of the asthenosphere to relax upon ice load changes. Here the rigidity of the lithosphere (He, see Fig. S8) and the upper mantle relaxation time (τ, see Fig. S8) are considered perturbation parameters in our ensemble.

The calving rate is calculated following the Von Mises stress criterion62, here rewritten in terms of strain rates56. The constant value κt from Eq. (73) of ref. 56 is here interpreted as a tunable parameter. The submarine melting around the ice sheet is computed through the anomaly method described in 63, such that

Bm(t) = Bref + κ ΔTocn(t)

where Bref is a reference melting rate (set as equal to 1 m/yr) and ΔTocn is the paleoclimatic ocean temperature anomaly around Greenland. κ is the ice-water heat flux coefficient (also called sensitivity factor), and it is treated as a tunable parameter. Melt below the ice shelves is considered here as the 10% of the grounding line melt.

Climate and ocean forcings
Ice sheet evolution is forced by a transient climate experiment over the full-deglaciation period, similar to ref. 30. Atmospheric temperature is defined as

Tatm(t) = Tref + ΔTatm(t)

where ΔTatm is the 22 kyr-long monthly temperature anomaly from the fully coupled atmosphere-ocean TraCE-21ka experiment corrected for the AMOC strength24. Precipitation changes over the Holocene are imposed as

P(t) = Pref * Pf(t)

where the past precipitation fraction Pf is calculated at sea level through this precipitation–temperature relationship:

Pf(t) = P(t) * exp(-fp * (Tatm - Tsealevel))

Tsealevel is the atmospheric temperature corrected at sea level, as Tatm + δ * γatm, where γatm is the atmospheric lapse rate and δ is the climatology-modelled surface elevation anomaly. P(t) is the 20-kyr spatially variable precipitation fraction reanalysis reconstructed from a new ice-core data assimilation approach25, and fp is the precipitation scaling factor for past temperature anomalies (similar to β of Eq. (8) of ref. 25). This factor should be defined as spatially variable to represent the heterogeneous precipitation-temperature relationship throughout the ice sheet25, however, we consider it as a tunable parameter only for the Northeast basin. Three different precipitation scenarios P(t) (low, moderate, high) representing the relative amount of accumulation during the deglaciation are generated. Here, we use the high precipitation scenario, which suggests a rapid and substantial increase in accumulation in the early Holocene. Pref has been processed to refer to the 1981–2010 climatology of ref. 51 as Tref.

Transient ocean temperature anomalies ΔTocn are taken from the TraCE-21ka main experiment26. Temperatures are averaged at 300–350 m in the ocean column. We correct the TraCE-21ka temperature anomaly for the Northeast marine basin with a homogeneous temperature signal that mimics the continuous advection of subsurface warmer Atlantic Water (AW) to the shelf from the early to the mid-Holocene, following the well-constrained evidence from sediment core analysis (33 and references therein). The corrected signal shows a positive subsurface ocean temperature anomaly with respect to glacial conditions already at 13.3 kyr BP, a time that represents the minimum age of presence of AW found in both troughs in the northeast31,32. ΔTocn is set as constant until 10 kyr BP, after which it declines slowly until 7.5 kyr BP. A subsequent cooldown of subsurface waters due to the reduced intrusion of AW in the shelf31–34 is then approximated by a constant low-temperature anomaly until almost the end of the Holocene (the oceanic temperature anomaly recovers to present values during the last 100 years, see Fig. S3). In this setup, maximum and minimum Holocene-present temperature anomalies are defined as +2 °C (early Holocene) and -1 °C (mid-late Holocene). The anomaly from the late glacial to the early Holocene corresponds well with the warming of ~3 K found from marine records of subsurface ocean temperature in Fram Strait64.

In an additional explorative test, the model is forced by one of the Greenland-wide 20 kyr-to-present annual air temperature scenarios reconstructed through the paleoclimatic assimilation technique presented in ref. 25 (S3 scenario). Here, annual temperature anomalies are converted into monthly fields by applying the seasonal variability from temperature data from ref. 24, having them smoothed below 50 yr through a low-pass Butterworth filter. The considered scenario has been reconstructed by assuming a spatially homogeneous but stronger precipitation seasonality. We choose this scenario among the others available because (1) it is one of only two scenarios that, by applying a monthly anomaly, allows a proper glacial expansion, and (2) it presents the sharpest early-Holocene transition, approximating best the temperature signal of ref. 24. Here, precipitation and oceanic temperatures are the same as in the main experiment.

Deglaciation simulations
A classic steady-state spin-up approach is used for model initialisation. The ice sheet is allowed to freely evolve under LGM (~20 kyr ago) conditions for 50 kyr to reach the thermomechanical equilibrium54. The model is then forced toward the present by the transient climatologies described above. Full deglaciation experiments are performed for a large ensemble of 3000 simulations where five basal friction parameters (fmid, fmin, q, z0, δ), three climate variables (c, b, fp), two isostasy-related parameters (He and τ), two frontal/submarine ablation terms (κt and κ) and one material property (Eshear) are perturbed using the Latin Hypercube Sampling approach (e.g., ref. 59) (see Table S1 for an overview of the perturbed parameters). Note that while q, z0, δ, c and b are set as uniform around Greenland, the friction coefficient parameter is changed only beneath the northern/southern/centre branches of the NEGIS, following the temporal scheme of Fig. S1. All other basins hold the same cf = 1. This approach allows us to directly evaluate the effect of a fast-flowing ice stream on the thinning in central Greenland. The ice-sheet topography is taken from BedMachine v362. The bedrock elevation is reduced by 1-σ beneath the main tail of both paleo and present NEGIS to constrain the direction of the flow during the retreat, where σ is the BedMachine-estimated elevation uncertainty (Fig. S17, left panel). This simple approach allows us to reproduce both the northern paleo NEGIS-like ice stream that was active during the Holocene and the complex present-day NEGIS velocity pattern after a full deglaciation simulation with unprecedented detail (Fig. 1). The geothermal heat flow is prescribed as in ref. 12 below the ice sheet and as in ref. 65 outside the continental borders.

Validation process
The large ensemble undergoes a validation test to rank the goodness of the model simulations and to statistically weigh the results. The approach follows a simple averaging skill-score method based on a Gaussian error distribution66. For each ensemble member, we calculate 13 individual scores based on paleo and present model-data misfits and then aggregate them into a single total score. However, contrary to ref. 66, here we do not account for data and structure model uncertainties in the misfit calculation. The misfits considered are: mean-square error (MSE) of modelled to observed present GrIS thickness (HPD, data from ref. 62), MSE of modelled to observed present GrIS surface velocity (UPD, data from ref. 67), absolute modelled-observed modern surface elevation difference at NGRIP (SNGRIPPD, data from ref. 68), MSE of modelled to reconstructed deglaciation elevation changes at NGRIP (SNGRIPDG, data from ref. 2), MSE of modelled to observed present GrIS ice-covered area (APD, data from ref. 62) and MSE of modelled to observed ice-covered area at the LGM (ALGM, data from ref. 5). The errors HPD, UPD and APD are additionally calculated for the northeast sector only, so that three further misfits are considered (HNE, UNE, and ANE). Moreover, we calculate the MSE of modelled-observed deglaciation timings computed and averaged (1) at four sediment core locations from the northeast continental shelf (TSED, data from refs. 31–34), 2) at three locations in the outer shelf calculated from 10Be exposures dates from morain/boulder samples (TOUTDG, data from ref. 14) and (3) at other six locations closer to the present NEGIS ice front (TINDG, data from ref. 14) (see Fig. 3). For moraine and boulder data, deglaciation occurs when the cell is modelled as partially or completely ice-free. For sediment core data, we consider a point as deglaciated if it is either ice-free or still presents a floating shelf, as the paleo environment interpretation of such data is uncertain. Errors in the deglaciation timing are here calculated with respect to the time interval of 10.2–9.6 kyr ago for ref. 34, 13.4–11.2 kyr ago for ref. 31, 10.9–10 kyr ago for ref. 33 and 14.0–13.3 kyr ago for ref. 32. The resulting three marine-moraine MSEs describe the lag in simulating the time of deglaciation of specific locations in the northeast sector during the early Holocene. Finally, we introduce a last binary control on the upstream retreat of 79°N glacier: simulations that retreated more than 70 km upstream will present a maximum error. Each misfit Mi,j, where i is the data class and j is the ensemble simulation, is then interpreted as a score Si,j following the equation66:

Si,j = exp(-Mi,j / median(Mi,j))

The simple multiplication of the 13 Si,j yields the simulation-dependent total score Sj:

Sj = Π Si,j

which is further normalised between 1 (best simulation) and 0 (worst simulation). We assume a good total score if it is higher than the product of the mean score for each misfit for best simulations, i.e. 1e-7 (Fig. S10). Thus, the score-weighted ensemble mean of a variable x is calculated as

x = (Σ xj * Sj) / (Σ Sj)

with its standard deviation as

σ = sqrt( (Σ (xj - x)^2 * Sj) / ( (N-1) / N * Σ Sj) )

where N is the total number of ensemble members.

To further reduce the uncertainty in our model results, we add another constraint that distinctly separates plausible from implausible model representations within the ensemble. Following this approach, simulations are considered valid as long as the absolute elevation difference between modelled and observed elevations is < 50 m (e.g., SNGRIPPD < 50 m) so that 563 simulations over 3000 are retained as valid. This additional constraint rewards simulations exhibiting an accurate fit to present-day elevation data, whose errors in the estimate are usually below 10 m. Thus, the ensemble statistics are calculated only for this subset of valid simulations. This means that we calculate the score-weighted mean and standard deviation of the plausible simulations only within the ensemble (N = 563). All simulations are a priori considered equally likely (i.e., given equal weight in the statistics).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Main comments\nOverall, this is a well-designed study that tackles an interesting feature of Greenland Ice Sheet history. The findings advance our understanding of the ice stream’s response to climate change. The methodology mostly appears sound, and the sensitivity experiments and validation process provide confidence in the results. However, there are some aspects of the study and manuscript that should be improved, which I outline below.\n1. Proposed mechanism for the drop in ice surface elevation\nFirst, the hypothesis that is being tested regarding the mechanism of the drop in ice surface elevation at NGRIP needs to be made explicit. Second, the glaciological basis should be more clearly outlined. There are several characteristics of ice sheet change, such as grounding line migration, ice surface elevation change at NGRIP, basal shear stress, ice surface velocity change, but what is the relative order of events? Taking the abstract sentence ‘Fast ice-stream flow caused by reduced basal shear stresses following the NEGIS formation’ as an example, does fast ice-stream flow not dictate the formation of NEGIS?\n2. Sensitivity to ocean forcing\nThe modelling experiments are heavily focused on atmospheric forcing, principally on air temperature, and there is limited consideration and discussion of ocean forcing. This may be justified in studies such as Briner et al. (2020) where Holocene retreat was mostly across a terrestrial margin, but NEGIS would have had a marine ice margin during its Holocene retreat that would be sensitive to oceanic forcing. An ice stream may respond differently dependent on whether retreat is primarily atmosphere or ocean or atmosphere driven, and while the simulations fit reconstructions fairly well, they could appear correct for the wrong reasons.\nAn ocean temperature anomaly derived from the TraCE-21ka climate model is applied as a forcing, and the magnitude of this forcing is varied in a sensitivity experiment. However, 1) How realistic is this forcing? Syring et al. (2020, Paleoceanography and Paleoclimatology) state that warm ocean inflow occurred near the mouth of NEGIS at 10-7.5 kyr ago, which corresponds to the timing of reconstructed retreat and ice surface elevation lowering at NGRIP. This study and potential driver should be discussed. 2) Does the model adequately allow for a sensitivity to ocean forcing? The model resolution is 8 km, yet the width of NEGIS trough (where Holocene retreat occurred) is 20 km, meaning just 2.5 model cells are used to capture the retreat of the grounding line through this important, overdeepened trough. Does sensitivity to ocean forcing increase with a finer model resolution?\n3. Clarity of manuscript\nThere are multiple parts of the manuscript where the meaning of the wording is ambiguous, preventing a clear understanding of the motivation, method and interpretation. Several parts of this text are identified below.\nFurther comments\n(unfortunately there are no line numbers to pin-point locations in the text)\nAbstract:\n- “both reconstructions and models” – Reconstructions from ice cores? Ice sheet models, climate models, or both?\n- “The hypothesis…” is unclear – What hypothesis? That NEGIS has been in a similar configuration to today since 8000 years ago?\n- “thinning gap” – What thinning gap? The wording needs to more clearly link to the previous sentences.\n- “early retreat in northeast Greenland” – here and elsewhere, ‘retreat’ needs more specific wording – retreat of the grounding line in northeast Greenland?\nMain\nPart 1:\n- “Punctual information” – What does this mean?\n- “surface elevation drop” – This wording may be misunderstood by a reader, and should be more descriptive as the first time it is mentioned. \"Abrupt lowering of the ice surface elevation by X over Y years\"?\n- “data-model mismatch” – It would be useful to first describe the timing and magnitude of this data-model mismatch.\n- “meltdown” – I suggest more common wording (e.g. retreat, demise).\n- “Northeast Greenland Ice Stream (NEGIS)” – It would be useful to include a map with the ice stream identified, perhaps with the locations of various reconstructions that are used in other figures.\n- “this aspect” – What aspect?\n- “inversion techniques are not valid…” – Explain why they are not valid, and how that is relevant to this study/your approach.\n- “Here we explore…” – This sentence could include “(see Methods for details)”, or similar, at the end.\n- “inland penetration of the NEGIS flow upon margin recession” – Unclear.\nPart 2:\n- “the large ice cap” – What ice cap? Are you referring to the Greenland Ice Sheet in the northeast region of the continental shelf?\n- “agreement with reconstructions…” – What aspect is in agreement with what reconstructions? The cited 2011 review paper shows that there are few constraints (at the time of publication) in this region.\n- “high Arctic temperatures” – needs to be more specific (e.g. air temperature, annual air temperature).\n- “trigger and drive” – What is meant by \"trigger\" vs \"drive\"?\n- “ice sheet continental border” – What does this mean? The present-day ice sheet margin/coastline?\n- “the NEGIS margin” – Needs to be more specific. The grounding line or the ice (shelf) margin?\n- “Such retreat history is congruent with…” – refer to the relevant figures where the simulations are compared to the data.\n- “ice stream penetrates the interior” – Unclear and confusing wording. Describe the glaciological characteristics.\n- “tail of the stream\" - What does this mean? The lower reaches of the ice stream? The terminus?\n- “In our simulations,…” - Why does the ice surface elevation increase prior to 9 kyr ago when air temp increases and retreat occurs? Why is there an offset between the timing of reconstructed and modelled ice surface elevation change?\nPart 3:\n- “potential relationship between…” - What is the argument here? That basal sliding state explains the NGRIP elevation change?\n- “This implies that a faster-sliding and…” - What is the chicken, what is the egg? The wording suggests that basal sliding affects the NGRIP elevation, but this relationship needs a clearer glaciological explanation. What if these characteristics are both effects of retreat, and basal sliding changes are an effect of dynamic thinning?\nPart 4:\n- This is clear at the NEGIS margin,…” - Is the air temperature shown in Fig.4b just of the NEGIS margin?\nPart 5:\n- “Our late-Holocene results are consistent with…” - What aspect is consistent?\n- “NEGIS continental intrusion” - Unclear.\nFigure 1:\n- It should be made clear at the start of the caption that the modelled data in blue is from this study.\n- Wording of “NEGIS ocean border” and “present margin was ice-free” is not clear.\n- “Annual temperature” - specify that this is air temperature.\nFigure 2:\n- Try for a colorblind friendly sequential color scheme, rather than divergent scheme.\n- The NEGIS basin should be identified in this figure, if not elsewhere.\n- “best run” - refer to Figure 3.\nFigure 3:\n- “coloured following their Holocene maximum thinning” - The total amount of thinning across the Holocene at NGRIP?\nFigure 4:\n- How was the dashed red curve determined? I assume that this curve is supposed to be the same as that in Fig. 2a in Larsen et al. 2018, but it does not seem to be exactly the same (e.g. differences in the variations in the slope of the curve, and retreat in last ~500 years).\nMethods:\n- “the upper mantle is uniformly set to relax in 3000 yr” - Why? Is this over 3000 years from the perturbation, by 3 kyr ago, or over the last 3 kyr?\n- “RMSE of simulated to reconstructed grounding line movement…” - Is this done using the approximated grounding line position time series from Larsen et al., or against the specific data points for different locations?\nData availability:\n- “All Yelmo simulations presented here are available upon request” - The authors should provide the code used to run the experiments, and at least the best fit simulation.", "reviewer_2": "The paper investigates a possible explanation for the significant ice sheet elevation drop (~200m) reconstructed during the last deglaciation and through much of the Holocene, from analyses of the NGRIP ice core. It undertakes an ensemble numerical modelling exercise where they assess the dynamics of the Northeast Greenland Ice Stream (NEGIS), and the role it can play in the drawdown of elevation from the ice divide. The model is forced with high temporal resolution (not defined) climate data, has reduced basal friction and high geothermal heat flux and provides a good representation of the velocity pattern of NEGIS. As the ice sheet margin migrates landward during the deglaciation so the onset zone of the NEGIS migrates farther back towards the ice divide and thus explains about 50% of the ice-core derived elevation change at NGRIP. Overall the findings of the paper support the general consensus of knowledge at the present time.\n1 – It is unclear to me that this inference can be so straightforwardly made. During the deglaciation NEGIS evolved from a shelf-based terminal zone some 300 km in length into a more topographically constrained terminal zone (esp 79°N and Zacharie Isstrøm). If in the future the margin steps back into the “basin” area then it is perhaps less clear how the dynamics of the ice stream will evolve. It seems that the boundary conditions may be very different. In addition, Fig 2 e, f, g show the ice margin inboard of the present day and the drawdown of ice at NGRIP continues until ~4 ka (and longer in the ice core reconstruction). There is no modelled recovery in ice elevation presumably associated with the ice margin advance as seen in the ice core reconstruction. So how well the model is capturing the detail of the drawdown is less clear than being able to explain at least half of the total drawdown. Additionally thinking of a forward look is where of course where the numerical modelling will be invaluable. However, without such modelling I think this inference is ill-founded.\n2 – are these numbers correct. It appears from Fig 2e that the ice margin at 8 ka is further inboard than the current grounding line location and retreat continues until at least 6 ka. The “further recession” seems to be rather more than ~100 ±20 m, which I would have thought to be rather below the model spatial resolution? Indeed, there is no information provided on the spatial resolution at which the model has been run?\n3 – it is unclear why this description is used? The surface velocities depicted in Fig 2 A-D appears to show two well-defined ice streams linked to the broad NEGIS catchment draining thought Westwind and Norske troughs? This fits well with mapping of the sea floor which shows evidence of glacial lineations (e.g. Arndt et al, 2017 http://dx.doi.org/10.1016/j.quascirev.2017.01.018; Olsen et al, 2020 - doi.org/10.5194/tc-14-4475-2020) which in such settings is often liked to ice stream. The fast ice flow seems well structured and topographically controlled by the troughs crossing the continental shelf?\n4 - The model fails to capture the duration and rate of NGRIP drawdown (see also point 1) and some comment on this is warranted. In light of this, the reduction in modelled surface velocities is likely to be a function of the “steady” NEGIS. The longer and lower rate drawdown reconstructed from NGRIP does not at first sight square with the modelled surface velocity reduction. Some comment on this is warranted.\n5 – The grounding line location of NEGIS characterised by the record of Larsen et al, 2018. It appears that in this paper the grounding line in the model (easily identified) is assumed equivalent to the ice margin location defined by Larsen et al, 2018. However I do not think they are the same thing. My understanding of the Larsen reconstruction is that the distance is measured from the present ice margin i.e. the edge of the ice shelf, and not the present grounding line. From the caption for his Fig 2b – “¹⁰Be ages with external uncertainties from outer coast (dark gray) and outside the Little Ice Age moraine (black), and ¹⁴C dates from raised marine deposits presently dammed by NEGIS up to 70 km upstream the present ice margin at NG”. I am also not clear how the Bennike and Weidick (2001) data have been interpreted in this paper i.e. what “Radiocarbon dates at the NEGIS ocean border” means. I assume these ages relate to the samples collected around the epi shelf lake Blåsø. These samples are taken proximal, +10 km beyond, the present grounding line. The samples include marine mammal bones (seal and whale) and shells and driftwood. The bones and shells are assumed to indicate that Blåsø is no longer a lake but is fully connected to the fjord with no ice shelf but possibly permanent/semi-permanent sea ice in the fjord. The driftwood indicates that there was no ice shelf at Blåsø and the fjord was of sea ice free in summer. Bennike and Weidick (2001) provide no estimate for a grounding line location just saying that “the open fjord may have extended as far west as the present grounding line, just south of the western end of Blåsø.” My reading of Figs 1, 2, 4 and 5S and the text indicates that this paper has assessed the modelled grounding line distance measured from the present location of the GL at 79°N. The validation process has calculated using the goodness of fit based on the “RMSE of simulated to reconstructed grounding line movement at 79°N (79_GL, data from Larsen et al., 2018)”. Larsen provides the location of the ice margin (actually the ice shelf margin) not the grounding line. If I am correct then the scoring matrix needs to be recalculated to and the results reassessed. If I am incorrect then the Figures need to be amended and how the ice margin record of Larsen et al, 2018 can be correlated to the grounding line of the model output.\n6 – given the potential importance of ice shelves to the flow of their up-stream ice streams, it is surprising that there is no mention in the model, nor in the modelling of an ice shelf? Does the model include an ice shelf? It does not seem so from the surface velocity show in Fig 2, or is the velocity only for grounded ice shown? Some explanation of how the model deals with ice shelves is warranted, what exactly is being shown in Fig 2 and the potential impact of ice shelf loss, for example their loss e.g. Zacharie Isstrom, on the flow of the NEGIS is needed.", "reviewer_3": "Dear authors,\nThank you for the opportunity to review this manuscript. In general, I found the manuscript enjoyable to read, the results exciting and novel, and the conclusions well-supported. This manuscript will be of interest to many in the glaciology and paleoclimatology communities. The results have important implications for projecting future sea-level rise and highlight an area of ice-sheet model development that has been elusive for decades (representing the flow of large ice-streams and in particular the NEGIS) but is clearly critical for both hindcasting and forecasting ice-sheet dynamic responses to a changing climate. This is an excellent manuscript and with attention to a few key adjustments outlined below, I think this manuscript is very suitable for publication in Nature Communications.\nI applaud the authors for their clever experimental design and novel approach to diagnosing the role of ice-stream dynamics in driving thinning in the ice-sheet interior. In my line edits, I have flagged where I think a statement is too strong or slightly over-interpreted, but in general the conclusions drawn here are well supported by a clear experimental design and a careful analysis.\nI have one overarching question about the framing in this paper – in some places (like the title) there is a geographic restriction to the NEGIS influence implied (thinning in central Greenland) but in other places (the abstract and the conclusion) the influence of fast ice-stream flow for the full ice sheet is gestured at. I think the authors could refine the framing of this article to be focused just on the NEGIS influence for this catchment alone (or north-central Greenland, or demonstrate exactly where for the ice sheet the NEGIS has a dominant dynamic effect on thinning through the Holocene). After all, just the NEGIS catchment amounts to something like 12% of the entire ice-sheet volume or 1.1 meters so it’s a significant portion of the whole Greenland ice mass and thus there is a clear case that we should care about how it behaves dynamically (e.g. Mouginot et al. 2015). Alternatively, if the authors want to maintain the ice-sheet wide perspective, I think they need to show that NEGIS impacts more than just the NGRIP site during the Holocene. Perhaps a map of where in Greenland the ensemble shows divergent surface elevation histories versus where the ensemble members are in greater agreement about the elevation history would demonstrate this nicely, and allow the authors to better support statements such as “Our findings hint at a wider and more dramatic influence of the fast ice-stream flow on the whole ice-sheet behavior than previously expected” and “Such a strong dynamics-induced destabilization of the ice divide warns us about a possible extreme response of the ice-sheet interior to future warming ice loss.”\nIn summary, this is a strong and clear manuscript that advances the field and I hope to see this exciting work published soon with a few slight modifications.\nREFERENCES\nMouginot, J., Rignot, E., Scheuchl, B., Fenty, I., Khazendar, A., Morlighem, M., Buzzi, A., and Paden, J., 2015, Fast retreat of Zachariæ Isstrøm, northeast Greenland: Science, v. 350, p. 1357–1361, doi:10.1126/science.aac7111.\nLINE EDITS/QUESTIONS\nPlease note – the PDFs I received did not have line numbers so the page numbers below refer to the page of the PDF.\nMain text.\nPage 1.\nABSTRACT. “Our findings hint at a wider and more dramatic influence of the fast ice-stream flow on the whole ice-sheet behavior than previously expected.” I do not think you have shown fast ice-stream flow impacts the whole-ice sheet behavior, but you have very clearly and directly demonstrated it for the northeast catchment, from the oceanic margin all the way to the ice divide. To me this is impact enough, but if you want to keep your statement as-is then I would like to see a clear demonstration that the impacts of NEGIS reach the whole ice-sheet.\nMAIN\n“Punctual” – I’m not sure what you mean by this – could you use another word? Punctual sounds like it could mean abrupt, discontinuous, on-time? Do you mean well-dated, or something else?\n“~300 km long” check this length for NEGIS? I thought it was closer to 700km.\nPage 2. GHF – should this be geothermal heat flow or geothermal heat flux?\n“Penetration” – suggest a different word. Propagation?\n“Such a strong dynamics-induced destabilization of the ice divide warns us about a possible extreme response of the ice-sheet interior to future warming ice loss.” I think this is too strongly worded and I don’t think “destabilization” is the right word either. What about “Such a far-inland reach of dynamic thinning in northeast Greenland during a major retreat phase suggests that such processes may also amplify the response of this sector to future warming and ice loss.”\n“Several tens of meters” based on Figure 1 I think you mean kilometers?\nPage 3 “Penetrating” again – could be replaced with “more inland-propagating” or something else? Or is there a reason why you prefer the term “penetrating” here – I may be missing something but to me, NEGIS is not “penetrating” the ice sheet but instead the dynamic controls on the ice-sheet position are allowing the fast-flow to propagate farther and farther upstream. Note – later on page 5 you use “continental intrusion” – I think if you pick one phrase that you think best captures the phenomenon you are talking about it, and use it throughout, that would be helpful for the reader.\n“indicating that any other ice stream across the ice sheet has no significant involvement in the NGRIP thinning.” – reword: indicating that no other ice stream across the ice sheet has significant involvement in the NGRIP thinning. What does “other ice stream” refer to here – may be confusing to some readers who are not familiar with the difference between glacier versus ice stream and whether some ice flow features we identify as “glaciers” today may be more aptly considering “ice streams” under paleo flow conditions. What if you reworded this to instead say “indicating that there is no catchment except for NEGIS where ice flow has significant involvement in the NGRIP thinning” or something similar? Also, can you really say this? If you did the experiment to reduce the bedrock by one sigma beneath all of the major flow features in Greenland, would none of them contribute to the NGRIP thinning? Maybe just say “No similar correlation is found for any other GrIS basins in our experiments (Fig. 2S).”\nSentence “although even higher GHF values…” does not make sense to me. Do you mean “ Although even higher GHF values have been suggested in northeast Greenland as attributed to the passage of the Icelandic hotspot (Rogozhina et al., 2016), a GHF of just 57 allows our model…” ?\nCheck units of GHF – should be mW per m2 ?\nSpell “MacGregor”\nThe paragraph starting “Runs with a very fast ice stream…” Could you add one or two sentences describing any secondary controls on the NGRIP thinning rates from your ensemble, and comment briefly on what that means in terms of the influence of ice-stream dynamics on the thinning? For example in Fig 3S – z0 plays a less important but still influential role on inland thinning rates.\nPage 4\n“Here, bedrock uplift is due to the sole GrIS mass loss” reword: “Here, bedrock uplift is solely due to GrIS mass loss.”\n“Yet, the GIA signal from the Greenland peripheral ice sheets deglaciation” – what ice sheets are you referring to here – you mention the Innuitian, is that what you mean? Why not just say that – “Greenland peripheral ice sheets” is confusing.\n“crucial for a correct retreat” is it possible to quantify this following your method – e.g. the median skill score for the ensemble using each of the forcings? This would provide some quantitative foundation for what you mean by “correct retreat”\nPage 5. \nThe final sentence is really important and I suggest some slight changes to wording to really drive your main point home. Your logic here makes sense to me but it’s a very strong statement with big implications so I think it could be important to guide readers through your thought process here.\nI don’t really understand what you mean by “destabilization.” I don’t think you have shown that the Holocene warming destabilized northeast Greenland – what you have shown is that dynamic thinning in central Greenland is largely attributable to the development of the proto-NEGIS during the early/mid Holocene. But couldn’t this still be the “steady-state” response of the ice sheet to the transient climate forcing? To me “destabilize” implies a runaway or positive-feedback process which I don’t think you have shown here. I would avoid the word destabilize and try to stick closer to the physical processes that you have studied in this final paragraph. I suggest something like: Our results imply that a future northeast margin retreat stemming from increasing temperatures might similarly provoke tens of meters of dynamic thinning in central Greenland. Given that the present-day NEGIS velocity pattern is not captured by most ice-sheet models (Goelzer et al., 2018), our results suggest that inland dynamic thinning [in north-central Greenland] is underestimated by future projections and Greenland’s contribution to future sea level rise may be greater than currently anticipated.\nFIGURES\nFigure 1(b), Figure 4\nHow is “distance from present position at 79ºN” calculated? Along the central flowline for the measured velocity? Consider using “79ºN Glacier” everywhere you refer to this (e.g. Figure 4 caption) so that it is easily differentiated from “79ºN” the geographic line of latitude – for example in Figure 1 caption (e) where is this climate reconstruction from – the outlet of 79ºN glacier?\nFigure 2\nPlease label the colorbar. Consider adding 50 m/a contour from observed velocities on top of the “present day” modelled velocity field – I think that would emphasize that your model does a great job at reproducing NEGIS.\nConsider adding to at least on panel the location of Agassiz ice cap as this is referred to throughout the text.\nFigure 3\nDelete “for valid runs” – modify sentence to read something like: “ Thinning rate and shear stress are averaged between the time of Holocene maximum elevation and the present; coloured points follow the colour palette of panel (a); runs that do not meet scoring criteria are shown as grey dots; correlation coefficient and p-value for scored runs are also reported. The best run of the ensemble is shown as a blue triangle.”\nFigure 4\nWhy is the spread for the yellow ensemble so much greater than for the blue ensemble? Also these ensembles seem to be diverging already by 12ka but it is difficult to see what the reconstructions are doing at 12ka in the inset. Is it possible to show more of the reconstructions (maybe just to 14ka so it matches the ensemble plot)? Also consider labelling the gray line in (b) directly on the plot as “Agassiz ice cap temperature” or similar – because this does not correspond to plot (a) I found it confusing but I could see from the text why you wanted to include it.\nMETHODS\nPage 10. “It is only one (with another)” what does this mean? Consider “It is one of only two scenarios…” ?\nPage 11.\nCan you add a figure to the supplement showing how the bed elevation was changed (minus one sigma) within the NEGIS boundaries?\nI think this sentence should read: “absolute modelled-observed modern surface elevation DIFFERENCE at NGRIP” ?\nHow robust are your estimates of the % of thinning attributable to NEGIS propagating inland depending on which of these misfits you use? Is the best-fit scenario the best fit across all of the different misfits, or does one or a few of these misfits dominate the results? Can you say something about why these parameters were chosen to calculate the misfit?\nHow did you choose the criteria for runs to discard as NGRIP_PD<40m ? Assuming the above edit is right, maybe also just drive it home: “…runs are considered valid as long as the absolute elevation difference between modelled and observed elevations is less than 40 meters (e.g., NGRIP_PD<40m), so that 122…”\nPage 12\nFix: Reference climate FROM Box\nPage 13 – Karlsson reference should come before Le Meur reference I think?\nSUPPLEMENT\nFigure S5 – Can you add the blue triangle to each of these panels for your best simulation?"}, "author_response": "We thank the three anonymous reviewers for the valuable comments and suggestions for our manuscript. Revising the paper has required a large effort, specifically due to an improvement in the methodology that required additional model calibration and tuning. We believe the new set of experiments is improved in many aspects and the results are more robust. Yet it confirms the physical outcomes pointed out in the previous version of our manuscript. We hope that the new manuscript (attached) will be well received by the reviewers.\n\nWe agree that some critical points needed to be clarified and others improved. We took advantage of these criticisms to redo all experiments with a newer version of the ice-sheet model Yelmo (v1.801, available on Github under tagged version “ngrip_thinning”), which features several improvements such as:\n1. New calving law based on von Mises stress (Lipscomb et al., 2019), which is more suitable for the narrow, elongated Greenland outlet glaciers (Choi et al., 2018);\n2. New treatment of ice-free and partially ice-free points (temperature, viscosity and strain rate are extrapolated to those points) to improve stability;\n3. Improved viscosity routine in the DIVA solver;\n4. Several minor bugs are fixed.\n\nIn this revised manuscript we also include major improvements in the methodology:\n1. Improvement in the oceanic forcing. As suggested by Reviewer #2, we improved the ocean temperature anomaly in the northeast, making it closer to recent reconstructions from marine sediment cores in the Arctic. The temperature signal now mimics the continuous advection of subsurface warmer Atlantic Water (AW) into the shelf from early to mid Holocene, following the well-constrained evidence from sediment core analysis (Lloyd et al., 2023, Hansen et al., 2022, Davies et al., 2022, Syring et al., 2020). Thus, ocean temperature now increases between 13300 and 10000 years ago to represent the minimum period of presence of AW found in both Westwind Trough and Norske Trough. Then it declines slowly until 7500 years ago. A subsequent cooldown of subsurface waters due to the reduced intrusion of AW in the shelf (Lloyd et al., 2023, Davies et al., 2022, Hansen et al., 2022) is then represented by a constant low temperature anomaly until almost the end of the Holocene (the oceanic temperature anomaly recovers to present values during the last 100 years, see Fig. S11).\n2. Improved configuration of the paleo ice stream based on recent findings. We characterise the past northeast fast flow regime as a system of a northern, now-extinct paleo ice stream branch and a southern present-day-like ice stream activated during the Holocene, following Franke et al., 2022. There is large uncertainty around the timing of activation/deactivation of the ice-stream system in the northeast and the processes that might have caused it. Yet, significant restructuring of ice streams’ shape and strength is more likely during the initial phase of a deglaciation due to the large ice-sheet geometry reorganisations upon margin retreat. To reproduce such change in configuration during the deglaciation, we start to deactivate the paleo NEGIS-like ice stream and activate the present-NEGIS at 8 kyr ago by changing the local basal friction. This is also supported by the recent hypothesis of the NEGIS showing a similar configuration to today since 8000 years ago (Gerber et al., 2021).\n3. Inclusion of additional perturbation parameters to better cover the ensemble phase space. Seven other perturbation parameters have been added to the simulation ensemble: 1. the ice-ocean heat flux parameter (i.e. κ of Eq. 7), to test the ice-sheet sensitivity to the ocean forcing during the deglaciation; 2. the enhancement factor Ef to explore different dynamic states of the ice sheet due to ice anisotropy; 3. the lithosphere relaxation time τ and 4. the effective elastic thickness of the lithosphere to investigate the role of the isostatic adjustment; 5. the calving rate parameter k of the von-Mises-like calving law (Lipscomb et al., 2019); 6. the scaling factor f for past temperature anomalies in the precipitation-temperature relationship (similar to β from eq. 8 of Badgeley et al., 2020) and 7. the additional scaling factor f to be applied to the NEGIS basal friction parameter f during the last part of the Holocene to reproduce the switch from the paleo-NEGIS to the present-NEGIS (see Figure 1). The large ensemble now includes a total of 3000 simulations.\n4. Inclusion of additional data used to evaluate the model: the ensemble is now evaluated against a total of thirteen individual scores based on paleo and present model-data misfits. In addition to previous misfits, we add now: RMSE of modelled to observed ice-covered area at the LGM (AREA_LGM, data from Lecavalier et al., 2014), and RMSE of modelled-observed deglaciation timings computed and averaged 1. at four sediment core locations from the northeast continental shelf (SEDIM_DEG, data from Hansen et al., 2022, Lloyd et al., 2023, Davies et al., 2022, Syring et al., 2020), 2. at three locations in the outer shelf as calculated from 10Be exposures dates from morain/boulder samples (OUT_MORAIN_DEG, data from Larsen et al., 2018) and 3. at other six locations closer to the present NEGIS margin (INN_MORAIN_DEG, data from Larsen et al., 2018). We also introduce a last binary control on the upstream retreat of 79N glacier: simulations that retreated more than 70 km upstream will present a maximum error and vice versa (see Methods for more details).\nSpecific reviewers’ questions/comments have been answered point by point below."}]

Decision letter