Linear regression indicates a weakly significant correlation between variability in our δ18O-shell record and the North Icelandic shelf regional marine radiocarbon reservoir ages (ΔR)21 (r=−0.39±0.29, P=0.056; Fig. 2 and Supplementary Figs 4 and 5). The ΔR data represent the regional difference from the global mean ocean radiocarbon reservoir age, which on the North Icelandic shelf also reflects the degree of entrainment of relatively old (with respect to carbon) polar AIW into the younger SPMW transported by the NIIC21,25. The coherence between the North Icelandic shelf ΔR21 and δ18O-shell series indicates that during periods characterized by relatively older (younger) ΔR values, due to a relative increase (decrease) in the proportion of AIW within the NIIC, the δ18O composition of the shell aragonite increases (decreases), indicative of an increase (decrease) in the density of the ambient seawater bathing the site.

It has previously been proposed that changes in the strength of the Atlantic Meridional Overturning Circulation are linked to variability in the degree of AIW entrainment into the NIIC and therefore the SPMW/AIW composition of the waters influencing the North Icelandic shelf21,31,32. Although instrumental oceanographic observations are currently of insufficient length to evaluate this suggestion, evidence derived from an ensemble of 17 Coupled Model Intercomparison Project phase 5 (CMIP5) unforced pre-industrial control simulations suggests that changes in SST and SSS on the North Icelandic shelf correlate significantly with inter-annual and multi-decadal scale variability in North Atlantic barotropic and overturning stream functions (Fig. 4 and Supplementary Fig. 1). The stability of these correlations, irrespective of the smoothing function applied, indicates that the SST and SSS (along with the associated seawater density) variability on the North Icelandic shelf are sensitive to changes in Atlantic circulation dynamics (Supplementary Fig. 1).

Figure 4: Model analyses of the coherence between North Icelandic SST and SSS variability and North Atlantic circulation. Correlation between annual North Icelandic SSTs against (a) barotropic stream function and (c) overturning stream function, and SSS against (b) barotropic stream function and (d) overturning stream function, in the North Atlantic derived from an ensemble CMIP5 preindustrial control simulation. The solid black lines show the time-mean stream functions. Full size image

To better constrain the uncertainties associated with the role that the ocean system plays within the wider climate system, it is necessary to assess quantitatively the degree to which ocean variability is related to external forcing mechanisms. To do this we compare the new δ18O-shell record to millennial climate forcing time series (volcanics, TSI and modelled combined forcings), NHSAT, derived from a composite dendrochronological record11, and Greenland air temperatures, derived from an ice core stack13.

Estimates of past TSI variability can be characterized into two distinct intervals, the observed period (∼AD 1750 to present) and the proxy-generated reconstruction period (years before ∼AD 1750)2,33. The proxy reconstruction of TSI, derived from radiogenic isotopes (10Be and 14C) is typically deficient in the high-frequency component of the solar variability spectrum (Supplementary Fig. 6). The observed TSI period however captures both the longer-term trends, such as the gradual increase in TSI since the end of the Maunder Minimum in the early AD 1700s and high-frequency variability such as the 11-year sunspot cycle2,33. We identify a significant correlation between the δ18O-shell series and TSI over the last 1,000 years (annual resolution r=0.27, P<0.05; decadal resolution r=−0.34, P<0.05). The correlation is significantly stronger over the observed period of AD 1750–1997 (r=−0.68 P<0.05), suggesting that TSI variability may explain up to 47% of the variance in seawater density over this period. The coherence spectra of the correlations show that the increase in correlation over the observed period is not likely to be due to low-frequency variability, as similar levels of coherence were identified across these portions of the spectrum for both the observed period and the entire millennium (Supplementary Fig. 6). However, over the observed period, significant correlations (r=−0.39, P<0.01) were identified within the high-frequency domains (using 50-year loess first-order high-pass-filtered data; Fig. 5d). This high-frequency component of the solar spectrum is unrepresented in the proxy reconstruction derived from radiogenic isotopes (Supplementary Fig. 7). Although these analyses and several previous studies (e.g. refs 16, 17) indicate that solar variability probably plays a role as a natural driver of marine variability, the presence of several large volcanic eruptions over the observed period and the lack of any reliable proxy for high-frequency solar variability over earlier parts of the last millennium complicates the interpretation of the role of solar variability in driving oceanographic variability.

