Greenland ice mass loss during the Younger Dryas driven by Atlantic Meridional Overturning Circulation feedbacks

Understanding feedbacks between the Greenland Ice Sheet (GrIS) and the Atlantic Meridional Overturning Circulation (AMOC) is crucial for reducing uncertainties over future sea level and ocean circulation change. Reconstructing past GrIS dynamics can extend the observational record and elucidate mechanisms that operate on multi-decadal timescales. We report a highly-constrained last glacial vertical profile of cosmogenic isotope exposure ages from Sermilik Fjord, a marine-terminating ice stream in the southeast sector of the GrIS. Our reconstruction reveals substantial ice-mass loss throughout the Younger Dryas (12.9-11.7 ka), a period of marked atmospheric and sea-surface cooling. Earth-system modelling reveals that southern GrIS marginal melt was likely driven by strengthening of the Irminger Current at depth due to a weakening of the AMOC during the Younger Dryas. This change in North Atlantic circulation appears to have drawn warm subsurface waters to southeast Greenland despite markedly cooler sea surface temperatures, enhancing thermal erosion at the grounding lines of palaeo ice-streams, supporting interpretation of regional marine-sediment cores. Given current rates of GrIS meltwater input into the North Atlantic and the vulnerability of major ice streams to water temperature changes at the grounding line, this mechanism has important implications for future AMOC changes and northern hemisphere heat transport.

In today's warming climate, the persistent North Atlantic cooling anomaly off the southeast coast of Greenland is thought to be caused by the accelerating input of meltwater from the Greenland Ice Sheet (GrIS) and considered an indicator of a weakening of the Atlantic Meridional Overturning Circulation (AMOC) 1,2 . Modelling studies suggest that increased freshwater input from Greenland has the potential to slow down AMOC in the future 3 , thus substantially weakening heat transport to the North Atlantic and potentially driving positive ice sheet-ocean feedbacks 4 .
Unfortunately, contemporary observations of the GrIS and North Atlantic are too short (<100 years) to inform on our understanding of long-term trajectories and feedbacks over multi-decadal to centennial timescales. Periods of large-scale rapid climate reorganisation in the recent geological past offer considerable potential for understanding ice-sheet-ocean dynamics, thereby extending the observational window. The Younger Dryas (YD, broadly coincident to Greenland Stadial 1, ~12.9-11.7 ka 5 ) was a period of sustained Northern Hemisphere atmospheric cooling of up to 10 °C 6,7 that is considered to have been initiated by freshwater input into the Arctic Ocean and/or North Atlantic from the expanded Laurentide Ice Sheet 8 . In contrast to widespread terrestrial glacier expansion in Europe during the YD 9 , ice streams and local glaciers in southern Greenland may have reacted in antiphase to atmospheric temperature changes 9 ; local terrestrial glaciers near Scorsby Sund, east Greenland, show evidence of Younger Dryas standstill 10   magnitude of thinning as well as the rate of terminus or lateral retreat. Here we present the results of new in situ cosmogenic isotope analysis of vertical transects which we combine with published lateral transects up Sermilik Fjord 21 and examine possible drivers of the dynamic changes we discover using transient experiments performed with an Earth system model.

Results
Cosmogenic isotope analysis. Glacial erratics and bedrock were sampled along two individual vertical transects (one from an island in the north of the fjord, and another on the east wall), between ~40-430 m above present day sea level in Sermilik Fjord, as well as two additional discrete sites at lower elevations (Fig. 1). Samples were analysed for the in situ cosmogenic nuclide 10 Be, and exposure ages calculated using the Northeast North American production rate (see Methods; Tables S1 and S2) 25 . The highest samples imply that the ice sheet surface was over 400 m a.s.l before 13 ka, inferring the presence of a fully-grounded ice stream at this time, given the depth of the fjord and the flotation point of ice.
To constrain the age-elevation history and produce a thinning model, we applied Bayesian modelling to the cosmogenic ages, in combination with two basal organic 14 C ages from elevated lakes on the eastern side of Sermilik 26 , using OxCal 4.1 27 (see Methods). This analysis reveals that ice stream drawdown was underway by 13.0 ± 0.4 ka (1 s.d.), and continued into the Holocene until 10.3 ± 0.2 ka (1 s.d.), prior to its retreat upstream, indicating the ice ungrounded and underwent surface lowering of ~376 m. This continued throughout the Younger Dryas, with ~280 m of the of the thinning occurring during this period ( Fig. 2A). The timing of this significant mass loss agrees with the implied early YD retreat of Jakobshavn Isbrae 11,14,28 as well as with the southeast Greenland freshwater spike at 12.8-12 ka 17 (Fig. 2B), and suggests that late YD retreat of the Kangerlussuaq Fjord ice stream (east Greenland, previously known as Kangerdlugssuaq Fjord, and distinct from Kangerlussuaq Fjord in west Greenland) 12 may have been preceded by a similar thinning. Crucially, the well-constrained trajectory we demonstrate here confirms previous suggestions of GrIs Younger Dryas ice retreat [9][10][11]14,28 . Our results appear to represent a coherent pattern of major ice stream downwasting and significant GrIS mass loss through a period of sustained atmospheric cooling in this important sector of the ice sheet ( Fig. 2A).

