Sea ice and millennial-scale climate variability in the Nordic seas 90 kyr ago to present

In the light of rapidly diminishing sea ice cover in the Arctic during the present atmospheric warming, it is imperative to study the distribution of sea ice in the past in relation to rapid climate change. Here we focus on glacial millennial-scale climatic events (Dansgaard/Oeschger events) using the sea ice proxy IP25 in combination with phytoplankton proxy data and quantification of diatom species in a record from the southeast Norwegian Sea. We demonstrate that expansion and retreat of sea ice varies consistently in pace with the rapid climate changes 90 kyr ago to present. Sea ice retreats abruptly at the start of warm interstadials, but spreads rapidly during cooling phases of the interstadials and becomes near perennial and perennial during cold stadials and Heinrich events, respectively. Low-salinity surface water and the sea ice edge spreads to the Greenland–Scotland Ridge, and during the largest Heinrich events, probably far into the Atlantic Ocean.

D ansgaard/Oeschger (D/O) events in Greenland ice cores consist of warm interstadial (IS) and cold stadial events 1 and are strongly imprinted in sediments from the northern North Atlantic region and Nordic seas 2,3 . In general the warming to insterstadial conditions was abrupt as seen in Greenland ice cores and marine records. The warm conditions were followed by gradual cooling called the insterstadial transitional cooling phase, and a rapid transition to cold stadial conditions. Larger and/or longer-lasting stadials correlate with North Atlantic Heinrich events (H-events) 2 , where numerous icebergs were released from the Laurentide ice sheet and melting over the North Atlantic region in the so-called Ruddiman belt 4,5 (Fig. 1). Even though D/O events have been extensively studied, changes in sea ice cover have only been inferred by indirect evidence for presence or absence of sea ice (for example, deposition patterns of ice-rafted debris, oxygen isotope records and palaeo-temperature reconstructions) (Supplementary Fig. 1 and Supplementary Table 1).
The Nordic seas are characterized by northward inflow of warm, saline Atlantic surface Water (North Atlantic Current, Northwest Atlantic Current, North Atlantic surface Water, Faroe Current) and southward outflow of cold Polar surface Water (East Greenland Current and East Icelandic Current) 6 ( Fig. 1). In the Fram Strait, the Atlantic Water continues its flow below the sea ice-covered Polar surface Water as an intermediate water mass 7 . In the central part of the Nordic seas cooling and sinking of the salty surface water during the winter months generate cold deep overflows over the Greenland-Scotland Ridge into the North Atlantic 6,7 . The inflow of Atlantic surface Water is the major source of heat to the Arctic and Nordic seas, and it is generally agreed that changes in ocean circulation and sea ice cover has played a major role in the control of past millennial-scale climate changes of the glacial D/O events 2,8,9 . The Atlantic surface Water is ice-free throughout the year, while the East Greenland Current is covered by drifting nearperennial sea ice. In the central parts of the Nordic seas, mixing of Atlantic Water with Polar Water forms the zone of Arctic surface Water, which is located between the Arctic and Polar fronts 6 ( Fig. 1). The Arctic surface Water is seasonally sea ice covered and comprises the marginal ice zone (MIZ). The location of the MIZ and the Arctic and the Polar fronts changes with the seasons and on inter-annual and longer-time scales 10 . In the East Greenland Current behind the Polar front productivity is very low, while intermediate to high productivity is found in the ice-free zone of Atlantic surface Water. The highest seasonal  productivity occurs at the frontal areas and in the MIZ 11,12 . The positions of the Arctic and Polar fronts and the degree of sea ice cover thus depend on the distribution of the major surface water masses in the Nordic seas. A recent study showed that in the Arctic Ocean, the flow of Atlantic Water has a direct impact on sea ice distribution 13 . Previous studies of a C 25 isoprenoid lipid (IP 25 ) synthesized mainly by diatoms have shown its potential as a valuable new proxy for the reconstruction of the presence of seasonal sea ice [14][15][16][17][18][19] . IP 25 reportedly is produced by a few sea ice diatom species within the genera Haslea and Pleurosigma 20 . Furthermore, IP 25 is a stable organic compound preserved in sediments as old as Late Miocene and Pliocene [21][22][23] , and can be used in environments where other micropalaeontological sea ice proxies are absent or disturbed by dissolution effects 14,24 . Open-ocean conditions have been successfully reconstructed by using phytoplankton-derived sterols as a proxy for increased surface productivity such as brassicasterol and dinosterol 15,25,26 . The calculated ratio IP 25 to brassicasterol and dinosterol, the so-called P B IP 25 and P D IP 25 , respectively, can be related to sea ice coverage on a scale from permanent sea ice to year-round open water. In both ends of the scale, the IP 25 proxy is practically absent 15,16,26 , but accompanied by either low or high brassicasterol and dinosterol values, respectively 26,27 .
In this study we present results of IP 25 , brassicasterol, dinosterol, total organic carbon (%TOC) and d 13 C org (as terrigenous/marine organic-matter proxy; see ref. 28) together with the calculated sea ice indicators P B IP 25 and P D IP 25 , in the interval 801-2 cm of sediment core JM11-FI-19PC from the SE Norwegian Sea (see Methods and Fig. 1 for core location). In addition, low-resolution counts and identification of diatom frustules have been performed (see Methods). The purpose of this study is to reconstruct sea ice cover in the past in relation to millennial-scale climate change during the last 90 kyr in medium to high resolution. We compare our records with previously published records from the Fram Strait and published records that have provided indirect evidence for sea ice distribution and ocean circulation on millennial timescale from the northern North Atlantic (Fig. 1, Supplementary Fig. 1 (Fig. 1). The northern Faroe margin is presently sea ice free year round, but historical records show that during the Little Ice Age sea ice drift reached to the Faroe Islands via the East Icelandic Current 29 .
Our study shows that the sea ice retreats abruptly at the start of warm Interstadials and spreads rapidly during cooling phases of the Interstadials, before it becomes near-perennial and perennial  during cold stadials and H-events, respectively. The distribution of sea ice correlates closely with climate and variations in ocean circulation, in flow of meltwater and stratification of the ocean surface.

Results
Correlation and age model. The age model for the last 90 kyr interval of sediment core JM11-FI-19PC is based on wellknown tephra layers (Figs 2 and 3 and Supplementary Table 2 (Fig. 2). The characteristic saw-tooth pattern of the D/O events in the ice cores is recognized in the MS of our core ( Fig. 2). High values correlate with insterstadial and minimum values with stadial (S) climate, respectively. The presented interval of JM11-FI-19PC comprises the interstadials IS21-IS1 (IS1 ¼ Bølling and Allerød interstadials 1 ) and H-events H8-H1 (refs 8,31), which is equivalent to the last 90 kyr b2k (before year AD 2000) (Figs 2 and 3). The age control of the core JM11-FI-19PC younger than 65 kyr b2k has been presented by Ezat et al. 30 , while the age model for the interval older than 65 kyr b2k is new (see Methods and Supplementary Information for details). Sample resolution in the studied interval 90 kyr ago to present (per 1-cm thick sample) varies between 35 and 320 years, while the sedimentation rates vary between 0.28 mm per year (Fig. 3).
IP 25 and other biomarkers. For the last 30,000 years, trends in our PIP 25 record are similar to core MSM5/5-712-2 from the western Svalbard margin (for core location see Fig. 1) with very similar IP 25 , P B IP 25 and P D IP 25 values 17 (Fig. 4). In the Svalbard record, maxima in sea ice cover occurred in marine isotope stage (MIS) 2 (between 30 and 17 kyr ago), while the sea ice cover was slightly reduced after 17 kyr ago in the warm Bølling and Allerød interstadials (IS1) and in the Holocene 17 (Fig. 4). A very similar development is seen at the northern Faroe margin and with the same events clearly marked. The similarity of the data indicates that the sea ice proxy results for the Faroe margin are overall reliable. The minor offsets in the timing of the IP 25 , to P B IP 25 and P D IP 25 might reflect uncertainties in the age models and the fact, that the Svalbard record is located in the direct flow of Atlantic Water, while JM11-FI-19PC is also influenced by the IFF (Fig. 1).
For the part older than 30 kyr, sea ice variability (reflected in the IP 25 , P B IP 25 and P D IP 25 records) and related shortterm changes in surface-water productivity (reflected in the brassicasterol and dinosterol records and in the flux patterns of diatoms) seem to follow very clearly the D/O events (Fig. 5).   Interstadials are characterized by absence or near absence of IP 25 , maximum content of brassicasterol and dinosterol, and low P B IP 25 and P D IP 25 values (Fig. 5). The (near) absence of IP 25 indicates absence of sea ice-related diatoms, while high brassicasterol and dinosterol values indicate favourable conditions for phytoplankton growth in general and higher primary   productivity at the surface ocean. Thalassiosira oestrupii (Ostenfeld) Hasle, an indicator of relatively warm temperatures 32,33 is present in higher amounts in IS21 (MIS 5a) and in the Holocene, coinciding with high values of brassicasterol and dinosterol (Fig. 5, see Methods for details on diatoms). The species occurs sporadic or is absent in interstadials of MIS 3 and 2 (IS17-IS3) indicating colder than modern conditions. Altogether, the biomarker signals indicate open-ocean conditions and relatively high surface water temperatures over the northern Faroe Islands margin, with the Arctic Front located north of the core position (Fig. 6a). Planktic foraminiferal species from nearby core ENAM93-21/MD95-2009 indicated inflow of Atlantic Water during interstadials of MIS 3 and MIS 2 3 in support of the interpretation based on diatom floras and organic biomarkers. Furthermore, the d 13 C org values are higher in interstadials, which point to a (relative) decrease of terrigeneous organic matter supply and/or an increase in marine organic matter input 28 . The latter option, that is, an increased influence of marine organic matter at times of higher water temperatures due to higher marine productivity at the core location, is supported by higher flux of diatoms and higher concentrations of brassicasterol and dinosterol (Figs 5 and 6a). The optimum conditions are followed by a cooling phase defined by gradual changes to medium or higher values of IP 25 , brassicasterol, and dinosterol, P B IP 25 and P D IP 25 (Fig. 5). The relatively higher P B IP 25 and P D IP 25 values together with the increase in brassicasterol and dinosterol values (most clearly seen in IS21, IS20, IS13-10 and IS8; marked by arrows in Fig. 5) indicate increasing seasonal sea ice cover and elevated surface productivity, respectively, typical of the marginal ice zone 11,12 (Fig. 6b). These findings suggest that the IFF had moved to a probably variable position south-east of the core location and that the area was in the zone of Arctic surface water, likely resulting in decrease of both atmospheric and surface-water temperatures (Fig. 6b). This phase correlates with an increase in ice rafting and an increase in cold-water planktic foraminiferal species as seen in nearby core ENAM93-21/ MD95-2009 (ref. 3).

Discussion
The cooling phase terminated in a phase of extended sea ice cover. This is associated with a maximum in IP 25 values (higher than the calculated mean value), an abrupt decrease in the concentration of brassicasterol, dinosterol and high P B IP 25 and P D IP 25 values (Fig. 5, light-blue colour bars). Generally supporting this, the sea ice-associated Fragilariopsis oceanica (Cleve) Hasle 33-35 is present in higher relative abundance in H-event 8, S16 and S1 (Fig. 5, see Methods for details on diatoms). The lack of diatoms in most of the identified phases of extended sea ice cover (Fig. 5) is in accordance with the low abundances of diatoms found in areas with near-perennial sea ice cover today 16,20 (for diatom preservation, see Methods). With cold surface conditions and extended sea ice cover as far south as the core location, the Polar Front migrated to a position close to or probably south of the core site (Fig. 6c). Conditions were probably similar to conditions currently observed in areas governed by the East Greenland Current with (very) low primary productivity 36,37 and a dense pack-ice cover. Rather low d 13 C org values might indicate decreased marine organic carbon flux due to the low productivity and/or an increased influence of terrigenous organic matter, probably being ice-rafted to the core site. The presence of near-perennial sea ice prevented heat exchange and the atmosphere was probably very cold. These last phases of the D/O events represent the stadial conditions of the Greenland ice cores.
The sea ice proxy data for the H-events 6, 4 and 1 (Figs 4  and 5) show a very prominent period of perennial or nearperennial sea ice cover (absent or medium IP 25 , respectively, absent or minimal brassicasterol and dinosterol, and maximum P B IP 25 and P D IP 25 values). In these intervals, the IP 25 values of zero cannot be interpreted as absence of sea ice. Instead, the combined signals likely represent a thick permanent sea ice cover and very cold temperatures 15,26 (see also discussion in ref. 17). As long as the sea ice cover is thin enough for sunlight to penetrate, sea ice diatoms synthesizing IP 25 can grow attached beneath the ice 14,16,20 , which was probably the case at the beginning and end of some of the stadial intervals, where high peaks in P B IP 25 and P D IP 25 can occur (Fig. 5). When the sea ice becomes permanent and too thick for sunlight to penetrate, the signal of IP 25 will drop to zero (Fig. 5). To allow for the perennial sea ice and extreme cold temperatures, Polar surface water most likely spread over the area north of the Faroe Islands and the Polar Front was located far more southerly than today and during the smaller stadials (Fig. 6d). The transition from stadial and H-events to interstadials is rapid and interpreted as a sudden transition from extended sea ice cover to open-ocean conditions, seen as a decrease in values of IP 25 , P B IP 25  Our data generally show a good correlation between climate and sea ice cover for MIS 3, as well as for the last 30,000 years (Figs 4 and 5). We demonstrate that the presence/absence of sea ice varies closely in pace with the different climatic phases of the D/O millennial-scale climate events (Fig. 5). The peak warm interstadials with no sea ice (Fig. 7a) (Supplementary Fig. 1 and Supplementary Table 1) resembled the modern conditions, which have been shown by numerous marine core studies from the northern North Atlantic and Nordic seas, and have generally been interpreted as a sign of strong flow of Atlantic surface water (Fig. 6a). The transitional cooling phase of the insterstadial with gradually expanding sea ice cover from northwest (Figs 6b and 7b) correlate with increase in ice rafting from icebergs, decreasing atmospheric temperatures 1 and an increasing amount of meltwater over a larger area of the northern North Atlantic region and Nordic seas. In the following stadial events iceberg rafting reached a maximum (Figs 6c and 7c) as also seen in nearby core ENAM93-21/MD95-2009 and other records from the Nordic seas and North Atlantic (Supplementary Fig. 1 (Fig. 5). Sea ice advanced and the sea ice cover became near-perennial in the case of smaller stadial events (Figs 6c and 7c) or perennial as during the colder H-events H6, H4 and H1 (Figs 6d and 7d). The latter events are the three strongest events during the last 90 kyr, probably due to orbital forcing 2,38 . All stadial and H-events in the North Atlantic and Nordic seas show dominance by the polar planktic foraminifera Neogloboquadrina pachyderma sinistral (s) and cold polar conditions (see references in Supplementary  Information; Supplementary Fig. 1 and Supplementary Table 1).
Sea ice was an active player in millennial climate change, in most cases probably enforcing trends already caused by the predominantly cold, glacial atmospheric conditions during the last glacial period. The peak warmth of the interstadials lasted only shortly 1 and was immediately followed by cooling and spreading of sea ice. The inflow of Atlantic surface Water to the core area diminished to the extent that deep-water formation became very slow and stopped and sea ice cover became perennial or near perennial. The ocean circulation in the Nordic seas was probably more similar to the system in the northern Fram Strait today, in our view the closest analogue to the situation during stadials and North Atlantic H-events 38 and with similar circulation patterns of water masses with warmer Atlantic water flowing at intermediate depth below Polar surface Water [38][39][40][41] . In other words, the present-day conditions of the Fram Strait moved far south into the Atlantic Ocean (Fig. 7c,d). The remarkable abrupt disappearance of sea ice at the end of stadials/H-events correlates with sudden renewed inflow of Atlantic surface Water, peak insterstadial warmth and probably onset of deep-water formation. The distribution of sea ice thus correlated closely with variations in ocean circulation and expanding-retreating ice sheets and variations in flow of meltwater and stratification of the ocean surface.

Methods
Core logging and sampling. The core handling, sampling and diatom flora analyses of sediment core JM11-FI-19PC (62°49.97 N; 03°52.03 W; 1179 m water depth; 1109 cm core recovery) were performed at the Department of Geology, UiT, Arctic University of Norway, Tromsø, Norway. Before opening the magnetic susceptibility (MS) was measured with a Bartington MS meter MS2 (loop; sample resolution: 1 cm). After opening of the core, but before sampling, it was scanned for light and heavy elements, such as potassium (K) or titanium (Ti) with an Avaatech XRF core scanner in a resolution of 0.5-1 cm (ref. 30).
Samples taken at 5-cm intervals in 1-cm thick slices were weighed, freeze-dried and weighed again before wet sieving over 63 and 100-mm sieves. Tephra particles were counted in the size fraction 4100 mm. The method for counting tephra particles followed the procedures described in Wastegård and Rasmussen 42 .
Subsamples at 5-cm intervals from the presented core section (801-2 cm) were used to determine TOC. For %TOC dried and powdered aliquots of the samples were treated with 10% hydrochloric acid (HCl), for carbonate removal, before being measured with a Leco CS-200.
Carbon isotopes of organic material. Carbon isotopes d 13 C org were measured every 5 cm in a powdered and carbonate-free aliquot of the samples, using a Finnigan MAT Delta-S mass spectrometer equipped with a FLASH elemental analyser and a CONFLO III gas mixing system for the online determination of the carbon isotopic composition. Measurements were conducted at the Quantification of diatom floras. The preparation method for diatom samples given in Koc¸et al. 43 was followed, but using hydrogen peroxide (H 2 O 2 37%; 10 drops) instead of nitric acid (HNO 3 65%) to remove existing organic matter around the diatom valves. Quantitative diatom slides were produced as described in KocĶ arpuz and Schrader 44 . A sample was considered to contain enough diatoms for quantification, when there were more than 40 diatom valves ( ¼ 20 frustules) within B350 fields of view (2 vertical transects) and a minimum of 300 diatom valves per sample were counted. A total of 52 species that previously have been used by Koc¸Karpuz and Schrader 44 for establishing a sea-surface-temperaturetransfer function (including additional species added by Andersen et al. 45 ) have been identified and counted on a Chaetoceros spp. free basis 44,46 . Diatom counts followed the procedures of Schrader and Gersonde 47 and taxonomic identification mainly followed those of Hustedt 48 57 .
For the interval between 78 and B15 kyr b2k numerous fragmented frustules were found. This might indicate that this interval was influenced by dissolution processes 58 , which additionally could have altered the species composition of the diatom flora, and/or, that the sediments were directly influenced by the presence of sea ice as a mechanical component, breaking diatom valves by shear forces 59 . Diatom data for this interval have thus to be considered with caution and our reconstruction of past sea ice cover is exclusively based on the biomarker proxies that are more resistant and preserved in this type of sediments 14,16,23 (for general aspects of preservation of biomarkers see refs 28,60). Between B15 kyr b2k and the present the preservation of diatom frustules is excellent.
It is important, when using the PIP 25 index to distinguish between different sea ice conditions, that coevally high amounts of both biomarkers (suggesting ice-edge conditions) as well as coevally low contents (suggesting permanent-like ice conditions) would give a similar or even the same PIP 25 value. Especially, for the latter situation of permanent sea ice conditions both biomarker concentrations may approach values around zero and the PIP 25 index may become indeterminable (or misleading). For a correct interpretation of the PIP 25 data the individual IP 25 and phytoplankton biomarker concentrations must be considered 15,26 . Recently, Smik et al. 63 introduced a HBI-III alkene as phytoplankton biomarker replacing the sterols in the PIP 25 calculation. This modified PIP 25 approach is based on biomarkers from the same group of compounds (that is, HBIs) with more similar diagenetic sensitivity, which is important for palaeo-sea ice reconstructions and comparison of records from different Arctic areas.
Radiocarbon dating. Eleven AMS-14 C dates of piston core JM11-FI-19PC have previously been presented by Ezat et al. 30 . They were performed on monospecific samples of the planktic foraminiferal species N. pachyderma s (Supplementary Table 3). Measurements were performed at the 14 CHRONO Centre for Climate, the Environment, and Chronology, at Queens University Belfast, Northern Ireland. For this study all 14 C ages were calibrated to calendar years using the CALIB Radiocarbon Calibration 7.0.2. software and the Marine13 data set (including 400 year correction for surface reservoir ages) 64,65 . To make the age scale comparable to the ice core timescale (GICC05) 50 years were added. All ages are in ice core years b2k (before year AD 2000).
Construction of the age model. The age model younger than 65 kyr is based on the age model published in Ezat et al. 30 . The age model from 90 to 65 kyr b2k was constructed by using the same approach and at the same time validating the age model of the younger sediments published by Ezat et al. 30 (Supplementary  Tables 2, 3 and 4). To further solidify the age model, piston core JM11-FI-19PC was not only correlated to the NGRIP ice core but also to the nearby compiled marine cores ENAM93-21 (ref. 3  Interstadial transitional cooling phase with variable sea ice cover and variable positions of AF (IFF) and PF, which is indicated by white arrows. The core location is always between AF and PF. (c) Stadial with extended sea ice cover. The core location is again between AF and PF, but markedly closer to PF (potentially PF can be south of the core location; see text). Both fronts are located further towards the south compared with interstadial and interstadial cooling conditions. (d) H-events (H6, H4 and H1) with perennial sea ice cover. Both Fronts, AF and PF, are now located south of the core location. Stars mark the location of the studied core JM11-FI-19PC, as well as the location of core MSM5/5-712-2 17 used for comparison; dashed white lines give approximate positions of fronts AF (IFF) and PF. The reconstructions of regional distribution of sea ice cover are based on published records together with the records of JM11-FI-19PC and MSM5/5-712-2 17 ( Supplementary Fig. 1 MS, common tephra layers and stable oxygen isotopes (d 18 O) measured in the planktic foraminiferal species N. pachyderma s ( Fig. 2; see also Supplementary Table 2). The ages in between the tie points and/or tephra layers were calculated by interpolation assuming a constant sedimentation rate (Fig. 3).
Data availability. Data referenced in this study are available in PANGAEA with the identifier DOI 'https://doi.org/10.1594/PANGAEA.859992'.