Figure 5: The δ18O-shell series and climate forcings. Comparison between (a) the inverted δ18O-shell anomalies (annually resolved and 50-year loess low-pass-filtered data shown by the pink and red lines, respectively) with (b) individual climate forcing indices (TSI and volcanics shown by the orange and green lines, respectively), and (c) combined natural forcings (solar and volcanics; green line) and combined total forcing (solar, volcanics, greenhouse gases, aerosols; black line). (d) Plot of the 50-year loess first-order high-pass-filtered TSI (black line) and inverted δ18O-shell anomalies (red line) over the period AD 1750–1997. (e) Timing of volcanic eruptions over the AD 1750–1997 period. The grey box highlights the 2 consecutive 30 periods following the Laki (1783) and Tambora (1815) volcanic eruptions, which probably dominate the climate signal over this period. Climate forcing data from Crowley2. Full size image

We examined the influence of volcanic aerosols injected into the stratosphere, a consequence of large explosive volcanic eruptions, on the δ18O-shell series by conducting a superposed epoch analysis (SEA34; Fig. 6). We find a significant positive shift in the δ18O-shell anomalies in the first (P<0.05) and second (P<0.1) year following volcanic eruptions. This result is consistent using two approaches to the SEA (see Supplementary Fig. 7). Although evidence from Swingedouw et al.35 suggest that volcanic eruptions can have an influence on the North Atlantic that spans several decades, as part of a bi-decadal oscillation, our analyses only indicates a significant response of the marine system in the two years following the largest eruptions of the last 1,000 years. The bi-decadal oscillation present in the SEA analyses of the δ18O-shell series, although similar in amplitude and timing to the previously suggested North Atlantic circulation variability associated with volcanic forcing, is not statistically significant35. These results indicate that seawater densities north of Iceland show a rapid response to volcanic forcing with an immediate decrease in seawater density over the two years following the eruption. The significant response of the δ18O-shell series to the volcanic forcings probably explains the drop in correlation observed between the δ18O-shell series and TSI over the period from 1783 to 1843. This period is characterized by three large eruptions, Laki in 1783, an unknown eruption in 1809 and Tambora in 1815, which probably drive the climate signal over this interval.

These analyses indicate that both the TSI and volcanic aerosol forcings are playing a major role as natural drivers of high-frequency (inter-annual to multi-decadal) marine variability on the North Icelandic shelf. However, the comparison between the δ18O-shell record and the combined climate forcing index (using TSI, volcanics, stratospheric aerosols and greenhouse gases; Fig. 5c and ref. 2) suggests that, even with the gradual 0.01% increase in TSI since the Maunder Minimum3, natural climate forcings alone cannot account for the significant change in δ18O-shell values over the industrial period.

Given the precisely dated nature of the δ18O-shell series and its sensitivity to sub-polar North Atlantic variability, these data present the first opportunity to assess quantitatively the previously untestable hypothesis that ocean circulation dynamics might have played an active role in driving the key climate transitions of the last 1,000 years. We have therefore examined the timing of variability in the absolutely dated δ18O-shell series relative to reconstruction of NHSAT, derived from a composite dendrochronological series11 and individual tree ring series, and Greenland air temperatures, derived from an ice core stack13, over the last 1,000 years (Fig. 7 and Supplementary Figs 9 and 10). Lead–lag correlation analysis, conducted between the normalized δ18O-shell anomalies and NHSAT and Greenland air temperatures, indicates that over the pre-industrial period (AD 953–1800) multi-decadal to centennial frequency atmospheric temperature variability lagged changes in the δ18O-shell seawater density anomalies by 40±30 years (peak correlation r=−0.27 and −0.29, P<0.05 against NHSAT and Greenland air temperatures, respectively). Such a temporal offset is beyond the dating uncertainty of these records. The lead of North Icelandic shelf hydrographic variability over both the NHSAT and Greenland air temperature series strongly indicates that ocean variability played an active role in driving the major pre-industrial climate variability of the past 1,000 years. The coincident variability observed in the Greenland air temperature and wider NHSAT records argues against the possibility that the ocean lead we observe is driven by regional changes in atmospheric circulation patterns.

