A 600 kyr reconstruction of deep Arctic seawater δ¹⁸O from benthic foraminiferal δ¹⁸O and ostracode Mg ∕ Ca paleothermometry

The oxygen isotopic composition of benthic foraminiferal tests (δ¹⁸O_b) is one of the pre-eminent tools for correlating marine sediments and interpreting past terrestrial ice volume and deep-ocean temperatures. Despite the prevalence of δ¹⁸O_b applications to marine sediment cores over the Quaternary, its use is limited in the Arctic Ocean because of low benthic foraminiferal abundances, challenges with constructing independent sediment core age models, and an apparent muted amplitude of Arctic δ¹⁸O_b variability compared to open-ocean records. Here we evaluate the controls on Arctic δ¹⁸O_b by using ostracode $\mathrm{Mg}/\mathrm{Ca}$ paleothermometry to generate a composite record of the δ¹⁸O of seawater (δ¹⁸O_sw) from 12 sediment cores in the intermediate to deep Arctic Ocean (700–2700 m) that covers the last 600 kyr based on biostratigraphy and orbitally tuned age models. Results show that Arctic δ¹⁸O_b was generally higher than open-ocean δ¹⁸O_b during interglacials but was generally equivalent to global reference records during glacial periods. The reduced glacial–interglacial Arctic δ¹⁸O_b range resulted in part from the opposing effect of temperature, with intermediate to deep Arctic warming during glacials counteracting the whole-ocean δ¹⁸O_sw increase from expanded terrestrial ice sheets. After removing the temperature effect from δ¹⁸O_b, we find that the intermediate to deep Arctic experienced large (≥1 ‰) variations in local δ¹⁸O_sw, with generally higher local δ¹⁸O_sw during interglacials and lower δ¹⁸O_sw during glacials. Both the magnitude and timing of low local δ¹⁸O_sw intervals are inconsistent with the recent proposal of freshwater intervals in the Arctic Ocean during past glaciations. Instead, we suggest that lower local δ¹⁸O_sw in the intermediate to deep Arctic Ocean during glaciations reflected weaker upper-ocean stratification and more efficient transport of low-δ¹⁸O_sw Arctic surface waters to depth by mixing and/or brine rejection.

1 Introduction

