Observations and modelling studies have shown that sea ice perturbations can extend beyond the period over which the external forcing was present9,11. 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)23 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 change24. Radiative forcing over the pre-industrial component of the last millennium (850 Common Era (CE) to 1849 CE) was dominated by volcanic and solar activity24. 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 δ18O 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 δ18O 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 work25. The first interval coincides with the Medieval Climate Anomaly, where there is widespread evidence for anomalous warmth in the region surrounding the North Atlantic26, and little evidence of sea ice in the vicinity of the Icelandic coast27 (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 δ18O eqil. convectively-mixed waters and the warm and saline low δ18O eqil. 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 δ18O 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 δ18O eqil. 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 CE27 (Fig. 2a) we propose that the bivalve site was close to the boundary between the low δ18O eqil. 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 δ18O eqil. Arctic Halocline waters at the bivalve site (Fig. 4b), the opposite δ18O response than occurred during the preceding interval. The most recent interval (1600–1850 CE) follows a large increase in recorded sea ice off North Iceland27 (Fig. 2a) and corresponds to the expression of the Little Ice Age in the North Atlantic region28. During this most recent interval, we propose that the increased climatological sea ice extent would have moved the boundary between the high δ18O eqil. Arctic Halocline water and the cold and relatively fresh low δ18O eqil. 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 δ18O equil. Polar Mixed Layer waters on the site, and again a negative δ18O versus sea ice relationship.

Figure 4 Diagrammatic explanation of the influence of variability in sea ice extent on δ18O equil. 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 δ18O equil. variability in response to multidecadal (i.e. small-scale) sea ice variability around these three climatological states. Full size image

Correlation between the PMIP3 multi-model-mean sea ice and the bivalve δ18O timeseries reinforces our interpretation of the δ18O 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 δ18O versus sea ice correlations (Fig. 2b). The shifts between 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 δ18O 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 watermass structure, occurs through eddies17. These fine-scale features will not be captured in low resolution model simulations.