Palaeo-productivity record from Norwegian Sea enables North Atlantic Oscillation (NAO) reconstruction for the last 8000 years

The North Atlantic Oscillation is the dominant atmospheric driver of North Atlantic climate variability with phases corresponding to droughts and cold spells in Europe. Here, we exploit a suggested anti-correlation of North Atlantic Oscillation-phase and north-eastern North Atlantic primary productivity by investigation of south-eastern Norwegian Sea sediment cores spanning the last 8000 years. Age model uncertainties between 2 and 13 years for the period 1992–1850 AD allows for the proxy to observational data calibration. Our data suggest that Ca/Fe core-scanning results reflect sedimentary CaCO3 variability in the region. Cross-correlating the Ca/Fe record with nearby phytoplankton counts and dissolved O2 data suggests that Ca/Fe can be used as a proxy for primary productivity variability in the region. Our data support an anti-correlation of primary productivity to the winter North Atlantic Oscillation index. Hence, we propose a sub-decadally resolved palaeo-North Atlantic Oscillation reconstruction based on an open-ocean record spanning the last 8000 years.


INTRODUCTION
The North Atlantic Oscillation (NAO) is the dominant driver of atmospheric variability in the North Atlantic region, determining the strength and position of the westerlies, which are the prevailing winds from the west across the North Atlantic. The westerlies strongly influence the surface ocean circulation in the region, affecting current speeds and mixing the surface layers 1,2 . The NAO index is defined as the sea-level pressure difference between the Icelandic Low and the Azores High, a positive (negative) phase generating warmer (colder) and wetter (drier) climate conditions in north-western Europe 3 . The NAO is most pronounced during the winter months (December-March) and varies on timescales ranging from days to centuries 4 . There are currently only a few published pre-instrumental NAO reconstructions based on an open-ocean sediment record 5,6 among the numerous reconstructions based on palaeo data [7][8][9][10][11][12][13] . This is likely due to challenges such as limited temporal resolution and lack of precise (sub-decadal) age control in marine records. However, if these challenges are overcome, marine records could offer insights from, e.g., key features of the North Atlantic circulation and enable the investigation of the important ocean-atmosphere link. Specifically, in regards to a non-stationary phenomenon like the NAO, the ocean-wide "catchment" size can be advantageous. A NAO reconstruction based on an open-ocean proxy record from the North Atlantic has, therefore, the potential to be an important advance on existing proxy-based reconstructions. Several studies have suggested and numerically tested the influence of the NAO on the dynamics of the phytoplankton productivity through air-sea fluxes and its impact on the circulation in the eastern North Atlantic 2, [14][15][16] .
Variability in biogenic calcium carbonate (CaCO 3 ) content in marine sediments has been used to infer Quaternary palaeoceanographic changes on glacial-interglacial timescales in the North Atlantic 17,18 . However, the relationship between sediment carbonate and sea-surface conditions is complex, and general statements on such relationships have proven to be difficult to demonstrate. The sources of biogenic CaCO 3 in deep-sea sediments are primarily coccolithophores and planktonic foraminifera, followed by pteropods and calcareous dinoflagellates 19,20 . In addition to variability in the productivity of these organisms, the biogenic CaCO 3 content in the sediments is diluted through the flux of terrigenous material and may also be affected by chemical dissolution 19,21 . The productivity of CaCO 3 secreting algae is essentially constrained by light, nutrients and temperature 22 . As such, productivity has a distinct seasonal cycle, peaking during spring and autumn phytoplankton blooms 23 , but also relays decadal to millennial, climatically driven, changes in the surface water conditions 2,24 .
This study investigates which factors influence CaCO 3 concentrations in the slope sediments on the southern Norwegian continental margin. This is done in order to address the following objectives: (i) can ITRAX-derived Ca/Fe variability be used as a proxy for changes in sea-surface primary productivity in the region and (ii) how this proxy may reflect known modes of oceanic and atmospheric variability in the region. To address these questions, two well-dated, high-resolution sediment cores covering the last 8000 years were analysed. These cores were collected at 850 and 960 m water depth on the south-eastern Norwegian margin (Fig.  1a). Both core stratigraphies overlap with the late 20th-century observational time series of sea-surface conditions and productivity, enabling the calibration of proxy to observational data.
Slope Current (NwASC), and the western branch as the Norwegian Atlantic Front Current (NwAFC). The NwASC is a barotropic shelf edge current, streaming over the Vøring Plateau 25 , while the NwAFC is topographically constrained and flows around the edge of the outer Vøring Plateau (Fig. 1a). The total transport of Atlantic surface water within the NwASC and NwAFC has been estimated to 6.8 Sv 26 , with salinities >35.0 PSU and temperatures of 5-10°C 27 . The shallow, northwards flowing Norwegian Coastal Current (NCC) covers and partly overlies the surface waters east of the NwASC (Fig. 1b). The surface sediments in water depths up to 700 m consists of sorted sand and gravel, reflecting high-energy conditions within the Atlantic inflow water. Between 700 and 1200 m water depth, corresponding to the depth, the presented cores were retrieved from, the surface sediments consist of pelitic mud 28 , reflecting the low-energy regime within the homohaline Norwegian Sea Atlantic Intermediate Waters (NSAIW, Fig. 1b). These pelitic muds, which are considered to be mainly derived from erosion of upper slope and shelf sediments 28,29 contaiñ 15-30% CaCO 3 in surface samples 28 , changing to foraminiferal ooze below 1200 m water depth 28 .
Cores GS13-182-01CC (63.64°N, 05.51°E), hereafter GS13, and P1-003 (63.76°N, 05.25°E) penetrate into the slide debris of the Storegga Slide ( Fig. 1a and Supplementary Fig. 1), which is dated to 7650 ± 250 14 C years BP 30 . The slide debris in the lowermost part of GS13 is overlain by~12 cm of sorted, very fine sands (77% >63 µm), with a sharp transition into marine mud at 1777.7 cm (Fig. 2f). The structureless marine mud does not show signs of bioturbation and appears uniform with, on average, 98.5% silt and clay (<63 µm). The remaining 1.5% of the material (>63 µm) consists, on average, to 85% of very fine sand (<125 µm). The age models of both cores are primarily based on Icelandic tephra, radiocarbon and 210 Pb dates and consist of an ensemble of 19,000 equally likely Bayesian age models (see "Methods"). The age model suggests that the accumulation rate in core GS13 is relatively constant from 0.16 cm per year between ca. 8000 to 2000 a BP, rising to 0.8 cm per year in the last 1000 years, and on average to 1 cm per year in the last 100 years (see "Methods" and details on chronological approach in Supplementary Note 1). The increase in accumulation rate appears to be unrelated to sediment-water content as this is relatively constant throughout the core (~50% before 1000 a BP, rising to 60% on average in for the youngest part). Chronological uncertainties are 2-13 years back to 1850 AD, and 17-190 for the section from 1850 AD to 8000 a BP (see detailed description in "Methods" and on the chronological approach in Supplementary Note 1).
To investigate, if bulk sediment CaCO 3 variability is reflected in the bulk elemental composition from an ITRAX X-ray fluorescence (XRF) core scanner, we focused on raw, bulk element XRF Ca counts and CaCO 3 measurements at selected depths (see Supplementary Note 2 and Supplementary Fig. 2). The XRF Ca counts were normalised with the counts of a lithogenic element (Fe) to account for closed-sum constrains of compositional data (see "Methods"). The raw Ca/Fe time series exhibits distinct changes to relatively high Ca concentration, centered on 5000, 3000, 2000 and 0 a BP (Fig. 2). Comparing the Ca/Fe ratio with the Ca counts normalised with the total kilo counts per second (kcps) does not reveal any major differences, suggesting no or little detectable terrigenous influence on the Ca variability (Fig. 2). Furthermore, small variations in the sub-annually resolved very fine sand fraction (Fig. 2f) do not correspond to the variability in the Ca/Fe ratio, indicating, that grain-size fluctuations related to current speeds are not principally driving the changes in relative Ca concentrations. Sortable silt data in the region were earlier interpreted to show lowerfrequency variability related to changes in NwASC current speed 31 . The grain-size data from GS13 are relatively uniform throughout the core, suggesting little direct current speed influence on the terrigenous flux at the core site but this cannot be entirely excluded.