Climate variability on orbital and suborbital timescales is amplified in polar regions as shown by changes in ocean temperature, sea-ice cover, deep-water formation, ecosystems, heat storage, and carbon cycling (Masson-Delmotte et al., 2006; Miller et al., 2010; Serreze and Barry, 2011; Cronin et al., 2017). Yet, the relationship between the Arctic Ocean and global climate remains poorly understood due in part to poor chronostratigraphic constraints of high-latitude marine orbital-scale records (O'Regan et al., 2008; Marzen et al., 2016; Wang et al., 2018) as well as challenges in applying traditional geochemical paleoceanographic proxies arising from an incomplete understanding of the biogeochemistry of Arctic Ocean systems. The global implications of a modern changing Arctic require a better understanding of past glacial–interglacial cycles in the Arctic and their link to orbital insolation cycles, atmospheric CO₂, and cryospheric processes.

The oxygen isotopic composition (${}^{\mathrm{18}}\mathrm{O}/{}^{\mathrm{16}}\mathrm{O}$, expressed as δ values in per mil [‰]) of benthic foraminifera (hereafter, δ¹⁸O_b) is a widely used tool in isotope stratigraphy and correlation of orbital cycles in open-ocean sediments (e.g., Lisiecki and Raymo, 2005; Ahn et al., 2017). Pairing δ¹⁸O_b with an independent paleotemperature estimate (for instance, from benthic microfossil $\mathrm{Mg}/\mathrm{Ca}$) potentially removes the temperature effect on δ¹⁸O_b, isolating the signature of seawater δ¹⁸O change from δ¹⁸O_b (Lear et al., 2000; Billups and Schrag, 2002; Elderfield et al., 2010, 2012).

The use of δ¹⁸O_b for stratigraphic and paleoclimate investigations in marginal seas such as the Arctic Ocean has significantly lagged its use in the open ocean. This reflects the interpretive challenges posed by marginal seas, where the histories of local temperature and δ¹⁸O_sw can greatly diverge from the global signal (e.g., Vergnaud-Grazzini et al., 1977; Cacho et al., 2000; Bauch et al., 2001; Sagawa et al., 2018). These issues are further compounded in the Arctic Ocean by variable microfossil preservation, poorly constrained chronologies, and low and variable sedimentation rates (e.g., Poirier et al., 2012; Alexanderson et al., 2014). As a result, δ¹⁸O_b has not been routinely measured in previous investigations of Arctic Ocean sediment cores.

Despite these challenges, recent work has highlighted the potential for δ¹⁸O_b applications to the Arctic Ocean. Mackensen and Nam (2014) concluded that live epifaunal benthic foraminifera calcify close to equilibrium with the oxygen isotopic composition of Arctic bottom waters, suggesting that δ¹⁸O_b faithfully records the δ¹⁸O_sw of Arctic bottom waters. Sediment cores from the central and western Arctic possess microfossil abundances sufficient for paleoceanographic reconstruction of glacial–interglacial cycles (Polyak et al., 2004, 2013; Löwemark et al., 2014; Cronin et al., 2014, 2017; Marzen et al., 2016; Wang et al., 2018). Cronin et al. (2019) presented a preliminary δ¹⁸O_b compilation from nine sediment cores obtained from the central and western Arctic Ocean spanning the last ∼600 kyr (Marine Isotope Stages (MISs) 1–13).

Here we interpret downcore δ¹⁸O_b compiled from 12 Arctic Ocean sediment cores alongside constraints on past Arctic bottom-water temperature (BWT) provided by magnesium to calcium ($\mathrm{Mg}/\mathrm{Ca}$) ratios in the calcite shells of benthic ostracodes (Dwyer et al., 1995; Cronin et al., 2012, 2017, Farmer et al., 2012). From this, we calculate a record of seawater δ¹⁸O (δ¹⁸O_sw) in the intermediate to deep Arctic Ocean. We interpret the δ¹⁸O_sw record in terms of the isotopic effects of global ice volume and drivers of Arctic hydrographic change over the last five glacial–interglacial cycles. Our results show that the Arctic Ocean experienced relatively large (∼1 ‰) local δ¹⁸O_sw variations coherent with glacial–interglacial cycles. This previously unrecognized local δ¹⁸O_sw variability helps to explain the disagreement between Arctic and open-ocean δ¹⁸O_b stratigraphies and provides a new dataset for testing hypotheses of Arctic cryosphere–ocean change over the late Quaternary.

2 Background

2.1 Controls on δ¹⁸O_b

Fossil benthic foraminifera δ¹⁸O_b values are controlled mainly by seawater temperature and local δ¹⁸O_sw, with the latter integrating contributions from changes in local hydrography and global terrestrial ice volume (Shackleton, 1967; Waelbroeck et al., 2002):

$$\begin{array}{}\text{(1)}& {\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\mathrm{b}}={\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\mathrm{T}}+{\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\mathrm{L}}+{\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\mathrm{IV}},\end{array}$$

where the subscripts T, L, and IV denote temperature, local hydrography, and ice volume, respectively.

The relationship between ambient seawater temperature (T in Eq. 1) and the fractionation of oxygen isotopes in biogenic calcite tests has been assessed in several paleoceanographic studies (e.g., Shackleton, 1974; Marchitto et al., 2014). At thermodynamic equilibrium, the δ¹⁸O of calcite increases by ∼0.20 ‰ to 0.25 ‰ for each degree Celsius (^∘C) of water cooling (Epstein et al., 1953; Shackleton, 1974; Kim and O'Neil, 1997; Matsumoto and Lynch-Stieglitz, 1999; Ravelo and Hillaire-Marcel, 2007). However, many species of foraminifera precipitate in isotopic disequilibrium with ambient seawater due to kinetic effects related to the hydration of CO₂ in seawater, metabolic effects (widely referred to as “vital effects”), and/or preservation effects related to shell dissolution (e.g., Duplessy et al., 2002; Hoogakker et al., 2010; Zeebe, 2014; Poirier et al., 2021). If these factors can be accounted for using a constant correction, as is typically applied for benthic foraminifera (e.g., Graham et al., 1981; Katz et al., 2003), isotopic calcification temperature can be expressed in a paleotemperature equation that links the δ¹⁸O of foraminiferal calcite to seawater temperature during precipitation and the δ¹⁸O of ambient seawater (Shackleton, 1974; Matsumoto and Lynch-Stieglitz, 1999; Ravelo and Hillaire-Marcel, 2007).

In addition to temperature, the growth and decay of ice sheets during Pleistocene glaciations drives changes to δ¹⁸O_sw that are expected to be globally uniform on the timescale of ocean mixing (∼1000 years) (IV in Eq. 1). Growth of ice sheets during glacial inception results in increased δ¹⁸O_sw and subsequently a greater incorporation of the heavier ¹⁸O into foraminiferal tests relative to ¹⁶O.

Finally, δ¹⁸O_sw can also reflect the input of waters with different δ¹⁸O signatures because of local hydrography (L in Eq. 1). Changes in the δ¹⁸O of seawater can result from local processes such as sea-ice formation and brine rejection, precipitation–evaporation balances, and processes related to advection. In intermediate to deep waters, the isotope signature of seawater is primarily a function of water mass flow and mixing (Ravelo and Hillaire-Marcel, 2007, and references therein; Lisiecki and Stern, 2016). The polar Arctic Ocean is uniquely affected by processes that could modify δ¹⁸O_sw. During sea-ice formation, the oxygen isotopic fractionation between ice and water is minimal (Lehmann and Siegenthaler, 1991), but the exclusion of salt tends to densify surface waters that have characteristically low δ¹⁸O in the Arctic Ocean (Fig. 1). Thus, brine rejection can lead to the sinking of low-δ¹⁸O surface waters into the intermediate and deep Arctic (Hillaire-Marcel and de Vernal, 2008; Meland et al., 2008; Stanford et al., 2011). Additionally, the episodic presence of large ice shelves (e.g., Mercer, 1970; Polyak et al., 2001; Jakobsson et al., 2010, 2016; Geibert et al., 2021), which putatively had extremely low δ¹⁸O compared to ocean waters (e.g., Spielhagen et al., 2022), would have further modified the δ¹⁸O of surface waters and potentially deep waters through mixing and/or brine formation. To isolate the effect of local hydrography, δ¹⁸O_sw records can be corrected for the influence of ice volume and are then referred to as “ice volume corrected” or δ¹⁸O_sw-ivc.

Figure 1 The oxygen isotopic composition of Arctic Ocean seawater and core locations. Measured δ¹⁸O_sw in samples from the surface (0–30 m; panel a) and intermediate to deep depths (700–2750 m; panel b) (data compiled in GLODAP2020 version 2; Olsen et al., 2020). Color bar for (a) and (b) is on the right; black box indicates inset of panel (c). (c) Locations of sediment cores used in this study and key geographic features. Basemap for (c) from the International Bathymetric Chart of the Arctic Ocean (IBCAO; Jakobsson et al., 2012). Panels (a) and (b) made in Ocean Data View (Schlitzer, 2022). Orange circles denote cores covering the last 50 kyr; white circles denote piston cores contributing to the 600 kyr composite record.

2.2 Study area: Arctic oceanography

Four primary water masses characterize the modern Arctic Ocean with characteristic water depth, temperature, salinity, and δ¹⁸O_sw ranges (Fig. 1): the Polar Surface Layer (PSL; ∼0–50 m, 0 to −2 ^∘C, ∼32 to 34 psu, 0 ‰ to $<-\mathrm{4}$ ‰), Atlantic Water (AW; ∼200 to 1000 m, ∼0 to 2 ^∘C, ∼35 psu, 0.21 ‰), Arctic Intermediate Water (AIW; 1000–2000 m, $\sim -\mathrm{0.5}$ to 0 ^∘C; 34.6 to 34.8 psu, 0.26 ‰), and Arctic Ocean Deep Water (AODW; 2000–4000 m, −0.9 to −0.95 ^∘C, 34.9–34.95 psu, 0.28 ‰) (Aagaard and Carmack, 1989; Anderson et al., 1994; Jones, 2001; Rudels et al., 2004; Mackensen and Nam, 2014). Warm AW, which enters the Arctic through the eastern Fram Strait and the Barents Sea, is separated from colder and fresher surface waters by a deep halocline (∼100–200 m). The AW flows along the Eurasian Basin continental shelf and splits into two branches at the Lomonosov Ridge, with the eastern branch flowing poleward along the eastern Lomonosov Ridge and the other branch flowing into the Amerasian Basin. As AW circulates through the Arctic Basin proper, it is transformed into what has been called the Arctic Circumpolar Boundary Current (Rudels et al., 1999). Eventually, the AW exits the Arctic through the Fram Strait on the Greenland Slope.

Our study reconstructs the properties of AIW and the upper AODW (700–2726 m) in the Amerasian Basin (including the Canada and Makarov basins; Fig. 1), which have similar δ¹⁸O_sw to, but are slightly (<0.5 ^∘C) warmer than, AIW and AODW in the Eurasian Basin. This temperature difference may reflect geothermal warming or sinking of plumes along the slope entraining warmer AL water (Timmermans and Garrett, 2006; see discussion in Rudels et al., 2015). As this modern temperature difference is small compared to the typical precision of paleotemperature reconstructions (±1.0 ^∘C; Farmer et al., 2012), we interpret our records as being representative of Arctic-wide intermediate to deep conditions.

3 Materials and methods

3.1 Materials

This study compiles measured benthic foraminifera oxygen isotopic compositions and ostracode $\mathrm{Mg}/\mathrm{Ca}$ ratios from eight Arctic piston cores (P1-92-AR-P30, P1-92-AR-P39, P1-92-AR-P40, P1-93-AR-P21, P1-92-AR-P23, P1-94-AR-P9, P1-92-AR-P15, and HLY0503-06JPC), three box cores (P1-94-AR-B8, P1-94-AR-B16, and P1-94-AR-B17), and one trigger core (HLY0503-18TC) (Cronin et al., 2012, 2017, 2019) (Table 1). All cores were recovered from the Northwind and Mendeleev ridges in the western Arctic Ocean at water depths ranging from 700 to 2726 m, except for HLY0503-18TC, which was recovered from the Lomonosov Ridge in the central Arctic (Table 1; Fig. 1c).

Table 1 Core locations, data availability, and age model sources for the combined Arctic δ¹⁸O_swreconstruction.

¹ Cronin et al. (2019). ² This study. ³ Farmer et al. (2023).

3.2 Geochemical analyses

3.2.1 Oxygen isotopes

We compiled oxygen isotope measurements from Poore et al. (1999a) and Cronin et al. (2019) made on the benthic foraminifer species Cassidulina teretis, Oridorsalis tener, and Cibicidoides wuellerstorfi (equivalent to Fontbotia wuellerstorfi; e.g., Wollenburg and Mackensen, 1998; Osterman et al., 1999) and briefly summarize the methods used. In the study of Cronin et al. (2019), foraminifera were brush-picked, and individual specimens were scored for visual preservation ranging from 1 (transparent) to 4 (visual signs of alteration). Each specimen was carefully selected to best avoid any evidence of alteration related to fine-scale dissolution, analogous to the preservation distinctions of Poirier et al. (2021). A minimum of ∼30 µg of foraminiferal calcite (two to eight specimens) was used for each stable isotope analysis. Analyses for cores AOS94-B8, AOS94-B16, and AOS94-B17 were performed on a Finnigan MAT252 with a Kiel device at Woods Hole Oceanographic Institution (WHOI) and published by Poore et al. (1999a). Analyses for HLY0503-18TC were performed on a Thermo Delta V+ with a Kiel IV device at the Lamont-Doherty Earth Observatory (LDEO). The remaining analyses were conducted at Rensselaer Polytechnic Institute (RPI) using an Isoprime dual-inlet isotope ratio mass spectrometer. LDEO and RPI analyses were first reported in Cronin et al. (2019). In these previous studies, all measurements are reported in delta notation relative to Vienna Pee Dee Belemnite (VPBD) corrected to the NBS19 standard, with an average analytical precision (1 SD) of ±0.05 (measurements from LDEO and WHOI) and (2 SD) of ±0.08 (measurements from RPI).

3.2.2 Ostracode $\mathrm{Mg}/\mathrm{Ca}$ paleothermometry

The $\mathrm{Mg}/\mathrm{Ca}$ of benthic ostracodes Krithe hunti and Krithe minima (Yasuhara et al., 2014) was previously measured via methods described in Cronin et al. (2012, 2017). Briefly, in these studies Krithe specimens were brush-picked under a binocular microscope and assigned a visual preservation index of 1 (transparent) to 7 (opaque white). Only adult specimens rating 1–5 on this index were measured for $\mathrm{Mg}/\mathrm{Ca}$. These specimens were soaked in an oxidative solution of 5 % NaOCl for 24 h to remove any reactive organic material, rinsed five times in high-purity deionized water under mild sonication to remove any residual inorganic surface material, dissolved in ultrapure 0.05 N nitric acid, and analyzed by atomic emission spectrometry on a Fisons Instruments Spectraspan 7 at Duke University using matrix-matched calibration standards from ultrapure, plasma-grade standards (SPEX brand). Unlike foraminiferal shells, which are also commonly used for marine $\mathrm{Mg}/\mathrm{Ca}$-BWT reconstructions, ostracodes produce a smooth, solid, chamber-free shell, making ostracodes far less susceptible to clay particle contamination than foraminiferal shells (Dwyer et al., 2002). Nevertheless, as a precaution, contamination-prone metals (Al, Fe, and Mn) were simultaneously monitored to ensure the absence of any signal for these constituents. External precision on $\mathrm{Mg}/\mathrm{Ca}$ ratios using this method yields a relative error (1 SD) of ≤1.5 % based on replicate analyses of an in-house limestone consistency standard prepared and analyzed along with each batch of ostracode specimens (Dwyer et al., 2002).

Krithe $\mathrm{Mg}/\mathrm{Ca}$ was converted to bottom-water temperature (BWT) using the linear least squares regional calibration of Farmer et al. (2012):

$$\begin{array}{}\text{(2)}& \mathrm{BWT}\left({}^{\circ }\mathrm{C}\right)=\mathrm{0.439}\cdot \left(\mathrm{Mg}/{\mathrm{Ca}}_{K.\mathit{hunti}}\right)-\mathrm{5.14}.\end{array}$$

Following observations of consistently elevated $\mathrm{Mg}/\mathrm{Ca}$ in K. minima relative to K. hunti, K. minima $\mathrm{Mg}/\mathrm{Ca}$ was scaled by 0.77 and converted to BWT using Eq. (2) (Cronin et al., 2012, 2017). The prediction error (1 SD) of Eq. (2) is ±1.0 ^∘C (Farmer et al., 2012).

The potential for non-thermal factors to bias Krithe $\mathrm{Mg}/\mathrm{Ca}$ has been addressed in detail in previous studies (Cronin et al., 2012, 2017; Farmer et al., 2012) and is summarized here. First, the calibration of Krithe $\mathrm{Mg}/\mathrm{Ca}$ to BWT shows no evidence of a carbonate ion effect when oxidative cleaning is performed (Farmer et al., 2012), in contrast to epifaunal benthic foraminifera (Elderfield et al., 2006). Although other studies note differences in ostracode $\mathrm{Mg}/\mathrm{Ca}$ when a reductive cleaning step is performed (Elmore et al., 2012; Gray et al., 2014), only oxidative cleaning was performed on the samples used in this study. Second, the high $\mathrm{Mg}/\mathrm{Ca}$ of glacial- and stadial-aged Arctic Krithe cannot be explained by coeval changes in sediment properties. Arctic sediments contain distinct layers of dolomite-rich ice-rafted debris (IRD) sourced from Paleozoic sediments in the Canadian Arctic (Bischof and Darby, 1997; Phillips and Grantz, 2001). However, intervals of elevated Krithe $\mathrm{Mg}/\mathrm{Ca}$ occur outside of these discrete dolomitic IRD layers (Cronin et al., 2017). Furthermore, even if dolomite dissolution did occur outside of these intervals, such a dissolution could not plausibly raise Krithe $\mathrm{Mg}/\mathrm{Ca}$. This is because dolomite dissolution releases Mg and Ca at a 1:1 ratio, thereby decreasing seawater $\mathrm{Mg}/\mathrm{Ca}$ from its current ratio ($\sim \mathrm{5}:\mathrm{1}$). Thus, if significant dolomite dissolution did occur outside of dolomite IRD layers to a sufficient degree to lower the $\mathrm{Mg}/\mathrm{Ca}$ of Arctic seawater, dolomite dissolution would be expected to lower, not raise, Krithe $\mathrm{Mg}/\mathrm{Ca}$ (Cronin et al., 2017).

3.3 Chronology

3.3.1 Box and trigger cores (last ∼50 kyr)

Chronologies for box cores AOS94-B8, AOS94-B16, and AOS94-B17 and the HLY0503-18TC trigger core are based on radiocarbon (¹⁴C) dating on Neogloboquadrina pachyderma (Poore et al., 1999a; Hanslik et al., 2010; Poirier et al., 2012). To generate consistent chronologies for these cores, radiocarbon-based age models were recalculated following the approach of Farmer et al. (2021, 2023). Briefly, ¹⁴C dates were recalibrated using Marine20 (Heaton et al., 2020) using the default marine reservoir age and without further reservoir correction (ΔR=0). The justification for using ΔR=0 is discussed in detail in Farmer et al. (2021). Age–depth models were created by Bayesian modeling in rBacon (Blaau and Christen, 2011) using a section thickness of 3 cm and an accumulation rate prior of 500 yr cm⁻¹, with all other settings as default. New ¹⁴C-based age–depth models for AOS94-B16 and HLY0503-18TC are provided in the Supplement.

3.3.2 Piston cores

The orbital-scale chronology for piston cores is unchanged from the approach detailed in Cronin et al. (2017, 2019). Briefly, an initial age model was created for each piston core using ¹⁴C dates in the upper ∼30 cm and two biostratigraphic datums. First, the benthic foraminifer Bulimina aculeata has a distinct double spike in abundance that occurs stratigraphically just below ¹⁴C-dated MIS 3 sediments, supporting its age assignment to Marine Isotope Stage (MIS) 5a (∼80 ka; Polyak et al., 2004; Cronin et al., 2014). Second, the planktonic foraminifer Turborotalita egelida exhibits a unique abundance peak in western Arctic Ocean sediment cores dated to MIS 11 (∼400 ka; Cronin et al., 2019; O'Regan et al., 2020). From these tie points, orbital-scale age models were established using cyclostratigraphy of benthic calcareous microfossil density (foraminifera and ostracodes), with intervals of high microfossil abundance aligned to interglacials (Marzen et al., 2016; Cronin et al., 2019).

Two sources of uncertainty in these age models concern the age assignment of B. aculeata to MIS 5a and the T. egelida abundance zone to MIS 11. Regarding the B. aculeata zone, alternative age models based on the extinction of excess ²³¹Pa and ²³⁰Th (∼140 and 300 ka, respectively) in Arctic sediments assign the B. aculeata zone to older ages of MIS 8/7 (Hillaire-Marcel et al., 2017). It should be noted that identifying the depth of excess ²³¹Pa and ²³⁰Th is not trivial in Arctic cores given stratigraphic and redox-driven alterations in excess ²³¹Pa and ²³⁰Th (Not and Hillaire-Marcel, 2010; Purcell et al., 2022). For instance, the causes of the absence of their excesses during some glacials and stadials are not understood and may be from dilution due to excess terrigenous inputs (Hillaire-Marcel et al., 2022) or changes in the composition of Arctic waters (Geibert et al., 2021). Additionally, it is not currently possible to unequivocally accept age models derived by excess ²³¹Pa and ²³⁰Th, as there is a major disagreement between these age models and those derived from the documented extinction and appearance of calcareous nannofossils (Jakobsson et al., 2001; O'Regan et al., 2020) and results from optically stimulated luminescence dating of quartz grains on both the central and southern Lomonosov Ridge (Jakobsson et al., 2003; West et al., 2021). Finally, the B. aculeata zone has only been dated by excess ²³¹Pa and ²³⁰Th in one location from the Lomonosov Ridge (PS87/030; Hillaire-Marcel et al., 2017) and not in any of our studied western Arctic cores. Therefore, we maintain the assignment of the B. aculeata zone to MIS 5a following previous studies (Polyak et al., 2004; Cronin et al., 2014) but encourage the application of excess ²³¹Pa and ²³⁰Th extinction to these well-studied western Arctic cores to explore the potential age revision of the B. aculeata zone.

Regarding the T. egelida abundance zone, a recent study comparing nannofossil stratigraphy from the Lomonosov Ridge (O'Regan et al., 2020) with microfossil biostratigraphy in a sediment core from the Alpha Ridge suggested that the T. egelida abundance peak could potentially be older, occurring in MIS 15 or 17 at these locations (Vermassen et al., 2021). It should be noted that all cores covering MIS 11 in our study are sourced from the Northwind and Mendeleev ridges in the western Arctic and not from the Alpha Ridge or Lomonosov Ridge where disparate ages for the T. egelida abundance peak were determined (Fig. 1c). As such, while the basin-scale synchroneity of this biostratigraphic marker is uncertain and requires further work to establish, an MIS 11 age remains the preferred interpretation for the T. egelida zone in the western Arctic. This age assignment is further supported by previous multiproxy reconstructions from these western Arctic cores that exhibit suborbital variability aligning with independently dated MIS 11 records from the Nordic Seas and North Atlantic Ocean (Cronin et al., 2019).

Taken together, consensus exists that microfossil-rich sediments characterize interglacials and interstadials in the Arctic Ocean (e.g., Poore et al., 1993; Spielhagen et al., 1997, 2004; Polyak et al., 2009, 2013; Hanslik et al., 2013; Marzen et al., 2016), but we acknowledge that the assignment of specific interglacials prior to MIS 5 in Arctic sediments is still uncertain and may require further stratigraphic refinement. Therefore, we focus our results and discussion primarily on comparing interglacial and glacial intervals and note that the exact interglacial assignment of features in our records may be subject to future revision.

3.4 Calculation of δ¹⁸O_sw

We used PSU Solver (Thirumalai et al., 2016) to calculate δ¹⁸O_sw and its probabilistic uncertainty (presented throughout at 32 % and 68 % 1 SD equivalent confidence intervals) derived from 1000 bootstrap Monte Carlo propagations of Krithe $\mathrm{Mg}/\mathrm{Ca}$ and δ¹⁸O_b errors. To perform this calculation, C. wuellerstorfi-equivalent δ¹⁸O_b was first converted to Uvigerina scale by adding 0.64 ‰ (Shackleton, 1974). Next, as our $\mathrm{Mg}/\mathrm{Ca}$ and δ¹⁸O_b data were measured in separate phases, Krithe $\mathrm{Mg}/\mathrm{Ca}$ and δ¹⁸O_b data had to be matched. Of a total of 352 δ¹⁸O_b measurements, 126 δ¹⁸O_b measurements had Krithe $\mathrm{Mg}/\mathrm{Ca}$ measured in the same sediment sample (Cronin et al., 2017). For the remaining unmatched samples, we interpolated the higher-resolution ostracode $\mathrm{Mg}/\mathrm{Ca}$ stack from Cronin et al. (2017) to estimate Krithe $\mathrm{Mg}/\mathrm{Ca}$ for each unmatched δ¹⁸O_b. Third, accounting for the analytical error in the $\mathrm{Mg}/\mathrm{Ca}$ measurement (±1.5 %, Sect. 3.2.2) led to an overly confident BWT prediction relative to the stated 1 SD calibration error of ±1.0 ^∘C (Farmer et al., 2012). To address this, we set the $\mathrm{Mg}/\mathrm{Ca}$ measurement uncertainty in PSU Solver to ±2.15 mmol mol⁻¹ such that BWT calculated by PSU Solver using Eq. (2) averaged ±1.0 ^∘C error at 32 % and 68 % (1 SD equivalent) confidence. Last, δ¹⁸O_sw was calculated using the linear approximation of Shackleton (1974):

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\mathrm{sw}}\left(\mathrm{‰}\mathrm{VSMOW}\right)& =\left({{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\mathrm{b}\text{-}\mathrm{uvi}}+\mathrm{0.27}\right)\\ & +\mathrm{0.25}\cdot (\mathrm{BWT}-\mathrm{16.9}),\end{array}\end{array}$$

where δ¹⁸O_b-uvi is the C. wuellerstorfi δ¹⁸O_b converted to Uvigerina scale versus VPBD.

Figure 2 Chronological constraints on Arctic Ocean sediments over the last 50 kyr and comparison of Arctic δ¹⁸O_b to δ¹⁸O_b from selected deep ocean sites. (a) δ¹⁸O_b data. White, light blue, medium blue, and dark blue diamonds indicate δ¹⁸O_b data points from AOS94-B8, AOS94-B16, AOS94-B17, and HLY-0503-18tc, respectively. Blue line indicates smoothed spline fit (created in MATLAB using the smoothing parameter p=0.5); shading indicates ±1 RMSE for spline smoothing function. Other δ¹⁸O_b records shown are the LR04 global stack (black line; Lisiecki and Raymo, 2005), CHN-82-24-4PC (red line; Boyle and Keigwin, 1985), ODP 1123 (orange line; Elderfield et al., 2010), and KNR0734-3PG (brown line; Curry et al., 1988); all records are shown as C. wuellerstorfi-equivalent δ¹⁸O_b values. (b) Abundance of planktonic foraminifera per gram of sediment for western Arctic sites B8 (gray), B16 (light blue), and B17 (medium blue) (left axis) and the central Arctic Site HLY0503-18TC (dark blue, right axis). (c) Radiocarbon dates on planktonic foraminifera; colors the same as in panel (a).

Extra-Arctic δ¹⁸O_sw reconstructions from DSDP Site 607 and ODP Site 1123 were also calculated with PSU Solver using Eq. (3) and BWT derived from the $\mathrm{Mg}/\mathrm{Ca}$-BWT calibrations appropriate for benthic foraminifera used at those locations (calculations from Ford and Raymo, 2020).

4 Results

The results are separated into two intervals: the last 50 kyr and the entire 600 kyr record. This distinction reflects the different cores used (box and trigger versus piston cores), chronological constraints (radiocarbon versus tuning based on benthic microfossil density), and data availability (greater in the last 50 kyr) between the two intervals. We first investigate the last 50 kyr to contextualize the longer, albeit discontinuous records.

4.1 Arctic δ¹⁸O_b and δ¹⁸O_sw: last 50 kyr

Radiocarbon-dated Arctic δ¹⁸O_b records from AOS94-B8, AOS94-B16, AOS94-B17, and HLY-0503-18TC are shown with representative extra-Arctic δ¹⁸O_b records in Fig. 2. A smoothed spline interpolation to the four Arctic δ¹⁸O_b records fits within ±0.17 ‰ (root mean square error), suggesting that different Arctic core locations record similar δ¹⁸O_b patterns for the last 50 kyr. However, this Arctic-wide δ¹⁸O_b is noticeably distinct from representative global ocean records, particularly since 20 ka. Between 50 and 20 ka, a period encompassing Marine Isotope Stages (MIS) 3 and 2, Arctic and global records all converge on δ¹⁸O_b values of 3.5 ‰ to 4.5 ‰ (Fig. 2a). Over the last 10 kyr (approximately the Holocene), Arctic δ¹⁸O_b averages 3.80±0.14 ‰ versus 2.68±0.10 ‰ for the global LR04 δ¹⁸O_b stack (Lisiecki and Raymo, 2005), 2.74±0.18 ‰ in the North Atlantic (CHN82-24; Boyle and Keigwin, 1985, 2.72±0.04 ‰ in the southwest Pacific Ocean (ODP 1123; Elderfield et al., 2010), and 2.66±0.25 ‰ in the equatorial Pacific (KNR0734-3PG; Curry et al., 1988) (all δ¹⁸O_b data reported on C. wuellerstorfi scale). As a result, whereas δ¹⁸O_b in open-ocean cores was 1 ‰ to 2 ‰ higher during MIS 3 and 2 compared to the Holocene, Arctic δ¹⁸O_b was only ∼0.4 ‰ to 0.8 ‰ higher during MIS 3 and 2 compared to the Holocene (Fig. 2a).

Based on planktonic foraminifer radiocarbon dates, previous studies have noted that sedimentation in the Arctic Ocean was extremely condensed or possibly absent during part of MIS 2, including the Last Glacial Maximum (LGM) and early deglaciation (Darby et al., 1997; Poore et al., 1999b; Polyak et al., 2009; Hanslik et al., 2010; Poirier et al., 2012). In the four Arctic cores used for our composite 50 kyr record, planktonic foraminifer abundance is notably reduced between ∼35 and ∼14 ka (Fig. 2b). Additionally, only eight radiocarbon dates fall within this interval, with no dates calibrated to between 21 and 13 ka (Fig. 2c). Given this, it is likely that the Arctic δ¹⁸O_b composite does not capture peak LGM conditions ca. 20 ka or the early deglaciation (Fig. 2a). However, this absence would not greatly affect our interpretation as there is only a modest (<0.5 ‰) δ¹⁸O_b difference between the peak LGM and late MIS 3 (∼35 ka) in both regional records and the LR04 δ¹⁸O_b stack (Fig. 2a). In summary, the potential absence of Arctic benthic foraminifera from peak LGM conditions would not compromise our ability to capture the majority of the glacial–interglacial δ¹⁸O_b change observed in extra-Arctic records because Arctic benthic foraminifera from late MIS 3 should capture near-LGM δ¹⁸O_b values.

Since both the Arctic and extra-Arctic δ¹⁸O_b should have experienced the same contribution from ice volume changes (e.g., the same Δδ¹⁸O_IV), differences between Arctic and extra-Arctic δ¹⁸O_b records must reflect different temperature (Δδ¹⁸O_T) and/or local (Δδ¹⁸O_L) effects (Eq. 1). Figure 3 deconvolves the contributions to δ¹⁸O_b in the Arctic and ODP Site 1123 in the southwest Pacific (Elderfield et al., 2010). Site 1123 was chosen for comparison because of its high-resolution δ¹⁸O_b and infaunal benthic foraminifera (Uvigerina) $\mathrm{Mg}/\mathrm{Ca}$-derived BWT records covering the last 50 kyr. At Site 1123, Uvigerina $\mathrm{Mg}/\mathrm{Ca}$ was converted to BWT using the preferred equation of Elderfield et al. (2010).

Figure 3 (a) δ¹⁸O_b, (b) $\mathrm{Mg}/\mathrm{Ca}$-derived BWT, (c) δ¹⁸O_sw, and (d) δ¹⁸O_sw-ivc (ice volume corrected) over the last 50 kyr from Arctic cores (diamonds are individual data points; blue lines are composited records) and ODP 1123 (orange circles are individual data points; orange lines are 1 kyr interpolations; Elderfield et al., 2010). Shading denotes propagated 1 SD equivalent errors (32 % to 68 % confidence intervals) calculated in PSU Solver (Thirumalai et al., 2016). Blue and orange circle in (c) and (d) denote the modern δ¹⁸O_sw in Arctic Deep Water and at ODP Site 1123, respectively (Mackensen and Nam, 2014; Adkins et al., 2002).

In addition to differences in δ¹⁸O_b (Fig. 3a), reconstructed BWT is also notably different between Site 1123 (Elderfield et al., 2010) and the Arctic (Cronin et al., 2012) for the period since 50 ka (Fig. 3b). During the interval encompassing late MIS 3 and the LGM, Arctic BWTs were higher (1 to 3 ^∘C) in comparison to Site 1123 (−1 to 1 ^∘C). However, this pattern reversed in the late deglaciation and over the Holocene, with ostracode $\mathrm{Mg}/\mathrm{Ca}$-derived BWTs of <1 ^∘C compared to benthic foraminifera $\mathrm{Mg}/\mathrm{Ca}$-derived BWTs at Site 1123 of ∼2 ^∘C. Holocene $\mathrm{Mg}/\mathrm{Ca}$-derived BWTs at both locations are consistent with modern bottom waters that are warmer at Site 1123 (1.3 ^∘C) compared to the intermediate and deep Arctic (−0.3 to −0.9 ^∘C).

Arctic and Site 1123 δ¹⁸O_sw records calculated from the respective δ¹⁸O_b and BWT histories are shown in Fig. 3c. Arctic and Site 1123 δ¹⁸O_sw were similar during the interval encompassing late MIS 3 and the LGM, with slightly higher δ¹⁸O_sw in the Arctic. As for BWT, Arctic and Site 1123 δ¹⁸O_sw diverged in the late deglacial or Early Holocene, with Arctic δ¹⁸O_sw elevated by up to 0.5 ‰ to 1 ‰ over Site 1123 in the Early Holocene. This is directionally consistent with, albeit slightly higher than, the modern ocean δ¹⁸O_sw elevation in the intermediate to deep Arctic Ocean (0.26 ‰, Fig. 1b; Mackensen and Nam, 2014) relative to Site 1123 (−0.10 ‰; Adkins et al., 2002) (orange and blue circles in Fig. 3c).

To investigate the local hydrographic contribution to δ¹⁸O_sw at each location, we subtracted the estimated mean ocean δ¹⁸O_sw change attributed to ice sheets (from Waelbroeck et al., 2002) from Arctic and Site 1123 δ¹⁸O_sw records; following convention, we refer to these as ice volume corrected (δ¹⁸O_sw-ivc). The δ¹⁸O_sw-ivc record of Site 1123 is roughly constant, with only a slight δ¹⁸O_sw-ivc elevation before the LGM (0 ‰ to 0.2 ‰) compared to after the LGM (−0.2 ‰ to 0 ‰) (Fig. 3d). In contrast, the Arctic record shows a δ¹⁸O_sw-ivc minimum (−0.2 ‰ to 0 ‰ ) during the LGM and early deglaciation, with higher δ¹⁸O_sw-ivc (0 ‰ to 0.4 ‰) during the Holocene and MIS 3. We note that the Arctic δ¹⁸O_sw-ivc minimum falls in the interval where foraminifer abundance and radiocarbon dating are particularly limited (Fig. 2a), and the presence of this δ¹⁸O_sw-ivc minimum should be confirmed with additional data. Nonetheless, these results suggest that, in contrast to its elevated δ¹⁸O_sw today, the LGM and early deglacial intermediate to deep Arctic contained seawater with an oxygen isotopic composition similar to the deep southwest Pacific Ocean.

4.2 Arctic δ¹⁸O_b and δ¹⁸O_sw: last 600 kyr

We next assess whether the >0.5 ‰ of δ¹⁸O_sw variation in the intermediate to deep Arctic Ocean over the last 50 kyr also characterized previous glacial cycles. One caveat to calculating Arctic δ¹⁸O_sw-ivc on these longer timescales is that the ice volume contribution to global δ¹⁸O_sw is weakly constrained over older glacial cycles in part because of a paucity of geologic sea-level markers (e.g., Hibbert et al., 2016). Therefore, instead of calculating Arctic δ¹⁸O_sw-ivc by correcting Arctic δ¹⁸O_sw for global ice volume as we have done for the last 50 kyr (Fig. 3d), we estimate its equivalence by calculating the δ¹⁸O_sw gradient between the Arctic composite and two deep ocean δ¹⁸O_sw reconstructions: one from ODP Site 1123 covering the last 600 kyr (Elderfield et al., 2012) and one from DSDP Site 607 in the North Atlantic covering 150 to 600 ka (Ford and Raymo, 2020). We assume that changes to the δ¹⁸O_sw gradient between the Arctic and these locations (expressed as Δδ¹⁸O_sw) are driven by changing local Arctic δ¹⁸O_sw.

Figure 4 (a) Individual sample Arctic δ¹⁸O_b (blue diamonds) and three-point running mean (blue line) compared to interpolated ODP Site 1123 (orange line; Elderfield et al., 2010, 2012) and global ocean LR04 stack (gray line; Lisiecki and Raymo, 2005). (b) Reconstructed bottom-water temperatures from Arctic ostracode $\mathrm{Mg}/\mathrm{Ca}$ (blue circles) and five-point running mean (blue line; Cronin et al., 2012, 2017) and Site 1123 benthic foraminiferal $\mathrm{Mg}/\mathrm{Ca}$ (orange line). (c) δ¹⁸O_sw reconstructed from the Arctic (blue diamonds are discrete data points; blue line is three-point running mean), Site 1123 (orange line; shading is 1 SD propagated error; data from Elderfield et al., 2012, and recalculated by Ford and Raymo, 2020), and Site 607 (purple circles; shading is 1 SD propagated error; Ford and Raymo, 2020). (d) Δδ¹⁸O_sw (Arctic – 1123 in orange; Arctic – 607 in blue) as a proxy for local Arctic δ¹⁸O_sw change. The modern Δδ¹⁸O_sw between the Arctic and Site 1123 is indicated by the dashed orange line. Gaps in the Arctic record illustrate intervals with limited microfossil abundance. Interglacial marine isotope stages are numbered red; glacial marine isotope stages are shaded gray and numbered; interglacial substages 7b and 13b are shaded blue.

First comparing δ¹⁸O_b, Arctic δ¹⁸O_b typically exhibits a lower range than δ¹⁸O_b at Site 1123 or the LR04 stack over the last 600 kyr, as previously noted for the last 50 kyr (Figs. 3a and 4a). Arctic δ¹⁸O_b was generally higher than Site 1123 or the LR04 stack during interglacials (MIS 1, 5, 7, 9, 11, and 13) (Fig. 4a), but the δ¹⁸O_b values for the three records were generally comparable during glacial periods (MIS 4, 8, and 14) with the caveat that peak glacial intervals may not be fully represented in the Arctic δ¹⁸O_b as was the case for MIS 2 (Fig. 2a). Arctic BWTs were similar to or slightly lower than Site 1123 during interglacials (MIS 5 and putatively MIS 9, 11, and 13 based on our age models, but curiously not during MIS 7) and were generally higher than Site 1123 during glacials (MIS 2, 4, 6, and 8). Note the general paucity of benthic foraminifera (and hence δ¹⁸O_b data) in Arctic sediments during peak glacial MIS 6, and putatively MIS 10 and 12 that leads to gaps in our records. The assignment of these gaps to peak glacials results from our age model construction because intervals of low microfossil density are attributed to glaciations (Sect. 3.3 and Marzen et al., 2016).

Comparing δ¹⁸O_sw, both sites 607 and 1123 show the recognized pattern of higher glacial and lower interglacial δ¹⁸O_sw from the waxing and waning of continental ice sheets (Fig. 4c; Elderfield et al., 2010; Ford and Raymo, 2020). However, this pattern is less distinct in Arctic δ¹⁸O_sw. The highest Arctic δ¹⁸O_sw of ∼1 ‰ is not unique to glacials; while such values are observed during MIS 2 and 4, they also occur during MIS 5, 7, and 9. Similarly, the lowest Arctic δ¹⁸O_sw of ∼0 ‰ to 0.3 ‰ is not unique to interglacials, occurring during MIS 1, 7, 11, and 13 but also MIS 6 (Fig. 4c). In general, Arctic δ¹⁸O_sw appears higher than at sites 607 and 1123 during most interglacials (MIS 1, 5, 9, 11, and 13) and lower during glacials (most notably MIS 6 and 14).

To quantitatively compare δ¹⁸O_sw between the intermediate to deep Arctic and the deep North Atlantic (Site 607) and southwest Pacific (Site 1123), we interpolated the δ¹⁸O_sw records to 3 kyr resolution and calculated δ¹⁸O_sw gradients (as ${\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\mathrm{sw}}={{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\mathrm{sw}\mathrm{\_}\mathrm{Arctic}}-{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\mathrm{sw}\mathrm{\_}\mathrm{607}\mathrm{or}\mathrm{1123}}$). To avoid gradient artifacts, Δδ¹⁸O_sw was not calculated over intervals where gaps ≥ 5 kyr existed in the Arctic δ¹⁸O_sw record, particularly during glacial maxima in MIS 6, 8, 10, and 12 (Fig. 4d). The δ¹⁸O_sw gradient between the Arctic and Site 1123 (orange line, Fig. 4d) averages 0.27±0.31 ‰ (1 SD, n=120) during interglacials over the last 600 kyr, which agrees with the modern δ¹⁸O_sw gradient of 0.36 ‰ (Fig. 3d). However, Δδ¹⁸O_sw is significantly lower (at α=0.05) during glacials (0.03±0.34 ‰, n=71) over the last 600 kyr (two-tailed heteroscedastic t test, P<0.001). The δ¹⁸O_sw gradient between the Arctic and Site 607 confirms these findings, with significantly higher (P=0.009) interglacial Δδ¹⁸O_sw (0.25±0.43 ‰, n=77) than glacial Δδ¹⁸O_sw (0.04±0.40 ‰, n=42). In addition to the broad glacial–interglacial differences, there are notable orbital and suborbital variations of up to 1 ‰ that are particularly well-resolved during MIS 6, 7, 11, and 13. Such variations appear consistent between gradients calculated relative to sites 607 and 1123, suggesting that the origin of these δ¹⁸O_sw variations lies within the intermediate to deep Arctic Ocean (Fig. 4d).

5 Discussion

Our combined ostracode $\mathrm{Mg}/\mathrm{Ca}$ and benthic foraminifera δ¹⁸O data indicate that the intermediate to deep Arctic Ocean experienced large (≥1 ‰) changes in local δ¹⁸O_sw relative to both the southwest Pacific Ocean and the North Atlantic over the last five glacial cycles (Fig. 4). The presence of large local δ¹⁸O_sw variations is further illustrated in Fig. 5. Arctic samples with paired ostracode $\mathrm{Mg}/\mathrm{Ca}$-derived BWT and δ¹⁸O_b exhibit a positive correlation between δ¹⁸O_b and BWT (Fig. 5a); this is the opposite of the trend expected based on inorganic calcite precipitation (Kim and O'Neil, 1997). In contrast, at Site 1123, lower δ¹⁸O_b corresponds to higher $\mathrm{Mg}/\mathrm{Ca}$-derived BWT in the same samples (Elderfield et al., 2010, 2012), following the trajectory expected from inorganic calcite precipitation (Fig. 5b). Thus, it appears that local δ¹⁸O_sw variations in the intermediate to deep Arctic Ocean were sufficiently large as to fully overwrite the influence of temperature on benthic foraminiferal δ¹⁸O.

Figure 5 Paired $\mathrm{Mg}/\mathrm{Ca}$-derived bottom-water temperature (^∘C) and equilibrium δ¹⁸O_b (calculated by subtracting 0.35 ‰ from C. wuellerstorfi δ¹⁸O_b; Mackensen and Nam, 2014) for the Arctic Ocean (a) and ODP Site 1123 (b). The black dashed line is the linear regression to the data; the red line is the equilibrium change in calcite δ¹⁸O as a function of temperature (Kim and O'Neil, 1997; Matsumoto and Lynch-Stieglitz, 1999).

In the last 50 kyr, the lowest Arctic δ¹⁸O_sw-ivc was reconstructed from 29 to 14 ka, encompassing the peak of the LGM and the early deglaciation (Fig. 3d) but with high uncertainty considering limited sedimentation and foraminifera in this interval (Fig. 2a). Identifying the phasing of previous low local δ¹⁸O_sw intervals is even more uncertain given chronological challenges and the absence of δ¹⁸O_b from previous glacial maxima. Nevertheless, we observe a recurring pattern of lower intermediate-to-deep-Arctic δ¹⁸O_sw relative to the deep North Atlantic or southwest Pacific during glacial periods over the last 600 kyr (Fig. 4d). This observation requires a source of ¹⁸O-depleted water to the intermediate and deep Arctic during glacial periods (and possibly deglaciations) that was absent during interglacials. We next evaluate several potential candidates for this glacial ¹⁸O-depleted water source.

5.1 Ice shelves and the “fresh Arctic” hypothesis: a critical evaluation

Based on numerical ice models and erosional features along Arctic ridges indicated in geophysical surveys, it has been argued that large ice shelves existed within the Arctic Ocean during at least some previous glacial maxima (Polyak et al., 2001; Jakobsson et al., 2010, 2016; Niessen et al., 2013; Nilsson et al., 2017; Gasson et al., 2018). Recently, Geibert et al. (2021) proposed that the Arctic Ocean became filled with freshwater, at least down to 2500 m water depth, during two intervals assigned to MIS 4 (62–70 ka) and 6 (131–151 ka). The evidence for this fresh Arctic hypothesis comes from the absence of excess ²³⁰Th in multiple Arctic sediment cores ranging from 1000 to 2700 m depth. Because ²³⁰Th is sourced from the radioactive decay of uranium, a conservative element in seawater, Geibert et al. (2021) interpreted the absence of excess ²³⁰Th as representing the absence of uranium in Arctic waters, which could be explained by the displacement of intermediate to deep Arctic seawater with freshwater in these glacial intervals. Additionally, the absence of excess ²³⁰Th corresponded to the absence of cosmogenic ¹⁰Be in several cores, suggesting that the absence of excess ²³⁰Th was coeval with intervals where the Arctic Ocean was shielded from the input of cosmogenic nuclides, potentially by circum-Arctic ice shelves (Geibert et al., 2021).

The fresh Arctic hypothesis (Geibert et al., 2021) has generated significant discussion about both the interpretation of excess ²³⁰Th and cosmogenic ¹⁰Be in Arctic sediments (Hillaire-Marcel et al., 2022; Geibert et al., 2022a) and the evidence supporting or refuting this hypothesis from outside the Arctic Ocean (Spielhagen et al., 2022; Geibert et al., 2022b). Importantly, any significant contribution of high-latitude, ¹⁸O-depleted freshwater should be readily traceable in Arctic δ¹⁸O_sw, provided that our signal carriers (benthic foraminifera and ostracodes) were present to trace intervals of freshwater input. Importantly, the foraminifera and ostracodes need not be coeval with the freshwater event. For comparison, during the last deglaciation, a flooding event from the Mackenzie River estimated to last <1000 years (Keigwin et al., 2018) appeared as a low-δ¹⁸O interval covering 5–10 cm of sediment in AOS94 box cores with similarly low sedimentation rates (Poore et al., 1999b). The amount of freshwater from the deglacial Mackenzie River flooding event would have been orders of magnitude less than would be required to freshen the entire Arctic Ocean. Thus, the Arctic-wide freshening events proposed by Geibert et al. (2021) would conceivably leave a lasting imprint on δ¹⁸O_sw that should be detectable even with the low-resolution, time-averaged sediments of our study.

We evaluate the coherence of the fresh Arctic hypothesis with our data through two approaches: an isotope mass balance and an analysis of the temporal trends in our records and the excess ²³⁰Th records of Geibert et al. (2021).

5.1.1 Oxygen isotope mass balance

We created a simple mixing model to calculate how Arctic δ¹⁸O_sw would change with various contributions of freshwater input from the melting of an ice shelf with δ¹⁸O_w of −20 ‰ and −40 ‰ (Fig. 6a). The addition of only 1 % of freshwater to the intermediate to deep Arctic Ocean would decrease Arctic δ¹⁸O_sw by 0.2 ‰ and 0.4 ‰ for a −20 ‰ and −40 ‰ freshwater source, respectively, with additional δ¹⁸O_sw decline scaling linearly with additional freshwater inputs (Fig. 6a).

Figure 6 Two-component mixing models for the δ¹⁸O_sw of Arctic Deep Water (a) by mixing with hypothetical ice melt of −20 ‰ (blue) and −40 ‰ (cyan) and (b) based on the modern δ¹⁸O_sw–salinity relationship of the deep Arctic Ocean (LeGrande and Schmidt, 2006). Gray boxes denote the observed range of local Arctic δ¹⁸O_sw (from Δδ¹⁸O_sw) in our records (Fig. 4d). Dashed arrows in (a) indicate the maximum ice melt contribution needed to explain the lowest Arctic δ¹⁸O_sw in our records.

Our Δδ¹⁸O_sw reconstruction (Fig. 4d) indicates that, over the last 600 kyr, the intermediate to deep Arctic Ocean experienced at most a −0.7 ‰ δ¹⁸O_sw deviation relative to Site 1123 when both benthic foraminifera and ostracodes were present in Arctic sediments. Were this Arctic δ¹⁸O_sw decline attributed solely to ice melt, this could be explained by at most a 5 % contribution of meltwater from a −20 ‰ ice shelf or a 2.5 % contribution of meltwater from a −40 ‰ ice shelf to the modern δ¹⁸O_sw of the intermediate to deep Arctic Ocean (Fig. 6a). Alternatively, considering the modern δ¹⁸O_sw–salinity relationship of the deep Arctic Ocean (LeGrande and Schmidt, 2006), a local δ¹⁸O_sw decline to −0.7 ‰ would correspond to a 1.1 unit salinity reduction (Fig. 6b). Thus, the Arctic-specific δ¹⁸O_sw variations observed in our 600 kyr records can be explained by only a small percentage contribution of freshwater and are not consistent with the large freshwater volumes proposed by Geibert et al. (2021).

It should be noted that our records exhibit gaps in δ¹⁸O_b data coverage due to absent benthic foraminifera and ostracodes, particularly during glacial maxima intervals ascribed to MIS 6, 10, and 12 in our age models. Such occurrence gaps could be interpreted as representing freshwater intervals (as in Geibert et al., 2021) that excluded marine benthic foraminifera. However, we discount that the existence of these gaps by themselves represents intervals of a fresh Arctic Ocean (at least down to 2500 m water depth) for two reasons. First, the absence of benthic foraminifera and ostracodes is not uniquely attributable to freshwater. Benthic microfossil abundance is controlled by numerous factors, including dissolution, ecological preferences such as food availability, and dilution with terrigenous material (e.g., Wollenburg and Kuhnt, 2000; Wollenburg et al., 2004; Scott et al., 2008; Cronin et al., 2013; Polyak et al., 2013; Lazar and Polyak, 2016). In addition, ostracodes inhabit both fresh and marine environments and experience large faunal transitions in many locations, yet no freshwater ostracodes have been observed from Arctic Ocean sediments (Cronin et al., 1995; Poirier et al., 2012). Instead, terrigenous material dilution provides an alternative explanation consistent with all observations, including low excess ²³⁰Th, low-to-absent ¹⁰Be, and limited microfossil abundance (e.g., Hillaire-Marcel et al., 2022; but see rebuttal by Geibert et al., 2022a).

Second, even when missing data from glacial maxima, our Arctic data proximal to these glacial maxima do not show extremely low-δ¹⁸O_sw values, as might be expected from residual freshwater and the time averaging inherent in our low sedimentation rate sites. While our records lack the chronological resolution to speak to millennial-scale events, it is noteworthy that we reconstruct similarly low-δ¹⁸O_sw values during interglacial substages MIS 7b and 13b and early during MIS 11, as we also do during MIS 6, 10, and 14 (Fig. 4d). Provided that this age assignment is correct, the occurrence of low-δ¹⁸O_sw intervals outside of glacial maxima rules out an explanation for their cause that is unique to glaciations. As above, these δ¹⁸O_sw values can be explained by a small (<5 %) contribution of freshwater to the deep Arctic and are not consistent with wholesale shifts to a freshwater Arctic Ocean on millennial or longer timescales.

5.1.2 Temporal trends in δ¹⁸O_sw and excess ²³⁰Th

Geibert et al. (2021) argued that two intervals of low excess ²³⁰Th observed across multiple Arctic Ocean sediment cores date to MIS 4 and 6 and that those intervals arose because freshwater filled the deep basins of the Arctic. We assess the comparability of our δ¹⁸O_sw reconstructions to excess ²³⁰Th data from Site PS2185-3/6 on the Lomonosov Ridge (Fig. 1c), which shows a similar excess ²³⁰Th pattern observed elsewhere in the Arctic (Geibert et al., 2021). To evaluate this argument, it must first be understood that the age assignment from Geibert et al. (2021) for these two intervals is not consistent with the more widely accepted age model published for this core (Spielhagen et al., 2004; see below) or other marine sediment records from this region of the Arctic. The longer-standing and more commonly cited age model would place the two low-excess ²³⁰Th intervals observed in PS2185-3/6 in a thick (70–80 cm) coarse-grained unit deposited during MIS 3 and 4 (Spielhagen et al., 2004; see below).

Figure 7 LR04 δ¹⁸O_b stack (a; Lisiecki and Raymo, 2005), local intermediate-to-deep-Arctic δ¹⁸O_sw from differencing Site 1123 (orange) and Site 607 (blue) (b; this study), central Arctic Site PS2185-3/6 metal oxide ε_Nd (c; Haley et al., 2008), excess ²³⁰Th (d; Geibert et al., 2021), and planktonic foraminifera abundance (e; Spielhagen et al., 2004). The age model for the PS2185-3/6 data (c–e) was based on aligning low excess ²³⁰Th to MIS 4 and 6 as proposed by Geibert et al. (2021). High ε_Nd (plotted downward) is argued to represent increased input of Siberian weathering products to the intermediate Arctic Ocean, presumably via brine formation on the Eurasian shelves (Haley et al., 2008). Gray shading indicates glacials MIS 2, 4, and 6; orange shading indicates interglacials MIS 1, 3, and 5a, 5c, and 5e.

Figure 8 As in Fig. 7 but with PS2185-3/6 data plotted on the age model of Spielhagen et al. (2004). LR04 δ¹⁸O_b stack (a; Lisiecki and Raymo, 2005), local intermediate-to-deep-Arctic δ¹⁸O_sw from differencing Site 1123 (orange) and Site 607 (blue) (b; this study), central Arctic Site PS2185-3/6 metal oxide ε_Nd (c; Haley et al., 2008), excess ²³⁰Th (d; Geibert et al., 2021), and planktonic foraminifera abundance (e; Spielhagen et al., 2004). High ε_Nd (plotted downward) is argued to represent increased input of Siberian weathering products to the intermediate Arctic Ocean, presumably via brine formation on the Eurasian shelves (Haley et al., 2008). Gray shading indicates glacials MIS 2, 4, and 6; orange shading indicates interglacials MIS 1, 3, and 5a, 5c, and 5e. Note the alignment of LR04 δ¹⁸O_b stack minima (a) with PS2185-3/6 planktonic foraminifera abundance maxima (e).

To evaluate this age assignment, we additionally present planktonic foraminifer abundance (Spielhagen et al., 2004), an indicator of glacial–interglacial changes in Arctic sedimentation (Fig. 2; Poore et al., 1993; Polyak et al., 2009; Marzen et al., 2016), and neodymium isotopes (ε_Nd) of sediment metal oxide coatings, an indicator of the contribution of weathering input from the Siberian Traps to Arctic intermediate waters via brine formation on the Eurasian shelves (Haley et al., 2008), from PS2185-3/6. We created an age model for PS2185-3/6 that aligns the two intervals of low excess ²³⁰Th with MIS 4 (62–70 ka) and MIS 6 (131–150 ka) as proposed by Geibert et al. (2016). The PS2185-3/6 ε_Nd, ²³⁰Th, and planktonic foraminifera abundance data are plotted on the Geibert et al. (2021) age model in Fig. 7. Although this age model aligns the excess ²³⁰Th minima to MIS 4 and 6 (Fig. 7d), it creates unrealistic interpretations for the other datasets. First, the Haley et al. (2008) ε_Nd record shows an extended period of more positive ε_Nd throughout MIS 4, 5, and 6 with this age model (Fig. 7c). Based on the interpretation of Haley et al. (2008), this would indicate increased Eurasian shelf brine input to the central Arctic during interglacial MIS 5, an unusual situation considering the absence of brines during the Holocene. This age model also removes any apparent glacial–interglacial variability from the ε_Nd data, leaving the low ε_Nd of the Holocene through MIS 3 as unique over the past 200 kyr. More concerningly, the Geibert et al. (2021) age model requires that planktonic foraminifera were completely absent from PS2185-3/6 during MIS 5 except for a small abundance (35 planktonic foraminifera per gram sediment) in MIS 5a (Fig. 7e). Both the general absence of MIS 5 planktonic foraminifera and the magnitude of this MIS 5a abundance “peak” sharply contrast with their prevalence during the Holocene and MIS 3 (>1000 planktonic foraminifera per gram sediment), which is independently constrained in time by ¹⁴C dating (Nørgaard-Pedersen et al., 1998). Considering the well-established link between interglacial periods and carbonate microfossil-rich sedimentation in the Arctic Ocean (Fig. 2b; e.g., Poore et al., 1993; Spielhagen et al., 1997, 2004; Polyak et al., 2009, 2013; Hanslik et al., 2013; Marzen et al., 2016; Vermassen et al., 2021), the Geibert et al. (2021) assignment of low-excess ²³⁰Th intervals in PS2185-3/6 with MIS 4 and 6 is unconvincing and disregards other important stratigraphic indicators that underpin the more widely accepted age model for this core (Spielhagen et al., 2004).

To illustrate this, Fig. 8 presents the PS2185-3/6 data from Fig. 7 on the original age model of Spielhagen et al. (2004) (which was also used by Haley et al., 2008). Notably, with this age model, maxima in planktonic foraminifera abundance align with the Holocene, MIS 3, and MIS 5 substages a, c, and e (compare Fig. 8e with Fig. 7e). Additionally, two intervals of higher ε_Nd correspond to MIS 4 and 6 on this age model (Fig. 8c), suggesting a greater contribution of Eurasian shelf-sourced weathering products to the intermediate to deep Arctic Ocean via brines during MIS 4 and 6 as originally proposed (Haley et al., 2008). With respect to excess ²³⁰Th, the Spielhagen et al. (2004) age model compresses (in time) the two intervals of low excess ²³⁰Th, with one approximately at the onset of and the second during MIS 4 (Fig. 8d). It is notable that our δ¹⁸O_sw reconstruction shows values indistinguishable from modern δ¹⁸O_sw during this interval (Fig. 8b), which presumably would rule out even a small freshwater contribution to the Arctic Ocean at this time. However, we caution that the composite Arctic δ¹⁸O_sw record does not have a sufficiently well-constrained chronology or high enough sedimentation rates to resolve submillennial-scale features like the two low-excess ²³⁰Th excursions in MIS 4 (Fig. 8d). Still, if the two low-excess ²³⁰Th excursions in MIS 4 are related to large freshwater inputs to the Arctic Ocean, these freshwater intervals must have dissipated sufficiently rapidly to allow for excess ²³⁰Th to accumulate within ∼1000 years of their occurrence and also to leave no apparent trace in Arctic δ¹⁸O_sw on multi-millennial timescales.

5.2 A mixing explanation for lower glacial Arctic δ¹⁸O_sw

Instead of freshwater accumulation, our preferred explanation is that the ¹⁸O-depleted water in the intermediate and deep Arctic Ocean during glacials reflected stronger mixing between the surface and deep Arctic Ocean. Today, relatively fresh waters in the Polar Surface Layer (PSL) possess a low δ¹⁸O_sw (of $\le -\mathrm{1}$ ‰) from the input of ¹⁸O-depleted water from rivers and across the Bering Strait (Fig. 1a; e.g., Bauch et al., 1995). These low-δ¹⁸O_sw waters are restricted to the surface by the strong density stratification between the PSL and underlying Atlantic Water that prevents downward mixing of the low-δ¹⁸O_sw waters to intermediate depths. Vertical mixing is weak in the Arctic Ocean today due to sea-ice cover and strong stratification (e.g., Rainville et al., 2011). While brine rejection during sea-ice formation on the shelves densifies low-δ¹⁸O_sw surface waters, the resultant waters are not sufficiently dense to penetrate the strong Arctic halocline and reach intermediate depths (Bauch and Bauch, 2001; Mackensen and Nam, 2014), except when the sinking waters entrain sediments and particulates (e.g., Rogge et al., 2023).

During the glacial periods, however, it is likely that the Arctic halocline was weaker due to reduced freshwater input to the Arctic Ocean. In both simple models and coupled ice–ocean circulation models, the strength of the Arctic halocline and the depth of the mixed layer vary as a function of net freshwater input to the Arctic Ocean (Nilsson and Walin, 2010; Pemberton and Nilsson, 2016). Today, the Arctic halocline is strong due to an excess of precipitation over evaporation, high rates of river discharge, and freshwater input across the Bering Strait. However, during glacial periods, precipitation over the Arctic Ocean was likely greatly reduced, rivers were ice dammed, and the Bering Strait was closed due to lower sea level. Together, these changes would have reduced freshwater input to the Arctic Ocean, weakened the halocline, and deepened the mixed layer relative to the Holocene (see discussion in Farmer et al., 2021).

Proxy evidence supporting this scenario over the last 50 kyr is shown in Fig. 9. Between approximately 29 and 14 ka, Arctic δ¹⁸O_sw-ivc and Δδ¹⁸O_sw were lower (Fig. 9a), nitrogen isotope ratios of organic matter bound within the planktonic foraminifera N. pachyderma (δ¹⁵N${}_{N.p.}$) were lower (Fig. 9b; Farmer et al., 2021, 2023), and BWTs reconstructed from Krithe $\mathrm{Mg}/\mathrm{Ca}$ were higher (Fig. 9c), but sediment leachate ε_Nd appeared unchanged (Fig. 9d; Haley et al., 2008) relative to the Holocene. The interpretations for these proxies coalesce on weaker Arctic stratification from 29 to 14 ka, although they are ambiguous about the role of brine formation. The lower δ¹⁵N${}_{N.p.}$ from 29 to 14 ka (Fig. 9b) indicates a decline in nitrate consumption by Arctic phytoplankton that most likely reflects increased nitrate supply to the Arctic mixed layer arising from weaker density stratification (Farmer et al., 2021). The increase in BWTs have been interpreted to indicate a deeper warm Atlantic inflow to the Arctic Ocean during MIS 2, which also has been explained as a response to a deepening of the mixed layer from weaker stratification (Cronin et al., 2012). This interpretation is additionally supported by a shift in ostracode faunal assemblages in our studied cores from Krithe- and Cytheropteron-dominated during the Holocene to Polycope-dominated from ∼29 to 14 ka (Poirier et al., 2012).

Figure 9 Low intermediate-to-deep-Arctic δ¹⁸O_sw during MIS 2 and records speaking to its cause. (a) Arctic δ¹⁸O_sw-ivc (blue line and shading) and Δδ¹⁸O_sw normalized to Site 1123 (orange line and diamonds). (b) N. pachyderma-bound δ¹⁵N from Lomonosov Ridge Site AOS94-B28 (Farmer et al., 2021; 2023). (c) Arctic BWT stack from Krithe $\mathrm{Mg}/\mathrm{Ca}$ (Cronin et al., 2012). (d) PS2185-3/6 metal oxide ε_Nd (Haley et al., 2008). Interpretive framework for each proxy is shown next to the record. The interval from 14 to 29 ka is marked by vertical dashed lines.

Weaker stratification during glacial periods could have contributed to lower intermediate-to-deep-Arctic δ¹⁸O_sw directly via enhanced vertical mixing or indirectly by allowing a greater fraction of PSL waters densified by brine rejection to penetrate to intermediate depths. Although of low resolution, the constant sediment leachate ε_Nd at Site PS2185-3/6 across the last ∼40 kyr (Fig. 9d) indicates that the weathering products of the Siberian Traps were not being delivered to the intermediate-depth Arctic Ocean over this time (Haley et al., 2008). At face value, this could be interpreted as a lack of contribution of brine-modified Eurasian shelf waters to the intermediate to deep Arctic. However, we note that the radiogenic Nd isotope signature (high ε_Nd) of weathering products is only observed in shelf areas adjacent to the Ob and Yenisey rivers today. In contrast, unradiogenic Nd isotope signatures (low ε_Nd) are observed associated with the Indigirka, Lena, and Khatanga rivers (Guo et al., 2004). Thus, while the Site PS2185-3/6 ε_Nd data rule out brine-derived surface waters carrying a radiogenic Nd isotope weathering signature into the intermediate Arctic during MIS 2, it is unclear whether the lower intermediate-to-deep-Arctic δ¹⁸O_sw reflected input of brine-derived surface waters formed in a shelf region with a different (and presumably lower) ε_Nd signature. This uncertainty also applies to previous glacial intervals where Arctic δ¹⁸O_sw and Nd isotopes are both available. For instance, during MIS 6 higher sediment leachate ε_Nd at Site PS2185-3/6 corresponds to low intermediate-to-deep-Arctic δ¹⁸O_sw; brines produced on the Eurasian shelves could explain both observations (Fig. 8b and c). By contrast, during MIS 4 higher sediment leachate ε_Nd at Site PS2185-3/6 aligns with δ¹⁸O_sw values that are indistinguishable from the modern ones (Fig. 8). In summary, Arctic Ocean intermediate waters have not exhibited a consistent pattern in δ¹⁸O_sw and ε_Nd over the last three glacial periods (MIS 2, 4, and 6).

5.3 Prospects for understanding Arctic δ¹⁸O_sw during glacial periods

That lower δ¹⁸O_sw in the intermediate to deep Arctic Ocean appears to be a recurrent feature of previous glacial periods (Fig. 4d) highlights the importance of understanding its cause. Yet, from the above discussion it is apparent that understanding the source(s) of lower-δ¹⁸O_sw waters in the intermediate to deep Arctic Ocean is hampered by a paucity of data even during the most recent glacial period. Here we outline several approaches that may help to shed light on past Arctic δ¹⁸O_sw variations.

5.3.1 Weathering input tracers

The Nd isotope record from PS2185-3/6 (Haley et al., 2008) is highly informative of the potential for brine contribution to the intermediate to deep Arctic through its ability to track weathering inputs from the Siberian Traps to the Eurasian shelf seas. However, this record is relatively low resolution (e.g., Figs. 8 and 9), and the potential for brine formation on different parts of the Arctic shelves to influence Nd (and other radiogenic) isotope compositions is poorly known. Expanding the spatiotemporal coverage of these weathering input proxies in the intermediate to deep Arctic Ocean will help elucidate where, when, and how shelf-derived materials were transported into the deep Arctic basins.

5.3.2 Upper-ocean δ¹⁸O_sw

If vertical mixing was strengthened during glacial periods in the Arctic, as we have suggested above, this may be traceable as a reduction in the water column (surface-to-deep) δ¹⁸O_sw gradient. To this end, generating surface ocean δ¹⁸O_sw reconstructions from paired planktonic foraminiferal $\mathrm{Mg}/\mathrm{Ca}$ and δ¹⁸O measurements in Arctic cores and their comparison with the deep Arctic records would be beneficial. However, we note that upper-ocean δ¹⁸O_sw records may have their own interpretive challenges related to the uncertainty in, and dynamic controls on, the depth habitat of the planktonic foraminifer N. pachyderma (e.g., Greco et al., 2019). Considering the steep δ¹⁸O_sw gradients in the upper Arctic water column (e.g., Bauch et al., 1995; Mackensen and Nam, 2014), constraining the paleo-depth habitat of N. pachyderma will be particularly important for interpreting δ¹⁸O_sw reconstructions from its geochemistry.

5.3.3 Reconstructions from outside of peak glacial intervals

The data gaps due to a lack of microfossils during some peak glacial intervals (e.g., Fig. 4) pose a significant challenge toward understanding the causes of lower δ¹⁸O_sw intervals. A potential workaround to this problem would be to apply additional proxy reconstructions that could speak to the cause of lower intermediate-to-deep-Arctic δ¹⁸O_sw over intervals with low δ¹⁸O_sw and high microfossil abundance. Examples include low Arctic δ¹⁸O_sw intervals during presumably “weaker” glacials (MIS 14), early in glacial inceptions or late in deglaciations (e.g., early MIS 11, 10, and 6), and “stadial”-type events (MIS 7b and 13b). Provided that these age assignments are correct, these intervals could be a useful target for better constraining the processes leading to Arctic δ¹⁸O_sw decline with multiproxy approaches.

6 Conclusions

We reconstructed the local δ¹⁸O_sw of the intermediate to deep Arctic Ocean over the last 600 kyr using measurements of oxygen isotopes in benthic foraminifera and $\mathrm{Mg}/\mathrm{Ca}$ in benthic ostracodes. Arctic δ¹⁸O_b shows a reduced glacial–interglacial amplitude compared to open-ocean δ¹⁸O_b records, with higher Arctic δ¹⁸O_b relative to the open ocean during interglacials and similar Arctic δ¹⁸O_b values to the open ocean during glacials. This reduced Arctic δ¹⁸O_b amplitude is not altogether surprising considering the unique glacial–interglacial temperature history of the intermediate to deep Arctic Ocean, with warmer glacials and cooler interglacials (Cronin et al., 2012, 2017). However, removing the temperature contribution to δ¹⁸O_b shows that the intermediate to deep Arctic Ocean also experienced relatively large (∼1 ‰) changes in local δ¹⁸O_sw across glacial cycles, with a tendency for lower δ¹⁸O_sw during glacials and higher δ¹⁸O_sw during interglacials.

In assessing potential explanations of lower intermediate-to-deep-Arctic δ¹⁸O_sw during glacials, we find that the magnitude of glacial δ¹⁸O_sw decline does not require dramatic inputs of freshwater to the Arctic Ocean. Moreover, we find that previous arguments for glacial freshwater events in the Arctic Ocean using ²³⁰Th excess are in part based on age models that are inconsistent with accepted patterns of Arctic sedimentology. A reinterpretation of these age models shows that intervals of low ²³⁰Th excess do not correspond with low intermediate-to-deep-Arctic δ¹⁸O_sw, raising uncertainty about the causal relationship between ¹⁸O-depleted freshwater input and low ²³⁰Th excess intervals. Instead, we find that intervals of lower intermediate-to-deep-Arctic δ¹⁸O_sw could be explained by reduced Arctic Ocean stratification during glacial periods, with stronger mixing and/or brine-enhanced transport of low-δ¹⁸O_sw surface waters to intermediate depths.

Ultimately, our results show that δ¹⁸O_b integrates complex histories of temperature and local δ¹⁸O_sw changes in the intermediate to deep Arctic Ocean. Over the glacial–interglacial cycles of the last 600 kyr, the δ¹⁸O influences of temperature, global ice volume, and local δ¹⁸O_sw on δ¹⁸O_b tend to act in opposing directions within the Arctic Ocean. During glaciations, larger continental ice volume worked to raise δ¹⁸O_b, but warming and lower local δ¹⁸O_sw in the intermediate to deep Arctic Ocean worked to lower Arctic δ¹⁸O_b. Consequently, the glacial–interglacial δ¹⁸O_b ranges in the Arctic Ocean are modest. Considering this modest δ¹⁸O_b range and the limited availability of benthic foraminifera in Arctic sediments, δ¹⁸O_b does not appear to hold the same promise as a chronostratigraphic tool in Arctic sediments as in the open ocean. Instead, paired δ¹⁸O_b and ostracode $\mathrm{Mg}/\mathrm{Ca}$ paleothermometry provides a new approach for investigating past Arctic δ¹⁸O_sw variations, which we suggest will serve as a sensitive integrator of oceanographic and cryosphere changes.