High-resolution seawater δ18O records, derived from coupled Mg/Ca and benthic δ18O analyses, can be used to evaluate how global ice volume changed during the mid-Pleistocene transition (MPT, ca. 1250–600 ka). However, such seawater δ18O records are also influenced by regional hydrographic signals (i.e., salinity) and changes in deep-ocean circulation across the MPT, making it difficult to isolate the timing and magnitude of the global ice volume change. To explore regional and global patterns in seawater δ18O records, we reconstruct seawater δ18O from coupled Mg/Ca and δ18O analyses of Uvigerina spp. at Ocean Drilling Program Site 1208 in the North Pacific Ocean. Comparison of individual seawater δ18O records suggests that deep-ocean circulation reorganized and the formation properties (i.e., salinity) of deep-ocean water masses changed at ca. 900 ka, likely related to the transition to marine-based ice sheets in Antarctica. We also find that an increase in ice volume likely accompanied the shift in glacial-interglacial periodicity observed in benthic carbonate δ18O across the MPT, with increases in ice volume observed during Marine Isotope Stages 22 and 16.

Even after decades of research, the mid-Pleistocene transition (MPT) remains a perplexing interval in Earth’s climate history (Shackleton and Opdyke, 1976; Raymo and Nisancioglu, 2003; Elderfield et al., 2012). A change from low-amplitude 41 k.y. glacial-interglacial cycles to larger-amplitude 100-k.y.-dominated glacial-interglacial cycles is observed in all foraminiferal δ18O records (as well as many other climate proxy records), although the timing of this transition varies by location. Some records suggest that the MPT was gradual (Ruddiman et al., 1989; Clark et al., 2006; Sosdian and Rosenthal, 2009), whereas others show a more abrupt change (Mudelsee and Schulz, 1997; Elderfield et al., 2012). In the spectral domain, a weak 100 k.y. signal is observed in records of sea-surface temperature and δ18O of benthic foraminifera (δ18Obenthic) as early as 1200 ka (Clark et al., 2006; McClymont et al., 2013). By Marine Isotope Stage (MIS) 22 (ca. 900 ka), the 100 k.y. cycle emerges as statistically significant, and at MIS 16 (ca. 600 ka), the transition is considered complete (Mudelsee and Schulz, 1997). Of course, all spectral methods are based on averaging a signal across a generally long time interval.

Another difficulty in evaluating the timing (and mechanisms) of the MPT is that δ18Obenthic values, our primary climate proxy, reflect not just the δ18O of seawater (δ18Osw, salinity and ice volume), but also bottom-water temperature and the circulation history of the water mass over a particular location (Clark et al., 2006; Ford et al., 2016). High-latitude conditions imprint temperature and salinity properties on water-mass parcels that sink and ultimately traverse the ocean’s depths (Talley, 2013). In polar regions, the salinity of sinking water masses is controlled by atmospheric processes (evaporation and precipitation), as well as by sea ice processes such as brine rejection (Brennan et al., 2013). To separate these effects, studies have used the Mg/Ca of benthic foraminifera to isolate the temperature contribution to δ18Obenthic (Sosdian and Rosenthal, 2009; Elderfield et al., 2012). However, changes in ocean circulation remain a confounding factor because associated water-mass salinity variations can also influence these records (Ford et al., 2016).

Here we present a new Mg/Ca-derived δ18Osw reconstruction from Ocean Drilling Program (ODP) Site 1208 in the North Pacific Ocean (Fig. 1). This location is at the end of the “global conveyor belt” and sits within Pacific Deep Water (PDW). PDW is an old, well-mixed water mass primarily composed of Southern Component Water (SCW) (Talley, 2013) and is also the most voluminous water mass globally (Worthington, 1981). We compare and combine our record with previously published Mg/Ca-derived δ18Osw reconstructions from ODP Site 1123 in the South Pacific (41.8°S, 171.5°W, 3290 m water depth; Elderfield et al., 2012) and Deep Sea Drilling Project (DSDP) Site 607 in the North Atlantic Ocean (41.0°N, 33.0°W, 3427 m water depth; Sosdian and Rosenthal, 2009; Ford et al., 2016) to generate a δ18Osw stack. We show that deep-ocean circulation reorganization and salinity played an increasingly important role in determining the density structure of the glacial deep ocean over the MPT. With this first δ18Osw stack across the MPT, we are able to provide an estimate of changes in global ice volume across this important climate transition.