Exploring the drivers of Younger Dryas melt.
Our observations raise the question: what could have driven substantial retreat of a major outlet glacier of the GrIS despite a marked drop in atmospheric temperature at the Younger Dryas? Rising sea levels can cause ice-stream instability and retreat; however, although global sea level was rising prior to and during the YD 29 , Greenland was undergoing substantial isostatic uplift following post-LGM ice-sheet retreat 30 , which may have proved an additional stabilising factor. Heightened seasonality (and therefore length of melt season) caused by the 11 ka peak in obliquity 31 may also have contributed, but the YD expansion of a local terrestrial ice cap in the far north of Greenland 32 implies that this effect was not strong enough to drive major ice stream retreat.
An alternative mechanism for negative mass balance could be reduction in precipitation in the accumulation zones of the GrIS as a result of sea ice expansion in the Northern Atlantic 33,34 . Reconstructions of Arctic sea ice at the time are limited and contradictory 35 , but modelling suggests that winter sea ice at 12 ka was expansive 36 (Fig. S3). The GISP2 ice core shows a marked decrease in accumulation throughout the YD 37,38 , consistent with the idea of precipitation starvation (Fig. S4), but equally low levels of precipitation occurred during the LGM when the GrIS was expanding 39 . Furthermore, significant Arctic sea ice expansion during the Little Ice Ageincluding across the GrIS precipitation source of the Labrador Sea 40 -apparently did nothing to inhibit glacial expansion 41 , which strongly suggests that sea-ice induced precipitation starvation alone could not have driven the YD mass loss.
Although sea surface temperatures in southern Greenland were markedly depressed during the Younger Dryas 42 , the contrasting mass trends between the marine-terminating ice streams of the GrIS 11 and terrestrial glaciers elsewhere in the Northern Hemisphere 9 suggests oceanographic changes may have played a key role. To investigate potential feedback mechanisms between the GrIS and the North Atlantic oceanic response to a weakened AMOC, we analyse model outputs from a transient deglacial experiment performed with the Earth system model LOVECLIM 43 , where a weakening of the AMOC during the YD is simulated by a 0.25 Sv freshwater flux applied in the Arctic Ocean (175°W-95°W, 67°N-83°N) between 13 and 12.2 ka B.P (see Methods). Despite a marked simulated North Atlantic sea surface temperature decline and sea ice expansion (Fig. S3) during the early part of the YD, a sustained sub-surface ocean warming is simulated at mid-depth along the entire south-eastern Greenland margin (Figs 3 and S5). Before the YD event, North Atlantic Deep Water (NADW) formation at the confluence of the Greenland and Norwegian Seas (75-80°N) drives a subsurface poleward advection of relatively warm water to that region, as well as the southward advection of cold and dense NADW along the eastern margin of Greenland (Figs 1 and 3). Meltwater input into the Arctic at the beginning of the YD leads to the shutdown of NADW production at 12.8 ka (Fig. S3), which weakens the southward advection of cold subsurface waters through the Denmark overflow, and the advection of warm water from the Norwegian Atlantic Current (NAC) to the Norwegian Sea due to a reduced meridional density gradient in the North Atlantic and changes in the wind stress curl. This significantly reduces the strong east-west temperature gradient at depth, supporting regional palaeoceanographic reconstructions 17 (Figs 2 and 3). This is a robust feature of North Atlantic Themohaline circulation in both earth system models of intermediate complexity and fully coupled atmosphere-ocean general circulation models 44 . Furthermore, within 100 years of the meltwater input, our results suggest the subsurface cyclonic circulation south of Iceland intensified, strengthening the Irminger Current and enhancing the subsurface advection of warm waters to the south of Greenland, as evidenced by marine micropalentological investigations from the Kangerlussuaq 17 (Figs 1, 3 and S5). As a consequence, ocean temperatures below ~350 m increase markedly from 14 ka to 12.5 ka, with positive ocean temperature anomalies of ~3 °C and ~0.4 °C southeast  10 Be exposure ages combined into a single event by C_Combine) and two radiocarbon ages from basal lake sediments (AAR-1542 and K-0638) 26 . Dark blue envelope to 1 s.d., pale blue to 2.s.d. Light grey cumulative probability functions show raw 10 Be exposure ages after outliers were omitted by C_Combine function (Table S2). Dark grey probability density functions show Bayesian modelled ages (Table S4) 55°W, 60-68°N) of Greenland respectively, a depth comparable to those of cross-shelf troughs of the major ice streams (Figs 1, 3 and S5).

