# Natural drivers of multidecadal Arctic sea ice variability over the last millennium

Paul Halloran, Ian Hall, Matthew B Menary, David Reynolds, James Scourse, James Screen, Alessio Bozzo, Nick Dunstone, Steven Phipps, Andrew Schurer, et al.

► To cite this version:

Paul Halloran, Ian Hall, Matthew B Menary, David Reynolds, James Scourse, et al., Natural drivers of multidecadal Arctic sea ice variability over the last millennium. Scientific Reports, 2020, 10 (1), pp.688. 10.1038/s41598-020-57472-2. hal-03147371

HAL Id: hal-03147371

https://hal.science/hal-03147371

Submitted on 21 Feb 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
OPEN

# Natural drivers of multidecadal
Arctic sea ice variability over the
last millennium

Paul R. Halloran¹*, Ian R. Hall², Matthew Menary³, David J. Reynolds²,⁴,
James D. Scourse⁵, James A. Screen⁶, Alessio Bozzo⁷, Nick Dunstone⁸, Steven Phipps⁹,
Andrew P. Schurer¹⁰, Tetsuo Sueyoshi¹¹,¹², Tianjun Zhou¹³ & Freya Garry¹

The climate varies due to human activity, natural climate cycles, and natural events external to the climate system. Understanding the different roles played by these drivers of variability is fundamental to predicting near-term climate change and changing extremes, and to attributing observed change to anthropogenic or natural factors. Natural drivers such as large explosive volcanic eruptions or multidecadal cycles in ocean circulation occur infrequently and are therefore poorly represented within the observational record. Here we turn to the first high-latitude annually-resolved and absolutely dated marine record spanning the last millennium, and the Paleoclimate Modelling Intercomparison Project (PMIP) Phase 3 Last Millennium climate model ensemble spanning the same time period, to examine the influence of natural climate drivers on Arctic sea ice. We show that bivalve oxygen isotope data are recording multidecadal Arctic sea ice variability and through the climate model ensemble demonstrate that external natural drivers explain up to third of this variability. Natural external forcing causes changes in sea-ice mediated export of freshwater into areas of active deep convection, affecting the strength of the Atlantic Meridional Overturning Circulation (AMOC) and thereby northward heat transport to the Arctic. This in turn leads to sustained anomalies in sea ice extent. The models capture these positive feedbacks, giving us improved confidence in their ability to simulate future sea ice in a rapidly evolving Arctic.

Sea ice is one of the most dynamic components of our climate system. Changes in sea ice amplify energy imbalance in the polar regions, and may have significant impacts on mid-latitude weather and extremes¹,². The Arctic sea ice decline observed in recent decades, while robustly linked to anthropogenic forcing³, includes a contribution from internal variability⁴ and volcanic activity⁵. While the short-term response of sea ice to volcanic forcing can be understood in terms of reduction of incoming radiation, atmospheric reorganisation and changes in surface ocean temperature⁶, which can be examined in the observational record⁷,⁸, the multidecadal to centennial response is likely to be mediated through large-scale ocean dynamics such as the inflow of warm water to the Arctic within the upper limb of the AMOC⁹⁻¹². Modelling studies have proposed a number of mechanisms through which volcanic activity can lead to decadal and longer changes in AMOC strength¹³⁻¹⁵, but agreement about how ocean circulation has changed on these timescales remains elusive. Understanding Arctic sea ice variability on multidecadal timescales therefore not only helps to better explain recently observed changes and potential predictability, but also to uncover the mechanisms of multidecadal variability in the wider climate system.

¹College of Life and Environmental Sciences, University of Exeter, Exeter, UK. ²School of Earth and Ocean Science, Cardiff University, Cardiff, CF10 3AT, UK. ³LOCEAN/IPSL, Sorbonne Universités (SU)-CNRS-IRD-MNHN, Paris, France. ⁴Laboratory of Tree Ring Research, The University of Arizona, Arizona, USA. ⁵College of Life and Environmental Sciences, University of Exeter, Penryn Campus, Treliever Road, Penryn, Cornwall, UK. ⁶College of Engineering, Mathematics and Physical Sciences, University of Exeter, Exeter, UK. ⁷Eumetsat, Darmstadt, Germany. ⁸Met Office Hadley Centre, FitzRoy Road, Exeter, EX1 3PB, UK. ⁹Institute for Marine and Antarctic Studies, University of Tasmania, Private Bag 129, Hobart, TAS, 7001, Australia. ¹⁰School of Geosciences, The University of Edinburgh, Edinburgh, UK. ¹¹National Institute of Polar Research, 10-3 Midori-cho, Tachikawa-city, 190-8518, Japan. ¹²Japan Agency for Marine-Earth Science and Technology, 3173-25, Showa-machi Kanazawa-ku, Yokohama City, 236-0001, Japan. ¹³LASG, Institute of Atmospheric Physics, Chinese Academy of Sciences, Beijing, 100029, China. *email: p.halloran@exeter.ac.uk
## Oxygen Isotope Variability Recorded from North Of Iceland