Sources of biogenic CaCO 3
The main sources of biogenic CaCO 3 in (hemi-)pelagic sediments are commonly foraminifera and coccolithophores with coccoliths often contributing 40-60%, although usually even more in higher-latitude  in the investigated cores, which is in line with previous studies suggesting coccolithophores as the primary source of CaCO 3 in the fine fraction of sediments in the Nordic seas 34 . The relative abundance of species Emiliania huxleyi in surface sediments in the Norwegian Sea is ca. 60%, while Coccolithus pelagicus accounts for~20-40% 35 . The contribution of C. pelagicus to the CaCO 3 flux to the sediments is approximately one order of magnitude higher than the contribution of E. huxleyi 36 due to its larger size. On the North Iceland Shelf, C. pelagicus accounts for >95% of the total coccolith carbonate mass throughout the last 8000 years 37 . We suggest that the Ca/Fe time series of the investigated cores primarily reflects CaCO 3 , which in turn primarily reflects the variability in the deposition of coccoliths, with the species C. pelagicus as the main contributor. Factors such as nutrients and irradiance have been suggested as the primary control on the distribution of C. pelagicus 35,37,38 . Within the yearly cycle, the highest concentrations of C. pelagicus are reported in July in the Norwegian Sea, while the state of preservation increasingly deteriorates throughout the season 36 . Given the evidence and the reasoning for the relation between Ca/Fe variability and coccolith abundance described above, we hypothesise that the variability in Ca/Fe can be used as a proxy for primary productivity of CaCO 3 secreting algae in the south-eastern Norwegian Sea. This relationship might, however, be distorted by the dissolution of the tests before or after they are enclosed in the sediments. Several tests run with scanning electron microscopy revealed no signs of dissolution like overgrowth or increased fragility on the CaCO 3 microfossil assemblages investigated in GS13 and P1-003 (Supplementary Figs. 3 and 4). This is in line with previous findings within the Nordic seas 39 . However, postdepositional dissolution will likely have a negligible effect in this case due to the extremely high sedimentation rates, which is underlined by the absence of signs of dissolution in the samples.

