Northwest Paciﬁc ice-rafted debris at 38 ◦ N reveals episodic ice-sheet change in late Quaternary Northeast Siberia

The ice-rafted-debris (IRD) record of the open Northwest Paciﬁc points towards the existence of substantial glacial ice on the Northeast Siberian coast during the late Quaternary. However, the scale and timing of glaciation and de-glaciation remains controversial due to the dearth of both onshore and offshore records. Existing IRD data suggests at least one event of dynamic and abrupt change during mid-late Marine Isotope Stage (MIS 3) which mimics the massive collapse of the Laurentide ice sheet during Heinrich Events. It is uncertain whether other events of this magnitude occurred during the late Quaternary. Here we present a ∼ 160,000 yr IRD series, planktic foraminiferal counts and an age model, derived from a benthic δ 18 O curve, radiocarbon dates and tephrochronology, from core ODP 1207A (37.79 ◦ N, 162.75 ◦ E), revealing the presence of low but episodic ﬂux of IRD. We conclude that glacial Northwest Paciﬁc icebergs spread further south than previously thought, with icebergs emanating from Northeast Siberia being transported to the transition region between the subpolar and subtropical waters, south of the subarctic front during at least the Quaternary’s last two glacial periods. The episodic nature of the 1207A IRD record during the last glacial, combined with coupled climate-iceberg modelling, suggests occasional times of much enhanced ice ﬂux from the Kamchatka-Koryak coast, with other potential sources on the Sea of Okhotsk coast. These ﬁndings support the hypothesis of a variable but extensive ice mass during the last glacial over Northeast Siberia, particularly early in the last glacial period, behaving independently of North American and Eurasian ice masses. In strong contrast, IRD was absent during much of the penultimate glacial Marine Isotope Stage (MIS) 6 suggesting the possibility of very different Northeast Siberian ice coverage between the last two glacial periods.  2020 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
The Northeast (NE) Siberian coast has been the source of significant iceberg flux to the Northwest (NW) Pacific throughout the entire Quaternary (Conolly and Ewing, 1970;Haug et al., 1995;St. John and Krissek, 1999). As the NW Pacific is the end-member of the global ocean circulation (Keigwin, 1998), significant environmental changes associated with large-scale iceberg flux in this sev, 1974;Glushkova, 2001). However, this last glacial history of NE Siberia relies on a limited number of terrestrial dating tie points (Glushkova, 2001;Solomina and Calkin, 2003;Braitseva et al., 2005), and most of the glacial landforms in the region have no absolute or relative chronological control.
A NW Pacific marine core record of NE Siberian glaciation is therefore vital in constraining the timing and scale of significant glacial change. The limited number of key ocean core datasets from the NW Pacific and the marginal seas of the Sea of Okhotsk and the Bering Sea show enhanced IRD flux during the penultimate glacial period (191,,000 yrs ago), around the onset of the last glacial period (∼120,000 yrs ago), and during MIS 3 (∼57,000-29,000 yrs ago) (Nürnberg et al., 2011;St. John and Krissek, 1999;Levitan et al., 2013). Nevertheless, the NW Pacific IRD does not show compelling evidence of synchronicity with either the Atlantic or the NE Pacific IRD records (Maier et al., 2018). This suggests that the last glacial NE Siberian ice masses were not dynamically coupled to the Laurentide and Cordilleran ice sheets, despite the latter showing linkage between freshwater and iceberg injections from the Cordilleran ice sheet and Northeast Pacific IRD and planktic (surface) δ 18 O data during Atlantic Heinrich stadials 4 (∼40-38,000 yrs ago), and 1 (17-15,000 yrs ago; Maier et al., 2018).
The highest concentrations of glacial IRD in the open NW Pacific are situated on the Meiji Drift ( Fig. 1) in a NE-SW band extending from the east coast of Kamchatka (Conolly and Ewing, 1970). However, IRD is present during a range of glacial periods at least as far south as 41 • N, off the coast of northern Japan. The southern known limit of Quaternary IRD marks the approximate position of the subarctic front (SAF) where the cold Oyashio Current and the warm Kuroshio Current converge (Thompson and Shackleton, 1980;Chiyonobu et al., 2012). This convergence zone between the NW Pacific's sub-polar and sub-tropical gyres (Fig. 1), while it will have moved equatorward during cold periods, is likely to be roughly fixed geographically, like the equivalent Gulf Stream detachment zone in the western North Atlantic.
This wide extent of IRD distribution supports the existence of a large glacial ice mass or masses reaching the NE Siberian coast (Kent et al., 1971; last glacial record summarised in Fig. 2a of Bigg et al., 2008). At some sites, observed IRD concentrations imitate the sub-millennial duration and magnitude of the catastrophic iceberg armada events discovered from N. Atlantic Heinrich layers (Hemming, 2004). This is particularly apparent during MIS 3 at ODP 883 (51.2 • N, 167.8 • E; Fig. 1; Bigg et al., 2008) where a sudden increase in IRD occurs ∼40,000 yrs ago. Another site, ODP 882 (50.2 • N, 167.4 • E) off the east coast of Kamchatka records similar increases in the siliciclastic fraction during late MIS 4 to mid MIS 3 (Haug et al., 1995) and it is assumed that due to its relatively close proximity, this increase relates to the same abrupt ice-rafting event observed in ODP 883. Given the high magnitude of this event, with IRD counts in the >150 µm size fraction of >2000 grains/g of bulk

