Fennoscandian freshwater control on Greenland hydroclimate shifts at the onset of the Younger Dryas

Sources and timing of freshwater forcing relative to hydroclimate shifts recorded in Greenland ice cores at the onset of Younger Dryas, ∼12,800 years ago, remain speculative. Here we show that progressive Fennoscandian Ice Sheet (FIS) melting 13,100–12,880 years ago generates a hydroclimate dipole with drier–colder conditions in Northern Europe and wetter–warmer conditions in Greenland. FIS melting culminates 12,880 years ago synchronously with the start of Greenland Stadial 1 and a large-scale hydroclimate transition lasting ∼180 years. Transient climate model simulations forced with FIS freshwater reproduce the initial hydroclimate dipole through sea-ice feedbacks in the Nordic Seas. The transition is attributed to the export of excess sea ice to the subpolar North Atlantic and a subsequent southward shift of the westerly winds. We suggest that North Atlantic hydroclimate sensitivity to FIS freshwater can explain the pace and sign of shifts recorded in Greenland at the climate transition into the Younger Dryas.

C dates for Hässeldala based on selected terrestrial plant remains and used to construct the composite age-depth model. Blue samples refer to dates transferred from Core 2 to Core 5.

Sample depth (cm)
Thickness error (cm) Sample ID Material analysed 14 C age (yr) 14 C error (1 sigma) Used in the age model     14 C dates constraining the age of the first drainage of the Baltic Ice Lake. Ages refer to lake isolations from the sea, which indicate timing of deglaciation west of the drainage outlet area at Mt. Billingen,near 11,000 radiocarbon years BP. Two dates refers to isolation owing to concomitant lowering of the Baltic Ice Lake in Blekinge. Dates were combined using the 'R-combine' function in OxCal4. 2 (ref. 4) after calibration with the IntCal13 calibration curve (ref. 2). Only ages considered reliable by the original study authors and directly related to lake isolations were used. The temporal correspondence between the supposed first drainage of the Baltic Ice Lake and the regional pollen stratigraphy is well constrained and thus offers an indirect time marker in other regional sedimentary sequences (see text). 14 C dates from other sequences used to constrain timing of deglaciation at Mt. Billingen through pollen stratigraphic evidence were listed and noted with the entry "Pollen time marker". One 14 C date constraining the timing of the Baltic Ice Lake water-level fall is also listed.

Supplementary Discussion
Sources of n-alkanes in sediments from Hässeldala Lake and suitability of δD

records for hydrological investigations
In order to use the δD composition of sedimentary n-alkanes as a proxy for hydroclimatic conditions (e.g. evapotranspiration), it is important to first assess the possible aquatic versus terrestrial origin of these compounds in their specific lake environment 57 . In our study we employed the δD composition of n-C 21 , and n-C 27 , n- This interpretation is also confirmed by geochemical data. For instance, C/N ratios in preferred to only use n-C 21 in order to minimize the potential bias from higher terrestrial plants, which also produce, to some extent, mid-chain n-alkanes.
Conversely, long-chain n-alkanes, like the n-C 27 , n-C 29 and n-C 31 alkanes, are main constituents of leaf waxes of higher terrestrial plants 21 . In Hässeldala sediments, we assigned these compounds to a mixture of leaves from Betula nana and Salix polaris.
The n-C 27 is generally synthesized by these two plant species 22  In summary, we argue that the short-chain and long-chain n-alkanes discussed here originate from separate aquatic and terrestrial sources, respectively. This is also supported by the lack of any significant correlation between the δD values of n-C 21 and n-C 27 , n-C 29 , n-C 31 alkanes, which gives us confidence that the ΔδD long-short represents the actual ΔδD terr-aq (ref. 11).
It should be noted that the δD records of long-chain n-alkanes exhibit a lower variability than the δD records of short-chain n-alkanes ( Supplementary Fig. 2) and that the ΔδD terr-aq appears thus mainly driven by changes in δD aq . We argue that the freshwater forcing in the North Sea caused marked declines in δD aq values at Hässeldala during the Late Allerød and the Younger Dryas pollen zone, respectively, but also caused a hampered moisture flux from the near marine source. The freshwater forcing thus indirectly induces drought conditions at Hässeldala.
Mechanistically, the δD terr responds to changes in evapotranspiration with less negative values relating to greater evapotranspiration under dry atmospheric conditions. Therefore, in our record the expected variability in δD terr predicted by negative shifts in δD aq values is offset by a concurrent rise in evapotranspirationand accordingly in ΔδD terr-aq values.
The effect of evaporative deuterium enrichment owing to enhanced evapotranspiration is particularly pronounced in shrub waxes as compared to other functional groups such as trees or grasses 24 due to the drought tolerance and high water-use efficiency in shrubs (e.g. ref. 25). However, the apparent low variability of δD terr relative to δD aq could have also been caused, to some extent, by admixture of older organic material, which resulted in a longer temporal integration of the stable isotopic composition of the n-alkane pool (see further discussion below).
Nonetheless, we consider the effect of enhanced evapotranspiration as the most likely primary factor subduing the variability of δD terr at Hässeldala. In fact, we observe that the overall pattern of the evapotranspiration history predicted by the