The Greenland-Iceland-Norway (GIN) Seas are characterised by a number of stratified water-masses with contrasting temperature and salinity properties. The position of these water masses is tightly coupled to that of the sea ice edge (Fig. 1b). The large area of land that drains into the Arctic, the net precipitative regime over the basin, and relatively fresh inflow through the Bering Strait results in a stratified low-salinity layer approximately 200 m thick sitting directly below the sea ice across much of the Arctic¹⁶. This low-salinity layer is made up of very fresh Polar Mixed Layer waters at the surface and the more saline Arctic Halocline waters below¹⁶ (Fig. 1a,b). Warm relatively saline waters enter the Arctic from the North Atlantic Current¹⁷, representing the upper limb of the AMOC (Fig. 1b). Where the North Atlantic Current waters meet the cold fresh waters, they are either entrained into the GIN Seas gyres, lose heat to the atmosphere and convectively mix with deeper Arctic waters¹⁸, or subduct beneath the fresher Arctic water-masses, where they fill much of the Arctic basin between about 200 and 800m below the surface¹⁷ (Fig. 1a,b). Within the region bounded by the winter sea ice edge, the Arctic Halocline waters are too fresh to mix convectively with the waters below¹⁶. The winter sea ice edge therefore marks the region where Atlantic waters subduct below the Halocline waters and, seaward of which, convection occurs (Fig. 1b).

The North Iceland Shelf, although it remains in communication with the subpolar North Atlantic at the surface via the Irminger Current¹⁸ (Fig. 1b), in the subsurface transects the water-bodies of the GIN Seas (Fig. 1c,d). The western edge of the North Iceland Shelf is connected along isopycnal surfaces with the Arctic Halocline and Polar Mixed Layer waters near the surface, and the subducted Atlantic waters near the shelf edge (Fig. 1c). In contrast, the eastern edge of the North Iceland Shelf is connected along isopycnal surfaces to the convectively mixed Atlantic/deep-Arctic waters (Fig. 1d). In the centre of this transect an annually-resolved stable oxygen isotope (δ¹⁸O) record has been constructed spanning the last millennium from aragonite samples drilled from shell growth bands belonging to an *Arctica islandica* bivalve chronology¹⁹,²⁰ (site location: 66°31.59′N, 18°11.74′W, 80 m water depth). The North Iceland bivalve δ¹⁸O timeseries, which is derived from aragonite precipitated in apparent equilibrium with the ambient seawater salinity and temperature²¹, displays marked multidecadal variability (Fig. 2b). Reynolds *et al*.²⁰, take advantage of the annual-resolution and absolute dating in this record to identify the relative timing of change occurring within the ocean and that occurring on land, as a way to ask whether the ocean is passively responding to atmospheric change over this interval, or whether ocean change led atmospheric change and therefore played an active role in the key climate transitions of the last millennium. These authors found that the change recorded in the bivalve timeseries tends to precede Northern Hemisphere air-temperature change, and therefore suggested that the ocean played an important role in driving atmospheric change. However, the specific drivers of the ocean change remain unresolved. While solar irradiance and volcanic aerosol forcing have been linked to the high-frequency variability in the bivalve δ¹⁸O timeseries, no direct coherent multidecadal relationship is found between the record and natural external forcing²⁰. Here we propose that the seawater δ¹⁸O recorded over time in bivalve shells from the North Iceland Shelf is faithfully documenting change in the relative influence of the different GIN Seas water masses at the site, as the sea ice extent varies in response to the integration and mediation of the external forcing by coupled climate system processes.

The expected isotopic composition of the bivalve aragonite can be predicted from thermodynamic equilibrium with the isotopic composition of the ambient water (δ¹⁸O<sub>water</sub>) at a given temperature, according to the empirical equation of Grossman and Ku²¹. Four distinct water-masses can be seen in δ¹⁸O<sub>water</sub>-space across the East-West GIN Seas transect, mapped out by calculating the δ¹⁸O of aragonite in equilibrium with seawater at observed temperatures and salinities²², described herein as the δ¹⁸O<sub>equil</sub>. (Fig. 3a) (see Methods).

## External Forcing of Sea Ice Extent