Site selection
A key factor in the uncertainty surrounding NE Siberian glaciation, in addition to the incomplete terrestrial record, is the unsuitability of late Quaternary NW Pacific marine sediments for radiocarbon dating in this carbonate-poor basin (Korff et al., 2016). However, on preliminary inspection ODP 1207A, situated on a bathymetric high of 3112 m depth at (37.79 • N, 162.75 • E), was found to contain abundant carbonate that had not been previously dated, and had not previously been inspected for IRD, possibly because it was considered to be beyond the southern limit of icerafting as delimited by Conolly and Ewing (1970) and Kent et al. (1971).
The upper 300 cm of ODP 1207A was selected for sampling. This section was largely well preserved, with few missing sections apart from a 13 cm gap between 137 and 149 cm. The chronology of the upper 300 cm was uncertain due to not having been previously dated, however the upper 300 cm of other NW Pacific cores in the Shatsky Rise region were typically accumulated over approximately the last 180 ka BP (e.g. Korff et al., 2016) to 200 ka BP (e.g. Yamane, 2003), therefore it was assumed that the upper 300 cm of ODP 1207A would cover a suitable time period for the study of sediments spanning at least the last glacial period. In what follows depths are given relative to the core top, and so is expressed as below sea floor (bsf).

Accelerator mass spectrometry 14 C dating
Around 20 mg of the planktic foraminifera species Neogloboquadrina (N.) pachyderma were picked from the >150 µm size fraction. Samples were cleaned by sonication in deionised water, followed by sonication in methanol. Pre-treatment of the foraminifera specimens followed the approach by Nadeau et al. (2001). Foram samples were weighed into septa seal vials and cleaned by adding 2 ml of 15% hydrogen peroxide and placed in an ultrasonic bath for 3m i n . The peroxide was removed with a syringe, the forams rinsed with Milli-Q ® water and then most of the water removed again with a syringe. While the forams were still damp, the septa seal vial was evacuated. A quantity of 1-2 mil of 85% phosphoric acid was added (depending on sample size) to hydrolyse the samples to CO 2 . The vials were placed on a heating block at 70 • C until the forams were completely dissolved. The CO 2 was converted to graphite on an iron catalyst using the hydrogen reduction method (Vogel et al., 1987) and the 14 C/ 12 C ratio and 13 C/ 12 C were measured by accelerator mass spectrometry (AMS) at the 14 CHRONO Centre, Queen's University Belfast. The sample 14 C/ 12 C ratio was background-corrected and normalised to the HOXII standard (SRM 4990C; National Institute of Standards and Technology). The radiocarbon ages were corrected for isotope fractionation using the AMS measured δ 13 C which accounts for both natural and machine fractionation. The radiocarbon age and one standard deviation were calculated using the Libby half-life of 5568 yrs following the conventions of Stuiver and Polach (1977). The ages were corrected for isotope fractionation using the AMS measured δ 13 C (not given) which accounts for both natural and machine fractionation. The background sample used for samples UBA-28906 and 28907 was Icelandic spar calcite. For UBA-26539 we used the 14 C/ 12 C ratio of forams processed at the same time from a depth expected to correspond to >60 ka BP. The CALIB8.2 Radiocarbon Calibration Program (Stuiver et al., 2020) with the Marine20 curve (Heaton et al., 2020) was used to calibrate the dates, and includes a global offset and a regional correction ( R) of 328 14 C yrs (Yoneda et al., 2007).