Ca/Fe measurements and primary productivity variability
To test if Ca/Fe can be used as a proxy for total primary productivity, the Ca/Fe data are statistically validated against available instrumental observations of primary productivity. Surface phytoplankton concentrations from the region, as available from the Continuous Plankton Recorder (CPR) survey since 1958 (see "Methods"), are the most direct observation of primary productivity in the region 40 . Available time series of primary productivity observations suggest that CaCO 3 secreting plankton (coccolithophores and foraminifera) co-vary with non-CaCO 3 secreting plankton (diatoms, dinoflagellates, radiolarians, etc.; Supplementary Fig. 5). Therefore, here, the record of CaCO 3 primary producers 41 is assumed to reflect the total primary production variability in the region (see "Methods"). The time period used for calibration is anomalous from the past due to possible human impact on the marine environment 42 , but the absolute relationship of the boundary conditions to the environment should, however, stay valid.
Cross-correlations suggest a strong link between the Ca/Fe record and the summer total phytoplankton concentrations of standard area B4 (Figs 1a and 3a-c). About 84% of all 19,000 age models (see "Methods") suggest a significant (P < 0.05) correlation, while almost 70% of all age models suggest a lag of 1 or 0 years (Fig. 3b). The median cross-correlation coefficient (R) of all significant fits is R = 0.54, with the best fits exceeding R = 0.80 ( Fig. 3c and see "Methods"). As expected from the relatively large age uncertainty in the age model input data in the youngest 10 years of the core (see "Methods"), the fit and the spread of age models deteriorate slightly between 1983 and 1993 (Fig. 3a). Notably, the fit of Ca/Fe and the autumn and spring data is significantly lower (median R = 0.37), even though the majority of age models still suggest a similarly sized lag, while there is no reasonable correlation with the winter concentrations (Supplementary Fig. 6).
Dissolved O 2 data are closely related to the strength of the primary production 43 . Hence, the cross-correlation script was additionally used to evaluate the Ca/Fe variability in relation to the seasonally resolved, dissolved O 2 data from Ocean Weather Ship Mike (OWSM; integrated over 1-50 m water depth, Fig. 1). Here, only the 10,000 age models with the best fit to standard area B4 were used (Fig. 3a-c). The best cross-correlations were computed with spring/early summer (MAMJJA) dissolved O 2 (median R = 0.40). The strong correlation with dissolved O 2 in spring/early summer coincides with the main production period, the phytoplankton blooms 23 and the highest concentrations of C. pelagicus in the region 36 .
These results give confidence in using the variability of Ca/Fe as a proxy for relative changes in total primary productivity in the south-eastern Norwegian Sea. This conclusion is supported by the trends of Ca/Fe closely resembling the temporal variability of phytoplankton mass measured in chlorophyll concentration in the eastern North Atlantic through the last century 24 (Fig. 4). A very similar trend over the last century is also seen in a recent Greenland ice core-based subarctic Atlantic primary productivity reconstruction 44 . The weak cross-correlation to other available areas with CPR phytoplankton concentration data (including the area linked to the NCC; Fig. 1a and Supplementary Fig. 6) supports the assumption that the Ca/Fe record is an open-ocean signal, reflecting the productivity from the inflow trajectory of the Norwegian Atlantic Current into the Norwegian Sea (Supplementary Fig. 6).