Observations and modelling studies have shown that sea ice perturbations can extend beyond the period over which the external forcing was present⁹,¹¹. This means that sequential forcing events can interfere with each other, obscuring the relationship between the forcing and the resultant sea ice timeseries. The complexity of the sea ice response to external forcing necessitates the use of climate models to aid interpretation of records of past sea ice change. Phase three of the Paleoclimate Model Intercomparison Project (PMIP3)²³ included simulations of the last millennium, forced by time varying land-use and greenhouse gas concentrations, orbital parameter changes, solar irradiance variability and volcanic aerosol concentration change²⁴. Radiative forcing over the pre-industrial component of the last millennium (850 Common Era (CE) to 1849 CE) was dominated by volcanic and solar activity²⁴. We find remarkable agreement between the PMIP3 multi-model mean GIN-Seas sea ice extent (mean of eight simulations, see Methods), and the annually resolved North Iceland Shelf δ¹⁸O recorded in the bivalve shells (Fig. 2b). Significant correlations of −0.58, +0.58 and −0.34 are seen respectively in the early (950–1200 CE), middle (1200–1600 CE) and late (1600–1849 CE) intervals of the study period over multidecadal timescales (see Methods). The early and late intervals present negative correlations, i.e. greater sea ice extent corresponding with a lower bivalve δ¹⁸O value (Fig. 2c,e), and the middle interval presents a positive correlation (Fig. 2d).

We propose that the three intervals identified in Fig. 2 correspond to three different hydrographic regimes where in the real ocean, but as discussed below not necessarily in models, multidecadal variability in the sea ice extent is superimposed on three different climatological sea ice states, and therefore water-mass geometries. The transition between regimes of Atlantic and Arctic water-mass influence at this site is supported by previous work²⁵. The first interval coincides with the Medieval Climate Anomaly, where there is widespread evidence for anomalous warmth in the region surrounding the North Atlantic²⁶, and little evidence of sea ice in the vicinity of the Icelandic coast²⁷ (Fig. 2a). Within this interval, similar to the situation at the present day, we propose that the boundary between the cool and moderately saline high δ¹⁸O<sub>eqil</sub>, convectively-mixed waters and the warm and saline low δ¹⁸O<sub>eqil</sub>, subducting Atlantic waters was close to the bivalve site. In this state, multidecadal variability in sea ice extent, i.e. relatively small changes around the climatological state for that period, will have varied the influence of these two water-masses over the δ¹⁸O values recorded by the bivalves (Fig. 4a). Under these conditions an increase in sea ice extent would push the subducting Atlantic waters eastward and reduce the δ¹⁸O<sub>eqil</sub>.
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

**Figure 1.** Oceanographic context. (a) Annual mean climatological (1955–2012 mean) temperature (colours) and salinity (contours) from a vertical section across 75°N (highlighted by the yellow bar), highlighting the different water-masses²². 75°N is chosen because it most clearly illustrated the water-mass structure. This water-mass structure is present throughout the open water of the GIN seas, but with the volume of Arctic-sources waters decreasing to the south. (b) Greenland-Iceland-Norway Seas schematic with bivalve site location (yellow circle), at 66°31.59′N, 18°11.74′W, 80 m water depth). c and d: Zonally averaged annual mean climatological (1955–2012 mean) temperature (contour colours, °C) and density (contour lines, kg.m⁻³) from the (c) West (averaged from 15.75–13.75°W) of Iceland and (d) East (averaged from 23.75–21.75°W) of Iceland²² as illustrated by the grey regions on the outline map of Iceland.
<figure><img src="image_9.png" /><figcaption>**Figure 2.** Correspondence of δ¹⁸O record and modelled Greenland-Iceland-Norway sea ice fraction. (a) Severity of ice off North Iceland coast calculated from documentary evidence of the duration of ice (total weeks) multiplied by the extent along the coast²⁷. (b) Area-averaged PMIP3 multi-model mean Greenland-Iceland-Norward sea sea ice fraction anomaly (blue) and bivalve Arctica islandica δ¹⁸O record (red), both timeseries are linearly detrended and sea ice fraction is normalised by its standard deviation. Between 1200 CE and 1600 CE the model sea ice timeseries is inverted (multiplied by -1). Note that following convention, the bivalve δ¹⁸O axis is inverted so that warming is represented by a positive change. Lower Panels: Correlation coefficients between the sea ice and δ¹⁸O timeseries for the intervals 950–1200 CE (c), 1200–1600 CE (d) and 1600–1849 CE (e) presented across a range of smoothing windows (see Methods). Stippling highlights significance (see Methods). Correlation coefficients are presented with smoothing windows of up to 20 years away from the 1:1 line to accommodate smoothing already introduced by calculating the multi-model mean. Maximum positive or negative correlations are -0.58, +0.58 and -0.34 for the intervals shown in panels c, d and e respectively, leading to a maximum variance explained in any phase (R²) of 34%.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 3.** Relationship of oceanography, temperature and salinity to δ¹⁸Oₐₔ₉. δ¹⁸Oₐₔ₉ is defined here as the δ¹⁸O of aragonite in equilibrium with seawater at observed temperatures and salinities. (a) Calculated δ¹⁸Oₐₔ₉ climatology from climatological temperature and salinity values²² from vertical section across the GIN Seas at 75°N (see Methods), with water-masses annotated. (b) Temperature (contour colours) and salinity (contour lines) from vertical section across 75°N within the GIN Seas Regional Climatology²². The solid black lines depict the boundaries between water masses. The dashed black lines highlight where the water stratifies in summer, leading seasonally to the warm surface layer evident in panel b.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 4.** Diagrammatic explanation of the influence of variability in sea ice extent on δ¹⁸Oₐₔ, at the bivalve site. The relationship between the GIN Seas water-masses and the bivalve site under the three different approximate climatological sea ice states proposed for the three intervals highlighted in Fig. 2, and the implications for δ¹⁸Oₐₔ, variability in response to multidecadal (i.e. small-scale) sea ice variability around these three climatological states.</figcaption></figure>