Benthic δ 18 O stratigraphy
Samples were chosen from the benthic foraminifera species Cibicidoides wuellerstorfi. Samples were picked from the >150 µm size fraction, of a consistent diameter of ∼500 µm. These were cleaned by sonication in de-ionised water, and then methanol. After various trials, the oxygen isotope composition (δ 18 O) for each sampled depth was taken to be the mean of analyses of three separate individuals, as this yielded a low standard deviation between replicates. The upper 130 cm was sampled at 1 cm intervals, while the deeper core, down to c. 260 cm was sampled every 5 cm, with some additional samples near depths recording big swings in values.
Approximately 60-100 micrograms of carbonate were used for isotope analysis using an IsoPrime dual inlet mass spectrometer plus Multiprep device at the British Geological Survey, Keyworth. Isotope values ( 13 C, 18 O) were reported as per mille ( ) deviations of the isotopic ratios ( 13 C/ 12 C, 18 O/ 16 O) calculated to the VPDB scale using a within-run laboratory standard calibrated against NBS-19. The Calcite-acid fractionation factor applied to the gas values was 1.00798. Due to the long run time of 21 hours a drift correction was applied across the run, calculated using the standards that bracket the samples. The Craig correction was also applied to account for 17 O. The average analytical reproducibility of the standard calcite (KCM) was 0.05 for δ 13 C and δ 18 O.

Volcanic glass counts and tephrochronology
Counts of volcanic glass were carried out in the >150 µm fraction alongside IRD counts, to search for evidence of potentially dateable tephra layers. A distinct tephra layer at 160-169 cm was geochemically fingerprinted with the help of electron microprobe analysis and laser ablation inductively coupled mass-spectrometry (LA-ICP-MS) of single glass shards. Major and selected volatile elements (including S, F and Cl) in glass shards from ODP 1207A (166-167 cm, within the visible tephra layer) were analyzed in 2015 at GEOMAR Helmholtz Centre for Ocean Research Kiel, using a JEOL JXA 8200 electron microprobe. All analyses were performed at 15 kV accelerating voltage with a beam diameter of 5µ m and a beam current of 6n A . A detailed account of analytical conditions can be found in Portnyagin et al. (2020), including internal standards and settings for spectrometers. For comparison with literature data, the analyses were normalized to an anhydrous basis (i.e. 100% total oxides). Trace element analyses were performed using laser ablation inductively coupled mass-spectrometry (LA-ICP-MS) at the Institute of Geosciences, Christian-Albrecht University of Kiel, Germany. The analyses were performed in 2019 using a quadrupole-based (QP) ICP-MS (Agilent 7900) and a Coherent Ge-oLas ArF 193 nm Excimer LA system. Analytical conditions and subsequent data reduction have been described by Portnyagin et al. (2020). The data on glass from ODP 1207A and reference materials used for EMP and LA-ICP-MS analyses are provided in Supplementary Information Table 1 (Appendix A).

IRD counts
Ice-rafted debris counts were carried out down the entire upper 300 cm of core ODP 1207A, at a high resolution of every 1c m . Each 1c m sample, weighing between approximately 2g and 5g was wet-sieved through a 150 µm sieve. The >150 µm size fraction was then closely inspected using a standard reflected light binocular microscope. Consistent with earlier studies, this 150 µm level was deemed to be the lower size threshold for IRD categorisation Hemming, 2004). Due to the expected low concentrations of IRD, the sample was not split before counts took place. Grains with the character of lithic material were categorised as IRD. Using a standard method of quantification (Hemming, 2004), IRD was expressed as grains per gram of the dry unsieved bulk sample. Selected grains were analysed further using a scanning electron microscope to obtain high-resolution images confirming ice-rafted features including conchoidal fractures, scour marks, and angularity (Immonen et al., 2014). Energy Dispersive Spectroscopy (EDS) was also carried out on selected grains to semiquantitatively analyse their major element composition.
In this work we claim that ODP 1207A contains the most southerly NW Pacific record of IRD. Previous work (Taylor et al., 1990;Bigg et al., 2008) had recorded IRD in ODP 792/793 at 32 • N, thus further south than ODP 1207A. However, on closer inspection in this study the "pebbles" in the ODP 792/793 record were found to be pumice, and were associated with ash events. The material could have been ice-rafted, but at this subtropical latitude, there are other more likely pumice entrainment methods, such as flotation. Therefore, this pumice is not categorised as IRD here.

Planktic foraminifera counts
We investigated the relationship between IRD and surface oceanographic shifts through high-resolution counts of planktic foraminifera. Following wet sieving, the >150 µm size fraction was split until a sample of approximately 300 grains was obtained. Counts of all planktic foraminiferal species were then carried out in the split sample.