DISCUSSION
A number of studies have suggested that annual primary productivity within the south-eastern Norwegian Sea anticorrelates to the polarity of the NAO 2,15,16,45 . This effect has been suggested to result from an increase in wind-driven mixing and the enhanced inflow of Atlantic surface waters into the Norwegian Sea through strengthened westerlies during the winter season 2,14,46,47 . These factors disturb the water mass stratification and result in a greater mixed layer depth (MLD) 48 . The shoaling of the MLD in spring is, next to sunlight, an important driver of the onset of the subpolar spring bloom. A greater MLD ultimately delays the start of the spring bloom, effectively shortening the growing season and reducing the magnitude of the phytoplankton bloom and with that the primary production 16 . This effect might be increased through stronger surface currents driving the spatial distribution of the bloom, potentially adding to the net reduction in primary productivity at the core site 49 . Thus, the positive (negative) phase of the NAO leads to decreased (increased) primary production within the south-eastern Norwegian Sea.
To examine a possible correlation between the suggested primary productivity proxy and the NAO, the Ca/Fe record of GS13 was cross-correlated to the instrumentally measured, seasonal NAO index of the years 1824-1992 AD 50 (Fig. 5 and Table 1). Crosscorrelating the data within the years with the best data quality (1950( -1992, the results suggest a significant (P < 0.05), robust anti-correlation with the winter NAO (median R = −0.48, best fit R = −0.72) in line with a previously published correlation between proxy-based data and the instrumental NAO record (e.g., Fig. 3 in Ortega et al. 10 ). The distribution of lags is relatively uniform, with L.W.M. Becker et al.
most age models suggesting a lag of 0 or 1 year (Fig. 5b). The covariability of the Ca/Fe record with the NAO summer, spring and autumn seasons shows a decreasing co-variability throughout the year (Table 1 and Supplementary Fig. 7). Apart from winter and intra-seasonal variability, other seasons therefore most likely influenced the signal as well. However, when extending the cross-correlation to the available instrumental data, i.e., until the year 1824 (Fig. 5d-f), the fit weakens significantly, despite of a similar lag size (median R = −0.20, best fit R = −0.44, see Table 1 for an overview). This reduction might be caused by several factors, e.g., chronological uncertainties, data quality in the Ca/Fe record or data quality degradation in the earlier instrumental NAO data. The change in the degree of correlation might also be a factor of the proposed non-stationarity of the NAO signal. However, the strong correlation between the pressure records from different stations of the southern NAO node 51 suggests that discrepancy between the instrumental and proxy records are stemming from the proxy data and not from the instrumental data. Taken together, these results and the suggested mechanistic link between primary productivity and NAO variability indicates that Ca/Fe in core GS13 does reflect NAO variability in the last ca. 150 years. However, the variance observed in the Ca/Fe record can not all be explained by the NAO variability. Apart from influencing factors discussed above such as geochronological uncertainty and sedimentological processes, another climate/weather regimes may explain some of the variability. Tests of the relationship to such regimes pointed to influences of e.g., the Atlantic Multidecadal Oscillation 52 (AMO; see Supplementary Note 3 and Supplementary Fig. 8) and the subpolar gyre index 53 (SPG; Supplementary Fig. 9) in the signal. The co-variability between the SPG and the Ca/Fe data appears to be strong and direct, but the cross-correlation is limited to a period of 42 years. In conclusion, based on these considerations, we suggest that our   . Only significant correlations (P < 0.05) are used. For orientation, the annually subsampled Ca/Fe on the median age model is also shown (blue, dotted). b Distribution of median cross-correlation coefficients (R-value, black dots) across 31 lead/ lags (negative values denote a lag of Ca/Fe relative to compared proxy) underlain by 5-95% quantiles (grey area). The mode of lag, % significant age models and the median R-value at this lag are given at the top. The histogram denotes at which lead/lag each age model has the highest R (n = 19,000). Only significant correlations (P < 0.05) are used. c Density plot of the distribution of R at a lag of 1 year (mode of the histogram in panel b). Significant correlations (P < 0.05, red) are underlain by all correlations (grey). d-f Same as panel a-c, but with mean spring-early summer (MAMJJA) dissolved O 2 data from OWSM integrated over 1-50 m water depth. Only the best-fitting age models of panel b (lag of −1 and 0, n = 10,000) are used.
Ca/Fe record may be used as a proxy for winter NAO variability back in time, beyond instrumental data.
The raw 8000 years, Ca/Fe record of GS13 (Fig. 6b) exhibits a non-stationary signal with a broad frequency range. The underlying lower-frequency component (periods >500 a) is removed from the Ca/Fe record in the following, as these frequencies may mainly be related to trends in the development of the sedimentary supply system, like the sea level and the geometry of water masses. Additionally, the lower-frequency component shares some similarities with the low-frequency components in the reconstructed total radiocarbon production rate (Fig. 6a), which indicates, among other things, solar activity changes 54 . The higher-frequency components (periods <500 a) of the Ca/Fe record show annual to (multi-)centennial-scale variability. This part of the record, hereafter NAO CaFe , is therefore used to reconstruct the NAO variability through the last 8000 years. To improve the readability of the NAO CaFe index, the data are bandpass filtered (10-500 a), when focusing on the complete record (Fig. 6e). To investigate only the high-frequency components of the NAO CaFe index, the data were additionally high-pass filtered (periods <100 a, Fig. 6f) and compared to existing reconstructions (Fig. 7).
The high-frequency components (periods <100 a) of NAO CaFe correspond relatively well with available NAO reconstructions covering the last 1000 years (Fig. 7). Cross-correlations run with NAO growth 13 ( Supplementary Fig. 10) and NAO med 7 (Supplementary Fig. 11) show a significant fit within the last 1000 years (see Table 1 for an overview). Specifically, NAO growth has a good fit within the last 300 years (median R = 0.55, best fit R = 0.70). Ca/Fe correlates in this case directly with a stalagmite growth index, NAO growth , as both are proposed to be anti-correlating to the NAO polarity. The fit decreases throughout the last 1000 years (median R = 0.28, best fit = 0.55) with a still reasonably well peak-to-peak correlation, with the exception of the period between 1200 and 1300 a AD (Supplementary Fig. 10). Before 1000 a AD, the fit decreases further, which could in part be explained by the increasing chronological uncertainties. Notably, NAO CaFe does not show a distinct prevailing positive NAO phase during the Medieval Warm Period (MWP), as seen in the NAO growth , in line with other studies 55 . The fit to NAO med is weaker than the fit to NAO growth and best within the last 300 years (median R = −0.27, best fit R = −0.45). The fit reduces further throughout the last 1000 years ( Supplementary Fig. 11, see Table 1 for an overview). The fit to NAO PC1 8 is of a similar size within the last 300 years, but the crosscorrelation analysis reveals a lag of 39 years. Additionally, NAO CaFe corresponds well with a record showing the likelihood (degree of belief) for the NAO being in a specific phase, based on a North Atlantic functional palaeoclimate network 56 . A sign-test between NAO CaFe (Fig. 7a) and the degree of belief record (Fig. 7b) of Franke et al. 56 showed an agreement of~50% in the last 1000 a (thresholds set to > ±0.001 in Ca/Fe and >0.5 in the degree of belief data). For this comparison, a 50-year running mean is calculated for both records. However, introducing a lag of 25 years into the degree of belief record relative to NAO CaFe increases the agreement between the records to >77%. The lag offset of the degree of belief data could solely be due to uncertainties within the model used to calculate the degree of belief data, with respect to the timing of changes 56 .
Beyond the last 1000 years, reconstructions become scarcer, with only three reconstructions extending beyond 2000 a BP 9,13,57 ( Supplementary Fig. 12). Difficulties in reconstructing a nonstationary signal like the NAO 1 become apparent when comparing the available long-term reconstructions with each other, resembling each other only during few periods ( Supplementary Fig. 12, and reviewed by Pinto et al. 4 ). This might again be an effect of a small catchment, which results in an incomplete NAO signal due to a shift in the centre of action of the NAO, as indicated by Lehner et al. 55 . However, the sampling site of GS13 is, as detailed above, influenced by the variations in surface currents of the North Atlantic and the dominant atmospheric modes of the region, giving confidence in recording a regional signal. Thus, it could be argued that the NAO CaFe record might be capable of recording a signal, which is less biased by the non-stationarity of the NAO.
In summary, the NAO CaFe record shows a relatively robust, statistically significant fit to the instrumentally measured NAO index. Beyond the instrumental time period, NAO CaFe has a good fit with NAO growth 13 and a weaker but significant fit to NAO med 7 and NAO PC1 8 (see Fig. 7c-e and Table 1 for an overview of correlation coefficients). Additionally, NAO CaFe mostly agrees with the NAO phase indicated by the North Atlantic palaeoclimate network model from Franke et al. 56 (Figs 6d and 7b). Extending beyond 1000 a BP, NAO CaFe agrees with some phases of the available reconstructions ( Supplementary Fig. 12). The available reconstructions are however shown to be to some degree inconsistent among each other 4,7,8 , which makes it challenging to assess direct correlations between the reconstructions. Taken together, the primary productivity variability within the southeastern Norwegian Sea reflected in the Ca/Fe record shows potential for a NAO index reconstruction throughout the last 8000 years, as it performs at least as well in correlating to the available instrumental NAO data compared to previous reconstructions (Table 1). Comparing the last 8000 years of high-frequency variability of NAO CaFe with known climate anomalies like the Little Ice Age (LIA), the MWP and Holocene Thermal Maximum (HTM) (Fig. 6) do not evidence any general changes in the amplitude of the NAO associated with these events. Studying the lowerfrequency variability in the Ca/Fe record, there appear to be some trends and phases which may be associated with larger-scale climatic changes. However, it is challenging to determine what part larger-scale changes in ocean current systems are playing and to what extend primary productivity changes reflect this lowfrequency variability. Future studies should focus on a detailed discussion of the Holocene NAO variability and its potential links to known climatic events and human history 58,59 .  . Only significant correlations (P < 0.05) are used. b Distribution of median cross-correlation coefficients (R-value, black dots) across 31 lead/lags (negative values denote a lag of Ca/Fe relative to compared proxy) underlain by 5-95% quantiles (grey area). The mode of lag, % significant age models and the median R-value at this lag are given at the top. The histogram denotes at which lead/lag each age model has the highest R (n = 19,000). Only significant correlations (P < 0.05) are used. c Density plot of the distribution of R at a lag of 0 year (mode of the histogram in panel b). Significant correlations (P < 0.05, red area) are underlain by all correlations (grey area). d-f As in panel a-c, but for 1824-1992 (bandpass filtered, periods 3-30 a).

Core material
The 1970-cm long CALYPSO piston core GS13 has raised 18.5-km southsoutheast of core P1-003 at 960 m water depth ( Supplementary Fig. 1), recovering 1777.7 cm of marine deposits. Core P1-003, collected at 850 m water depth, is a spliced archive, consisting of a 42-cm long MultiCore (MC) and a 690-cm long SelanticCore (SC) 60 . The top of the SC section was spliced with the MC section at 32.05 cm, in total recovering 670.3 cm of marine deposits. All core depths presented are on merged and void corrected scales.

Geochemical analysis
Cores GS13 and P1-003 were analysed for their bulk elemental composition on an ITRAX X-ray fluorescence (XRF) core scanner, using  Fig. 6 The complete NAO CaFe record. a Total radiocarbon production rate (grey line, inverted) 54 , overlain by low-passed version (red line, inverted, cut-off 500 a). b The raw Ca/Fe data (black line) overlaid with the low-passed version (red line, cut-off 500 a). c Bandpass filtered Ca/ Fe data (10-500 a periods) demonstrating the changing amplitude of the signal through time, lowest amplitudes around the Holocene Thermal Maximum (HTM) and the Medieval Warm Period (MWP). d Degree of belief, that the NAO is in a specific phase (50-year mean; indicated at the start of the window) based on a North Atlantic functional palaeoclimate network 56 . Grey areas indicate where less than 66% of all ensemble members agreed on the same NAO phase (modified from Franke et al. 56 ). e The NAO CaFe index consisting of annually resampled, bandpass filtered (cut-off 10-500 a) XRF Ca/Fe data, shown in an inverted, two-coloured plot. f The NAO CaFe index (annually resampled, bandpass filtered, cut-off 10-100 a, inverted, two-coloured plot). Known, sustained warm periods are marked in light red and respective cold periods in light blue. HTM Holocene Thermal Maximum, RWP Roman Warm Period, MWP Medieval Warm Period, LIA W Norway Little Ice Age western Norway.
the Molybdenum X-ray tube at room temperature. GS13 and P1-003SC were scanned with a resolution of 500 µm and a counting time of 10 s (100 s in P1-003 MC). The present work focuses on the XRF data recorded on GS13, whereas the XRF data from P1-003 61 are primarily used for age model refinement (Figs. 2 and 8). The XRF core-scanning technique produces semi-quantitative elemental counts rather than elemental concentrations. Grain-size variations can influence the elemental counts 62 , but the homogenous grain-size distribution (about 2% >63 µm, Fig. 2) indicates little grain-size bias of the XRF scanning results. To account for closed-sum constraints 62 , the Ca counts were normalised with the respective Fe counts (Fig. 2). Normalisation with other lithogenic elements like Ti or K yielded very similar results ( Supplementary Fig. 13). The log transformation of XRF ratios has been suggested to overcome closed-sum constraints 63 . However, once applied to GS13, the time series were practically identical, indicating a normal distribution of the XRF data. For simplification, the log transformation has thus not been applied. Additionally, interstitial water content and water films forming below the plastic foil can also impact XRF scan results 64 . For this study, the cores were immediately scanned after splitting the cores, hence reducing potential water-film issues. The sediment-water content increases throughout the core linearly from 50 to 60%, suggesting only little influence on the XRF scan results. The comparison of the down-core Ca intensities with the Cl intensities confirmed this assessment ( Supplementary Fig. 13).

Chronology
To establish chronological control of GS13, the planktonic foraminiferal species Neogloboquadrina incompta and Globigerina bulloides were picked in the 125-150 µm fraction for accelerator mass spectrometry (AMS) radiocarbon dating at INSTAAR, the University of Colorado at Boulder. The  Fig. 6f. b Degree of belief, that the NAO is in a specific phase (50-year mean, indicated at the start of the window) based on a North Atlantic functional palaeoclimate network 56 . Grey areas indicate where less than 66% of all ensemble members agreed on the same NAO phase (modified from ref. 56 ). c As in (a), but with inverted stalagmite growth rate, NAO growth 13 . d As in (a), but with NAO med index 7 . e As in panel a, but with NAO PC1 index 8 . For a full comparison throughout the last 8000 years, the reader is referred to Supplementary Fig. 12. radiocarbon dates were converted to calendar years, using the local delta R 60,65 and calibration curve Marine13 66 (Supplementary Table 1). Furthermore, based on the counting of basaltic-intermediate and rhyolitic type of volcanic shards in the 63-125 µm fraction, three Icelandic tephra were identified (Hekla 1947, Katla 1918 and Askja 1875, Fig. 8c-e). Finally, continuous 2-cm wide samples were taken in the top 50 cm of core GS13 for 210 Pb/ 137 Cs dating, analysed at Eawag, ETH Zürich. With this chronological information, an initial age model was established for GS13 using the Bayesian age modelling package "rbacon" (v. 2.3.6) 67 within R. To further decrease the chronological uncertainty beyond 1850 AD, the age model was refined with a detailed correlation of the XRF Ca/Fe signal (Fig. 2) to the welldated, nearby core P1-003 60,65 (Supplementary Fig. 14 and Supplementary  Table 2). All chronological information was combined in a final age model run, calculating the final age model (Fig. 8a). The two-sigma uncertainty of the age model is~2-13 (average 5) years within the last 150 years, which increases to~2-24 (average 7) years within the last 200 years. Before the year 1800 AD and until 0 AD/BC, the two-sigma uncertainty of the age model is~17-110 (average 52) years and~58-190 (average 100) years for the remainder of the core. For more details on the chronology see Supplementary Note 1.

Analysis of subsamples
Subsamples in 0.5 cm slices (~13 cm 3 ) were taken at 5 cm interval throughout core GS13, except for the top 2 m of the core, which were sampled continuously. The subsamples were each treated with 35% hydrogen peroxide (H 2 O 2 ), shaken for 48 h at 150 rpm and subsequently wet sieved with a sieve mesh size of 63, 125, 150 and 1000 µm (Fig. 2). Subsamples (43 samples in 0.5-1-cm slices) from core P1-003 were used to measure absolute element concentrations of CaCO 3 , calcium oxide (CaO) and iron oxide (Fe 2 O 3 ) with the sequential spectrometer Philips 1404 X-ray, using wavelength dispersive X-ray Fluorescence 61 . Furthermore, the total inorganic carbon content (TIC) of selected subsamples from P1-003 61 was measured with the multi EA 4000 for Macro Elemental Analysis (Analytik Jena). The percentage of CaCO 3 was calculated from the TIC content assuming that all evolved CO 2 was derived from the dissolution of CaCO 3 , using the molecular weight ratio of CaCO 3 to carbon (8.33), where~85% of the samples are averages of replicate measurements ( Supplementary Fig.  2). Internal standards were included in every run in a ratio of 3:10, indicating a measurement uncertainty of ca. 3%.

Statistical analysis
The time series were exclusively downsampled during re-sampling, using linear interpolation, with ca. 20 points per year in the youngest sections and ca. 4 points per year in the oldest sections of GS13. Similarly, R was used to compute the annual mean values of the ensemble of age models produced during the age modelling process. This was done to enable the computation of cross-correlations in respect to chronological uncertainties. The cross-correlation script takes advantage of the Bayesian age model output, which produced~19,000 equally probable age models, covering the total spread of age uncertainty for each depth in the core (for details see description of the chronological approach in Supplementary Note 1). The script computes about 19,000 cross-correlations, allowing for leads/ lags of up to 15 years, resulting in about 600,000 possible combinations. The leads/lags were increased to 50 if longer (>100 years) time series were cross-correlated. The resulting matrix was then further analysed, statistically verified and corrected for autocorrelation with the "cor.test" function. This function uses Pearson's product-moment correlation coefficient and tests the significance (P < 0.05). The result is corrected for the reduced number of degrees of freedom of autocorrelated time series using the autocorrelation parameter (AR1) of both time series, following Bartlett 68 .  All annual re-sampling and filtering of the acquired data were performed with the "stats" package (v. 4.0.1) in R 69 . All plotting was performed with packages "ggplot2" and "cowplot" in R 69 .
Observational productivity data The median chronology calculated 1993 AD as the core top age, hence there is an overlap of~40 years between available instrumentally measured time series and proxy data. The CPR survey has been recording and counting various plankton species in surface waters at up to daily intervals, across the North Atlantic since 1958, integrating the results in standard areas 40 . The most proximal standard area with sufficient valid data is B4, upstream of GS13 (Fig. 1a). For this study, the sum of total dinoflagellate and total diatom counts were used as an approximation of primary productivity changes 41 . These "total phytoplankton" data are therefore primarily based on SiO 2 -producing algae. Valid coccolithophore counts are only available between 1999 and 2016. During the time of overlap, these parameters appear to covariate (r = 0.81, summer in standard area B4, Supplementary Fig. 5), suggesting that the variability in the "total phytoplankton" concentrations in this area can also be used as a proxy for the productivity variability of CaCO 3 -secreting algae. At OWSM (Fig. 1a) monthly measurements of temperature and dissolved oxygen were carried out from 1948 to 2008 at depths between 0 and 2000 m 46,48 . The data from OWSM for this study were quality checked and integrated over 1-50 m water depth, avoiding potentially erroneous surface measurements 46 , while covering the depth of highest coccolithophore production 70 . Additionally, as the surface waters on the core site might occasionally be mixed with water from the NCC (Fig. 1b) 71 , integrating the data down to 50 m water depth would ensure an Atlantic, open-ocean signal. Seasonal means of each year in CPR and OWSM data were calculated from monthly mean data, from winter to autumn, averaging over months DJFM, MAMJ, JJAS and SOND, respectively. The seasonal mean was discarded if data of >1 month were not available.

DATA AVAILABILITY
The data analysed during this study are available from the PANGAEA ® data repository.