experienced at the bivalve site, and vice versa. Following a step increase in the documented incidence of sea ice off the coast of North Iceland around 1200 CE²⁷ (Fig. 2a) we propose that the bivalve site was close to the boundary between the low δ¹⁸Oₐₔ, subducting Atlantic waters and the cold and moderately saline Arctic Halocline waters. In this climatological state, multidecadal variability leading to an increase in sea ice extent would have increased the influence of low δ¹⁸Oₐₔ. Arctic Halocline waters at the bivalve site (Fig. 4b), the opposite δ¹⁸O response than occurred during the preceding interval. The most recent interval (1600–1850 CE) follows a large increase in recorded sea ice off North Iceland²⁷ (Fig. 2a) and corresponds to the expression of the Little Ice Age in the North Atlantic region²⁸. During this most recent interval, we propose that the increased climatological sea ice extent would have moved the boundary between the high δ¹⁸Oₐₔ. Arctic Halocline water and the cold and relatively fresh low δ¹⁸Oₐₔ. Polar Mixed Layer Waters to the site (Fig. 4c). In its 1600–1850 CE configuration, a multidecadal sea ice expansion would result in an increased influence of the low δ¹⁸Oₐₔ. Polar Mixed Layer waters on the site, and again a negative δ¹⁸O versus sea ice relationship.

Correlation between the PMIP3 multi-model-mean sea ice and the bivalve δ¹⁸O timeseries reinforces our interpretation of the δ¹⁸O record as a proxy for multidecadal sea ice variability and suggest that the models are successful in representing the dynamical feedbacks allowing persistent sea ice change. The correlation also indicates that a significant component of the multidecadal variability in sea ice extent is externally forced, with the maximum variance explained being 34% (Fig. 2c–e). The phase of internal variability within the individual model simulations will be independent from each other, and so will largely cancel when averaged together (Fig. S1). The remaining variability will typify the externally forced signal.

While the PMIP3 models appear to correctly simulate the multidecadal sea ice variability, they do not capture the shifts between positive and negative bivalve δ¹⁸O versus sea ice correlations (Fig. 2b). The shifts between
<figure><img src="image_12.png" /><figcaption>**Figure 5.** Proposed physical mechanism. Differences in sea ice fraction (a), surface ocean temperature (b), surface ocean salinity (c), surface ocean density (d) and normalised Atlantic Meridional Overturning stream function (e) composites from years with high and low sea ice extent within the multi-model mean sea ice timeseries, but excluding any years within the decade following a volcanic eruption (see Methods). Greyed out regions represent where fewer than 2/3rd of the model agreed on the sign of change.</figcaption></figure>

correlation regimes appear to arise from the interaction of stepwise sea ice extent advances from the Medieval Climate Anomaly into the Little Ice Age (Fig. 2a) with the detailed horizontal water-mass structure of the GIN Seas (Figs. 1 and 3a). The relatively coarse resolution PMIP3 models do not capture the observed water-mass structure of the GIN Seas (Fig. S2) and therefore, would not be expected to capture the switch in sign of the correlation between bivalve $\delta^{18}$O and sea ice timeseries. The transport of heat and freshwater from the Arctic coast into the interior of the basin to develop and sustain its stratification and therefore water-mass structure, occurs through eddies¹⁷. These fine-scale features will not be captured in low resolution model simulations.