Effect of vegetation and local hydrology on ΔδD terr-aq
The ΔδD terr-aq values are directly related to the isotopic fractionation between plant source water and leaf wax lipids, which is referred to as net or apparent fractionation 74 . The net fractionation values for C3 grasses (no C4 plants were present in southern Sweden) are generally lower than those for C3 shrubs 24 Fig. 3). This observation is corroborated by other independent local and regional vegetation reconstructions (e.g. ref. 5,12,14) and implies that changes in vegetation cover can be dismissed as a potential contributing factor to the shifts in δD terr values recorded at Hässeldala shortly prior to the Younger Dryas pollen zone.
We also dismiss the possibility that the trend observed in the δD aq record was caused by changes in net fractionation values due to a shift in the aquatic plant community producing n-C 21 alkanes. Vegetation reconstructions from nearby lakes 13 indicate that the aquatic plant community remained generally unchanged throughout the Late Glacial. This is also supported by aquatic environmental studies from the same sedimentary core at Hässeldala 27,29 , which show no significant change in the aquatic habitat structure throughout the Late Allerød and the first half of Younger Dryas pollen zones.
Moreover, we can also dismiss the contribution of local hydrologic changes as contributing factors driving the δD aq shifts at Hässeldala during this interval. For instance, under dry climatic conditions, such as during the Younger Dryas pollen zone (Fig. 2), evaporative deuterium enrichment on lake water δD may have occurred. However, this effect would have worked in the opposite direction as the one observed (δD depletion). Thus, if evaporation had played a role, it would mean that the signal recorded in both the δD aq and ΔδD terr-aq would underestimate the hydrological conditions. Similarly, this can be argued for the relatively dry phase that occurred during the Late Allerød pollen zone (Fig. 2), which has been identified in other independent proxy records from the same sedimentary core 27 .