Discussion
We propose that ice-mass loss from southern Greenland across the YD was driven by entrainment of warm sub-surface water in the cross-shelf troughs, leading to thermal erosion at the grounding line of southern GrIS marine-terminating outlets such as Sermilik and Kangerlussuaq 12 , consistent with the hypothesis that ocean circulation was more important than air temperature change in this region during the YD 15,17,19 . Increased stratification in the northeast Atlantic following the cessation of NADW formation increased the marine reservoir age in this region 18 , and the anomalously low marine reservoir age observed southeast of Greenland throughout the YD 16 could be explained by such subsurface current reorganisations, supporting the interpretation of isotopic  records from marine sediment cores from the North Atlantic off the SE GrIs 19 . Furthermore, micropalentological evidence suggests that recirculating intermediate depth North Atlantic water, advected via the Irminger Current, may have entered the Kangerlussuaq Trough during the YD chronozone, providing direct evidence for the subsurface current reorganisation we propose. Our model results suggest that the initial input of YD glacial meltwater into the North Atlantic from the Arctic and the Laurentide Ice Sheet triggered a positive ice sheet-ocean feedback, whereby a weakened AMOC -and the associated subsurface warming off the south Greenland shelf -helped sustain freshwater input into the North Atlantic, thus driving and potentially extending the Northern Hemisphere YD cooling. In our model experiments we do not add meltwater from the GrIs, however, as noted in other studies this may have enhanced this feedback mechanism along the eastern margin of the GrIs 17,19 .
Such a mechanism could have significant implications for future ice-sheet dynamics. A number of major outlets of the present-day southern GrIS have grounding lines at depths of 500-800 m b.s.l., including Helheim (which feeds into Sermilik Fjord), Kangerlussuaq (east Greenland) and Jakobshavn Isbrae 45 . The subsurface current reorganisation highlighted by our modelling study impacts depths of 150-900 m, which encompasses the grounding line depths near the margins of all these GrIS outlets 45 . If the contemporary observed AMOC weakening exacerbated by Greenland meltwater 1 continues, the mechanism highlighted here may lead to warmer Atlantic Water being drawn towards southern Greenland in the future, significantly increasing subsurface ocean temperatures at the grounding lines of major outlets. This could drive enhanced thermal erosion of outlet glaciers, and further freshwater input into the North Atlantic over and above that caused by current atmospheric temperature changes. Given that these glaciers are already projected to contribute 19-30 mm to sea-level rise by 2200 (ref. 46 ), it is vital that this feedback is included in ice-sheet modelling studies to constrain projections of future sea-level rise and northern hemisphere heat transport.
Our study reveals substantial Greenland ice stream collapse in the Younger Dryas during a period of considerable atmospheric and sea-surface cooling in the North Atlantic. Our finding supports the interpretation of marine sediment cores off the SE GrIS 17,19 , and suggest that ice sheet mass loss was driven by sub-surface warming driven by altered deep ocean convection in the Norwegian Sea and a subsequently weakened AMOC. This positive feedback mechanism could have significant implications for present-day GrIS/AMOC interactions 1 and highlights complex and non-linear ice-sheet-atmosphere-ocean dynamics that must be incorporated into future projections.