## Physical Mechanism

Multi-model-mean composites of Northern Hemisphere sea ice (see Methods) indicate that the variability constrained by the bivalve $\delta^{18}$O record is representative of change in the GIN, Barents and Kara Seas (Fig. 5a). Consistent with previous studies⁹,¹¹ sea ice expansion coincides with a reduction in AMOC strength (Fig. 5e). The reduced heat transport to the Arctic via the AMOC, in response to both reduced AMOC strength and negative temperature anomalies, allows the expanded sea ice state to persist⁹,¹¹, potentially reflecting a positive feedback between sea ice and the AMOC⁹,¹¹,²⁹. While there is substantial model evidence that natural external forcing drives AMOC variability⁴,¹⁵,³⁰,³¹, the mechanisms and timescales of variability proposed by these studies differ¹³. Here multidecadal sea ice variability provides a proxy for AMOC variability, and agreement between modelled and reconstructed sea ice variability indicates that the PMIP3 models capture the feedback sustaining the sea ice anomalies on multidecadal timescales. Within this ensemble, it appears that the persistence of enhanced sea ice after external volcanic or solar forcing occurs in response to an AMOC-driven Atlantic-Arctic sea surface temperature cooling (Fig. 5a,b,e). A latitudinally coherent AMOC reduction is likely to occur in response to a density-driven reduction in high-latitude convection. A reduction in density in the sinking regions of the GIN Seas, northern Irminger Sea and perhaps Labrador Sea (Fig. 5d) can be largely attributed to salinity reductions in these regions (Fig. 5c). The reduction of salinity in the regions of active convective mixing is likely to result from
the import of freshwater supplied by melting sea ice⁹,¹¹ coming from the central Arctic and Barents Sea (Fig. 5c). The mechanism described here aligns with the suggestion that the reduction in heat-transport to the GIN Seas is dominated by a freshwater driven AMOC reduction²⁹ rather than a sea ice driven Subpolar Atlantic cooling¹¹. In summary, we find that initial sea ice perturbations occurring in response to the volcanic and solar forcing prescribed within the PMIP3 models²⁴ appear to be sustained through export of freshwater in the form of sea ice to regions of active deep convection. Reduced convection leads to reduced AMOC strength, and thereby, to reduced heat transport to the sea ice edge via the upper limb of the AMOC and ultimately to expansion of sea ice cover.

Accurate near-term prediction, attribution of climatic events, and quantification of weather extremes relies on the ability of climate models to simulate internally and externally-driven variability on a range of timescales. Here we have presented a marine proxy-tested model reconstruction of multidecadal Arctic sea ice variability over the last millennium and identified an important contribution of natural external forcing to this variability. We have shown that climate models are capable of capturing the timing of this variability, despite the forcing being integrated by, and the response arising through, complex coupled atmosphere, sea ice and ocean processes. These results give us confidence that models are capable of simulating the oceanographic feedbacks which are likely to occur in response to anthropogenically driven sea ice decline.

## Methods

**Calculation of seawater $\delta^{18}\text{O}_{\text{equil}}$.** Annually-averaged temperature and salinity data from the NOAA/ NCEI GIN Seas Regional Climatology ²² was converted to a synthetic $\delta^{18}\text{O}$ ($\delta^{18}\text{O}_{\text{equil}}$) value using Eqs. 1 and 2:

$$^{18}\text{Os}_{\text{equil}} = 0.55 * S - 18.98 \quad (1)$$

where the slope and intercept values are those identified for the North Atlantic in LeGrande and Schmidt³²

$$^{18}\text{O}_{\text{equil}} = \frac{21.8 - T}{4.69} + (^{18}\text{Os}_{\text{equil}} - 0.27) \quad (2)$$

following Reynolds et al.²⁰ and Grossman and Ku²¹. Where S is the salinity and T is temperature in °C.

**Calculation of multi-model mean Arctic sea ice.** The first ensemble member of all PMIP3 last millennium model simulations with the required variables, and compliant netcdf files were used, other than that of HadCM3³³, which was assumed to be superseded by the more recent Met Office Hadley Center model (HadGEM2-ES) simulations. The model output used was from CCSM4³⁴, CSIRO-Mk3L-1-2³⁵, GISS-E2-R³⁶, HadGEM2-ES³⁷,³⁸, MIROC-ESM³⁹,⁴⁰, MPI-ESM-P⁴¹, MRI-CGCM3⁴² and bcc-csm1-1⁴³.