Figure 7: Lead-lag correlation analysis between the inverted δ18O-shell anomalies and NHSATs. NHSATs derived from a composite dendrochronological series11 and Greenland air temperatures derived from a stack of Greenland ice cores13. (a) Ten-year and (b) 50-year loess first-order low-pass-filtered δ18O-shell anomalies (red lines), Northern Hemisphere air temperatures (green lines) and Greenland air temperatures (blue lines). (c–h) Lead–lag correlation plots calculated between the δ18O-shell series and Northern Hemisphere air temperatures (green lines), the δ18O-shell series and Greenland air temperatures (blue lines), and Northern Hemisphere air temperatures and Greenland air temperatures (grey lines). The correlations are calculated over three periods. (c,f) Correlations calculated over the entire millennium. (d,g) Correlations calculated over the pre-industrial period (AD 1000–1800). (e,h) Correlations calculated over the industrial period (AD 1800–1997). Data used in the correlations in plots c–e were 10-year loess first-order low-pass filtered, whereas the data used to calculate the correlations in plots f–h were 50-year loess first-order low-pass filtered. The dashed black and red lines in plots (c–h) represent the respective 90 and 95% significance levels calculated using 1,000 Monte Carlo simulations using the Ebisuzaki method37. Full size image

The results of the lead–lag analysis of the reconstructions are further validated using the available ensemble of CMIP5 unforced pre-industrial control simulations (Supplementary Fig. 12). The CMIP5 preindustrial control simulations contain no external forcing; thus, all climate variability is internally generated9. Lead–lag analysis of the modelled North Icelandic shelf seawater density and NHSAT show that, consistent with the proxy records, the ocean variability leads the atmosphere by 10±10 years (Supplementary Fig. 12)18,19,35. The offset in the magnitude of ocean–atmosphere lag between model and real-world data are likely to be a result of differences in the detail of the surface and meridional circulations between models, and between the models and reality.

In comparison with the preindustrial interval, the modern trends (AD 1800–2000) in the NHSAT records relative to the δ18O-shell anomalies, coupled with the lead–lag analyses, suggests a reversal of the pre-industrial ocean–atmosphere coupling over the industrial period (AD 1800–2000) with changes in ocean dynamics more closely coupled in time, or lagging, increases in atmospheric surface air temperatures (Fig. 7). The close coupling identified by the proxy records over the industrial period is in agreement with the timing of coupling observed between instrumental NHSAT and sub-polar North Atlantic SSTs over the period AD 1880–2015 (Supplementary Fig. 11). The shift in coupling between the pre-industrial and industrial time periods is likely to be driven by a change in the balance between top-down and bottom-up forcings19.

The lead–lag correlations, coupled with the correlation of the δ18O-shell data with volcanic and high-frequency TSI forcing, suggest that at higher frequencies (<20 years) the ocean–atmosphere systems is closely coupled with periodic volcanic eruptions and high-frequency solar variability (for example, the sun spot cycle) driving a rapid response in both the atmosphere and ocean systems35. The lead–lag analyses suggest that at lower frequency domains (multi-decadal to centennial) marine variability is playing an active role in driving atmospheric variability. Such a divergence in the lagged response of the ocean and atmosphere systems at different frequency domains is detectable in the modern instrumental observations. Examination of a suite of instrumental SST and NHSAT’s indicates that over the twentieth century the atmosphere leads the ocean system at lower frequencies (multi-decadal; see Supplementary Fig. 11), probably driven by the faster response of the atmosphere to the warming influence of greenhouse gases. However, examination of the high-frequency variability (20-year high-pass-filtered data; Supplementary Fig. 11) implies a tighter coupling with high-frequency forcings/feedback mechanisms causing synchronous variability across both the ocean and atmosphere systems. These data suggest that during the pre-industrial period internal variability and feedback mechanisms within the North Atlantic substantially mediate the response of the climate system to top–down forcing (TSI, volcanic, atmospheric aerosols and greenhouse gases). Although, over the industrial period these data imply that internal oceanic mediation of top–down forcing has been overcome by the rate and nature of the NHSAT increase forced by increasing greenhouse gas concentrations.

These findings have implications for the interpretation of the modern climate system, as they highlight the problems that can result from using short modern oceanographic instrumental observations as a representative baseline of natural climate variability. Such records, although capturing a component of natural variability, are probably dominated by the anthropogenic signal. Furthermore, our data demonstrate shortcomings in methodologies that attempt to extend the modern instrumental ocean observations using absolutely dated records derived from terrestrial proxy networks (for example, see ref. 36). With the reconstructed reversal in the timing of ocean–atmosphere coupling, these results question the ability of terrestrial archives alone to reconstruct ocean variability and highlight the need to use independent absolutely dated marine archives to characterize the role ocean dynamics play in the global climate system.

The δ18O-shell record presented here provides the first precisely dated annually resolved record of past marine variability that spans the entire last 1,000 years. It provides a long-term baseline for the state of the North Atlantic coupled climate system and demonstrates that the role played by the ocean in naturally forced climate variability should be a key focus if we are to realise the societally crucial step forward in near-term climate prediction.