Iceberg modelling
The dominance of quartz in the IRD record suggests that IRD deposited at site 1207A may have originated from geological provinces in the Koryak Highlands (Fedorov et al., 2011) or western Kamchatka (Nürnberg et al., 2011). We investigate the most likely iceberg sources for IRD at Site 1207 using the Fine Resolution Greenland and Labrador Sea (FRUGAL) intermediate complexity coupled climate-ocean-iceberg model (Levine and Bigg, 2008).
We also investigate the range of iceberg sources to the wider NW Pacific basin by exploring the most likely iceberg sources for IRD supplied to Sites 580 and 883 (Fig. 1), where IRD has been recorded (Bigg et al., 2008, St. John andKrissek, 1999).

Model specifications and methodology
The FRUGAL (Fine ResolUtion Greenland And Labrador) intermediate complexity global climate model includes coupling between ocean, radiative-advective atmospheric and advectivethermodynamic iceberg trajectory and sea-ice models (Levine and Bigg, 2008). FRUGAL has been used for a range of palaeoclimate studies, including for the last glacial period (Fedorov et al., 2011;Bigg et al., 2012), the penultimate glaciation (Green et al., 2010) and the Early Pleistocene (Rea et al., 2018). It uses a curvilinear grid with the North Pole displaced to Greenland, giving an enhanced resolution in the Arctic and North Atlantic (Wadley and Bigg, 2000). This study uses a fine resolution grid of 182 × 211 cells, which is equivalent to approximately 2 • longitude by 1.5 • latitude in the Southern Hemisphere, but with a resolution of 20 km around the Greenland coast. The details of the ocean and sea-ice model configuration are given in a similar resolution, but present day, study (Wilton et al., 2015). The atmospheric part of FRUGAL is a simple radiative-advective atmosphere, which is an adaptation of the energy and moisture balance model of the UViC Earth Systems Model (Fanning and Weaver, 1996) that allows for advection of water vapour. The atmospheric model uses a monthly varying glacial wind stress, which has no feedback from the sea surface temperature (SST) field (Levine and Bigg, 2008). The basic topography and bathymetry of the model has been previously described (Levine and Bigg, 2008), but here we also take into account the finer resolution of the model grid in the shaping of local bathymetry for important global passages, shallow seas and the Kurile and Aleutian island chains in the North Pacific. The iceberg module has a range of both dynamic and thermodynamic components (Levine and Bigg, 2008), although in this study there is no feedback between the iceberg model and the rest of the FRU-GAL system as interest focuses on the iceberg trajectories only. Icebergs may roll over (Wagner et al., 2017) and grounded icebergs are allowed to melt instantaneously (Levine and Bigg, 2008). Model icebergs are divided into ten different size classes (Levine and Bigg, 2008) ranging from 0.491 to 492 × 10 9 kg in mass, based on observations of present day Arctic and Southern Ocean icebergs, excluding giant icebergs. Each model berg is assigned a scale factor appropriate to that size of berg from its specific seed site. Icebergs were seeded at locations where geomorphic evidence suggests that glaciers reached the coast during the late Quaternary, including on the Koryak coast, the eastern Kamchatka coast, southwest Kamchatka coast, the northern coast of the Sea of Okhotsk (Barr and Clark, 2012) and the main Kurile islands (Gorshkov, 1970), between Japan and Kamchatka. In the base simulation, the number of icebergs seeded at each location is proportional to the width of the palaeo-ice stream termini (Supplementary Information Table 2; Appendix A). This means that this simulation shows the maximal possible extent of iceberg distribution, whereas realworld iceberg calving would not occur simultaneously across the region at continuous maximum discharge. Note that as little is known about glaciation on the Kuriles, small, valley glacier-size fluxes were given to the Kurile seed points. Summing the mass of each berg, multiplied by its scale factor, gives the average ice discharge per quarter year expected from that seed site (Supplementary Information Table 2; Appendix A).