Model output was interpolated onto a regular 1 by 1 degree grid using the Climate Data Operator remapbil function⁴⁴, and the sea ice fraction extracted for the region 45W-25E, 0N-90N and between the years 950CE and 1850CE, and averaged from monthly to yearly means. The area-weighted average timeseries was calculated for each model, this was linearly detrended, then the multi-model mean timeseries calculated.

**Correlation between multi-model mean sea ice and bivalve $\delta^{18}\text{O}$ timeseries.** The correlation coefficient between the multi-model-mean sea ice timeseries and the bivalve $\delta^{18}\text{O}$ timeseries presented in Fig. 2c-e is calculated under rolling boxcar smoothing windows ranging from 1 to 50 years. The correlation coefficient is calculated for the two timeseries under independent smoothing windows to accommodate the fact that the calculation of a multi-model-mean is itself smoothing the sea ice timeseries. Autocorrelation associated with smoothing of the timeseries means that it is not appropriate to assign significance from the p-value of the correlations. Significance is therefore defined here as where the correlation coefficient exceeds (or is more negative than, in the case of a negative correlation) the correlation coefficient of at least 95% of 100 randomly generated timeseries of the same length which have undergone the same smoothing. The synthetic timeseries were generated with a uniform distribution between 0 and 1. It should be noted that this significance test is not directly analogous to a p-value, and the precise outcome is conditional on the spectral characteristics of the synthetic timeseries. The timing of the partitioning of the modelled sea ice timeseries into three periods (Fig. 1b) was undertaken based upon the timing of step changes in the observation-based index of ice severity off North Iceland (Fig. 2a).

**Composite anomalies.** The first ensemble member of all PMIP3 last millennium model simulations with the required variables, and compliant netcdf files were used, other than that of HadCM3³³, which was assumed to be superseded by the more recent Met Office Hadley Center model (HadGEM2-ES) simulations. For sea ice and salinity analysis these models were CCSM4, CSIRO-Mk3L-1-2, HadGEM2-ES, MIROC-ESM, MPI-ESM-P, MRI-CGCM3, bcc-csm1-1, and for AMOC analysis (calculated from ocean meridional velocities) these were bcc-csm1-1, CCSM4, CSIRO-Mk3L-1-2, FGOALS-gl⁴⁵, MIROC-ESM, MPI-ESM-P and HadGEM2-ES.

Model output was interpolated onto a regular 1 by 1 degree grid using the Climate Data Operator remapbil function⁴⁴, and the calculated Atlantic stream function data interpolated onto common depth levels using the Python scipy RectBivariateSpline function. High and low sea ice years were defined as where the multi-model-mean sea ice timeseries (normalised by its standard deviation) was greater than 1 standard deviation and less than −1 standard deviation from the mean respectively, but excluding the 10 years following a volcanic eruption with a global total stratospheric sulfate aerosol injection greater than 15Tg⁴⁶. A subset of each 3D (latitude-longitude-time, or for AMOC analysis, latitude-depth-time) model field for each variable was selected, containing data from only the high and low sea ice years respectively. These fields were averaged along the time dimension, and the difference for each model between the high and low sea ice composites calculated. The multi-model mean was calculated from these composite difference fields.
Received: 9 October 2019; Accepted: 21 December 2019;

Published online: 20 January 2020

References

1. Screen, J. A. et al. Consistency and discrepancy in the atmospheric response to Arctic sea-ice loss across climate models. *Nat. Geosci.*, https://doi.org/10.1038/s41561-018-0059-y (2018).

2. Francis, J. A., Vavrus, S. J. & Cohen, J. Amplified Arctic warming and mid-latitude weather: new perspectives on emerging connections. *Wiley Interdiscip. Rev. Clim. Chang.*, https://doi.org/10.1002/wcc.474 (2017).

3. Min, S. K., Zhang, X., Zwiers, F. W. & Agnew, T. Human influence on Arctic sea ice detectable from early 1990s onwards. *Geophys. Res. Lett.*, https://doi.org/10.1029/2008GL035725 (2008).

4. Ding, Q. et al. Influence of high-latitude atmospheric circulation changes on summertime Arctic sea ice. *Nat. Clim. Chang.*, https://doi.org/10.1038/nclimate3241 (2017).

5. Rosenblum, E. & Eisenman, I. Faster Arctic sea ice retreat in CMIP5 than in CMIP3 due to volcanoes. *J. Clim.*, https://doi.org/10.1175/JCLI-D-16-0391.1 (2016).

6. Stenchikov, G. et al. Volcanic signals in oceans. *J. Geophys. Res. Atmos.*, https://doi.org/10.1029/2008JD011673 (2009).