Age of terrestrial long-chain n-alkanes
One important consideration when interpreting isotopic records from terrestrial nalkanes is the possibility of smoothing and/or a delay in the recorded signal caused by pre-aging on land and time needed for transport to and deposition in the lake basin. However, most of the plant waxes are expected to have arrived either directly after abrasion from the leaves, or simply from leaf litter in autumn. Earlier studies on lakes with small catchment areas have indicated that terrestrial n-alkanes, or other soil-derived compounds, typically have radiocarbon ages that are less than a few hundred years older than those of the sediments in which they are embedded (e.g. ref. 30,31,32).
Older leaf wax ages have been observed in areas with higher erosion rates where a large stock of ancient soil carbon and/or large watersheds contributes carbon inputs (e.g. ref. 33,34). For smaller catchments the most viable explanation is that a small fraction of the total carbon pool has a significantly lower 14 C content, thereby exerting a large effect on the average 14 C content, while the vast majority has been deposited within a few years after biosynthesis. Under such a scenario, the stable isotopic composition of the total pool will remain relatively unaffected by admixture of older organic material, especially if their stable isotopic compositions are similar.
Hässeldala has a small catchment where transport time from the catchment to the basin is expected to have been short. Furthermore, from the age-depth relationship we observe no change in accumulation that could indicate relatively higher erosion rates in the catchment area during the period under study (Supplementary Fig. 7).
The only potential indication of an increase in erosional input of older organic material from land occurs at 12800 yr BP, when the total organic content and the abundance of n-alkanes rise ( Supplementary Fig. 1, 2). This, however, takes place 300 years after the δD aq values start to decline and the ΔδD terr-aq values start to increase ( Supplementary Fig. 2). Therefore, we dismiss the contribution of admixing of older n-alkanes on δD terr during Greenland Interstadial 1a.
In summary, we argue that the δD terr signal is virtually synchronous with that of the aquatic lipids. Even if there has been some smoothing and attenuation of the δD terr , this would not make any significant difference until the start of the regional Younger Dryas, as the δD terr shows significant changes only after the onset of the Younger Dryas ( Supplementary Fig. 2). A shift of a few hundred years of the δD terr record during the Younger Dryas and Holocene, together with a rise in amplitude to account for any smoothing that may have occurred, would also not change our conclusions on the late Allerød pollen zone hydroclimate shifts recorded at Hässeldala. This would also make the observed contrast between δD terr and δD aq even larger during the Younger Dryas.
In addition, we note that the δD terr values generally change in phase with δD aq and in particular at the regional Allerød-Younger Dryas pollen boundary ( Supplementary   Fig. 2), suggesting that the δD terr record is in fact synchronous with the other proxy records from Hässeldala.

Comparison with δD records from Meerfelder Maar
By comparing our records with a recently published high-resolution isotope reconstruction from Meerfelder Maar, Western Germany 35 , we can better decipher the European-scale regional hydroclimate expressions around the onset of Greenland Stadial 1 and the regional Younger Dryas pollen zone.
We  Fig. 4). This pattern is also evident in the ΔδD terr-aq profiles, which suggests drier conditions at Hässeldala as opposed to more humid conditions at Meerfelder Maar.
The observed opposite patterns of hydroclimate conditions are in agreement with the simulated sea-level pressure field under FIS freshwater forcing presented in

The first drainage of the Baltic Ice Lake
The deglaciation history of the Fennoscandian Ice Sheet (FIS) in southern Sweden is well known and documented (e.g. ref. 14,37,38). During the Late Allerød pollen zone, at ~13000 yr BP, the ice front was located at the southern margin of the southcentral Swedish low-land area. Here, owing to the isostatic depression, the land surface was considerably lower than the ice-free lowermost threshold region further south in the Öresund Strait. Since the Baltic Ice Lake (BIL) was dammed up 5-10 m a.s.l. prior to ~13000 yr BP, rapid deglaciation of this southern margin at the water divide near Mt. Billingen during the Late Allerød pollen zone generated a spillway system that connected the BIL to the sea in the west. As such, the rapid retreat of the FIS resulted in a 5-10 m lowering of the BIL also referred to as the first drainage of the BIL 37 .
Although this drainage hypothesis has been debated, new reconstructions now suggest that a catastrophic outflow of freshwater actually took place near Mt.
Billingen 39 . Unfortunately, these reconstructions lack precise chronological constraints that allow to confidently determining the age of the drainage event.
Shore displacement curves for the Mt. Billingen area provide a detailed framework for the timing of the deglaciation near the outlet and the freshwater connection to the sea (e.g. ref. 8,9). However, a proper effort to systematically calibrate the published radiocarbon dates has never been undertaken, especially since the advent of Bayesian inference to chronological modelling.

Air parcel back trajectories
The trajectory data is an output from the trajectory FLEXPART 40 . FLEXPART is widely used to trace transport of moisture and other chemical species (e.g. ref. 41,42).
Settings and specifications of the simulation used in this study can be found in Viste and Sorteberg 43 . The model was run in a global domain mode filled with 1 million particles and integrated for the time period 1998-2008. Three dimensional wind, temperature and specific humidity from ERA-interim were used to drive the model.
To trace the sources of air masses and moisture at Hässeldala, we have defined a target region of 1° by 1° centred over the site (Fig. 1a). We identified the particles arriving to the target region every six hours and traced the particles backward in time for six days to locate their origin. We only traced particles from within the troposphere, assuming that the moisture content of particles above this layer is negligible 44 . The position, mass and specific humidity of these particles were saved during back tracking. To map the starting points of the particles and their path to the target region, the data of the traced particles was re-gridded into 2° by 2° grid boxes.
By doing this, the analysis provides the flux of air and water masses that pass through each grid box to reach the target region.