Tephrochronology
The volcanic ash component was dominated by bubble wall glass shards ∼150-500 µm in diameter, with isolated pumice fragments up to 1000 µm in diameter. The ash in ODP 1207A between 160 and 170 cm was deemed to be a primary airfall deposit due to the homogenous geochemical population of the glass and the presence of a visible tephra layer ( Supplementary Information Fig. 1, Appendix A; Griggs et al., 2014), and was therefore considered suitable for tephrochronological analysis. The glass composition obtained permitted the explicit identification of the ODP 1207A ash as the geochemically unique Aso-4 marker tephra layer from the Aso caldera on Kyushu, Japan ( Fig. 2; Aoki, 2008;Kimura et al., 2015;Derkachev et al., 2016;Albert et al., 2019). The Aso-4 tephra has been found in numerous cores in the Sea of Japan, Okhotsk Sea and NW Pacific as well as in the Lake Suigetsu sedimentary record, and was independently dated at 86.8-87.3 ka based on the oxygen isotopic stratigraphy of core MD01-2421 (off central Japan; Aoki, 2008) and at 86.4 ± 1.1k a  (2008), which is also within the range calculated by Albert et al. (2019).
3.1.2. 14 C dates Data for the five 14 C dates are shown in Table 1. These dates suggest that most of the Holocene sediments are not present in ODP 1207A, and may have been lost during drilling. The four younger dates shown in Table 1 were used in the final chronology from ODP 1207A. These dates were obtained from depths between 9 and 51 cm. Most of the core deeper than 51 cm was concluded to be beyond the age range of 14 C analysis (Reimer et al., 2013). The oldest date of 47,495 cal. yrs BP at 150-151 cm, was not included in the final chronology because it was inconsistent with the age of the Aso-4 tephra by >40 ka and likely to have been contaminated by its proximity to the core hiatus just above. The two measurements of 34-35 and 38-39 cm are time-reversed, but close together physically and temporally. It is likely that local bioturbation has mixed sediment vertically over several cm.

Age model construction
The age model for ODP 1207A was derived from a synthesis of the dating techniques described in this section. Sedimentation rates were calculated from each radiocarbon, tephra and MIS boundary age point, the latter being derived from the benthic δ 18 O curve. The MIS boundaries were determined by visually comparing the δ 18 O signal with other NW Pacific curves (Thompson and Shackleton, 1980;Kiefer et al., 2001), deriving the ages of MIS boundaries from the LR04 global benthic stack of Lisiecki and Raymo (2005). A further age tie-point was derived from counts of the planktic foraminifera N. pachyderma (sinistral coiling, as described in Darling et al., 2006), where a known (Thompson and Shackleton, 1980;Yamane, 2003) NW Pacific temperature change occurs as shown in ODP 1207A by a rapid reduction in the abundance of this cold water species from around 66% of total planktic foraminifera to around 20%. This temperature change has been dated as the boundary between the penultimate glacial period during MIS 6 and the last interglacial during MIS 5E, and thus serves as a useful age marker for the MIS 6-5 boundary in ODP 1207A. Using these chronostratigraphic indicators we determine that the upper 300 cm of ODP 1207A extends from ∼175 ka BP, in the penultimate glacial MIS 6 to 11 ka BP in the early Holocene. This age model, with the oxygen isotope data, is shown on the core summary Fig. 3.

Composition of IRD in ODP 1207A
Analysis by SEM and EDS confirmed the majority of lithic grains in ODP 1207A as quartz. Two mafic dropstones were also identified. Some quartz grains exhibited high relief and angularity, while others were rounded (Fig. 4). We hypothesise that rounded grains ( Fig. 4: Images D and E) may have been transported supraglacially, thus avoiding intense glacial weathering. The morphology of quartz grains at Site 1207 is similar to those with a known glacial origin at Site 580 (41.6 • N, 154 • E; St. John and Krissek, 1999; Fig. 1;   Fig. 4: images G-I). Both sites contain populations of rounded and angular grains, demonstrating that glacially-derived quartz does not necessarily require an angular or fractured morphology to be categorized as glacial in origin. Rather, grain size could be regarded as the more important criteria to be definitely regarded as icerafted (Hemming, 2004).

Late quaternary iceberg flux to ODP 1207A
The IRD flux to ODP 1207A is characterised by generally low levels of <3g r a i n s / g of the bulk sample (Fig. 3), similar to IRD concentrations at some sites at the same latitude in the NE Atlantic during Heinrich Events 1 and 4 (Llave et al., 2006;Eynaud et al., 2009). However, the temporal pattern of IRD is distinctly episodic, with multi-millennial periods of no ice rafting, preceding and following short periods of IRD deposition. The highest concentrations of IRD occur in late MIS 3, with three peaks of similar magnitude (Fig. 3). These peak IRD concentrations occur ∼ over 36-43 kyr BP, approximately at the same time as the large Heinrich-style IRD layer in ODP 883 further north . There are slightly enhanced levels of IRD around the depth of the Aso-4 tephra (∼87,000 yrs ago). IRD counts at this depth reveal quartz and a small number of dark lithic fragments, distinct from the layer of volcanic glass shards and pumice, suggesting a glacial, rather than volcanic origin. All other glacial stages between MIS 2 and the MIS 6-5 transition are characterised by similarly low levels of episodic IRD flux. However, IRD is absent for most of MIS 6, with ice-rafting to Site 1207 during this glacial period only occurring towards the termination and transition to interglacial MIS 5e, with no IRD present older than ∼145 ka BP.