7. Driscoll, S., Bozzo, A., Gray, L. J., Robock, A. & Stenchikov, G. Coupled Model Intercomparison Project 5 (CMIP5) simulations of climate following volcanic eruptions. *J. Geophys. Res. Atmos.* **117** (2012).

8. Robock, A. & Mao, J. Winter warming from large volcanic eruptions. *Geophys. Res. Lett.* **19**, 2405–2408 (1992).

9. Slawinska, J. & Robock, A. Impact of volcanic eruptions on decadal to centennial fluctuations of Arctic sea ice extent during the Last Millennium and on initiation of the Little Ice Age. *J. Clim.*, https://doi.org/10.1175/JCLI-D-16-0498.1 (2018).

10. Swingedouw, D. et al. Bidecadal North Atlantic ocean circulation variability controlled by timing of volcanic eruptions. *Nat. Commun.* **6**, 6545 (2015).

11. Zhong, Y. et al. Centennial-scale climate change from decadally-paced explosive volcanism: A coupled sea ice-ocean mechanism. *Clim. Dyn.*, https://doi.org/10.1007/s00382-010-0967-z (2011).

12. Mahajan, S., Zhang, R. & Delworth, T. L. Impact of the Atlantic Meridional Overturning Circulation (AMOC) on Arctic Surface Air Temperature and Sea-Ice Variability. *J. Clim.* **24**, 110624114117014 (2011).

13. Swingedouw, D. et al. Impact of explosive volcanic eruptions on the main climate variability modes. *Global and Planetary Change*, https://doi.org/10.1016/j.gloplacha.2017.01.006 (2017).

14. Otterå, O. H., Bentsen, M., Drange, H. & Suo, L. External forcing as a metronome for Atlantic multidecadal variability. *Nat. Geosci.*, https://doi.org/10.1038/ngeo955 (2010).

15. Mignot, J., Khodri, M., Frankignoul, C. & Servonnat, J. Volcanic impact on the Atlantic Ocean over the last millennium. *Clim. Past*, https://doi.org/10.5194/cp-7-1439-2011 (2011).

16. Rudels, B. Arctic Ocean circulation and variability - Advection and external forcing encounter constraints and local processes. *Ocean Sci.*, https://doi.org/10.5194/os-8-261-2012 (2012).

17. Spall, M. A. On the Circulation of Atlantic Water in the Arctic Ocean. *J. Phys. Oceanogr.*, https://doi.org/10.1175/JPO-D-13-079.1 (2013).

18. Brakstad, A., Våge, K., Håvik, L. & Moore, G. W. K. Water Mass Transformation in the Greenland Sea during the Period 1986–2016. *J. Phys. Oceanogr.*, https://doi.org/10.1175/jpo-d-17-0273.1 (2018).

19. Butler, P. G., Wanamaker, A. D., Scourse, J. D., Richardson, C. A. & Reynolds, D. J. Variability of marine climate on the North Icelandic Shelf in a 1357-year proxy archive based on growth increments in the bivalve Arctica islandica. *Palaeogeogr. Palaeoclimatol. Palaeoecol.* **373**, 141–151 (2013).

20. Reynolds, D. J. et al. Annually resolved North Atlantic marine climate over the last millennium. *Nat. Commun.*, https://doi.org/10.1038/ncomms13502 (2016).

21. Grossman, E. L. & Ku, T.-L. Oxygen and carbon isotope fractionation in biogenic aragonite: Temperature effects. *Chem. Geol.* **59**, 59–74 (1986).

22. Seidov, D. et al. Greenland-Iceland-Norwegian Seas Regional Climatology version 2. Regional Climatology Team, NOAA/NCEI (2018).

23. Braconnot, P. et al. The Paleoclimate Modeling Intercomparison Project contribution to CMIP5. CLIVAR Exch. (2011).

24. Schmidt, G. A. et al. Climate forcing reconstructions for use in PMIP simulations of the Last Millennium (v1.1). *Geosci. Model Dev.* **5**, 185–191 (2012).

25. Wanamaker, A. D. et al. Surface changes in the North Atlantic meridional overturning circulation during the last millennium. *Nat. Commun.* **3**, 899 (2012).

26. Diaz, H. F. et al. Spatial and temporal characteristics of climate in medieval times revisited. *Bull. Am. Meteorol. Soc.*, https://doi.org/10.1175/BAMS-D-10-05003.1 (2011).

27. Schell, I. I. pd. The Ice off Iceland and the Climates During the last 1200 Years, Approximately. *Geogr. Ann.* **43**, 354–362 (1961).

28. Ahmed, M. et al. Continental-scale temperature variability during the past two millennia. *Nat. Geosci.*, https://doi.org/10.1038/ ngeo1797 (2013).