ODP Site 1208 is a single-hole drilling location on Shatsky Rise (36.1°N, 158.2°E, 3346 m water depth). Sampling varies from 1.4 to 3.7 k.y. resolution. For oxygen isotope analyses, between two and four Uvigerina spp. shells were picked from the >250 μm size fraction of each sample and were analyzed on an Elementar Isoprime 100 stable isotope ratio mass spectrometer with dual inlet at Lamont-Doherty Earth Observatory (Palisades, New York, USA). Long-term precision is 0.06‰. Using the δ18Obenthic record, an age model was generated using the HMM-Stack Matlab code (Lin et al., 2014; Butcher et al., 2017), which aligns to the Prob-stack (a probabilistic Pliocene–Pleistocene stack of benthic δ18O using a profile hidden Markov model; Ahn et al., 2017). Prob-stack includes 180 δ18Obenthic globally distributed records and is an updated version of the LR04 stack and time scale (Lisiecki and Raymo, 2005). Prob-stack is similar in form to LR04 over the interval studied here.

For Mg/Ca analyses, between six and 12 Uvigerina spp. shells from each sample were analyzed on Thermo Scientific iCAP-Q inductively coupled plasma mass spectrometers (ICP-MS) at Lamont-Doherty Earth Observatory and at Rutgers University (New Jersey, USA) (see the GSA Data Repository1). Mg/Ca values were converted to bottom-water temperature (BWT) using a modified version of the Uvigerina spp. paleotemperature equation (Elderfield et al., 2010, 2012; Ford et al., 2016) to account for the reductive cleaning step used in this study.

The δ18Obenthic and BWTs were then used to derive δ18Osw using PSU Solver computational toolkit, which uses Monte Carlo simulations to constrain uncertainty (Thirumalai et al., 2016). Our δ18Osw stack was generated by (1) interpolating each δ18Osw record to an even 3 k.y. interval (across the 338–1450 ka interval where all three records have ≤ 3 k.y. sample resolution), (2) bootstrapping the records to estimate error, and (3) averaging the interpolated means and bootstrapped error estimates (see the Data Repository). The age models for DSDP Site 607 and ODP Site 1123 were not modified because their age models were tied to LR04 and there is little difference between LR04 and Prob-stack in this interval.

The δ18Obenthic record from ODP Site 1208 (Fig. 2B) is similar in trend and structure to Prob-stack (Fig. 2A). At Site 1208, there is a gradual increase in glacial δ18Obenthic values from the early Pleistocene at ca. 1400 to ca. 600 ka when there is a stepped increase in δ18Obenthic values at MIS 16. Glacial BWTs show near-freezing values over most the record in the North Pacific (Fig. 2C; see the Data Repository), similar to the South Pacific (Elderfield et al., 2012). However, at Site 1208, glacial MIS 22 shows a strong warming and ∼1‰ transient increase in glacial δ18Osw values relative to the early MPT (Fig. 2D). The large δ18Osw increase at MIS 22 is not repeated in MIS 20 or MIS 18. We do not see changes in preservation or contamination that might affect our records (see the Data Repository). Because this positive δ18Osw excursion is associated with warm BWTs as indicated by the Mg/Ca data (Figs. 2C and 2D), this implies that either a large ice-volume growth event happened at MIS 22 and/or the warm bottom-water mass was also characterized by higher salinity.

The ODP Site 1208 and DSDP Site 607 δ18Osw histories are not consistent with Elderfield et al.’s (2012) conclusion that an abrupt, large increase in Antarctic ice volume occurred at MIS 22 and recurred during subsequent glacial periods, particularly MIS 20 and 18 (Figs. 3A–3C). The observation of a transient increase in δ18Osw at Site 1208 during MIS 22 (this study) and a gradual increase in δ18Osw throughout the MPT at Site 607 (Sosdian and Rosenthal, 2009; Ford et al., 2016) suggests that hydrographic changes also occurred over the MPT. In other words, if global ice volume grew during MIS 22 and persisted during subsequent glacials as Elderfield et al. (2012) suggested, the stepped glacial increase in δ18Osw observed at Site 1123 should certainly be present at both Pacific locations. However, the δ18Osw increase at MIS 22 is absent during the subsequent glacial MIS 20 and MIS 18 at Site 1208. As an alternative interpretation, we suggest that the δ18Osw increase at Site 1123 at MIS 22 and subsequent glacials can be attributed to changes in deep-ocean circulation, salinity stratification, and an increase global ice volume.