Relationship between planktic foraminifera and IRD
The four most common species are shown in Fig. 3. N. pachyderma shows a strong preference for cold sea surface temperatures (SST), while N. incompta shows a preference for temperate waters (Chiyonobu et al., 2012). Together these two species are indicators of SST (Thompson and Shackleton, 1980;Chiyonobu et al., 2012). Fig. 3 demonstrates the different temperature preferences between N. pachyderma and N. incompta. This is particularly apparent during MIS 6 when N. pachyderma dominates the assemblage, suggesting cold conditions, but also during MIS 5, 3 and 2 when N. incompta is dominant, suggesting more temperate SST. The abundances of Globigerina (G.) bulloides, a mixed-layer-dwelling (<10 m below sea surface), and Globorotalia (Gbr.) inflata, which prefers conditions at the base of the thermocline (∼1000 m) are largely dependent on nutrient supply at these two respective depths (Channell et al., 2013). In Fig. 3 these two, while always present, are often anticorrelated, suggesting significant variation in the depth of peak nutrient supply over time.

Modelling the sources and distribution of late quaternary IRD in the NW Pacific
An 1100 yr spin-up simulation was run, using orbital parameter and atmospheric CO 2 levels from 21 kyr BP (as in previ-  Table 1). The blue date, marked by the yellow bar denotes the Aso-4 tephra. Dates in black were extrapolated from the δ 18 O stratigraphy, 14 C and tephra ages, as described in section 3.1. Marine Isotope Stages are marked by alternating grey bars. Approximate timings of Atlantic Heinrich Events are shown by H1-H6 (Hemming, 2004). The absence of an IRD bar is a zero count, as 1 cm interval IRD recording extended to 300 cm. ous glacial runs (Levine and Bigg, 2008;Bigg et al., 2012)). By the end of this simulation, all large-scale ocean and atmospheric fields had reached an essential state of equilibrium. Fig. 1 shows a schematic of the mean glacial surface ocean circulation produced by the model, in which the icebergs are advected and melted. The Kuroshio/Oyashio convergence zone appears off southern Japan, with the Oyashio Current being fed from the Sea of Okhotsk and the Bering Sea to the north. The subpolar gyre then swings north in the mid-western Pacific, while further east an outflow from the Gulf of Alaska, past the Aleutians, joins the main eastward current flow. Icebergs will enter the Oyashio Current at various places along its southward journey and be able to penetrate different distances southwards and eastwards before they melt.
To investigate these possible iceberg trajectories several decadal-long experiments were then carried out where icebergs were seeded quarterly around the North Pacific and their lifetimes followed. The two simulations shown here are: i) maximal discharges from all seed points; and ii) discharges restricted to just the five smallest size classes of icebergs (Levine and Bigg, 2008). The iceberg density fields for experiments i) -a maximal flux scenario -and ii) -a "normal" flux scenario -are shown in Fig. 5, with the origin of the icebergs near each of the three core locations of Fig. 1 shown in Fig. 6. Note that for the maximal experiment (i) site 1207 is clearly within the main modelled iceberg alley. However, in the more typical situation of simulation ii) few icebergs reach the vicinity of 1207, consistent with the episodic nature of IRD presence at this location.
Icebergs seeded on the eastern Kamchatka coast account for almost 50% of all model icebergs reaching Site 1207 (Fig. 6), which is consistent with geomorphological evidence showing that most Kamchatkan coastal end-moraines are on the eastern side of the peninsula, particularly in the northern Sredinny Range (Barr and Clark, 2012). However, the Koryak-and-west Kamchatka coasts account for 46% of icebergs reaching 1207A, which is consistent with the prevalence of quartz IRD at this site (Fedorov et al., 2011;Nürnberg et al., 2011). Over 50% of all icebergs reaching Site 883 are also calved from the E. Kamchatka coast. Icebergs reaching Site 580 originate exclusively from E. Kamchatka and W. Kamchatka, despite this site's relatively close proximity to the Kurile Islands. Model icebergs seeded from the N American coast are limited to the NE Pacific gyre and do not reach any NW Pacific sites (Fig. 5). The similar proportions of modelled iceberg origins to Sites 1207 and 883 (Fig. 6) support our hypothesis that the enhanced late-MIS 3 iceberg flux to both sites is a marker of the same iceberg-flux events.

