# Precession and atmospheric CO2 modulated variability of sea ice in the central Okhotsk Sea since 130,000 years ago

Li Lo, Simon Belt, Julie Lattaud, Tobias Friedrich, Christian Zeeden, Stefan Schouten, Lukas Smik, Axel Timmermann, Patricia Cabedo-Sanz, Jyh-Jaan Huang, et al.

► To cite this version:

Li Lo, Simon Belt, Julie Lattaud, Tobias Friedrich, Christian Zeeden, et al. Precession and atmospheric CO2 modulated variability of sea ice in the central Okhotsk Sea since 130,000 years ago. Earth and Planetary Science Letters, 2018, 488, pp.36-45. 10.1016/j.epsl.2018.02.005. hal-02304989

HAL Id: hal-02304989

https://hal.science/hal-02304989v1

Submitted on 19 Jun 2025

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
## Highlights

* The first orbital timescale proxy-model sea ice-sea surface temperature records from the northwestern subarctic Pacific Ocean.

* Strong precession forcing controlled sea ice variations are modulated by greenhouse gas radiative forcing.

* Sea ice remained free in the central Okhotsk Sea during MIS 5e due to high greenhouse gas radiative forcing.
Precession and atmospheric CO₂ modulated variability of sea ice in the central
Okhotsk Sea since 130,000 years ago

Li Lo<sup>a,b,*</sup>, Simon T. Belt<sup>c</sup>, Julie Lattaud<sup>d</sup>, Tobias Friedrich<sup>e</sup>, Christian Zeeden<sup>f</sup>, Stefan

Schouten<sup>d,g</sup>, Lukas Smik<sup>c</sup>, Axel Timmermann<sup>e,h</sup>, Patricia Cabedo-Sanz<sup>c</sup>, Jyh-Jaan

Huang<sup>i,#</sup>, Liping Zhou<sup>j</sup>, Tsong-Hua Ou<sup>k</sup>, Yuan-Ping Chang<sup>l</sup>, Liang-Chi Wang<sup>m</sup>, Yu-

Min Chou<sup>n</sup>, Chuan-Chou Shen<sup>i</sup>, Min-Te Chen<sup>o</sup>, Kuo-Yen Wei<sup>i</sup>, Sheng-Rong Song<sup>i</sup>,

Tien-Hsi Fang<sup>o</sup>, Sergey A. Gorbarenko<sup>p</sup>, Wei-Lung Wang<sup>q</sup>, Teh-Quei Lee<sup>r</sup>, Henry

Elderfield<sup>b, †</sup>, David A. Hodell<sup>b</sup>

<sup>a</sup>State Key Laboratory of Isotope Geochemistry, Guangzhou Institute of Geochemistry, Chinese Academy of Sciences, Guangzhou 510640, China.

<sup>b</sup>Godwin Laboratory for Palaeoclimate Research, Department of Earth Sciences, University of Cambridge, Downing Street, Cambridge CB2 3EQ, UK.

<sup>c</sup>Biogeochemical Research Centre, School of Geography, Earth and Environmental Sciences, Plymouth University, Plymouth PL4 8AA, UK.

<sup>d</sup>NIOZ, Royal Netherlands Institute for Sea Research, Department of Marine Microbiology and Biogeochemistry, and Utrecht University, PO Box 59, 1790 AB Den Burg, The Netherlands.

<sup>e</sup>International Pacific Research Center, School of Ocean and Earth Science and Technology, University of Hawaii at Monoa, Honolulu, Hawaii 96822, USA.

<sup>f</sup>IMCCE, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Universités, UPMC Univ Paris 06, Univ Lille, 75014 Paris, France.

<sup>g</sup>Faculty of Geosciences, Department of Earth Sciences, Utrecht University, PO Box 80.021, 3508 TA, The Netherlands.

<sup>h</sup>Institute of Basic Science, Center for Climate Physics, Pusan National University, Busan 46241, Republic of Korea.

<sup>i</sup>Department of Geosciences, National Taiwan University, Taipei 10617, Taiwan ROC.

<sup>j</sup>Laboratory for Earth Surface Processes, Department of Geography, and Institute of Ocean Research, Peking University, Beijing 100871, China.

<sup>k</sup>Institute of Applied Mechanics, National Taiwan University, Taipei 10617, Taiwan ROC.

<sup>l</sup>Department of Oceanography, National Sun Yat-sen University, Kaohsiung 80424, Taiwan ROC.

<sup>m</sup>Collection Management Department, National Taiwan Museum, Taipei 10046, Taiwan ROC.

<sup>n</sup>Department of Ocean Science and Engineering, Southern University of Science and Technology, Shenzhen, Gongdong 518055, China.

<sup>o</sup>Institute of Applied Geosciences, National Taiwan Ocean University, Keelung 20224, Taiwan ROC.

<sup>p</sup>V.I. Il'ichev Pacific Oceanological Institute, Far East Branch Russian Academy of Science, Vladivostok 690041, Russia.

<sup>q</sup>Department of Biology, National Changhua University of Education, Changhua 50074, Taiwan ROC.

<sup>r</sup>Institute of Earth Sciences, Academia Sinica, Taipei 11529, Taiwan ROC.

*To whom correspondence should be addressed
Email address: lilo@gig.ac.cn (Li Lo)

#now at Institute of Geology, University of Innsbruck, Innrain 52, A-6020, Innsbruck, Austria

†Deceased

Revised and resubmitted to *Earth and Planetary Science Letters*
January 30, 2018
51 ABSTRACT

Recent reduction in high-latitude sea ice extent demonstrates that sea ice is highly sensitive to external and internal radiative forcings. In order to better understand sea ice system responses to external orbital forcing and internal oscillations on orbital timescales, here we reconstruct changes in sea ice extent and summer sea surface temperature (SSST) over the past 130,000 years in the central Okhotsk Sea. We applied novel organic geochemical proxies of sea ice (IP<sub>25</sub>), SSST (TEX<sup>L86</sup>) and open water marine productivity (a tri-unsaturated highly branched isoprenoid and biogenic opal) to marine sediment core MD01-2414 (53°11.77′N, 149°34.80′E, water depth 1123 m). To complement the proxy data, we also carried out transient Earth system model simulations and sensitivity tests to identify contributions of different climatic forcing factors. Our results show that the central Okhotsk Sea was ice-free during Marine Isotope Stage (MIS) 5e and the early-mid Holocene, but experienced variable sea ice cover during MIS 2-4, consistent with intervals of relatively high and low SSST, respectively. Our data also show that the sea ice extent was governed by precession-dominated insolation changes during intervals of atmospheric CO<sub>2</sub> concentrations ranging from 190 to 260 ppm. However, the proxy record and the model simulation data show that the central Okhotsk Sea was near ice-free regardless of insolation forcing throughout the penultimate interglacial, and during the Holocene, when atmospheric CO<sub>2</sub> was above ~260 ppm. Past sea ice conditions in the central Okhotsk Sea were therefore strongly modulated by both orbital-driven insolation and CO<sub>2</sub>-induced radiative forcing during the past glacial/interglacial cycle.

73

74 **Keywords:** Okhotsk Sea, Sea ice, Insolation, Greenhouse gases, Precession cycle

75

76 **Highlights:**

1. The first orbital timescale proxy-model sea ice-sea surface temperature records from the northwestern subarctic Pacific Ocean.

2. Strong precession forcing controlled sea ice variations are modulated by greenhouse gas radiative forcing.

3. Sea ice remained free in the central Okhotsk Sea during MIS 5e due to high greenhouse gas radiative forcing.

83
## 1. Introduction

Sea ice is one of the crucial components of the Earth's climate system, in part, due to its high reflectance or albedo, which influences energy budgets at both high and low latitudes (Serreze et al., 2016; Turner et al., 2016). In addition, brine rejection during seasonal sea ice formation supplies dense and well ventilated water to global deep/intermediate water circulation, while sea ice melt in spring causes stratification between near-surface and deeper water masses. Furthermore, the area bound by the retreating sea ice margin during spring and summer (the so-called marginal ice zone) represents a region of significant open water (pelagic) productivity.