Sediment core alignment
The isotopic data presented in this study refer to Core 5 (ref. 45), whereas the chironomid-based summer temperature reconstruction was obtained on a parallel core, HP4. Core HP4 was obtained in September 2005 close to where previously studied sediment cores had been sampled 12,46 , using a Russian corer (10 cm diameter, 1 m length). HP4 includes a longer section of the basal unit compared with Core 5 and thus covers a slightly longer time period.
To correlate HP4 to Core 5 we used a Monte Carlo algorithm for core-to-core stratigraphic alignment of proxy records 29 (Supplementary Fig. 5) and then interpolated the reconstructed temperatures to Core 5's sampling resolution using a spline interpolation (the chironomid-based temperature record from HP4 had a relatively higher temporal resolution than the data sets of Core 5). The method involves nonlinear deformation of the entirety of one record onto a reference sequence and makes use of a Monte Carlo technique to deliver the optimal correlation. The approach was designed to correlate proxy records from nearby sedimentary cores and to mutually transfer information between cores retrieved from within the same small depositional basin. We applied our method on loss-onignition (LOI) data from both cores. This is an ideal proxy for the application of the matching procedure since LOI time series exhibit the same high-resolution and identifiable stratigraphic pattern of variability in all the cores that have so far been obtained and studied at Hässeldala (e.g ref. 12). This is also due to the relatively small size of the site (approximately 20 m 2 ).
Using the same procedure, Core 2, which had previously been studied for TOC, pollen and tephra analysis, and which had been dated by 29 14 C dates 46 , was correlated to Core 5 ( Supplementary Fig. 6). In this case, the cores were correlated using the TOC data after transforming LOI values of Core 5 into TOC as described in Muschitiello et al. 29 . The alignment between these two cores enables transferring a number of the earlier published 14 C dates 46 to Core 5 and thus allows us to improve the published chronology of Core 5 (ref. 47), providing the framework for the agedepth model presented in this study and described in the following section.