Methods
Cosmogenic sampling strategy. Rock samples were collected for dating by in situ cosmogenic isotope analysis from ice-free areas along Sermilik Fjord, southeast Greenland during fieldwork in June/July 2009. Sample sites were chosen to span the greatest possible range within the time available for fieldwork, both geographically and attitudinally. Geological maps were consulted prior to sampling to ascertain ideal sampling areas for the quartz-rich rock types needed for 10 Be measurement. Accessibility was also a key control when selecting sample areas; sites were accessed by boat and foot, so terrain and local conditions had to be considered. Where possible, samples were taken in bedrock-erratic pairs or triplets (i.e. bedrock-erratic or bedrock-erratic-erratic), to allow for any inconsistencies in glacial history to be detected. Samples were selected based on vegetation cover, potential for snow cover (samples in raised and exposed areas, less prone to snow accumulation were preferred) and their stability (e.g. erratic cobbles were sampled from areas outside of drift accumulation, and with weathering pedestals beneath them; erratics boulders were not sampled if they seemed precariously perched or liable to down-slope movement; bedrock was not taken from areas of higher weathering). Erratic samples were taken either as whole cobbles (up to ~20 cm) or as subsample of larger erratic boulders (up to 1 m), in which case the top surface was taken. Bedrock samples were taken from glacially streamlined surfaces. Samples were taken up to 12.5 cm thick, although thinner samples were preferred where possible. Samples from lower elevations were above the local marine limit (~40 m a.s.l. (ref. 47 )). Samples were taken along two vertical transects, from an island in the northern sector of the fjord (Island transect; samples SF0915-30) and from the east wall of the fjord (East wall transect; samples SF0936-55), as well as from two additional sites (West wall; SF0962-64 and Mouth; SF0903-04); see Fig. 1.
Samples were taken by hand using a four-pound lump hammer and chisel; all samples were labelled and photographed before being securely wrapped and relabelled, and detailed field notes on lithology, vegetation/ snow cover and local geomorphology were made. Inclination measurements of the surrounding topography were measured to ascertain topographic shielding, using a handheld clinometer, every thirty degrees on the azimuth, with additional measurements taken if this missed significant features. Latitude, longitude and elevation readings were taken using a handheld GPS device (accurate to ±10 m in elevation).