Major changes in the extent and thickness of sea ice across the Arctic-subarctic and Antarctic regions during the last three to four decades have been revealed by direct and remote sensing observations (Serreze et al., 2016; Turner et al., 2016). In the Arctic, the average reduction rate in September sea ice extent, 13.4% per decade during 1981-2010, is higher than predicted by most model simulations (Serreze et al., 2016). This discrepancy emphasizes the limited knowledge of the high-latitude climate system to radiative perturbations and natural variability. In addition, a slow overall increase in annual mean sea ice extent in the Southern Ocean (ca. 1.6% per decade during 1979-2013), raises further questions regarding the relationships between anthropogenic radiative forcings, temperature, winds, and sea ice cover (Turner et al., 2016).

The Okhotsk Sea in the subarctic Pacific Ocean has responded rapidly to recent climate change (Mesquita et al., 2011; Kashiwase et al., 2014) and has been shown previously to be dynamically connected to Northern Hemispheric cooling events (e.g. Heinrich and Dansgaard-Oeschger (D/O) events) during the last glacial period (Ono et al., 2005; Sakamoto et al., 2005; Harada et al., 2008; Max et al., 2012; 2014). For example, Harada et al. (2008) demonstrated a strong link between a $U^{K'37}$-derived sea surface temperature (SST) record and Greenland Ice Sheet Project 2 ice core D/O events back to ~120 ka in the southwestern Okhotsk Sea. Further, sedimentary magnetic mineral and composition data have been used to reconstruct sea ice derived ice rafted debris (IRD) in the same region (Sakamoto et al., 2005) with larger fluctuations during last glacial period (Marine Isotope Stage, MIS 2-4) compared to interglacial periods (Holocene and MIS 5) for the southwestern Okhotsk Sea. Also within the study region, freshwater input from the Amur River and polar atmospheric dynamics are potential
candidates responsible for controlling sea ice dynamics in the southwestern Okhotsk Sea (Sakamoto et al., 2005; Harada et al., 2008). Previous low resolution sea ice reconstructions based on the Arctic sea ice biomarker proxy IP<sub>25</sub> (Ice Proxy with 25 carbon atoms, Belt et al., 2007; Belt and Müller, 2013), when combined with compiled SST data, have shown that sea ice extent variations in the central-west subarctic Pacific Ocean are tightly link to Atlantic meridional overturning circulation and atmospheric circulation between North Atlantic and North Pacific during the last termination (Max et al. 2012; 2014). More recently, Méheust et al. (2016) used IP<sub>25</sub> and other geochemical proxies to show that sea ice expanded significantly during Heinrich event 1 (H1) and the Younger Dryas (YD) in the western Bering Sea, in contrast to the Bølling-Allerød (B/A) and early Holocene, which experienced low/absent sea ice.

Despite the importance of its location and role in the subarctic sea ice system, very few studies have been reported from the central region of the Okhotsk Sea (Liu et al., 2006; Wang and Wang, 2008; Chou et al., 2011) and no detailed sea ice reconstructions have been conducted on orbital timescales. As such, the roles of insolation (external) and greenhouse gas radiative forcing (internal) on sea ice variation remain poorly understood. Therefore, long-term reconstruction of sea ice extent in the central Okhotsk Sea would be especially informative in understanding the interaction of sea ice and external/internal climatic forcings.

In this study, we reconstruct variations in sea ice extent and summer SST (SSST) using organic geochemical proxies and compare the proxy-derived records with Earth system modelling results in order to identify controlling mechanisms and interactions between sea ice and atmospheric-oceanic forcings in the Okhotsk Sea during the past 130,000 years. Our results reveal a strong precession control and a potential greenhouse gas induced radiative forcing threshold on sea ice variations since the penultimate peak interglacial period.

## 2. Regional setting

The Okhotsk Sea represents the southernmost region of contemporary seasonal sea ice formation in the Northern Hemisphere and has experienced a large decline rate of 11.4% per decade in sea ice extent during the past three decades (Kashiwase et al., 2014). The regional current system is affected by the north- and south-ward flowing
152 West Kamchatka Current (WKC), East Sakhalin Current (ESC), salty-warm Soya
153 Warm Current (SWC) and freshwater input from the Amur River (Fig. 1). The SSST is
154 in the range 5–13°C, while salinity varies between 31.5 to 33.2‰, being influenced
155 mostly by the Amur River discharge (Luchin et al., 2009). Sea ice easily forms on the
156 shallow continental shelves in the north-west Okhotsk Sea and is also influenced by the
157 large fresh water input by the Amur River. As a result, the Okhotsk Sea is the
158 southernmost region of subarctic sea ice distribution in the world (Kimura and
159 Wakatsuchi, 2004; Nishioka et al., 2014).

160

161 **3. Material and methods**

162 **3.1 Sediment core MD01-2414 and surface sediments**

163 Marine sediment core MD01-2414 (53°11.77′N, 149°34.80′E, water depth 1123
164 m, total sediment length 52.76 m, Chou et al., 2011, Supplementary Fig. 1) was drilled
165 from the Deryugin Basin during the circum-Pacific initiative cruise in 2001 as part of
166 the IMAGES project (Fig. 1). The upper 235-cm segment was not suitable for XRF
167 scanning due to its high water content (Liu et al., 2006). The main sediment
168 composition is terrestrial detritus (from sand to silty clay) with diatom and rare
169 calcareous (nannofossil and foraminifera) oozes (Liu et al., 2006; Wang and Wang,
170 2008; Chou et al., 2011). The ash layers and core gaps resulting from the coring process
171 were eliminated to prevent bias when performing sedimentation rate calculations (Chou
172 et al., 2011). In this study, we used sediment of the upper 700 cm to generate sea ice
173 and SSST histories covering the past 130,000-yrs. To complement the long-term
174 records, we also examined biomarker content in a number of surface sediments from
175 the region (Supplementary Fig. 2, Supplementary Table 3).

176

177 **3.2 Sediment core XRF scanning**

178 Non-destructive X-ray fluorescence (XRF) scanning was performed by the
179 ITRAX @ COX company. Continuous downcore measurements of elemental variations
180 were done in the ITRAX-XRF Core Scanner Laboratory, Department of Geosciences,
181 National Taiwan University (Huang et al., 2016). A U-channel of core MD01-2414 was
182 scanned using the 3 kW Mo source. The XRF measurements were analyzed at 30 kV,
183 24 mA, 2 mm resolution with a 30 second exposure time. The original XRF spectra
184 were processed by the Q-Spec software provided by COX Analytical Systems to obtain
185 element peak areas in counts.
186

187 3.3 *Age model*

188 The MD01-2414 age model was established based on accelerator mass
189 spectrometry radiocarbon (AMS ¹⁴C) dates and XRF data correlation to the global
190 composite benthic foraminiferal oxygen isotope curve (LR04, Fig. 3H; Lisiecki and
191 Raymo, 2005). Samples of the planktonic foraminifera *Neogloboquadrina pachyderma*
192 (>125 µm, *sinistral*) were picked from five depths (33, 113, 143, 170, and 210 cm).
193 AMS ¹⁴C dates were obtained by the College of Urban and Environmental Sciences,
194 Peking University. AMS ¹⁴C dates were calibrated to calendar ages using CALIB 7.1
195 software (Reimer et al., 2013) with a reservoir age calculated from the Marine
196 Reservoir Correction Database (http://calib.qub.ac.uk/marine/). Four sites near the
197 Okhotsk Sea were selected (Reimer et al., 2013). The calculated weighted mean ΔR
198 value was 450 ± 90 years. Radiocarbon results are listed in Supplementary Table 1.

199

200 ITRAX data show that Ba and Ti are negatively correlated with each other through
201 the whole study section (r² = 0.58, Fig. 2). Ba and Ti intensities vary from 340-680
202 counts and 1100-4800 counts, respectively. ITRAX data have also been confirmed by
203 traditional XRF measurements with lower time resolution (Liu et al., 2006). The log
204 (Ba/Ti) ratio was correlated with LR04 for the upper 7-m and shows a good correlation
205 with other physical parameters, including magnetic susceptibility (Chou et al., 2011),
206 color reflectance (Bassinot and Chen, 2002) and coarse fraction (C.F., this study; Fig.
207 2). We assumed the log (Ba/Ti) ratio to represent a biological/terrestrial input ratio for
208 our study site and this ratio varies with glacial/interglacial (G/IG) sea level changes due
209 to global ice volume variations. During glacial periods, the west and north parts of the
210 continental shelf were exposed in the Okhotsk Sea and more terrestrial sediment was
211 transported to the central basin. On the other hand, during intervals of sea level rise, the
212 central Okhotsk Sea became a region associate with a major biogenic bloom. Such
213 biogenic/terrestrial distributions are supported by a previous surface sediment study
214 (Strakhov et al., 1961). The similarity of the log(Ba/Ti) ratio of Site MD01-2414 to the
215 LR04 is shown in Fig. 3H, and it persists for older intervals (Supplementary Fig. 1).
216 The LR04 time is based on comparison to the Imbrie and Imbrie (1980) ice model, and
217 the fit for the last glacial cycle is excellent between the model, LR04, and the log(Ba/Ti)
218 ratio. The Imbrie and Imbrie (1980) model is delayed relative to insolation, and the
219 Okhotsk Sea can be assumed to have reacted to insolation no more than large Northern
220 Hemisphere ice sheets. A shift by half a precession cycle (or more) is, therefore, not considered realistic. Thus, we correlate log(Ba/Ti) directly to LR04, with no lags, and the magnetic susceptibility and major paleo-polar reversal events support this age model for the whole sediment section back to 1.55 million years ago (Supplementary Fig. 1).

225
226 The average sedimentation rate in core MD01-2414 is 9.5 cm kyr⁻¹ in the upper
227 235 cm, which indicates some piston stretching during the coring process. For the rest
228 of the sections, sedimentation rates are 2–4 cm kyr⁻¹, similar to that found for nearby
229 sites (Gorbarenko et al., 2010; Nürnberg et al., 2011).

230
231 3.4 Sea ice and phytoplankton biomarkers

232 Biomarker analysis was concentrated on the sea ice proxy IP25 (Belt et al., 2007; 233 Belt and Müller, 2013) and a further highly branched isoprenoid (HBI) lipid associated 234 with certain pelagic diatoms (HBI III; Belt et al., 2015). In the absence of a detailed 235 surface sediment-based calibration, we refrained from using PIP25 data and, instead, 236 focus on the distributions of the individual biomarkers. Part of the interpretation 237 described here is based around the absence of IP25 in certain sediment horizons. In some 238 previous investigations, absent IP25 has been interpreted in terms of representing either 239 ice-free or permanent ice cover, although this is likely an over-simplification, not least 240 as IP25 has been reported in regions of near permanent ice cover in the central Arctic 241 Ocean (Xiao et al., 2015). To distinguish between these two extreme conditions of sea 242 ice cover, here we use a combination of absent IP25 alongside HBI III and biogenic opal 243 data (Liu et al., 2006), the latter parameter having been used previously in the 244 neighboring western Bering Sea to identify intervals of higher productivity associated 245 with open water settings (Méheust et al., 2016). More recently, the combined IP25, HBI 246 III and biogenic opal approach was used to reveal changes in sea ice dynamics across 247 the Mid-Pleistocene Transition, also from the Bering Sea (Detlef et al., 2018). Further 248 indications of ice-free settings have been inferred from accompanying SST data, also 249 described herein.

250
251 IP25 and HBI III concentration data were obtained using methods described
252 previously (e.g. Belt et al., 2015). Briefly, 2–3 g of freeze dried sediment material was
253 extracted (dichloromethane/methanol; 3 x 12 mL; 2:1 v/v) by ultrasonication (15-min)
and centrifugation (2500 rpm; 1 min) following addition of internal standards (9-octylheptadec-8-ene, 9-OHD, 10 µL; 10 µg mL⁻¹) for quantification purposes. Dried (nitrogen) total organic extracts were re-dissolved in hexane (ca. 1 mL) and purified using column chromatography (silica), with IP₂₅ and HBI III (hexane; 6 mL) collected as a single fraction. Non-polar lipid fractions were further separated into saturated and unsaturated hydrocarbons using glass pipettes containing silver ion solid phase extraction material (Supelco Discovery® Ag-Ion). Saturated hydrocarbons were eluted with hexane (1 mL), while unsaturated hydrocarbons (including IP₂₅ and HBI III) were eluted with acetone (2 mL). All fractions were analyzed using gas chromatography–mass spectrometry (GC–MS) and operating conditions were as described previously (Belt et al., 2015). Mass spectrometric analyses were carried out either in total ion current or single ion monitoring mode. Identification of individual lipids was achieved by comparison of their characteristic GC retention times and mass spectra with those of reference compounds. Lipid quantification was achieved by dividing peak area integrations of selected ions (m/z 350 (IP₂₅); 346 (HBI III)) by those of the internal standard (m/z 350 (9-OHD)) in single ion monitoring mode, with these ratios then normalized according to their respective instrumental response factors and sediment masses (Belt et al., 2015). Analytical reproducibility was monitored by co-analyzing homogenized sediment material with known biomarker abundance (every 14–16 sediment samples extracted, analytical error 7%, n = 6). In surface sediment samples from the Okhotsk Sea, the occurrence of IP₂₅ and HBI III reliably reflect modern seasonal sea ice extent (Fig. 1, and Supplementary Fig. 2). Similar concentrations of HBI III in Holocene and MIS 5e indicate that the preservation is good throughout the study period (Fig. 3C).

## 3.5 TEX<sup>L</sup><sub>86</sub> summer sea surface temperature proxy.

Sediment material (1–10 g) was freeze-dried and homogenized by mortar and pestle. The sediments were extracted by Dionex accelerated solvent extraction (DIONEX ASE 200) using a mixture of dichloromethane (DCM)/methanol (MeOH) (9:1, v/v) at a temperature of 100°C and a pressure of 7.6 x 10⁶ Pa. The extracts were separated by Al₂O₃ column chromatography using hexane/DCM (9:1, v/v), hexane/DCM (1:1, v/v) and DCM/MeOH (1:1, v/v) as subsequent eluents. The polar fraction (DCM/MeOH) was dried under N₂, dissolved in hexane/isopropanol (99:1,
287 v/v), and filtered using a 0.4 µm Polytetrafluoroethylene filter prior to injection as 288 described by Hopmans et al. (2016).

289
290 Glycerol dialkyl glycerol tetraethers (GDGTs) were analyzed using an Agilent 1260
291 ultra high performance liquid chromatography (UHPLC) coupled to a 6130 quadrupole
292 mass selective detector in selected ion monitoring mode. Separation was achieved on
293 two UHPLC silica columns (unbonded ethylene bridged hybrid hydrophilic interaction
294 chromatography columns, 2.1 x 150 mm, 1.7 µm; waters) in series, fitted with a 2.1 x
295 5 mm pre-column of the same material (waters) and maintained at 30°C. GDGTs were
296 eluted isocratically using a two solvent gradient system (solvent A=hexane; solvent B=
297 hexane:isopropanol (9:1, v/v)). Initial conditions employed 18% solvent B for 25-min
298 followed by a linear gradient to 35% B over the next 25-min, then a linear gradient to
299 100% B over 30-min. Flow rate was 0.2 ml min⁻¹, resulting in a maximum back
300 pressure of 230 bar for this chromatographic system. Total run time was 90-min with a
301 20-min re-equilibration. Source settings were identical to Schouten et al. (2007).
302 Detection was achieved using atmospheric pressure positive ion chemical ionization
303 mass spectrometry analysis of the eluent. Conditions were: nebulizer pressure 60 psi,
304 vaporizer temperature 400°C, drying gas (N₂), flow 6 L min⁻¹, temperature 200°C,
305 capillary voltage -3 kV, corona 5 µA (~3.2 kV). GDGTs were detected via single ion
306 monitoring of their [M + H]⁺ ions (m/z = 1022, 1036, 1050, 1292, 1296, 1298, 1300,
307 1302) and quantified by integration of the peak areas.
308
309
310 The TEXL₈₆ proxy was proposed as a proxy to estimate SSST in the central
311 Okhotsk Sea (Kim et al., 2010) and equations are listed as below:
312
313
314
315

309
310 The TEXL₈₆ proxy was proposed as a proxy to estimate SSST in the central
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000

$$
311 \quad \mathrm{TEX}^{L}_{86} = \log\left(\frac{[\mathrm{GDGT-2}]}{[\mathrm{GDGT-1}]+[\mathrm{GDGT-2}]+[\mathrm{GDGT-3}]}\right)
$$

$$
312 \quad \text{SSST-TEX}^{L}_{86} = 46.9 + 67.5 \times \text{TEX}^{L}_{86}
$$

313
314

314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000

Cross spectral analysis was done using AnalySeries (Paillard et al., 1996) and its
implemented Blackman-Tukey method using a Bartlett window. Cyclic variations in
the IP<sub>25</sub> and SSST data were identified (Supplementary Table 2; Fig. 5; Taner, 1992;
Meyers, 2014). IP<sub>25</sub> lags precession by ca. 6-kyr (and lead Northern Hemisphere (NH)
summer (June-July-August) insolation by a similar amount, Supplementary Table 2).
Precession filters represent Taner filters using cut-off frequencies of 0.041 and 0.054
and a roll-off rate of $10^{54}$ (Fig. 5; Taner, 1992; Meyers, 2014). Note that the limited length of the time series limits precise phase statements.

## 3.7 Model simulations

To study the time-evolving aspects of orbital and greenhouse gas-driven climate change during glacial cycles, we conducted transient numerical modelling experiments using the earth system model LOVECLIM (Goose et al., 2010). LOVECLIM is a coupled ocean-atmosphere-sea ice-vegetation model. The ocean-sea ice component of LOVECLIM consists of a free-surface Ocean General Circulation Model with a 3° x 3° horizontal resolution coupled to a dynamic-thermodynamic sea-ice model (Fichefet and Morales Maqueda, 1997). The ocean and the sea-ice model use identical horizontal grids. The sea-ice model allows for the presence of open water areas (leads and polynyas) within an ice-covered ocean grid box. The fraction of a grid box covered by ice (sea-ice index) is computed from the heat budget of the open water area in this grid cell (Fichefet and Morales Maqueda, 1997). The sea-ice index is typically given as a number between zero and one with zero referring to ice-free conditions and one to a fully ice-covered grid cell.

The transient simulations were forced by time-dependent boundary conditions for orbital parameters, atmospheric greenhouse gas concentrations, NH ice sheet-orography and albedo following the methodology described in Timmermann et al. (2014). The time- and latitude-dependent orbital forcing was calculated according to Berger (1978). Atmospheric greenhouse gas concentrations (GHG) are prescribed according to reconstructions from EPICA Dome C for CO₂ (Lüthi et al., 2008) as well as CH₄ and N₂O (EPICA community members, 2004). Orbital forcing and atmospheric GHG concentrations are updated every model year. The ice sheets in the NH are prescribed according to Ganopolski and Calov (2011). The forcing is applied with an acceleration factor of 5 which compresses 784,000 years into 156,000 model years. This acceleration factor is appropriate for quickly equilibrating surface variables. The model simulation presented here is an updated version of the one presented in Timmermann et al. (2014) and uses a higher climate sensitivity (~4 K per CO₂-doubling). As a result, glacial-interglacial surface temperature amplitudes are simulated more realistically (Timmermann and Friedrich, 2016).
In addition to the full-forcing simulation described above, four sensitivity simulations were designed to elucidate the individual contributions by GHGs, NH ice sheets and orbital parameters to glacial-interglacial climate change. The first sensitivity simulation used transient forcing as described above but constant preindustrial atmospheric GHG concentrations. The “GHG effect” was then calculated as the difference between the simulation using the full forcing and this simulation. The second sensitivity simulation used transient forcing as described above but constant preindustrial NH ice sheets (extent and albedo). The “NH ice sheet effect” was calculated as the difference between the full-forcing simulation and this simulation.

The two remaining simulations were designed to study the role of orbital forcing under warm and cold climate respectively. For both simulations, transient orbital parameters are used. However, one simulation was run under constant pre-industrial atmospheric CO₂ concentration of 280 ppm whereas the second simulation uses a constant atmospheric CO₂ concentration of 200 ppm. It should be noted, however, that the forcing contributions cannot be deconvoluted entirely. Since variations in atmospheric GHG concentrations and NH ice sheets are actually feedbacks of the earth system, their temporal evolutions are ultimately driven by changes in orbital forcing.

# 4. Results and Discussion

## 4.1 Subarctic Pacific summer sea surface temperature and sea ice variations

Organic geochemical proxy analysis shows dramatic G/IG cycles in both sea ice extent and SSST (Fig. 3A, E). The absence of $IP_{25}$ during MIS 5e (130-118 ka, Fig. 3B) and the majority of the Holocene (10-3 ka, Fig. 4), coupled with highest $TEX^{L}_{86}$-derived SSST (10.7-13.6 °C, Fig. 3E) and open water productivity indicators (HBI III and biogenic opal; Fig. 3C, D), provide clear indication of ice-free conditions during these two interglacials. $TEX^{L}_{86}$-based SSST values for the late Holocene section of MD01-2414 (12.8-13.7 °C) are similar to that of the modern temperature of 13°C (Fig. 4) (Luchin et al., 2009; Seki et al., 2014), while the absence of $IP_{25}$ reflects that the study site was located beyond the position of recent average sea ice extent (Fig. 1; Supplementary Fig. 2). The modelled November sea ice index is also lowest during late MIS 5e and during the Holocene, while SSST data extracted from the simulations also reflect the trend of the $TEX^{L}_{86}$-derived SSST during the past 130 ka, albeit with a slight shift in absolute values (Supplementary Fig. 3). In contrast, the variable occurrence of
IP<sub>25</sub> and substantially lower SSST (3–4 °C) and concentrations of open water productivity indicators during 117-30 ka indicates that seasonal (at least) sea ice cover prevailed during most of the MIS 5d-5a and last glacial periods (Fig. 3A, E), supported further by the observation of maxima in the modelled November sea ice index (Fig. 3B).

394

## 4.2 Precession-cycle control of Okhotsk sea ice change

When examined in more detail, the proxy and model-based sea ice extent in the central Okhotsk Sea during 117-30 ka shows a cycle of sea ice expansion and retreat that follows generally the same trend as local autumn insolation (53°N September–November, SON). In contrast, no satisfactory trend was identified between proxy/model data and other seasonal insolation (i.e. spring/summer/autumn). The general occurrence of IP<sub>25</sub> throughout this interval is, however, punctuated by four intervals where IP<sub>25</sub> is absent. For three of these (i.e. MIS 5c (~97 ka), 5a (~75 ka) and early MIS 3 (~55 ka)), HBI III is also relatively low (but quantifiable) and coincides with low-intermediate values of biogenic opal and SSST (Fig. 3). We interpret these findings as indicative of extremely low or absent sea ice at these times, with generally lower SST and productivity, at least in comparison with the relatively warm and productive MIS 5e and Holocene. In support of this conclusion, relatively low (November) sea ice extent is also observed in the modelled data (Figs. 3B, 4B). Although the boundary conditions associated with IP<sub>25</sub> presence are not currently fully understood, the production of this biomarker by certain sympagic diatoms during the spring bloom (Brown et al., 2011) likely requires that sea ice cover extends at least beyond the winter months in order for this biomarker to provide a positive signature in underlying sediments. Further, IP<sub>25</sub> concentration in surface sediments from other regions is generally positively related to the extent of the overlying spring sea ice cover, even though a lower limit threshold has not, as yet, been identified. Indeed, the majority of the previously reported IP<sub>25</sub> data show lowest concentrations (bordering on limits of quantification in some cases) for subarctic locations characterized by low (e.g. <20%) spring sea ice extent, including the Okhotsk Sea (Supplementary Fig. 2), the neighboring Bering Sea (Méheust et al., 2013) and the Barents Sea in the eastern subarctic (Köseoğlu et al., 2017). As such, it is feasible that some regions experiencing extremely low spring sea ice concentration (e.g. <10%), or sea ice of relatively short seasonality, may not be amenable for study using the IP<sub>25</sub> method. Further, not all sea
ice diatoms produce IP<sub>25</sub>, even if the known producers are considered widespread (Brown et al., 2014). This potentially explains, in part, the absence of IP<sub>25</sub> in the Holocene record of MD01-2414 despite previous reports of some sea ice diatoms during this interval from the same core (Wang and Wang, 2008).

An alternative interpretation of the biomarker data for the three intervals MIS 5c, 5a and early MIS 3 is that permanent sea ice cover prevailed at these times. However, we consider this suggestion unlikely for the following reasons. First, we consider this end-member interpretation for absent IP<sub>25</sub> to be less robust than that of ice-free conditions, especially since IP<sub>25</sub> has been identified in some regions of the central Arctic Ocean experiencing near permanent sea ice cover (Xiao et al., 2015). As such, there are a number of reasons why the interpretation of IP<sub>25</sub> absence for sea ice reconstruction purposes is more challenging than that of its presence. Second, improvements to understanding the IP<sub>25</sub> absence scenario can generally be made through parallel measurement of proxies characteristic of open water conditions, including some phytoplankton biomarkers (e.g. Müller et al., 2009, Belt and Müller, 2013; Belt et al., 2015; Hoff et al., 2016; Méheust et al., 2016). Here, non-zero concentrations of phytoplankton-derived HBI III (and moderate biogenic opal content) are observed during MIS 5c (~97 ka), 5a (~75 ka) and early MIS 3 (~55 ka) (Supplementary Fig. 4), which is difficult to rationalize in terms of permanent ice conditions. In fact, despite their relatively low values, the individual HBI III concentrations during these intervals are all higher than the mean value for the interval 117-30 ka (ca. 0.05 ng g<sup>-1</sup>), which is characterized mainly by variable seasonal sea ice cover (Supplementary Fig. 4). Such an observation is, therefore, more consistent with intervals of low/absent sea ice during MIS 5c, 5a and 3, and with generally lower SSST and productivity compared to MIS 5e and the Holocene (vide supra), each of which are bracketed by longer intervals of more extensive sea ice extent. Third, we note that although HBI III was present in all surface sediments from across the Okhotsk Sea (Supplementary Fig. 2), its concentration was highly dependent on location (Supplementary Table 3). However, the HBI III concentrations in two surface sediments from sites in the central Okhotsk Sea (i.e. near MD01-2414) that experience low sea ice cover in modern times, and where IP<sub>25</sub> was also absent (i.e. stations 55-19-2 and 55-24-2; Supplementary Fig. 2 and Supplementary Table 3), are similar (ca. 0.2-0.3 ng g<sup>-1</sup>) to those found for intervals MIS 5c (~97 ka), 5a (~75 ka) and early MIS 3 (~55 ka) in MD01-2414. These core-top
data thus provide some potential context to those found for during MIS 5c, 5a and 3 in MD01-2414.

On the other hand, absent IP<sub>25</sub> and low HBI III during certain intervals in MIS 2-4 may, instead, reflect increased biomarker degradation downcore. While this may be a contributing factor, the general occurrence of both IP<sub>25</sub> and HBI throughout MIS 2-4 and the subsequent large increase in HBI III concentration during MIS 5e to values similar to those seen in the late Holocene, suggests that any climatic influence likely exceeds that resulting from diagenesis (Stein et al., 2016). Higher concentrations of IP<sub>25</sub> and HBI III in older sedimentary sequences (relative to those found in younger sections of the same cores) have also been reported in previous studies (e.g. Belt et al., 2015), where climatic influences were considered to be the controlling factors of biomarker distributions. As such, we conclude that the three intervals of absent IP<sub>25</sub> during MIS 5c (~97 ka), 5a (~75 ka) and early MIS 3 (~55 ka) represent periods of low/absent sea ice (as also seen in the modelled data), but otherwise relatively cool sea surface conditions and low productivity. These intervals coincide with precession minima and maximum insolation during autumn (Fig. 3B), which contrasts intervals of IP<sub>25</sub> occurrence when insolation was lower.

In contrast, we suggest that the absence of sea ice (IP<sub>25</sub>) and open water (HBI III) biomarkers at ~30 ka is more consistent with the presence of perennial (or near-perennial) sea ice cover at this time, based on related findings and interpretations in previous studies from other regions (e.g. Müller et al., 2009; Hoff et al., 2016). Although our observation is based on only one data point, and the related SSST and biogenic opal data were not obtained from exactly the same sediment horizons (albeit both were relatively low), our interpretation is consistent with earlier records of lithological, magnetic, and floral assemblages from the Okhotsk Sea (Sakamoto et al., 2005; Khim et al., 2012; Nürnberg et al., 2011) at ~30 ka and from the neighboring Bering Sea (Max et al., 2014) based on IP<sub>25</sub> and other geochemical data. Spatial and time differences within the Okhotsk Sea could be attributable to the different sources of IRD sediment source dynamics (Nürnberg et al., 2011) and further combined organic proxy, sedimentological, and inorganic geochemical analyses need to be conducted in different areas to evaluate the source dynamics during the past G/IG cycles in the Okhotsk Sea.
491

492

493

494

495

496

497

498

499

500

501

502

503

504

505

506

507

508

509

510

511

512

513

514

515

516

517

518

519

520

521

522

523

The variability in IP<sub>25</sub> from 120 to 30 ka, therefore, likely reflects reversible transitions between very low (IP<sub>25</sub> absent), seasonal (IP<sub>25</sub> present and variable), and possibly near-perennial (IP<sub>25</sub> absent) sea ice cover, with a cyclicity trend identified in the 20-kyr precession forcing (Fig. 5; Supplementary Table 2). These transitions are supported by previous and related observations in the northern Fram Strait during the last 22 ka (Müller et al., 2009) and Heinrich events in the Nordic Seas over the last 90 ka (Holf et al., 2016).

Overall, the apparent local autumn insolation control on the sea ice extent in the Okhotsk Sea suggests sensitivity to external (insolation) forcing. Further, since current sea ice formation in the Okhotsk Sea generally occurs somewhat later than autumn (typically in December), a further outcome of these findings is that any positive identification of seasonal sea ice in the proxy (i.e. IP<sub>25</sub>) and modeled data likely implies intervals of earlier freeze-up, at least compared to modern times. Consistent with this, although IP<sub>25</sub> is produced by certain diatoms during the spring, its sedimentary occurrence normally reflects regions of autumn sea ice presence (see Belt and Müller, 2013 for a review). Precession filters also align well between IP<sub>25</sub> data and the model-derived November sea ice index (Fig. 5). However, the limited length and resolution of the time series make quantification of (cross-)spectral analyses challenging, and results rely mainly on the clear cyclic expression of proxy data as observed, filtered, and compared to autumn insolation (Figs. 5, 6, 8).

Freshwater input from the Amur River has also been proposed as a potential controlling factor of sea ice formation in the Okhotsk Sea in modern and last G/IG cycles (Sakamoto et al., 2005; Harada et al., 2008). Further, it has been suggested that there may have been a strong coupling between the East Asian summer monsoon (EASM) and the amount of freshwater input on millennial-orbital timescales. However, there is no clear relationship between the composite speleothem δ<sup>18</sup>O record reported by Cheng et al. (2016) and our IP<sub>25</sub> record from the central Okhotsk Sea during the past 130 ka (Fig. 3A, G), even though both records exhibit a strong precession cycle. There is also no clear link between periods of sea ice expansion (high IP<sub>25</sub>) in the central Okhotsk Sea (this study) and intervals of northward shifting of East Asian
monsoons/Intertropical Convergence Zone (ITCZ) precipitation during the past 130 ka
(Cheng et al., 2016).

During the last termination, however, there is a near synchronous (within the error
of current age model) increase in IP₂₅ with a period of stronger EASM during the B/A
period (~15 ka). SSST and biogenic opal increased slightly earlier and HBI III
concentration increased coeval with IP₂₅. These proxy data suggest rapid fluctuations
between open water and seasonal sea ice covered conditions, with high surface
productivity during the B/A period, similar to previous studies in the subarctic Pacific
Ocean (Max et al., 2012). The relationship between EASM and sea ice variations on
millennial-orbital timescales, however, still needs further observational and physical
simulation studies.

### 4.3 CO₂ control for precessionally-paced sea ice

The observation that sea ice free conditions existed throughout MIS 5e, even when
autumn insolation reached a minimum, together with near-perennial sea ice cover
despite only relatively low insolation ~30 ka, suggests that forcing(s) other than
insolation controlled sea ice dynamics during these intervals. Coincidence of ice-free
conditions with atmospheric CO₂ levels in excess of 260 ppm (Figs. 3F, 4F) during both
MIS 5e and the Holocene suggests that such concentrations potentially represent an
important threshold that influences sea ice retreat in the central Okhotsk Sea. Similarly,
maximum sea ice extent in both proxy and simulated records are only attained during
periods of lowest CO₂ (180–200 ppm) in combination with relatively low insolation
(~30 ka, Figs. 3, 5, 6). In contrast, in the low-to-medium CO₂ range (190–260 ppm),
which is present in MIS 3, 5a and 5c, ice-free conditions are paced by precession and
controlled by maximum autumn insolation (Fig. 6). Our proxy-model results in the
central Okhotsk Sea therefore suggest that the major controlling factors of seasonal sea
ice extent are from both external orbital-driven insolation and internal CO₂
concentration forcings on orbital timescales (Fig. 8). On the other hand, the absence of
any clear relationship between sea ice and atmospheric CH₄ (not shown) suggests that
methane was not a significant factor during these interglacials, at least. This is
potentially due to the fact that the total radiative energy variations from atmospheric
CH₄ is believed to be less than 5% of the total radiative energy of greenhouse gases in
recent G/IG cycles (Lo et al., 2017).
4.4 Orbital pacing and CO₂ threshold mechanisms revealed by transient model simulations

The model sensitivity runs deciphering the contributions of different climatic forcings (orbital, CO₂, and NH ice sheet) show that sea ice variations can be regarded mainly as a superposition of the CO₂ effect and orbital forcing, with a small contribution from Northern Hemispheric ice sheets (Fig. 6B). During the sea ice-free MIS 5e, (130–117 ka), autumn insolation varied by nearly 70 W m⁻², yet the November sea ice index varied by less than 0.1 unit, suggesting that orbital forcing was muted during this period of high atmospheric CO₂ concentration (Figs. 6B, 8).

For the sensitivity simulation for orbital forcing under a warm climate (red squares in Fig. 7), a linear regression between insolation and the sea ice index yields a slope of -0.02/ W m⁻², pointing to a relatively weak impact of insolation changes under high atmospheric CO₂ conditions (CO₂ = 280 ppm), even during periods of large insolation changes, such as during MIS 5e (Figs. 6B, 8). However, the impact of orbitally-driven insolation changes on sea ice becomes larger for the lower atmospheric CO₂ condition run (CO₂ = 200 ppm, -0.03/ W m⁻², triangles in Fig. 7).

Under medium CO₂ concentrations during MIS 5a-5d and the last glacial period, insolation forcing and the CO₂ forcing are both important. Thus, precession-paced insolation maxima result in brief episodes of sea ice minima (Fig. 6B), although when the amplitude of the insolation forcing is weak, CO₂ variations appear to be the main driver of the simulated sea ice index during MIS 3 and 4. The maximum in simulated sea ice coverage coincides with minimum CO₂ values of ~190 ppm at ~30 ka. Local SON insolation increased towards the last glacial maximum, resulting in slightly decreased sea ice extent when compared to ~30 ka.

5. Conclusions

In this study, we reconstructed SSST and sea ice variations in the central Okhotsk Sea over the last 130 ka by combination of novel organic geochemical proxies including TEXL₈₆ and IP₂₅. To reveal the physical mechanisms responsible for the near-surface oceanographic changes, numerical simulation data and sensitivity runs were also performed. Our geochemical proxy analyses and transient simulation data show that the
precession pacing by local autumn insolation was a major control over variation in sea ice extent during the past 130 ka. Greenhouse gas (mainly CO₂) induced atmospheric radiative forcing acts as a further threshold for sea ice absence when CO₂ concentration exceeded ~260 ppm and for extensive sea ice cover when CO₂ concentration fell below 190 ppm. The dominant driver of, respectively, sea ice free conditions during MIS 5e and perennial sea ice conditions ~30 ka is thus atmospheric radiative forcing. Therefore, we suggest a combined orbital and greenhouse gas control over sea ice variations in the central Okhotsk Sea during the past 130 ka.

Near-surface oceanographic conditions (i.e. sea ice and SSST) may also have been influenced by the extent of fresh water input from the Amur River. However, further regional studies including seasonal temperature and more reliable freshwater records before the role of this potential forcing can be fully understood.

## Acknowledgements

This project was supported by grants of the Ministry of Science and Technology (MOST), Taiwan ROC (10-2119-M-002-003, 104-2917-I-564-046, 105-2119-M-002-001), State Key Laboratory of Isotope Geochemistry, Guangzhou Institute of Geochemistry, Chinese Academy of Sciences start-up fund (SKLaBIG-QD-16-04) and National Natural Science Foundation (41773116). MD01-2414 sediment samples were generously provided by Taiwan Ocean Research Institute (TORI). We thank the University of Plymouth for financial support. S.S. and J.L. were supported by European Research Council (ERC) under the European Union's Seventh Framework Program (FP7/2007-2013) ERC grant agreement [339206] and by the Netherlands Earth System Centre (NESSC) funded by the Dutch Ministry of OCW. Finally, we thank two anonymous reviewers and the Editor for providing extremely positive feedback on the original version of this manuscript, with further suggestions that improved its clarity in final form.

## References

Bassinot, F., Chen, M.-T., 2002. Physical properties of sediment core MD01-2414. doi:10.1594/PANGAEA.80047.

Belt, S.T., Massé, G., Rowland, S.J., Poulin, M., Michel, C., LeBlanc, B., 2007. A novel chemical fossil of palaeo sea ice: $IP_{25}$. Org. Geochem. 38, 16-27.

Belt, S.T., Cabedo-Sanz, P., Smik, L., Navarro-Rodriguez, A., Berben, S.M.P., Knies, J., Husum, K., 2015. Identification of paleo Arctic winter sea ice limits and the
628 marginal ice zone: optimised biomarker-based reconstructions of late Quaternary
629 Arctic sea ice. Earth Planet. Sci. Lett. 431, 127-139.

630 Belt, S.T., Müller, J., 2013. The Arctic sea ice biomarker IP25: a review of current
631 understanding, recommendations for future research and applications in palaeo sea
632 ice reconstructions. Quat. Sci. Rev. 79, 9-25.

633 Berger, A., 1978. Long-term variations of daily insolation and quaternary climate
634 change. J. Atmos. Sci. 35, 2362-2367.

635 Brown, T.A., Belt, S.T., Tatarek, A., Mundy, C.J., 2014. Source identification of the
636 Arctic sea ice proxy IP25. Nat. Commun. 5, 4197.

637 Cheng, H., Edwards, R.L., Sinha, A., Spötl, C., Yi, L., Chen, S., Kelly, M., Kathayat,
638 G., Wang, X., Li, X., Kong, X., Wang, Y., Ning, Y., Zhang, H., 2016. The Asian
639 monsoon over the past 640,000 years and ice age terminations. Science 534, 640-
640 646.

641 Chou, Y.-M., Lee, T.-Q., Song, S.-R., Chen, K.-J., 2011. Magnetostratigraphy of marine
642 sediment core MD01-2414 from Okhotsk Sea and its paleoenvironmental
643 implications. Mar. Geol. 284, 149-157.

644 Detlef, H., Belt, S.T., Sosdian, S.M., Smik, L., Lear, C.H., Hall, I.R., Cabedo-Sanz, P.,
645 Husum, K., Kender, S., 2018. Changes in sea ice dynamics in the eastern Bering
646 Sea during the Mid-Pleistocene climate transition. Nat. Commun. (in press).

647 EPICA community members, 2004. Eight glacial cycles from an Antarctic ice core.
648 Nature 429, 623-628.

649 Ganopolski, A., Calov, R., 2011. The role of orbital forcing, carbon dioxide and regolith
650 in 100 kyr glacial cycles. Clim. Past 7, 1415-1425.

651 Goosse, H., Brovkin, V., Fichefet, T., Haarsma, R., Huybrechts, P., Jongma, J.,
652 Mouchet, A., Selten, F., Barriat, P.-Y., Campin, J.-M., Deleersnijder, E.,
653 Driesschaert, E., Goelzer, H., Janssens, I., Loutre, M.-F., Morales Maqueda, M.A.,
654 Opsteegh, T., Mathieu, P.-P., Munhoven, G., Pettersson, E.J., Renssen, H., Roche,
655 D.M., Schaeffer, M., Tartinville, B., Timmermann, A., Weber, S.L., 2010.
656 Description of the Earth system model of intermediate complexity LOVECLIM
657 version 1.2. Geosci. Model Dev 3, 603-633.

658 Gorbarenko, S.A., Harada, N., Malakhov, M.I., Vasilenko, Y.P., Bosin, A.A., Goldberg,
659 E.L., 2010. Orbital and millennial-scale environmental and sedimentological
660 changes in the Okhotsk Sea during the last 350 kyr. Global Planet. Change 72, 79-
661 85.

662 Harada, N., Sato, M., Sakamoto, T., 2008. Freshwater impacts recorded in
663 tetraunsaturated alkenones and alkenone sea surface temperatures from the
664 Okhotsk Sea across millennial-scale cycles. Paleoceanography 23, PA3201,
665 doi:10.1029/2006PA001410.

666 Hoff, U., Rasmussen, T.L., Stein, R., Ezat, M.M., Fahl, K., 2016. Sea ice and millennial-
667 scale climate variability in the Nordic seas 90 kyr ago to present. Nat. Commun. 7,
668 12247, doi:10.1038/ncomms12247.

669 Hopmans, E.C., Schouten, S., Sinninghe Damsté, J.-S., 2016. The improved
670 chromatography on GDGT-based paleoproxies. Org. Geochem. 93, 1-6.

671 Huang, J.-J., Löwemark, L., Chang, Q., Lin, T.-Y., Chen, H.-F., Song, S.-R., Wei, K.-Y.,
672 2016. Choosing optimal exposure times for XRF core-scanning: Suggestions based
673 on the analysis of geological reference materials. Geochem. Geophy. Geosy. 17,
674 1558-1566.

675 Imbrie, J., Imbrie, J.Z., 1980. Modeling the climatic response to orbital variations.
676 Science 207, 943-953.

677 Kashiwase, H., Ohshima, K.I., Nihashi, S., 2014. Long-term variation in sea ice
678 production and its relation to the intermediate water in the Sea of Okhotsk. Prog.

679 Oceanogr. 126, 21-32.

680 Khim, B.-K., Sakamoto, T., Harada, N., 2012 Reconstruction of surface water

681 conditions in the central region of the Okhotsk Sea during the last 180 kyrs. Deep

682 Sea Res. II 61-64, 63-72.

683 Kim, J.-H., van der Meer, J., Schouten, S., Helmke, P., Willmott, V., Sangiorgi, F., Koç,

684 N., Hopmans, E.C., Sinninghe-Damsté, J.S., 2010. New indices and calibrations

685 derived from the distribution of crenarchaeal isoprenoid tetraether lipids:

686 Implications for past sea surface temperature reconstructions. Geochim.

687 Cosmochim. Ac. 74, 4639-4654.

688 Kimura, N., Wakatsuchi, M., 2004. Increase and decrease of sea ice area in the Sea of

689 Okhotsk: Ice production in coastal polynyas and dynamic thickening in

690 convergence zones. J. Geophys. Res. 109, C9, doi:10.1029/2003JC001901

691 Köseoğlu, D., Belt, S.T., Smik, L., Yao, H., Panieri, G., Knies, J., 2018.

692 Complementary biomarker-based methods for characterising Arctic sea ice

693 conditions: A case study comparison between multivariate analysis and the PIP<sub>25</sub>

694 index. Geochim. Cosmochim. Ac. 222, 406-420.

695 doi.org/10.1016/j.gca.2017.11.001.

696 Laskar, J., Robutel, P., Joutel, F., Gastineau, M., Correia, A.C.M., Levard, B., 2004. A

697 long-term numerical solution for the insolation quantities of the Earth. Astron.

698 Astrophys. 428, 261-285.

699 Lisiecki, L.E., Raymo, M.E., 2005. A Pliocene-Pleistocene stack of 57 globally

700 distributed benthic $\delta^{18}$O records. Paleooceanography 20, PA1003.

701 Liu, Y.-J., Song, S.-R., Lee, T.-Q., Lee, M.-Y., Chen, Y.-L., Chen, H.-F., 2006.

702 Mineralogical and geochemical changes in the sediments of the Okhotsk Sea

703 during deglacial periods in the past 500 kyrs. Global Planet. Change 53, 47-57.

704 Lo, L., Chang, S.-P., Wei, K.-Y., Lee, S.-Y., Ou, T.-H., Chen, Y.-C., Chuang, C.-K., Mii,

705 H.-S., Burr, G.S., Chen, M.-T., Tung, Y.-H., Tsai, M.-C., Hodell, D.A., Shen, C.-

706 C., 2017. Nonlinear climatic sensitivity to greenhouse gases over the past 4

707 glacial/interglacial cycles. Sci. Rep. 7, 4626, doi: 10.1038/s41598-017-04031-x.

708 Luchin, V., Kruts, A., Sokolov, V., Rostov, V., Rudykh, N., Perunova, T., Zolotukhin,

709 E., Pischalnik, V., Romeiko, L., Hramushin, V., Shustin, V., Udens, Y., Baranova,

710 O., Smolyar, I., Tarosh, E., 2009. Climatic Atlas of the North Pacific Seas: Bering

711 Sea, Sea of Okhotsk, and Sea of Japan, NOAA Atlas NESDIS 67, U.S. Gov.

712 Printing Office, Wash., D.C, p. 380 (CD Disc).

713 Lüthi, D., Floch, M.L., Bereiter, B., Blunier, T., Barnola, J.-M., Siegenthaler, U.,

714 Raynaud, D., Jouzel, J., Fischer, H., Kawamura, K., Stocker, T.F., 2008. High-

715 resolution carbon dioxide concentration record 650,000-800,000 years before

716 present. Nature 453, 379-382.

717 Max, L., Riethdorf, J.-R., Tiedemann, R., Smirnova, M., Lembke-Jene, L., Fahl, K.,

718 Nürnberg, D., Matul, A., Mollenhauer, G., 2012. Sea surface temperature

719 variability and sea-ice extent in the subarctic northwest Pacific during the past

720 15,000 years. Paleooceanography 27, PA3213, doi:10.1029/2012PA002292.

721 Max, L., Tiedemann, R., Fahl, K., Nürnberg, D., Riethdorf, J.-R., 2014. Rapid shifts in

722 subarctic Pacific climate between 138 and 70 ka. Geology 42, 899-902.

723 Méheust, M., Fahl, K., Stein, R., 2013. Variability in modern sea surface temperature,

724 sea ice and terrigenous input in the sub-polar North Pacific and Bering Sea:

725 Reconstruction from biomarker data. Org. Geochem. 57, 54-64.

726 Méheust, M., Stein, R., Fahl, K., Max, L., Riethdorf, J.-R., 2016. High-resolution IP<sub>25</sub>-

727 based reconstruction of sea-ice variability in the western North Pacific and Bering
728 Sea during the past 18,000 years. Geo-Marine Lett. 36, 101-111.

729 Mesquita, M.D., Hodges, K.I., Atkinson, D.E., Bader, J., 2011. Sea-ice anomalies in the Sea of Okhotsk and the relationship with storm tracks in the Northern Hemisphere during winter. Tellus A 63, 312-323.

732 Meyers, S.R., 2014. Astrochron: An R Package for Astrochronology Version 0.6.5.

733 Müller, J., Massé, G., Stein, R., Belt, S.T., 2009. Variability of sea-ice conditions in the Fram Strait over the past 30,000 years. Nat. Geosci. 2, 772-776.

735 Nishioka, J., Yasuda, I., Liu, H., Nakatsuka, T., Volkov, Y.N., 2014. Biogeochemical and physical processes in the Sea of Okhotsk and the linkage to the Pacific Ocean. Prog. Oceanogr. 126, 1-7.

738 Nürnberg, D., Dethleff, D., Tiedemann, R., Kaiser, A., Gorbarenko, S.A., 2011. Okhotsk Sea ice coverage and Kamchatka glaciation over the last 350 ka- Evidence from ice-rafted debris and planktonic δ¹⁸O. Palaeogr. Palaeoclimatol. Palaeoecol. 310, 191-205.

742 Ono, A., Takahashi, K., Katsuki, K., Okazaki, Y., Sakamoto, T., 2005. The Dansgaard-Oeschger cycles discovered in the upstream source region of the North Pacific Intermediate Water formation. Geophys. Res. Lett. 32, L11607, doi:10.1029/2004GL022260.

746 Paillard, D., Labeyrie, L., Yiou, P., 1996. Macintosh program performs time-series analysis, EOS Trans. 77, 379.

748 Reimer, P. J., Bard, E., Bayliss, A., Beck, J. W., Blackwell, P. G., Bronk Ramsey, C., Buck, C. E., Cheng, H., Edwards, R. L., Friedrich, M., Grootes, P. M., Guilderson, T. P., Haflidason, H., Hajdas, I., Hatté, C., Heaton, T. J., Hoffmann, D. L., Hogg, A. G., Hughen, K. A., Kaiser, K. F., Kromer, B., Manning, S. W., Niu, M., Reimer, R. W., Richards, D. A., Scott, E. M., Southon, J. R., Staff, R. A., Turney, C. S. M., and van der Plicht, J., 2013. IntCal13 and Marine13 radiocarbon age calibration curves 0-50,000 years cal BP. Radiocarbon, 55,1869-1887, https://doi.org/10.2458/azu_js_rc.55.16947, 2013.

756 Sakamoto, T., Ikehara, M., Aoki, K., Iijima, K., Kimura, N., Nakatsuka, T., Wakatsuchi, W., 2005. Ice-rated debris (IRD)-based sea-ice expansion events during the past 100 kyrs in the Okhotsk Sea. Deep-Sea Res. II 52, 2275-2301.

759 Schlitzer, R., 2017. Ocean Data View, odv.awi.de.

760 Schouten, S., Huguet, C., Hopmans, E.C., Kienhuis, M., Sinninghe Damsté, J.-S., 2007. Analytical Methodology for TEX₈₆ paleothermometry by high-performance liquid chromatography/atmospheric pressure chemical ionization-mass spectrometry. Anal. Chem. 79, 2940-2944.

764 Seki, O., Bendle, J.A., Harada, N., Kobayashi, M., Sawada, K., Moossen, H., Inglis, G.N., Nagao, S., Sakamoto, T., 2014. Assessment and calibration of TEX₈₆ paleothermometry in the Sea of Okhotsk and sub-polar North Pacific region: Implications for paleoceanography. Prog. Oceanogr. 126, 254-266.

768 Serreze, M.C., Stroeve, J., Barrett, A.P., Boisvert, L.N., 2016. Summer atmospheric circulation anomalies over the Arctic Ocean and their influence on September sea ice extent. A cautionary tale. J. Geophys. Res. Atmos. 121, 1463-11485.

771 Stein, R., Fahl, K., Schreck, M., Knorr, G., Forwick, M., Gebhardt, C., Jensen, L., Kaminski, M., Kopf, A., Matthiessen, J., Jokat, W., Lohmann, G., 2016. Evidence for ice-free summers in the late Miocene central Arctic Ocean. Nat. Commun. 7, 11148, doi:10.1038/ncomms11148.

775 Strakhov, N.M., Bezrukov, P.L., Yablokov, V.S., 1961. Recent precipitates in Seas and Oceans. Press of the Academy of Sciences USSR, Moscow, 644 pp. (in Russian).

777 Taner, M.T., 1992. In: Attributes revisted (Technical Report, Rock Solid Images, Inc),
778 http://www.rocksolidimages.com/attributes-revisited/#_Toc328470897.

779 Timmermann, A., Friedrich, T., Timm, O.E., Chikamoto, M.O., Abe-Ouchi, A., Ganopolski, A., 2014. Modeling obliquity and CO₂ effects on Southern Hemisphere climate during the past 408 ka. J. Clim. 27, 1863-1875.

780 Timmermann, A., Friedrich, T., 2016. Late Pleistocene climate drivers of early human migration. Nature 538, 92-95.

781 Turner, J., Scott Hosking, J., Marshall, G.J., Phillips, T., Bracegirdle, T.J., 2016. Antarctic sea ice increase consistent with intrinsic variability of the Amundsen Sea Low. Clim. Dyn. 46, 2391-2402.

782 Veres, D., Bazin, L., Landais, A., Toyé Mahamadou Kele, H., Lemieux-Dudon, B., Parrenin, F., Martinerie, P., Blayo, E., Blunier, T., Capron, E., Chappellaz, Rasmussen, S.O., Severi, M., Svensson, A., Vinther, B., Wolff, E.W., 2013. The Antarctic ice core chronology (AICC2012): an optimized multi-parameter and multi-site dating approach for the last 120 thousand years. Clim. Past 9, 1733-1748.

783 Wang, W.-L., Wang, L.-C., 2008. Reconstruction of oceanographic changes based on the diatom records of the central Okhotsk Sea over the last 500000 years. Terr. Atmos. Ocean. Sci. 19, 403-411.

784 Xiao, X., Fahl, K., Müller, J., Stein, R., 2015. Sea-ice distribution in the modern Arctic Ocean: Biomarker records from trans-Arctic Ocean surface sediments. Geochim. Cosmochim. Ac. 155, 16-29.

785
786
787
788
789
790
791
792
793
794
795
796
797
798
799

800 **Fig. 1.** Site location. Core MD01-2414 (53°11.77′N, 149°34.80′E, water depth 1123 m) is located in the Okhotsk Sea. The blue arrows indicate surface water circulation (WKC: West Kamchatka Current, ESC: East Sakhalin Current, and SWC: Soya Warm Current) and the dashed black line denotes the average position of modern seasonal sea ice extent, during the months of November to June. This map was generated with Ocean Data View (GMT) version 5 (Schlitzer, 2017).

806
807 **Fig. 2.** Age model of Site MD01-2414 for the upper 7 m. (A) ITRAX scanned log(Ba/Ti) ratio, (B) Ba, (C) Ti, (D) magnetic susceptibility in log scale (Chou et al., 2011), (E) b* (Bassinot and Chen, 2002), (F) coarse fraction (C.F., this study) in weight percent (this study), and (G) global benthic composite oxygen isotope curve (black curve, Lisiecki and Raymo, 2005) and sedimentation rate (orange). Dark yellow line represents 3-point running average. Green crosses and triangles represent radiocarbon dates and age model tie points, respectively.
808
809
810
811
812
813

814
815 **Fig. 3.** Geochemical proxy results from Site MD01-2414 during the past 130 ka. (E) SSST derived from TEX<sup>L86</sup> in red and 53°N June-July-August (JJA) insolation in pink (Lasakar et al., 2004). (A) IP<sub>25</sub> concentration in gray. (B) model-derived Okhotsk Sea November sea ice index, and 53°N September-October-November (SON) insolation in
orange (Lasakar et al., 2004) (C) HBI III (green) and concentration. (D) MD01-2414 opal content (dark yellow, Liu et al., 2006). (F) Atmospheric CO₂ (Lüthi et al., 2008) in AICC2012 timescale (Veres et al., 2013). (E) SSST derived from TEXL86 in red and 53°N June-July-August (JJA) insolation in pink (Laskar et al., 2004). (G) Compiled Chinese speleothem δ¹⁸O record (Cheng et al., 2016). (H) XRF scanned log(Ba/Ti) in purple and LR04 benthic oxygen isotope curve in black (Lisiecki and Raymo, 2005). Semi-quantitative sea ice conditions from: sea ice retreat (white), sea ice expand (gray), and possible perennial sea ice (blue with question mark). Slash bar in (E) represents critical CO₂ threshold from 260 to 280 ppm. Marine isotope stages (MIS) are listed in Fig. 3H.

**Fig. 4.** Geochemical proxy results from Site MD01-2414 during the past 29 ka. Captions and color bars are the same as Fig. 3. Black triangles represent AMS ¹⁴C dates with 2σ error bars in Supplementary Table 1. YD, B/A, and H1 represent Younger Dryas, Bølling-Allerød (yellow bar), and Heinrich stadial 1 periods, respectively.

**Fig. 5.** Precession filters comparison between model and proxy data. Data (all black) of the modelled November Sea ice index (A), the IP₂₅ (B) and precession (C), together with respective precession filters (all red, Taner filters using cut-off frequencies (1 kyr⁻¹) of 0.041 and 0.054 and a roll-off rate of 10⁵⁴ (Taner, 1992; Meyers, 2014) plotted versus time.

**Fig. 6.** Individual contributions by variations in GHG, orbital parameters and NH ice sheets. (A) 53°N SON insolation (Laskar et al., 2004). (B) Transient model results for Okhotsk Sea November sea ice index anomaly with respect to 1-kyr for full-forcing (black solid), GHG-effect (red dashed), orbital effect (green dashed) and effect of NH ice sheet (purple dashed). Background shading and MIS numbers are the same as Fig. 3.

**Fig. 7.** Scatter plot of transient model results for Okhotsk Sea November sea-ice index versus 53°N SON insolation (W m⁻²) (Laskar et al., 2004) for full-forcing (circles) and orbital effect for cold and warm climates (squares and triangles). Corresponding atmospheric CO₂ concentrations (ppm) are shaded. The orbital-effect simulations use
852 only transient orbital parameters and constant atmospheric CO₂ concentration of either
853 280 ppm (squares) or 200 ppm (triangles) respectively. Please see Section 3.7 (Model
854 simulations) for details regarding the model simulations.

855
856 **Fig. 8.** Schematic figure of the of CO₂ level and insolation forcings on sea ice variability
857 in the central Okhotsk Sea during the past 130 ka.

858
<figure><img src="image_2.png" /><figcaption>Figure 1</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 2</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 3</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 6</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 8</figcaption></figure>

882

883

884