Neodymium (Pena and Goldstein, 2014) and carbon records (Raymo et al., 1997; Venz et al., 1999; Lisiecki, 2014; Lear et al., 2016; Farmer et al., 2019) show that Atlantic Northern Component Water (NCW) decreased in spatial extent ca. 900 ka. During the early Pleistocene, NCW dominated the intermediate and deep Atlantic basin during both glacials and interglacials (Pena and Goldstein, 2014). At the “900 ka event,” during a time of anomalously low Southern Hemisphere insolation, Antarctic sea ice expanded and SCW flooded the South Atlantic at intermediate and deep ocean depths (Kemp et al., 2010; Pena and Goldstein, 2014). After MIS 22, carbon and neodymium isotopes show that NCW has been prominent during interglacials at intermediate and deep ocean depths, while SCW has dominated the deep ocean during glacials (Raymo et al., 1997; Pena and Goldstein, 2014).

What processes could have contributed to deep-ocean circulation and water-mass property changes at MIS 22? One possibility is that a transition from terrestrial to marine-based ice sheets in Antarctica (Raymo et al., 2006) around MIS 22 may have altered the distribution of deep-ocean water masses and water-mass properties in their formation areas. The transition from a largely terrestrial-based ice sheet to marine-based ice shelf would have also promoted the formation of coastal polynyas, like those in the Ross and Weddell Seas (Morales Maqueda et al., 2004), which increase salinity via brine rejection from sea-ice formation. Additionally, changes in sea ice in the Bering Sea (Stroynowski et al., 2017; Detlef et al., 2018; Kender et al., 2018) show northern high-latitude transformations that may have influenced NCW production in the Atlantic. Although sea-ice formation itself does not appreciably alter the δ18Osw value (Craig and Gordon, 1965), a salinity increase influences stratification and deep-ocean water mass formation.

These ocean-circulation and deep-ocean water mass formation changes influenced the regional 18Osw records at each site. For DSDP Site 607, a δ18Osw increase may be masked by corrosive, negative-δ18Osw SCW (LeGrande and Schmidt, 2006) penetrating further into the North Atlantic as shown from the carbonate saturation-state record which suggests that deep-ocean carbon storage increased during glacials after ca. 900 ka (Lear et al., 2016). On the Chatham Rise, ODP Site 1123 is sensitive to changes in water-mass boundary conditions (Shipboard Scientific Party, 1999). A switch to more NCW flowing over Site 1123 after the 900 ka event and later glacials may explain the persistent positive δ18Osw recorded there. In the North Pacific, Site 1208 has a large δ18Osw excursion due to warm BWTs during MIS 22 (and MIS 16), possibly related to North Pacific stratification (McClymont et al., 2008; Knudson and Ravelo, 2015). After the 900 ka event, δ18Osw is more positive during MIS 20 and 18 than prior to the transition, suggesting a modest increase in ice volume, though not to the same magnitude as originally interpreted by Elderfield et al. (2012).

Records at ODP Sites 1208 and 1123 and DSDP Site 607 each have different δ18Osw trends (Figs. 3B–3D), making it impossible to interpret them individually as recording purely a global ice volume signal; therefore, we stack these three high-resolution records and create a δ18Osw stack (Fig. 4C). We interpret the δ18Osw stack as primarily reflecting global ice volume; this assumes that the hydrographic changes are “averaged out” and that changes in atmospheric circulation (McClymont and Rosell-Melé, 2005) and ice sheet δ18O composition (Winnick and Caves, 2015) do not alter the overall δ18Osw composition in deep-ocean water mass formation regions. Although other δ18Osw records exist, they are discontinuous and do not have good interglacial-glacial resolution, so they are not included here. Similar to early δ18Obenthic stacks over the MPT (e.g., Williams et al., 1988; Raymo et al., 1990; Shackleton et al., 1995), our δ18Osw stack has its limitations (i.e., spatial coverage), and additional records will improve interpretation robustness (i.e., Lisiecki and Raymo, 2005; Ahn et al., 2017).