29. Lehner, F., Born, A., Raible, C. C. & Stocker, T. F. Amplified inception of European little Ice Age by sea ice-ocean-atmosphere feedbacks. *J. Clim.*, https://doi.org/10.1175/JCLI-D-12-00690.1 (2013).

30. Swingedouw, D. et al. Bidecadal North Atlantic ocean circulation variability controlled by timing of volcanic eruptions. *Nat. Commun.*, https://doi.org/10.1038/ncomms7545 (2015).

31. Menary, M. B. & Scaife, A. a. Naturally forced multidecadal variability of the Atlantic meridional overturning circulation. *Clim. Dyn.* **42**, 1347–1362 (2014).

32. LeGrande, A. N. & Schmidt, G. A. Global gridded data set of the oxygen isotopic composition in seawater. *Geophys. Res. Lett.* **33** (2006).

33. Gordon, C. et al. The simulation of SST, sea ice extents and ocean heat transports in a version of the Hadley Centre coupled model without flux adjustments. *Clim. Dyn.* **16**, 147–168 (2000).

34. Landrum, L. et al. Last millennium climate and its variability in CCSM4. *J. Clim.* **26**, 1085–1111 (2013).

35. Phipps, S. J. et al. The CSIRO Mk3L climate system model version 1.0 - Part 1: Description and evaluation. *Geosci. Model Dev.* **4**, 483–509 (2011).

36. Schmidt, G. et al. Configuration and assessment of the GISS ModelE2 contributions to the CMIP5 archive. *J. Adv. Model. Earth Syst.* **6**, 141–184 (2014).

37. Collins, W. J. et al. Development and evaluation of an Earth-system model - HadGEM2. *Geosci. Model Dev.* **4**, 1051–1075 (2011).

38. Jones, C. D. et al. The HadGEM2-ES implementation of CMIP5 centennial simulations. *Geosci. Model Dev.* **4** (2011).

39. Sueyoshi, T. et al. Set-up of the PMIP3 paleoclimate experiments conducted using an Earth system model, MIROC-ESM. *Geosci. Model Dev.* **6**, 819–836 (2013).

40. Watanabe, S. et al. MIROC-ESM: model description and basic results of CMIP5-20c3m experiments. *Geosci. Model Dev. Discuss.* **4**, 1063–1128 (2011).

41. Giorgetta, M. A. et al. Climate and carbon cycle changes from 1850 to 2100 in MPI-ESM simulations for the coupled model intercomparison project phase 5. *J. Adv. Model. Earth Syst.* **5**, 572–597 (2013).
42. Yukimoto, S. et al. A new global climate model of the Meteorological Research Institute: MRI-CGCM3 - model description and basic performance. *J. Meteorol. Soc. Japan. Ser. II* **90A**, 23–64 (2012).

43. Xin, X.-G., Wu, T.-W. & Zhang, J. Introduction of CMIP5 Experiments Carried out with the Climate System Models of Beijing Climate Center. *Adv. Clim. Chang. Res.* **4**, 41–49 (2013).

44. Schulzweida, U. CDO User Guide (Version 1.9.6), https://doi.org/10.5281/zenodo.2558193 (2019).

45. Li, L. et al. The flexible global ocean-atmosphere-land system model, Grid-point Version 2: FGOALS-g2. *Adv. Atmos. Sci.* **30**, 543–560 (2013).

46. Gao, C., Robock, A. & Ammann, C. Volcanic forcing of climate over the past 1500 years: An improved ice core-based index for climate models. *J. Geophys. Res. Atmos.* **113**, 1–15 (2008).

**Acknowledgements**

This research was supported by the UK Natural Environmental Research Council grants NE/N001435/1, NE/N001176/1 and NE/N018486/1. J.A.S. was supported by The Leverhulme Trust (PLP-2015-215). S.P. was supported by the Australian Research Council's Special Research Initiative for the Antarctic Gateway Partnership (Project ID SR140300001). M.M. was supported by the EPICE project funded by the European Union's Horizon 2020 programme, Grant Agreement number 789445.

**Author contributions**

P.R.H., D.J.R., I.R.H. and J.D.S. conceived and designed the analysis. D.J.R., A.B., S.P., A.S., T.S. and T.Z. contributed data. P.R.H. performed the analysis and led the study. All authors contributed to the writing of the paper.

**Competing interests**

The authors declare no competing interests.

**Additional information**

Supplementary information is available for this paper at https://doi.org/10.1038/s41598-020-57472-2.

Correspondence and requests for materials should be addressed to P.R.H.

Reprints and permissions information is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_13.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2020