Implications for NE Siberian glaciations during the last glacial period
The results of our marine core analysis and modelling add compelling evidence to support the hypothesis that during the late Quaternary, glaciation in NE Siberia was significant enough to deposit an IRD belt ranging from the Kamchatkan coast to the Shatsky Rise in the deep-sea NW Pacific, >1500 km from the nearest possible source of icebergs. This IRD belt is similar in scale to the 'Ruddiman Belt' of the N Atlantic (Andrews and Voelker, 2018), and the approximate southern limit of ice-rafting in the NW Pacific is shown here to be significantly further south than previously thought (Appendix B), with episodic iceberg drift southward of the SAF. It is also noteworthy that Fig. 3 shows very little evidence of temporal coupling between the NW Pacific and cold N Atlantic climatic events (marked by Heinrich Events H1-H6), unlike in the Northeast Pacific (Maier et al., 2018). We therefore also suggest that glacial change in NE Siberia occurs with different rhythms to those around the North Atlantic.
The episodic supply of icebergs to ODP 1207A, preceding and following multi-millennial periods of no iceberg supply, suggests that only the largest iceberg flux events were capable of generating sufficient bergs to reach this distal site. Therefore, ODP 1207A's IRD record serves as an important indicator of episodic large-scale glacial change or collapse on the east Siberian margin during the last, and penultimate, glacial periods. Maximum ice-rafting during MIS 3 around 40 ka BP in the 883 and 1207 IRD record also supports the hypothesis that the maximum glacial ice extent in eastern Siberia occurred ∼20,000 yrs earlier than the global LGM . We argue that the abrupt intensification of icerafting to ODP 1207A around 40 ka BP is the same ice-rafting signal observed in the dramatic Heinrich-style IRD peak in ODP 883 . This event has a basin-wide expression and can be considered a key event in the late Quaternary history of NE Siberian glaciation.
With the confirmation of an extensive IRD belt in the NW Pacific, and heightened IRD flux during late MIS 3, we strengthen the argument that a large glacial ice mass existed in NE Siberia during the last glacial period . Its collapse from its maximal state provides at least a partial explanation for the ∼ 20 m anomalous rise in global sea level around 40kyr BP (Siddall et al., 2003). Additionally, large-scale iceberg discharge into the NW Pacific, similar in magnitude to the N. Atlantic Heinrich Events, would have injected large volumes of freshwater into the ocean. There is strong and temporally persistent evidence for this in alkenone-derived salinity estimates between 40 ka and 25 ka BP in the SW Sea of Okhotsk (Harada et al., 2008), and some evidence for freshwater releases reaching the open ocean in the lower δ 18 O of ODP 883's (Kiefer et al., 2001) planktic δ 18 O record over a similar timeframe. As observed in the N. Atlantic (Hemming, 2004) and NE Pacific (Maier et al., 2018), such changes have the potential to significantly alter the regional and global climate, and further work should be carried out to investigate the climatic impact of sudden freshwater influx to the NW Pacific. This is particularly important given the modelling evidence for the possibility of enhanced modes of intermediate (Gong et al., 2019) to deep water (Okazaki et al., 2010;Green et al., 2011) formation in the northern Pacific.

Implications for NE Siberian glaciation in MIS6
Our analysis also suggests the size and behaviour of NE Siberian ice masses differ significantly between glacial periods. The longest period of no ice-rafting over the past >160,000 yrs occurs in MIS 6. According to our age model, this prolonged period of IRD absence lasts at least 30,000 yrs between at least ∼175 ka BP and ∼145 ka BP. Reduced IRD supply to the NW Pacific during midlate MIS 6 is also observed to the north-west of ODP 1207A at DSDP 580 ( Fig. 1; St. John and Krissek, 1999). IRD absence at ODP 1207A is followed by a period of ice-rafting towards the Termination of MIS 6 which coincides with the coldest SST. Nürnberg et al. (2011) also observe an increase in ice-rafting in late MIS 6 from West Kamchatka into the Sea of Okhotsk, suggesting that icebergs reaching ODP 1207A during this time may have originated from the Sea of Okhotsk coast. This is consistent with the evidence from Fig. 6 that icebergs from glacial West Kamchatka can readily reach site 1207. Alternatively, late MIS 6 ice-rafting may have been regionally synchronous, meaning the late MIS 6 IRD reaching site 1207 may have had another origin. Future, more detailed, analysis of the mineralogy of these late MIS 6 IRD fragments is needed to confirm their origin.