10
Be cosmogenic isotope extraction. Beryllium purification chemistry was undertaken at the University of Exeter Cosmogenic Isotope Laboratories. Chemistry followed protocols adapted from Ivy-Ochs 48 and Ditchburn and Whitehead 49 , modified by C. Fogwill and T. Barrows/J. Stone, respectively. Samples were then ignited to Be oxide, mixed with Nb, loaded into copper cathodes, and 10 Be/ 9 Be measured by accelerator mass spectrometry (AMS) at Scottish Universities Environmental Research Centre (SUERC). AMS measurements were normalized to the NIST 27900 standard (NIST SRM4325 standard material with an assumed isotope ratio of 2.79 × 10 −11 ). Procedural blank ratios were 5.16 × 10 −15 , 5.48 × 10 −15 , 3.65 × 10 −15 and 1.83 × 10 −15 . Blank ratios were 0.5-10.7% of sample ratios (mean: 3.4%). The isotopic ratio of the procedural blanks for each batch of samples was subtracted from the sample ratio, with errors in sample and background propagated in quadrature. See Table S1 for sample details.  51,52 ) and Be isotope ratio standardization of Nishiizumi et al. 53 . Age calculation assumed standard pressure, zero erosion and sample density of 2.62 g.cm 2 . Samples were scaled using the time-dependent Lal/Stone scaling scheme 54,55 , following similar studies and the recent analysis by CRONUS-Earth which shows the Lal scheme to provide a better fit to calibration date than the other, neutron-monitor based, schemes 56 ; using other scaling schemes produces exposure ages which differ by 1-4% (~120-480 yrs). For our analysis, we use the production rate from the North American (NENA) calibration set (reference production rate at SLHL: 3.87 ± 0.19 atoms g −1 a −1 using the Lal/Stone time-dependent scaling scheme) 25 , in line with other studies in the area, but also because of the robustness of this calibration set for studies in Greenland 55,57 . NENA exposure ages for all samples are shown in Figure S1 and Table S2. Samples were also calculated using the global production rate -which is increasingly being shown as inaccurate in higher latitudes owing to its reliance on low-latitude, high-altitude calibration sites -which produced results systematically 11.37% older than those using the NENA production rate. Other regional production rates were also used: Northern Norway 58 (mean = 2.6%/318 yr older than NENA); Western Norway 59 (mean = 4.9%/560 yr younger than NENA); and Baffin Island/Arctic 60 (mean = 1.5%/173 yr younger than NENA); see Table S3. We believe that the NENA is the most appropriate for our study site; however, the NENA, Northern Norway and Baffin Bay rates all produce ages within the range of the internal errors. The Western Norway rate is based on data from sites at ~60°N, significantly to the south of Sermilik Fjord. All regional production rates produce ages spanning the YD, meaning our interpretation of the results is not affected by the choice. Density was assumed to be 2.62 g cm −3 , and we made no correction for inheritance, vegetation or seasonal snow cover, all of which are assumed to be zero or minimal owing to our sampling methods and geomorphic assessment. We assume an erosion rate of zero; using an erosion rate of 2 × 10 −4 cm yr −1 increases ages by ~2% 61 .
Data analysis. Samples were taken in two vertical transects along the fjord (Island transect, samples SF0915-30, and East Wall transect, samples SF0936-55) as well as two additional discrete sampling sites from lower elevations (West Wall site, samples SF0962-64, and Mouth site, samples SF0903/04); comparing the box and whisker plots of the dates from both transects (Fig. S2) shows they have a similar spread of ages, suggesting that the timing and rate of thinning were indistinguishable along the fjord. As such, we are able to combine the samples to form one thinning profile of the ice stream. To test for internal consistency of ages and account for anomalous exposure ages (owing to inheritance, shielding etc.) we first analysed the dataset using the C_Combine function in the chronological software OxCal 4.1 (ref. 62 ); this performs a chi-square test on bedrock-erratic pairs/triplets, highlighting and discarding the outlier in those pairs/triplets that failed. The function then combines the remaining individual samples within each pair/triplet into one event using a cumulative probability density function. This process led to the exclusion of five samples (Table S2): four bedrock, of which two were anomalously old (SF0924/SF0941) and two anomalously young (SF0916/SF0947); and one erratic, which was anomalously old (SF0903). The 17 samples remaining were re-analysed, together with two radiocarbon dates from basal organic lake sediments on the east wall of Sermilik fjord ( Fig. 1 (ref. 26 )) using the Bayesian modelling of the T_Scaled analysis function within OxCal 4.1 (ref. 27 ). Using Bayes theorem with a Poisson process deposition model, the algorithms employed sample possible solutions with a probability that is the product of the prior and posterior likelihood probabilities. Taking into account the deposition model (i.e. those samples at higher altitudes must be the oldest and vice versa) and the actual 10 Be measurements, the posterior probability densities quantify the most likely age distributions; the outlier option was used to detect ages that fall outside the calibration model for each group and, where necessary, down-weight their contribution to the final age estimates. In this analysis, none of the remaining samples were considered outliers. This analysis generates the thinning profile of the ice stream shown in Fig. 2A. The Bayesian approach, together with the inclusion of the two radiocarbon dates from basal organic lake sediments, allows the errors on the age envelope to be wrapped down from the typically high errors of the 10 Be exposure ages, particularly at the lower end of the profile. Of the seventeen 10 Be exposure ages included in the model, only five were from unpaired events; of these five, only two were bedrock, one of which (SF0925) was younger than its discarded erratic pair, and therefore not experiencing inheritance. Although the remaining unpaired bedrock sample (SF0915, whose erratic pair was anomalously young when compared to other samples of similar heights, and therefore has the potential of inherited 10 Be) is used to help constrain the upper end of the profile, numerous other samples at this elevation and below (e.g. SF0936/37 (bedrock-erratic), SF0943 (erratic), SF0920/21 (bedrock-erratic)) have 10 Be exposure age envelopes well within or prior to the onset of the YD. As such, we can be confident that our thinning trajectory is not significantly skewed towards older ages by any potential 10 Be inheritance within bedrock samples. The ages resulting from the Bayesian modelling in OxCal are shown in Table S4. Our analyses show that deglaciation of the highest sample (423 m a.s.l.) commenced at 13.0 ± 0.4 ka (1.s.d.) (13.05 ± 0.75 ka to 2 s.d.) and that the ice stream thinned to an elevation of 47 m a.s.l. by 10.3 ± 0.2 ka to 1.s.d. (10.95 ± 0.65 ka to 2.s.d.). This is equivalent to an average modelled thinning rate of ~0.14 m yr −1 , or ~0.23 m yr −1 during the Younger Dryas alone (i.e. that due only to non-atmospheric drivers). This is on a comparable order of magnitude to the thinning seen in the fastest-flowing outlet glaciers of Greenland today (~0.84 m yr −1 in glaciers flowing faster than 100 m yr −1 (ref. 63 )); it is however important to note the difference in observational timescales when making these comparisons.

Climate modelling. A transient experiment of the last deglaciation was performed with the Earth
System model LOVECLIM (ref. 43 ). LOVECLIM comprises an ocean general circulation model and a thermodynamic-dynamic sea ice model, with an horizontal resolution of 3° × 3° and 20 vertical levels, coupled to a spectral T21quasi-geostrophic atmospheric model 64 . Starting from the background conditions of the Last Glacial Maximum, the model is forced with the time-varying evolution of solar insolation 31 , ice-sheet topography 65 , high latitude albedo and atmospheric CO 2 (ref. 66