The δ18Osw stack shows gradual and abrupt increases in ice volume over the MPT. Through the early Pleistocene, there was a gradual increase in ice volume (Fig. 4C). At MIS 22, there was an abrupt increase in ice volume, likely in Antarctica due to the transition from land-based to marine-based ice sheets. At MIS 16, there was another large increase in ice volume, likely in the Northern Hemisphere as evidenced by the first appearance of ice-rafted debris from the Laurentide Ice Sheet via the Hudson Strait (Hodell et al., 2008).

If we use the relationship that a 0.1‰ change in δ18Osw equals a 10 m change in sea level (Fairbanks, 1989), this equates to an ∼45 m additional drop in sea level over the MPT (calculated as the sea-level change from MIS 36 to MIS 12; Fig. 4C). MISs 21–13 are “weak” interglacials consistent with the lower pCO2 values (Fig. 4B) measured in ice cores (Lüthi et al., 2008) and boron isotope proxies (Hönisch et al., 2009; Chalk et al., 2017). Sea-level estimates from the δ18Osw stack (Fig. 4C), the Mediterranean sea-level curve of Rohling et al. (2014), and a model deconvolution of δ18Obenthic records (De Boer et al., 2010) (Fig. 4D) have similar amplitude and basic structure over the MPT, although the timing and intensity of sea-level change differs between records. For instance, MIS 16 does not appear as a large glacial in the Mediterranean sea-level curve but does in our δ18Osw stack and the δ18Obenthic deconvolution. These differences are probably largely attributable to the different assumptions entailed by these various methods of reconstructing past sea level.

The 100 k.y. cycle gains statistical significance at MIS 22 and dominates the spectral signal in a diverse set of records by MIS 16 (Clark et al., 2006; McClymont et al., 2013); both of these intervals are associated with large changes in our δ18Osw stack and thus ice volume. Wavelet analysis of the δ18Osw stack shows a similar pattern (see the Data Repository). Mechanistically, what caused the sensitivity of the climate system to the 100 k.y. eccentricity cycle to change? It was likely some combination of atmospheric CO2 drawdown and change in ice sheet dynamics (Raymo et al., 1997; Tabor and Poulsen, 2016; Farmer et al., 2019; Willeit et al., 2019). At MIS 22, ocean carbon storage increased (Lear et al., 2016; Farmer et al., 2019) and deep-ocean circulation reorganized (Pena and Goldstein, 2014). Additionally, the modest increase in Antarctic ice volume shown in our δ18Osw stack and the hypothesized transition to fully marine-based ice sheet margins (Raymo et al., 2006) may have initiated sensitivity to the 100 k.y. cycle. The growth of large Northern Hemisphere ice sheets at MIS 16 (Hodell et al., 2008) may have been supported by regolith removal (Clark and Pollard, 1998) and strengthened the 100 k.y. cycle.

This work greatly benefited from discussions with K. Thirumalai, Y. Rosenthal, S. Crowhurst, J. Farmer, L. Haynes, D. Hodell, B. Hönisch, S. Woodard, and three anonymous reviewers. We thank E. Arnold, P. Rumford, and the International Ocean Discovery Program core repository for providing samples. We are grateful to C. Chang, A. Dial, K. Esswein, M. Greaves, W. Huang, B. Lindsey, A. Piotrowski, and S. Severmann for laboratory support. This research used the Queen Mary University of London (QMUL) Apocrita high-performance computing facility, supported by the QMUL IT Services Research group. This work was funded by the U.S. National Science Foundation (OCE-1436014; HLF and MER) with support from the Natural Environment Research Council (NE/N015045/1; HLF), Lamont Climate Center (HLF), and the Vetlesen Foundation (MER).

1GSA Data Repository item 2020037, appendix and data set, is available online at http://www.geosociety.org/datarepository/2020/, or on request from [email protected].