Implications for NW Pacific ocean climate over MIS 2-6
The analyses for site 1207 can also contribute to information about climate variability in the NW Pacific during the last two glacial periods. Changes in the relative abundance of the subpolar N. pachyderma are likely to be linked to latitudinal shifts in the position of the SAF relative to Site 1207 (Fig. 3). As observed in other cores situated in the Shatsky Rise region of the NW Pacific (Thompson and Shackleton, 1980;Yamane, 2003), peak abundances of N. pachyderma at Site 1207 occur during MIS 6 (>65% of total planktonic foram population) and to a lesser extent during mid-late MIS 2 (>25% of total planktonic foram population). During these periods, and particularly MIS 6, the SAF was south of Site 1207, and cold subpolar waters were dominant. This corroborates the lipid palaeothermometer study of Jonas et al. (2017) who argue that reduced SST in the Kuroshio-Oyashio interfrontal zone during MIS 6 and MIS 2 are related to weakening of the subtropical gyre, thus permitting the southward influx of polar and subpolar waters as far south as the coast of southern Japan.
During late MIS 6 when there was maximal abundance of N. pachyderma, there was a concurrent onset of ice-rafting to ODP 1207A, suggesting that there may have been a relationship between colder SST and iceberg activity on the NE Siberian coast. However, earlier than late MIS 6, N. pachyderma abundance was relatively high, between ∼40 and 65% of total planktic foraminifera, but IRD was absent for a prolonged period during this time. Indeed, our data, notably, show no strong correlation between the timing of iceberg events and large-scale NW Pacific surface oceanographic variability during most glacial and interglacial stages throughout the late Quaternary.
The sub-surface foraminifer record may also provide information about deeper ocean conditions. G. bulloides is most abundant during MIS 6 and MIS 3, suggesting increased nutrient supply during these periods, which may be due to increased upwelling and possible influx of subarctic waters, the latter of which was likely the case during MIS 6 as recorded in the proliferation of N. pachyderma during this period. Gbr. inflata increases in abundance during early and late MIS 5, and during MIS 2, which may be due to increased nutrient supply at depth and a signal of increased deep-intermediate water convection somewhere in the North Pacific, most likely in the Bering Sea (VanLaningham et al., 2009).

Conclusions
The new ice-rafting record presented here demonstrates that NE Siberian icebergs were episodically transported further south than previously recorded over the past ∼160 ka. Furthermore, ice-rafting events which reached ODP 1207A appear to be independent of North Atlantic Heinrich Events and Northeast Pacific climatic changes. We argue that the enhanced IRD levels ∼40 ka BP, present in both ODP 883 and ODP 1207A originate from the same glacial collapse event, or series of events. This hypothesis is supported by iceberg model simulations which find that most icebergs reaching these two sites were likely to have been calved in similar coastal regions. The discovery of IRD in the mid-latitude NW Pacific suggests the need for a re-evaluation of the scale and importance of late Quaternary regional iceberg input, and its implications for the subarctic North Pacific's contribution to the atmospheric CO 2 levels (Gray et al., 2018).
The most striking result from this analysis is the strong contrast in ice-rafting patterns between the LGM and MIS 6. We argue that this suggests differing regional glacial dynamics from one glacial period to another. Thus, during MIS6, while there is evidence of iceberg supply to other NW Pacific sites during its earlier phase (St. John and Krissek, 1999;Nürnberg et al., 2011), the absence of IRD at ODP 1207A until its end phase of Termination II suggests that ice-rafting was not of sufficient intensity for most of MIS6 to transport bergs to this distal site, despite the planktic evidence (Fig. 3) of the ocean being colder. NE Siberian ice extent may therefore have been very different in each of the last two glacial periods, with implications for ascribing reasons for global sea level palaeo-reconstructions. If the last glacial period saw more extensive ice there much earlier terrestrial evidence will have been removed. More high resolution studies of NW Pacific cores extending beyond MIS5 are needed to resolve this question.

CRediT authorship contribution statement
APM, GRB, JDM and MR designed the study while APM and GRB wrote the manuscript. APM performed IRD and ash counts and picked and prepared benthic foraminifera for stable isotope analysis and radiocarbon dating, MR and APM picked ash samples for microprobe analysis. GRB designed model experiments and performed simulations. HB picked the planktic foraminifera with support from APM. MJL performed stable isotope analysis and PJR performed radiocarbon dating. VP and MP performed microprobe and LA-ICP-MS analyses of ash samples and tephra identification. The age model was constructed by APM and GRB. All authors contributed to the final version of the manuscript.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
The IRD, ash, planktic foraminifera and benthic δ18O data will be made available at https://pangaea .de/. The model data and code can be accessed by contacting GRB.