Bayesian age-depth modelling
By matching the TOC records of Core 2 and Core 5, we were able to transfer additional AMS 14  To dismiss spurious dates from the combined 14 C age data set, we first modelled the radiocarbon dating sequences of Core 2 and Core 5, separately. The 14 C dates for Core 2 and Core 5 were modelled multiple times prescribing different prior information parameters. At this stage we detected four potential age reversals on Core 2 and three on Core 5, which were rejected as probable outlying dates.
We then ran a series of preliminary age models using the composite series of 14 C dates in order to identify additional outlying dates. To do so, we approximated the error associated with each radiocarbon measurement to a normal-Gaussian probability distribution. This approach forces Bacon to become less conciliatory towards outlying ages and helped us to detect two additional outliers.
Finally, a series of 40 out of the 49 initial 14 C dates was employed to generate the final model (Supplementary Fig. 7 The age-depth model shows a solid structure and appears to be robust during the regional Late Allerød pollen zone and at the transition into the regional Younger Dryas pollen zone ( Supplementary Fig. 7). The reliability of the age model is further confirmed by the good agreement between the age of the Allerød-Younger Dryas pollen zone transition (12680 yr BP) at Hässeldala and other regional pollenstratigraphic records in southern Sweden 48 .
The age-model results obtained with Bacon were compared to output generated in OxCal4.2 (ref. 4). We run two age models using reasonably 'relaxed' k parameters of 0.5 and 1, respectively. The parameter defines the stiffness of the model upon the dating sequence when reproducing the sedimentation process. Therefore, the relatively low k values adopted constrain the radiocarbon sequence with a sizable degree of flexibility. Furthermore, we have not prescribed any lithological boundary, which in turn made the model even more flexible. In both cases the age results were equivalent to those yielded by Bacon but with relatively narrower age uncertainties.

Synchronization of NGRIP and IntCal13 time scales
The

Biomarker identification, quantification and hydrogen isotope analysis
The saturated hydrocarbon fraction was analysed by gas chromatography -mass spectrometry (GC-MS) for identification and quantification. GC-MS analysis was performed on a Shimadzu GCMS-QP2010 Ultra equipped with an AOC-20i auto sampler and a split-splitless injector operated in splitless mode. Volatile components were separated using a Zebron ZB-5HT Inferno GC column (30 m length, 0.25 mm inner diameter, 0.25 m film thickness). The GC oven was programmed to heat with 20 °C min -1 starting at 60°C, to 180°C, and subsequently the temperature was ramped to 320 °C at 4 °C min -1 , followed by an isothermal phase of 30 min. Helium was used as a carrier gas with a continuous flow rate of 1 ml min -1 . The ion source of the MS operating system was set to 200 °C and the ionization energy at 70 eV. C 18 to C 33 n-alkanes were identified based on mass spectra from the literature and retention times. Concentration of individual compounds was based on the comparison of peak areas relative to that of an internal standard (squalane) that was added to the samples before total lipid extraction.
Hydrogen isotope ratios were determined by gas chromatography-isotope ratio monitoring-mass spectrometry (GC-IRMS) with a Thermo Finnigan Delta XL mass spectrometer coupled to a Trace Ultra GC via a GC Isolink Conflo IV system. A standard mixture of n-alkanes with known isotopic composition (reference mixture A4, provided by Arndt Schimmelmann, Indiana University, USA) was run several times daily to calibrate the CO 2 reference gas against which the samples were measured. All analyses were performed in triplicate and results are reported as the mean, and expressed relative to the VSMOW scale.
Concentrations and hydrogen isotope ratios of n-C 21 , n-C 27 , n-C 29 , n-C 31 alkanes discussed in this study are presented in Supplementary Figure 1 and 2, respectively.

Chironomid-inferred quantitative temperature reconstruction
The core HP4 was sampled every 2 cm for chironomid analyses. Over 100 head capsules were counted for each sample, where possible. In most samples, 1-2 g of sediment was sufficient to obtain over 100 head capsules. Here, volumes processed

Removal of changing ice-volume and isotopic fractionation effects on sedimentary δD aq
To investigate the temporal evolution of 'non-amount' effects 60  The result was then removed from the δD aq record using the equation (1): (1) The effect of elevation changes in Blekinge was inferred from a glacial rebound model 64 Fig. 8).
We underscore that the predicted precipitation water isotope ratios generated with OIPC are based on a statistical model that uses as independent variables only latitudes and elevation 3 (Supplementary Fig. 9).
A potential pitfall of our approach is that all the 14 C dates used here are based on bulk sediment. In fact, the isolation of the lake basins was associated with sparse or no vegetation after ice sheet retreat and thus no plant macrofossils were available for radiocarbon dating 8,9 . Radiocarbon dating on bulk material may be biased owing to contamination of bulk sediment by different sources of carbon, which can carry significantly younger/older 14 C signatures 68 . Nonetheless, this drawback can be circumvented since a well-established stratigraphic relationship exists between the isolation of lake basins west of Mt. Billingen and specific traits featured in the regional pollen stratigraphies 7 . Therefore, the opening of the outlet area at Mt.
Billingen can be indirectly and independently dated in other regional sedimentary records where plant macro remains are available for 14 C measurements, thus providing a countercheck on the validity of the age results yielded by bulk dates.
Moreover, the available 14 C date from the Arkona Basin 10 is based on plant macro remains and it therefore provides a direct test of the reliability of the bulk dates. The isolation of several lake basins west of Mt. Billingen is thought to have occurred during the later stage of the Allerød regional pollen zone, corresponding with a phase of distinct increase in NAP values caused by increasing Artemisia pollen percentages 7 . This phase has been radiocarbon dated in other sedimentary sequences using plant macro remains and the inferred age is in excellent agreement with that obtained from bulk material elsewhere (Supplementary Table 2).

Transient climate model simulation
The location and the rate of the meltwater discharge have substantial uncertainties