Evolution of mean ocean temperature in Marine Isotope Stage 4

Deglaciations are characterized by relatively fast and near-synchronous changes in ice sheet volume, ocean temperature, and atmospheric greenhouse gas concentrations, but glacial inception occurs more gradually. Understanding the evolution of ice sheet, ocean, and atmosphere conditions from interglacial to glacial maximum provides insight into the interplay of these components of the climate system. Using noble gas measurements in ancient ice samples, we reconstruct mean ocean temperature (MOT) from 74 to 59.7 ka, covering the Marine Isotope Stage (MIS) 5a–4 boundary, MIS 4, and part of the MIS 4–3 transition. Comparing this MOT reconstruction to previously published MOT reconstructions from the last and penultimate deglaciation, we find that the majority of the last interglacial–glacial ocean cooling must have occurred within MIS 5. MOT reached equally cold conditions in MIS 4 as in MIS 2 (−2.7 ± 0.3 ^∘C relative to the Holocene, −0.1 ± 0.3 ^∘C relative to MIS 2). Using a carbon cycle model to quantify the CO₂ solubility pump, we show that ocean cooling can explain most of the CO₂ drawdown (32 ± 4 of 40 ppm) across MIS 5. Comparing MOT to contemporaneous records of benthic δ¹⁸O, we find that ocean cooling can also explain the majority of the δ¹⁸O increase across MIS 5 (0.7 ‰ of 1.3 ‰). The timing of ocean warming and cooling in the record and the comparison to coeval Antarctic isotope data suggest an intimate link between ocean heat content, Southern Hemisphere high-latitude climate, and ocean circulation on orbital and millennial timescales.

1 Introduction

The classical view of Pleistocene glacial cycles is a slow buildup of ice sheets followed by rapid disintegration (Abe-ouchi et al., 2013; Emiliani, 1955; Hays et al., 1976; Imbrie et al., 1993). However, the glacial inception – the transition from interglacial to glacial maximum – also involves global cooling, large-scale changes in ocean circulation, and carbon cycle reorganization that may not coincide with the gradual pacing of ice sheet growth. The last glacial inception was punctuated by a rapid global cooling at 70 ka at the Marine Isotope Stage (MIS) 5a–4 boundary (Lisiecki and Raymo, 2005). During this period, nearly half of the interglacial–glacial drawdown of atmospheric CO₂ occurred over roughly 4 kyr (Ahn and Brook, 2008). This transition also brought extensive global cooling, buildup of polar ice sheets, and changes in deep ocean circulation (Adkins, 2013; Bereiter et al., 2012; Cutler et al., 2003; Yu et al., 2016). The mechanisms behind these rapid changes are not yet fully understood.

Multiple lines of oceanographic evidence (Adkins, 2013; Piotrowski, 2005; Thornalley et al., 2013; Yu et al., 2016) suggest that the MIS 5a–4 boundary marks the transition from the interglacial to glacial mode of ocean circulation. MIS 4 (like MIS 2) is characterized by cold conditions in both hemispheres and by the near absence of millennial-scale variability (Fig. 1). Sea surface temperature records for MIS 4 and MIS 2 (Kohfeld and Chase, 2017; Snyder, 2016) indicate that these two intervals were comparably cold on the global scale, although the spatial distribution of temperature may have differed (Kohfeld and Chase, 2017). While similarities exist between the MIS 4 and MIS 2 intervals, there are notable differences. Atmospheric CO₂ was ∼20 ppm lower in MIS 2 compared with MIS 4. Northern Hemisphere ice sheets and total ice volume were not as extensive as they were during MIS 2 (Cutler et al., 2003), but MIS 4 conditions included a greater extent of many glacier systems across the globe (Doughty et al., 2021; Schaefer et al., 2015). Understanding how and why conditions in MIS 4 and MIS 2 differed provides important context for the evolution of climate conditions during glacial inception.

Figure 1 Climate records of the last glacial cycle. Panel (a) displays the mean ocean temperature (MOT) anomalies relative to the early Holocene (11–10 ka; Baggenstos et al., 2019; Bereiter et al., 2018a; Shackleton et al., 2019, 2020; this study). Shading shows the 1σ confidence envelope of MOT data from a spline with a 2500-year cutoff period and bootstrapping. For 0–25 ka, the spline includes all published records from this interval, and it ends at 25 ka to exclude EDC MOT data within the bubble–clathrate transition zone. Also shown are (b) the global benthic δ¹⁸O stack (Lisiecki and Stern, 2016), (c) EDC δ²H (Jouzel et al., 2007), (d) the global average surface temperature anomaly from present (Snyder, 2016), (e) the CO₂ composite record (Bereiter et al., 2015), (f) relative (light blue) (Grant et al., 2012) and eustatic (royal blue) (Lambeck et al., 2014) sea level, (g) NGRIP δ¹⁸O (Andersen et al., 2004), and (h) summer solstice insolation at 65^∘ N (Berger and Loutre, 1991). Dotted lines show boundaries between Marine Isotope Stages (MIS) from Lisiecki and Raymo (2005). Gray shading shows the interval of the mean ocean temperature (MOT) record presented in this study. Black bars at the top of the figure show the intervals used to define MIS 4 (this study) and MIS 2 (Bereiter et al., 2018a) MOT.

One powerful indicator of global climate is the mean ocean temperature (MOT), which can be reconstructed using atmospheric noble gas ratios in ice-core-trapped air (Headly and Severinghaus, 2007). The total inventory of krypton and xenon in the ocean–atmosphere system is fixed, and the portion of the total that is dissolved in the global oceans depends on the MOT, as solubility of these heavy noble gases is strongly temperature dependent (Ritz et al., 2011). Ice-core-trapped air's $\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$, $\mathrm{Xe}/{\mathrm{N}}_{\mathrm{2}}$, and $\mathrm{Xe}/\mathrm{Kr}$ reflect the fraction of the noble gas inventory not dissolved in the ocean, which allows MOT at that time to be reconstructed. High-resolution reconstructions of MOT have been limited to the last two glacial terminations (Baggenstos et al., 2019; Bereiter et al., 2018a; Shackleton et al., 2019, 2020) but have provided unique insight into the interplay of key climate variables. In addition to the long-term warming across these deglaciations, millennial-scale variations in MOT are observed, which are also seen in Antarctic isotope records (Masson-Delmotte et al., 2010), and correspond to changes in Atlantic Meridional Overturning Circulation (AMOC) (Mcmanus et al., 2004). These deglacial features of MOT suggest an intriguing link between ocean circulation and ocean heat content. However, it is unclear if this link is unique to terminations or also applies to Dansgaard–Oeschger (DO) events (Dansgaard et al., 1982), which are millennial-scale climate oscillations that are thought to be linked to AMOC variability within glacial intervals (Lynch-Stieglitz, 2017; Stocker and Johnsen, 2003).

Here, we reconstruct MOT from 74 to 59.7 ka, covering the MIS 5a–4 transition, MIS 4, and part of the MIS 4–3 transition. The new record serves several purposes. First, it allows for a direct MIS 2–MIS 4 comparison, to assess their relative climate and ocean states. Second, comparison of MOT to benthic δ¹⁸O changes from the onset of the last interglacial (MIS 5e) to MIS 4 and MIS 2 provides insight into the temporal evolution of ocean temperature and ice volume changes over the glacial cycle. Third, it allows us to test if the link between changes in ocean circulation and heat content exists during DO event 19 (DO19) at 72.1 ka. Last, using a simple carbon cycle model (Bauska et al., 2016), we estimate the contribution of whole-ocean cooling to the decrease in atmospheric CO₂ across MIS 5 and the MIS 5a–4 boundary due to the solubility pump.

2 Methods

2.1 Site description and ice core measurements

Ice samples were obtained by drilling a shallow (20 m, 0.24 m diameter) ice core at Taylor Glacier, Antarctica, a blue ice area located in the McMurdo Dry Valleys (Baggenstos et al., 2017). The core contains ice spanning gas ages from ∼58 ka near the surface to 74 ka at 20 m depth (Menking et al., 2019). We excluded samples above 4 m depth to avoid alteration/contamination due to near-surface thermal fractures (Baggenstos et al., 2017). A total of 56 samples (including 11 replicate samples from identical depths) from Taylor Glacier were measured, with an average sample weight of 806 g and a mean record temporal resolution of 330 years. In addition, four WAIS (West Antarctic Ice Sheet) Divide samples from late MIS 4 (∼66–64 ka) were measured to replicate the Taylor Glacier results using samples from a different ice core. All ice core samples were analyzed for $\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$, $\mathrm{Xe}/{\mathrm{N}}_{\mathrm{2}}$, and $\mathrm{Xe}/\mathrm{Kr}$ using the method described by Bereiter et al. (2018b). The average of the three noble gas ratios was used to determine the final MOT following procedures in Shackleton et al. (2019). For brevity, we will refer to the MOT reconstructed from measured noble gases as “MOT data” in this work.

2.2 Taylor Glacier age model

We apply the ice core age model of Menking et al. (2019) with slight modifications for the MOT reconstruction. The age model was developed by matching measured variations in CH₄ and δ¹⁸O_atm in the Taylor Glacier ice core to deep ice core records on the Antarctic Ice Core Chronology 2012 (AICC2012) timescale (Veres et al., 2013). Tie points were manually selected, and noble gas sample ages were determined from linear interpolation between tie points. For this study, we selected tie points from the higher-resolution North Greenland Ice Core Project (NGRIP; rather than EPICA Dronning Maud Land – EDML) CH₄ record on AICC2012, as well as three additional tie points from the EDML CO₂ record, also on AICC2012 (Table 1). Tie point uncertainties are reported relative to AICC2012 and do not include age uncertainty of the AICC2012 chronology itself. Tie point uncertainties have a minimal impact on the interpretation of the record.

Table 1 Tie points used in this study. Taylor Glacier CH₄, δ¹⁸O_atm, and CO₂ measurements are tied to preexisting records of CH₄ (Baumgartner et al., 2014), δ¹⁸O_atm (Capron et al., 2010), and CO₂ (Bereiter et al., 2012) from well-dated ice cores on the AICC2012 (Veres et al., 2013) chronology.

2.3 Fractionation corrections of the noble gas ratios

The noble gas ratios measured in ice cores must be corrected for fractionation that occurs within the firn, which alters the noble gas ratios from their original atmospheric values (Headly and Severinghaus, 2007). We apply the correction approach of Shackleton et al. (2019), which uses a linear least-squares method to solve and correct for gravitational (Schwander, 1989) and thermal (Severinghaus et al., 1998) fractionations using measured isotope ratios of argon, nitrogen, and krypton. Further details on the applied fractionation corrections are included in Appendix A.

Corrections are more robust when calculating relative MOT change, rather than absolute MOT values, because errors in the fractionation corrections may produce a systematic offset in the corrected noble gas ratios, whereas the relative changes in these ratios are minimally influenced (Shackleton et al., 2020 and Appendix A). Therefore, we report MOT relative to Holocene MOT measured in the same ice core with the same applied method of fractionation correction. For the Taylor Glacier samples, we compare our data to five early Holocene (10.6 ka) replicate Taylor Glacier samples from Shackleton et al. (2020). WAIS Divide samples are reported relative to the average of Holocene samples from 11 to 10 ka (n=4) (Bereiter et al., 2018a). WAIS Divide (Bereiter et al., 2018a) and EPICA Dome C (EDC) (Baggenstos et al., 2019) records both suggest that the entire Holocene was a very stable interval for MOT (respective 1σ standard deviations of 0.2 and 0.1 ^∘C for all Holocene samples), so the particular choice of reference Holocene interval has minimal impact on the reported results.

2.4 MOT calculations from corrected noble gas ratios

We employ the box model of Bereiter et al. (2018a), which calculates the MOT anomaly relative to the modern ocean from the firn-corrected noble gas ratios. Parameterizations of the MOT box model applied in this study are detailed in Baggenstos et al. (2019). In addition, we use the recently published xenon and krypton solubilities from Jenkins et al. (2019). The box model requires input of sea level (Grant et al., 2012) to account for changes in the oceanic reservoir of xenon, krypton, and nitrogen that are unrelated to ocean temperature change. This includes changes in ocean volume, salinity, and sea surface pressure (Headly and Severinghaus, 2007). For Holocene MOT reference data and the MIS 2 data against which the record is compared, we use the sea level record of Lambeck et al. (2014) in the MOT box model. We also reevaluate the WAIS Divide Holocene and MIS 2 MOT record (Bereiter et al., 2018a), applying the same box model parameterizations as applied in this study (and in Shackleton et al., 2020) for a consistent comparison.

2.5 Carbon cycle model calculations of the solubility pump

To estimate the effect of a cooling ocean on the atmospheric CO₂ concentration via the solubility pump, we use a simple carbon cycle model (Bauska et al., 2016) to run a forward scenario of prescribed ocean temperature change from MOT constraints. The model consists of 14 boxes, representing the surface oceans (6 boxes), intermediate oceans (2 boxes), and deep oceans (3 boxes), a well-mixed atmosphere (1 box), and a terrestrial biosphere (2 boxes). The model simulates thermohaline circulation and mixing, air–sea gas exchange, export production, sediment burial/CaCO₃ compensation, and exchange of carbon between the atmosphere and terrestrial biosphere. More details on the model can be found in Appendix B.

2.6 Error analysis

The error on our MOT reconstruction is estimated by propagating all known uncertainties with 10 000 Monte Carlo simulations of our data. Sources of uncertainty include the analytical uncertainties for the noble gas ratios as well as the isotope ratios used for firn fractionation corrections. Additional uncertainties include the age uncertainty for the Taylor Glacier tie points (Table 1) and temporal and analytical uncertainties in the sea level curve. The calculated uncertainty for an individual sample is 0.2 ^∘C (1σ), and the pooled standard deviation of the 11 replicate samples is 0.3 ^∘C. Systematic errors (such as changes in ocean saturation state) may cause us to underestimate total uncertainty in our record.

To produce the splined MOT record, we use the 10 000 Monte Carlo iterations of the dataset and randomly sample the 56 individual MOT points via bootstrapping using the “randsample” MATLAB function with replacement. We then fit each of the 10 000 time series using a spline with a 2500-year cut-off period using the “csaps” MATLAB function and average the resulting splines to produce a final, smoothed version of our MOT record including uncertainty estimates.

3 Results

3.1 MIS 5a–4 boundary

During the rapid drawdown in atmospheric CO₂ (∼72–68 ka, Fig. 2), we observe mean ocean cooling in two phases, with an overall net cooling of 0.9 ± 0.3 ^∘C (1σ). In the first phase (72–70 ka), MOT decreased by 0.7 ± 0.3 ^∘C over roughly 2 kyr, coincident with Antarctic cooling and Greenland Interstadial 19 (GI19). In the second phase (70–68 ka), MOT stabilized and then decreased by a further 0.2 ± 0.3 ^∘C, reaching a minimum around 67.5 ka.

Figure 2 Mean ocean temperature (MOT) anomalies relative to the Holocene versus key climate variables. Panel (a) shows MOT data from Taylor Glacier (turquoise) and WAIS Divide (blue). Crosses indicate individual Taylor Glacier MOT data, and shading shows the 1σ confidence envelope of the Taylor Glacier data from a spline with a 2500-year cutoff period and bootstrapping. Solid blue points show WAIS Divide data corrected as described in Sect. 2 (with 1σ error bars), and open circles show the MOT results if the firn corrections detailed in Bereiter et al. (2018a) are applied. Also shown are (b) EDC δ²H (Jouzel et al., 2007) corrected for changes in seawater δ²H (see Appendix B), (c) CO₂ from EDML (diamonds) (Bereiter et al., 2012) and Taylor Glacier (points) (Menking et al., 2019) on AICC2012, (d) the relative sea level record (Grant et al., 2012), and (e) NGRIP δ¹⁸O_ice (Andersen et al., 2004) on AICC2012. Gray panels show warm Greenland intervals (interstadials), and white panels indicate cold Greenland intervals (stadials). The black bar at the top of the figure shows the time interval used to calculate Marine Isotope Stage 4 MOT.

3.2 Comparison of MOT in MIS 4 and MIS 2

Here, we do not use the intervals identified and defined by benthic δ¹⁸O to compare MOT in MIS 4 and MIS 2, as the alignment of ice core and sediment records is uncertain, particularly in MIS 4. Instead, we define MIS 4 as the interval in which CO₂ and Antarctic temperature remain low and stable (70.3–63.7 ka, or Greenland Stadial 19 and Greenland Interstadial 18). For Taylor Glacier samples, we compare MIS 4 samples to five replicate MOT samples from MIS 2 (19.9 ka). For WAIS Divide samples, we compare the measured MIS 4 samples to all available, previously published (Bereiter et al., 2018a) MOT data from MIS 2 (24–18 ka) with the fractionation corrections and MOT box model parameterizations used in this study applied. The difference in WAIS Divide MOT results for the full MIS 2 interval (n=11) versus 20–19 ka (n=4) is less than 0.01 ^∘C, so the difference in the selected intervals to define MIS 2 for each core should not affect the MIS 4–2 comparison.

The Taylor Glacier data show that MIS 4 MOT was statistically indistinguishable from MIS 2 (−0.1 ± 0.3 ^∘C relative to MIS 2 or −2.7 ± 0.3 ^∘C relative to the Holocene). The four WAIS Divide MIS 4 samples cover a narrower interval (66.2–63.9 ka) but are consistent with the Taylor Glacier results (+0.1 ± 0.5 ^∘C relative to MIS 2 or −2.6 ± 0.5 ^∘C relative to the Holocene). However, the WAIS data show more scatter. If we instead correct the WAIS Divide data for thermal fractionation using a firn model (Buizert et al., 2015), as in Bereiter et al. (2018a), and compare the results to MIS 2 data using this method of fractionation correction, we find that the WAIS Divide MIS 4 data are slightly less scattered (Fig. 2). With these corrections applied, the MIS 4 interval in WAIS Divide is 0.2 ^∘C warmer than in MIS 2 (−2.5 ± 0.3 ^∘C relative to the Holocene). While all results are indistinguishable within error, they emphasize the importance of future work developing further understanding of firn air processes and their influence on MOT results. A detailed discussion of the choice in fractionation corrections and their effect on calculated MOT is included in Appendix A.

3.3 MIS 4–3 transition

While our record may not capture the full transition into MIS 3, we find that there was a substantial increase in MOT towards the end of MIS 4. By 59.7 ka, MOT had reached levels comparable to the MOT peak at the end of MIS 5a at 72 ka (Fig. 2). Because our record does not contain a clear leveling of MOT, it is uncertain if or by how much MIS 3 MOT exceeded levels found at the end MIS 5a.

4 Discussion

4.1 Coevolution of MOT, benthic δ¹⁸O, CO₂, and Antarctic temperature during the last glacial inception

While the MOT proxy was developed over a decade ago (Headly and Severinghaus, 2007), only in the last few years have high-resolution MOT records become available (Baggenstos et al., 2019; Bereiter et al., 2018a; Shackleton et al., 2019, 2020). With the additional data from this study, we take the opportunity to review available MOT records and their relation to other key climate variables with a particular emphasis on the glacial inception (Fig. 3). While the available MOT records do not cover the entire last glacial cycle, we may still gain insight into climate evolution during (de)glaciations by comparing contemporaneous MOT, benthic δ¹⁸O, CO₂, and Antarctic δ²H.

Figure 3 ΔMOT plotted against coeval (a) δ¹⁸O_benth (Lisiecki and Stern, 2016), (b) atmospheric CO₂ (Bereiter et al., 2015), and (c) EDC δ²H (Jouzel et al., 2007) corrected for changes in seawater δ²H. The δ¹⁸O_benth, CO₂, and EDC δ²H records were linearly interpolated in order to plot them against contemporaneous MOT. Additionally, the EDC δ²H was smoothed using a Gaussian filter with a 500-year window to remove high-frequency variability. The gray arrow in panel (a) shows the ΔMOT–δ¹⁸O scaling from (Shackleton, 1974). The gray arrow in panel (b) shows the ΔMOT–CO₂ relationship for the solubility pump from the carbon cycle model. Filled symbols include data from the last interglacial through glacial maximum (129–18 ka) to highlight the glacial inception, while open symbols indicate data from the last and penultimate deglaciations and the Holocene. Diamonds indicate MOT data constructed from the EDC record (Baggenstos et al., 2019; Shackleton et al., 2020), circles show data from WAIS Divide (Bereiter et al., 2018a, and this study), and stars show MOT data from Taylor Glacier (Shackleton et al., 2019, 2020, and this study). The color of data indicates the age.

4.1.1 Evolving control of ocean temperature and ice sheet volume on benthic δ¹⁸O

The link between ocean temperature and benthic foraminiferal δ¹⁸O (Fig. 3a) has long been recognized (Emiliani, 1955; Shackleton, 1974). While MOT represents volume-averaged ocean temperature, the intermediate and deep ocean make up the majority of total ocean volume. The benthic δ¹⁸O record (Lisiecki and Stern, 2016) shown in Figs. 1 and 3 contains stacked records from intermediate and deep sites, and (when binned into ocean regions) covers approximately 70 % of the total ocean volume. Thus, changes in MOT should be largely reflected in temperature-driven changes in this δ¹⁸O record. The scaling between MOT and δ¹⁸O for ocean temperature change at 3.5 ^∘C (Holocene/modern MOT, or ΔMOT = 0) from Shackleton (1974) (0.26 ‰ ${}^{\circ }{\mathrm{C}}^{-\mathrm{1}}$) is denoted by the gray arrow in Fig. 3a. While the temperature dependence of δ¹⁸O from Shackleton (1974) is quadratic, it is effectively linear in the temperature range of the plotted MOT data, i.e., $\mathrm{d}{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}/\mathrm{d}T$ varies by less than 6 % within the ΔMOT range shown in Fig. 3.

The other primary control on benthic δ¹⁸O is ice volume. Considering the temporal evolution of δ¹⁸O and MOT, it is possible to gain insight into the relative controls on δ¹⁸O during the intervals where δ¹⁸O and MOT data are available. Applying the benthic δ¹⁸O temperature sensitivity from Shackleton (1974), we find that the ocean temperature anomaly during MIS 4 accounts for 0.7 ‰ of the 1.3 ‰ δ¹⁸O anomaly relative to Holocene/modern benthic δ¹⁸O, implying that the remaining 0.6 ‰ is due to enhanced ice sheet volume. For comparison, the MIS 2–Holocene benthic δ¹⁸O change is 1.7 ‰. Considering the MOT–δ¹⁸O relationship for late MIS 5a/MIS 4 (light green in Fig. 3), late MIS 3 (cyan), and MIS 2 (light blue), there is some variability in MOT within these intervals, but average MOT across the intervals remains essentially unchanged. However, there is a clear, long-term increase in δ¹⁸O across these intervals. The similarity in MOT between MIS 4 and MIS 2 suggests that the more positive benthic δ¹⁸O during the latter stage is caused by a greater global ice volume. Taken together, these observations are consistent with previous studies (Cutler et al., 2003; Shakun et al., 2015; Waelbroeck et al., 2002) suggesting that ocean cooling outpaced Northern Hemisphere ice sheet growth in the last glacial inception. This decoupling of ocean cooling and ice sheet growth may be an important clue for future investigation of the mechanism of glacial cycles.

4.1.2 Early role of ocean cooling in atmospheric CO₂ drawdown

Here, we discuss two separate drawdowns of CO₂ during the glacial inception, each of which was approximately 40 ppm. The first occurred from MIS 5e to MIS 5a, and the second took place from MIS 5a to MIS 4. Using a carbon cycle model (Bauska et al., 2016), we estimate that the 3.1 ± 0.4 ^∘C MOT decrease between the onset of MIS 5e (129 ka) and the end of MIS 5a (72 ka) accounted for 32 ± 4 of the ∼40 ppm CO₂ reduction that occurred across this interval. We emphasize that the available MOT data spans 9 kyr at the onset and 2 kyr at the end of the long (∼57 kyr) MIS 5 interval, so our insight into the role of the solubility pump on CO₂ variations within MIS 5 is limited. However, our MOT data suggest a dominant role of ocean cooling on the ∼40 ppm CO₂ drawdown that occurred across MIS 5e–5a. Much of this drawdown is focused on the MIS 5e–5d transition around 115 ka, making this period a high priority for future detailed ice core MOT reconstruction.

During the MIS 5a–4 transition, we estimate that the reconstructed MOT net decrease of 0.9 ^∘C led to a CO₂ drawdown of 9 ± 3 ppm by solubility, which is a relatively small but not insignificant fraction of the ∼40 ppm drawdown that occurred over the full interval (Fig. 4). A comparison of the Taylor Glacier records of MOT and CO₂ over the MIS 5a–4 transition (∼72–68 ka) suggests that while both MOT and CO₂ decreased during this transition, the overall trends appear distinct in their shapes. While the rate of CO₂ decrease was relatively constant over the full transition, MOT decreases more rapidly during the first half of the transition (−0.41 ± 0.09 ^∘C kyr⁻¹ over GI19, 72.1–70.3 ka) than in the second half (−0.19 ± 0.07 ^∘C kyr⁻¹, 70.3–67.5 ka). This duality in trends over the MIS 5a–4 transition has been observed in other proxy records (Barker and Diz, 2014) and may provide important insight into the evolving controls on atmospheric CO₂ during this interval. While ocean cooling may explain roughly one-third of the CO₂ drawdown in the first half of the transition, significant carbon cycle reorganization is required to explain the majority of the atmospheric CO₂ decrease, particularly in the second half of the transition.

Figure 4 Results from a carbon cycle model estimating the magnitude of CO₂ drawdown due to mean ocean temperature cooling. Inputs to the model sea surface temperature changes (MOST, orange trace, c) are scaled to ice core δ²H data corrected for seawater δ²H changes (blue trace, a, see Appendix B). The sea surface temperature changes are transmitted to the deep ocean via circulation and mixing in the model, causing the mean ocean temperature (MOT, orange trace, d) to evolve through time. The modeled MOT history agrees well with the existing (but limited) ice core MOT data (blue traces, d). Ocean salinity also evolves in the model and is scaled to sea level data (Appendix B). The evolution of CO₂ in the model (orange trace, b) is only due to changes in ocean solubility. The modeled CO₂ history is compared to ice core CO₂ records (red markers) (Bereiter et al., 2015). Model results within 120–74 ka should be interpreted with caution, as MOT data do not exist for validation.

Again, we emphasize that available MOT records do not cover the full glacial cycle, and substantial gaps in the data exist for MIS 5 and MIS 3. However, it is notable that the coevolution of CO₂ and MOT over the last glacial inception shows a strikingly similar trend to that of benthic δ¹⁸O and MOT (Fig. 3); ocean cooling appears to play a dominant role in the net CO₂ decrease across MIS 5, but the long-term trend of CO₂ drawdown across MIS 4–2 does not correspond with ocean cooling, as MOT had reached levels comparable to MIS 2 by MIS 4. We speculate that, within the last glacial inception, the MIS 5a–4 boundary marks a distinct decoupling in MOT and CO₂ trends. This is consistent with the hypothesis presented by Adkins (2013): that the MIS 5a–4 boundary marks a transition between interglacial and glacial modes of ocean circulation and shifts the controls on carbon uptake from primarily temperature-driven solubility to circulation-driven storage, for example, via reduced ventilation of abyssal waters that allows respired carbon to accumulate there.

Our data allow us to put new constraints on the role of the solubility pump in atmospheric CO₂ variations across the studied intervals. Out of the full ∼80 ppm CO₂ from MIS 5e to MIS 4, our modeling suggests that MOT changes can explain 41 ± 4 ppm. These estimates of the solubility pump agree well with the canonical 10 ppm ${}^{\circ }{\mathrm{C}}^{-\mathrm{1}}$, (Williams and Follows, 2011). However, our MOT data provide no information on the spatial distribution of ocean temperature change, which a recent study suggests plays an important role in modulating the strength of the solubility pump via changes in ocean saturation state (Khatiwala et al., 2019). The referenced study found that changes in air–sea disequilibrium between interglacial and glacial ocean conditions enhanced the solubility pump by ∼60 % during the last glacial maximum. If such disequilibrium effects are also relevant for the timescales and periods considered here, our solubility-driven estimates from the carbon cycle model simulations may be considered a lower bound. In particular, if the enhanced disequilibrium effect is linked to the onset of the glacial mode of ocean circulation at the MIS 5a–4 transition, the solubility pump may play a larger role there than suggested from our simplified carbon cycle modeling.

We note that changes in ocean saturation state may also influence the noble gases and, thus, MOT estimates. However, the more rapid equilibration timescales of the noble gases versus CO₂ (a few weeks versus roughly a year) means that they are quite insensitive to disequilibrium (Ritz et al., 2011). However, given the recent improvements in analytical precision of the MOT technique, glacial–interglacial changes in noble gas disequilibrium merit future investigation.

4.1.3 Strong correlation between MOT and Antarctic climate on orbital and millennial timescales

As highlighted in this and several other MOT studies (Bereiter et al., 2018a; Shackleton et al., 2019, 2020), one of the most striking features of MOT records is their strong correlation to Antarctic water isotope records (Fig. 3c). For the MOT data from this study, we find a lower correlation between MOT and EDC δ²H (n=56, r²=0.59) than between all available MOT records (n=243, r²=0.94). However, MOT and δ²H data for this interval cover a relatively narrow range compared with other records, resulting in a lower signal-to-noise ratio, and thus may explain the lower correlation. To test this hypothesis, we use the pooled standard deviation of replicate MOT samples (0.3 ^∘C) as a predictor of random noise in the MOT record to estimate the expected correlation between δ²H and MOT if we assume they are perfectly correlated (r²=1). Under these assumptions, we would predict r² values of 0.58 ± 0.09 and 0.93 ± 0.01 for the MIS 4 subsamples and all MOT samples, respectively, which is consistent with the observed values. It is remarkable that the MOT–δ²H scaling is similar on millennial and orbital timescales, given that climate dynamics on these two timescales are likely to be different. Multiple explanations may be given for the strong correlation.

If there is indeed a causal relationship between MOT and Antarctic temperature, causality could plausibly run in either direction. First, it has been suggested that Southern Hemisphere high-latitude temperature, for which Antarctic δ²H is a proxy, provides a control on MOT (Bereiter et al., 2018a). Given that a large fraction of the global ocean interior is ventilated in the Southern Ocean (Johnson, 2008), processes acting in the Southern Ocean around Antarctica are likely to be important in setting the MOT. The temperature of deep waters formed in the Southern Ocean, as well as the rate at which they form, is probably linked to Southern Hemisphere high-latitude climate, providing a pathway to control MOT variations (Bereiter et al., 2018a).

Second, it is possible that causality runs in the opposite direction, with MOT being a strong control on Antarctic δ²H. In their modeling study, Pedro et al. (2018) proposed a mechanism linking MOT to Antarctic temperature on millennial timescales, as part of their effort to provide a more thorough dynamical framework for the bipolar seesaw. Briefly, during weakened AMOC intervals, ocean warming centered in the intermediate-depth North Atlantic is spread throughout the ocean basins via Kelvin and Rossby waves, which cannot cross the Antarctic Circumpolar Current. The enhanced temperature gradient across the Antarctic Circumpolar Current drives poleward ocean and atmospheric eddy heat fluxes, which are amplified by sea ice reduction and the ice–albedo feedback. The net result is a strong warming of the Antarctic continent. In this view, it is feasible that the MOT controls Antarctic temperature, via variations in Southern Ocean poleward eddy heat transport and sea ice feedbacks.

Finally, MOT and Antarctic temperature need not be causally linked; the tight correlation between them may reflect a shared dependence on a third variable that is most likely AMOC variability. It is well established that Antarctic temperature responds to AMOC variations via the bipolar seesaw mechanism (Stocker and Johnsen, 2003). Likewise, AMOC variations and their associated changes in oceanic heat loss to the Arctic atmosphere have been shown to influence MOT in model simulations (Galbraith et al., 2016; Pedro et al., 2018). Thus, it is conceivable that both variables respond to AMOC variations without the necessity for a direct causal link between them.

Here, we remain agnostic as to which of these three explanations is the correct one. Such a determination would require detailed modeling studies that are beyond the scope of the present work. However, our record demonstrates that Antarctic temperature and MOT covary on millennial timescales during DO19, suggesting that their link is not unique to deglaciations and is a general feature of the climate system.

4.2 The cold and stable MIS 4 interval

Given that global ice volume was greater (Cutler et al., 2003) and CO₂ concentrations were lower in MIS 2 than in MIS 4, the comparably cold conditions during these intervals suggested in this study and by recent (Doughty et al., 2021) and previous (Kohfeld and Chase, 2017; Snyder, 2016) work is somewhat puzzling. All else being equal, the 20 ppm higher atmospheric CO₂ in MIS 4 would lead to 0.5 ^∘C warmer global average surface temperatures than in MIS 2, assuming a climate sensitivity of 3.5 ^∘C. Other forcing may be required to resolve this conundrum. While changes in planetary albedo are often assumed to scale with ice volume, this may not be appropriate when comparing MIS 4 and MIS 2. Doughty et al. (2021) suggest that while ice volume was lower, glacier extent may have been greater in MIS 4 than in MIS 2. This may have led to a higher planetary albedo in MIS 4 than in MIS 2. However, the authors suggest that the greater glacial extent was due to the cold conditions in MIS 4, so there is some circularity to this argument.

Another notable difference between MIS 4 and MIS 2 exists for the Saharan region. Proxy records of the Sahara suggest that MIS 4 was a uniquely arid interval within the last glacial cycle (Castañeda et al., 2009; Skonieczny et al., 2019; Tierney et al., 2017), while regions of the eastern Sahara and Sinai Desert during MIS 2 may have been wetter than today (Hamdan and Brook, 2015), suggesting greener Saharan conditions in MIS 2 compared with MIS 4. Climate simulations of a greener Sahara suggest globally warmer temperatures due to lower albedo and higher atmospheric moisture (Tabor et al., 2020). We speculate that more arid Saharan conditions in MIS 4 may have contributed additional cooling of MIS 4 relative to MIS 2 and, in part, countered the warming effect of higher CO₂ during this interval. Of course, this new MOT record cannot shed light on the conditions of the low-latitude hydrosphere. However, it adds to a growing body of work suggesting that, despite the smaller ice volume and higher atmospheric CO₂, MIS 4 was comparably cold to MIS 2. We believe that this conundrum merits further investigation and may be a valuable target for forthcoming climate modeling efforts.

In addition to their cold temperatures, MIS 4 and MIS 2 also share an absence of millennial-scale variability. Theories explaining the apparent lack of bipolar seesaw behavior during very cold periods (such as MIS 4 and MIS 2) have invoked mechanisms related to thresholds in ice volume (McManus et al., 1999) and Southern Ocean temperature (Buizert and Schmittner, 2015). While ice volume during MIS 2 exceeds that of MIS 4, MOT during MIS 2 and MIS 4 indicate equally cold ocean conditions. This supports the idea that thresholds in ocean temperature, rather than global ice volume, may determine the presence or absence of millennial-scale variability within a glacial.

5 Conclusions and future outlook

Our record adds to the growing number of MOT reconstructions and provides unique insight into the climate conditions of MIS 4 within the context of the last glacial inception. The MOT record shows comparably cold and stable conditions in MIS 4 and MIS 2. As demonstrated in previous studies, MOT and Antarctic isotope records are remarkably correlated and covary on millennial timescales during DO19, providing the first evidence of the connection between MOT, Antarctic temperature, and inferred AMOC variability outside of deglaciations. Comparisons of coeval benthic δ¹⁸O, CO₂, and MOT show that while MOT reaches a minimum in the last glacial cycle by MIS 4, CO₂ and δ¹⁸O do not achieve their glacial extrema until MIS 2; ocean cooling outpaced ice sheet growth and atmospheric CO₂ drawdown in the last glacial inception.

Using a carbon cycle model and ocean temperature constraints provided from our record, we demonstrate that ocean cooling played a major role in the early (MIS 5) stages of atmospheric CO₂ drawdown in the glacial inception (32 of ∼40 ppm), a moderate role in the first half of the MIS 5a–4 transition (7 of ∼20 ppm), a minor role in the second half of the 5a–4 transition (2 of ∼20 ppm), and no measurable role in the ∼20 ppm decrease between MIS 4 and MIS 2. This suggests an evolving control on atmospheric CO₂ during the glacial inception in which the solubility pump initially dominates but plays a progressively smaller role as the glacial inception progresses. MOT reconstruction of the entire MIS 5 interval would provide valuable insight into this apparent trend.

Studies comparing the CO₂ (Bereiter et al., 2012) and Atlantic Western Boundary Undercurrent (Thornalley et al., 2013) response to millennial-scale variability during MIS 5 and MIS 3 suggest that the changes in ocean circulation at the MIS 5–4 boundary altered the nature of abrupt climate change between these two intervals. Comparison of the MOT response to DO cycles within MIS 5 and MIS 3 may prove useful in understanding the nature of this transition.

This study demonstrates that it is possible to capture MOT changes during the larger of the millennial-scale DO events using the noble gas ratio technique. However, comparison of the MOT records between smaller DO events will push the current analytical limits of this method. While improvements in analytical precision will benefit future studies, an improved understanding of gas fractionation processes within the ice and firn, and the mechanisms of air–sea gas exchange will be critical to accurate interpretation of ice core MOT records.

Appendix A: Comparison of MOT results using different methods of fractionation corrections and between $\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$, $\mathrm{Xe}/{\mathrm{N}}_{\mathrm{2}}$, and $\mathrm{Xe}/\mathrm{Kr}$

Gases are trapped in bubbles in ice during the process of firnification, as snow compacts and densifies into firn and eventually glacial ice. This process is gradual, occurring on timescales on the order of hundreds to thousands of years. During this time, the low permeability of the firn restricts bulk air motion but allows for air in the open pores to exchange with the overlying atmosphere and throughout the firn column primarily through molecular diffusion. This mechanism of air transport allows for processes such as gravitational settling (Schwander, 1989), thermal diffusion (Severinghaus et al., 1998), and kinetic fractionation (Birner et al., 2018; Buizert and Severinghaus, 2016; Kawamura et al., 2013) to alter $\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$, $\mathrm{Xe}/{\mathrm{N}}_{\mathrm{2}}$, and $\mathrm{Xe}/\mathrm{Kr}$ from their atmospheric compositions before bubble close-off. Correction of $\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$, $\mathrm{Xe}/{\mathrm{N}}_{\mathrm{2}}$, and $\mathrm{Xe}/\mathrm{Kr}$ for these processes may be done with output from a firn model and/or from measurements of isotope ratios of inert gases (such as argon, nitrogen, krypton, and xenon), which are also influenced by these processes but are unchanging in the atmosphere. Argon isotope ratios are a slight exception, due to the gradual degassing of ⁴⁰Ar from the solid earth (Bender et al., 2008). With the known rate of change in atmospheric ⁴⁰Ar and age of the samples, a small (<0.005 ‰) correction is applied to measured ${\mathit{\delta }}^{\mathrm{40}/\mathrm{38}}$Ar and ${\mathit{\delta }}^{\mathrm{40}/\mathrm{36}}$Ar. In the case of this study, we measure isotope ratios of argon (${\mathit{\delta }}^{\mathrm{40}/\mathrm{36}}$Ar, ${\mathit{\delta }}^{\mathrm{40}/\mathrm{38}}$Ar), nitrogen (δ¹⁵N–N₂), and krypton (${\mathit{\delta }}^{\mathrm{86}/\mathrm{82}}$Kr, ${\mathit{\delta }}^{\mathrm{86}/\mathrm{83}}$Kr, ${\mathit{\delta }}^{\mathrm{86}/\mathrm{84}}$Kr) for the applied fractionation corrections.

The approach of fractionation correction in this study differs slightly from that of Bereiter et al. (2018a), which uses previously published firn model output from Buizert et al. (2015) to correct $\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$, $\mathrm{Xe}/{\mathrm{N}}_{\mathrm{2}}$, and $\mathrm{Xe}/\mathrm{Kr}$ for thermal fractionation. To correct $\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$, $\mathrm{Xe}/{\mathrm{N}}_{\mathrm{2}}$, and $\mathrm{Xe}/\mathrm{Kr}$ for gravitational fractionation, the Bereiter et al. (2018a) study uses measured ${\mathit{\delta }}^{\mathrm{40}/\mathrm{36}}$Ar, which is first corrected for thermal fractionation using this same firn model output. While in the case of the WAIS Divide ice core, previously published firn model output was readily available for this purpose, no such model output exists for Taylor Glacier, which is why this method of fractionation correction was not considered in this study.

In this study, as in Shackleton et al. (2019), the corrections for gravitational and thermal fractionation are done with measured inert gas isotope ratios. The reason for this choice over other considered methods (see below) is that (1) it gives the best agreement in Taylor Glacier MOT results between replicate samples for the Holocene and MIS 2 (Shackleton et al., 2020), and (2) it gives the best results for calculated MOT in firn air and surface ice samples from a wide range of site conditions (Shackleton, 2019). The average magnitude of the fractionation corrections for $\mathit{\delta }\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$, $\mathit{\delta }\mathrm{Xe}/{\mathrm{N}}_{\mathrm{2}}$, and $\mathit{\delta }\mathrm{Xe}/\mathrm{Kr}$ were 5.0 ‰, 8.9 ‰, and 3.9 ‰, respectively, for the Taylor Glacier samples in this study. The fractionation-corrected, atmospheric values found for $\mathit{\delta }\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$, $\mathit{\delta }\mathrm{Xe}/{\mathrm{N}}_{\mathrm{2}}$, and $\mathit{\delta }\mathrm{Xe}/\mathrm{Kr}$ were −1.4 ‰, −3.6 ‰, and −2.1 ‰, respectively. While the fractionation corrections are large relative to the atmospheric signal, they are dominated by gravitational fractionation, which is physically well understood and constrained by the measured isotope ratios. For context, the average magnitude of fractionation correction for the WAIS Divide samples in this study are 22.7 ‰, 42.5 ‰, and 19.3 ‰ for $\mathit{\delta }\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$, $\mathit{\delta }\mathrm{Xe}/{\mathrm{N}}_{\mathrm{2}}$, and $\mathit{\delta }\mathrm{Xe}/\mathrm{Kr}$, respectively. The large difference in the magnitude of the applied fractionation correction between ice cores but good agreement in corrected noble gas ratios and resulting MOT gives us some confidence in our ability to apply these corrections.

To assess how differing methods of fractionation correction may impact the MOT results, we apply multiple corrections following Shackleton et al. (2020). Figure A1 shows the MIS 4 Taylor Glacier MOT data (average of the $\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$, $\mathrm{Xe}/{\mathrm{N}}_{\mathrm{2}}$, and $\mathrm{Xe}/\mathrm{Kr}$ results) with fractionation corrections for gravitational fractionation, gravitational and thermal fractionation, and gravitational and kinetic fractionation. For a detailed explanation of these methods of fractionation correction, see the supporting information of Shackleton et al. (2020). Results are compared between these differing methods of fractionation correction when (i) they are not calculated relative to a reference interval, (ii) they are calculated relative to Holocene MOT, and (iii) they are calculated relative to MIS 2 MOT. As previously shown, MOT reported relative to a reference interval in the same core is more robust to the method of fractionation correction than when no reference interval is used (Shackleton et al., 2020). However, even when the MIS 4 MOT data are referenced to Holocene MOT, there appears to be a small but systematic offset in the MOT results using different methods of fractionation correction. If the MIS 4 data are calculated relative to MIS 2 data from the same ice core, the offset is reduced.

Figure A1 Comparison of mean ocean temperature (MOT) anomalies between three methods of fractionation correction. Results are for the average of the three MOT proxies ($\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$, $\mathrm{Xe}/{\mathrm{N}}_{\mathrm{2}}$, and $\mathrm{Xe}/\mathrm{Kr}$). Panel (a) shows results if the noble gas ratios are corrected for fractionation and no reference interval is used. Panel (b) shows the results when the noble gas ratios are reported relative to Holocene data, using the same method of fractionation correction. Panel (c) is the same as panel (b) but MOT is reported relative to MIS 2. Individual MOT data from Taylor Glacier are shown as filled symbols and the 1σ confidence envelope from a spline with a 2500-year cutoff period, and bootstrapping is shown using shading. WAIS Divide data are shown as open symbols.

Comparison of the MOT results from $\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$, $\mathrm{Xe}/{\mathrm{N}}_{\mathrm{2}}$, and $\mathrm{Xe}/\mathrm{Kr}$ when normalized to Holocene versus MIS 2 MOT data show a similar phenomenon to the observed offset in results between differing fractionation corrections (Fig. A2). The offset between the MOT results for the three noble gas ratios is present, regardless of the fractionation correction applied; if the MOT results from the individual ratios are reported relative to those from the Holocene, there is a small offset between the MOT results from $\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$, $\mathrm{Xe}/{\mathrm{N}}_{\mathrm{2}}$, and $\mathrm{Xe}/\mathrm{Kr}$. However, if the MOT difference is calculated between MIS 2 and our MIS 4 data, the offset diminishes.

Figure A2 Mean ocean temperature (MOT) anomalies calculated from $\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$, $\mathrm{Xe}/{\mathrm{N}}_{\mathrm{2}}$, and $\mathrm{Xe}/\mathrm{Kr}$. Noble gas ratios are corrected for gravitational and thermal fractionation (Shackleton et al., 2019, and this study). Panel (a) shows the calculated MOT anomaly when no reference interval is used. Panel (b) shows the MOT anomaly relative to Holocene MOT results. Panel (c) shows the MOT anomaly relative to MIS 2. Individual MOT data from Taylor Glacier are shown using filled symbols and the 1σ confidence envelope from a spline with a 2500-year cutoff period, and bootstrapping is shown using shading. WAIS Divide data are shown using open symbols.

The observed patterns are consistent with systematic uncertainties in the fractionation corrections applied to the noble gas ratios. Differences in site conditions, such as temperature, accumulation, and atmospheric circulation, can lead to differences in firn column height, temperature profile, and dynamics of gas transport and mixing within the firn. These have implications for gravitational, thermal, and kinetic fractionation of $\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$, $\mathrm{Xe}/{\mathrm{N}}_{\mathrm{2}}$, and $\mathrm{Xe}/\mathrm{Kr}$. While the fractionation corrections should account and correct for these changes, a systematic error in these corrections or the presence of an additional fractionating process that has not been accounted for may result in systematic error that varies with site condition. This would result in a similar magnitude of systematic error under similar firn conditions. Thus, the systematic differences in MOT results using different fractionation correction methods or between the three noble gas ratios are largest (up to 0.3 ± 0.3 ^∘C between correction methods and 0.4 ± 0.3 ^∘C between individual noble gas ratios) when comparing results between glacial and interglacial intervals but are minimal (up to 0.1 ± 0.3 ^∘C between correction methods and <0.01 ± 0.3 ^∘C between noble gas ratios) when comparing MOT results between the MIS 4 and MIS 2 glacial intervals.

However, a systematic error in fractionation correction may not be the only explanation for the offset in MOT results between $\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$, $\mathrm{Xe}/{\mathrm{N}}_{\mathrm{2}}$, and $\mathrm{Xe}/\mathrm{Kr}$. Processes that decouple atmospheric noble gas exchange from ocean heat exchange may also introduce systematic error in MOT reconstructions and may affect the krypton, xenon, and nitrogen to different degrees, resulting in differences in MOT reconstructed from $\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$, $\mathrm{Xe}/{\mathrm{N}}_{\mathrm{2}}$, and $\mathrm{Xe}/\mathrm{Kr}$. If this were the cause of the observed offset in MOT results between the three noble gas ratios, we would predict that the offset would be consistent between ice cores.

Considering the MIS 4, MIS 2, and Holocene data from the WAIS Divide record, the relatively sparse and somewhat noisier data make it difficult to discern any trends. However, if anything, the relative offset in the three noble gas proxies is the opposite of that observed for Taylor Glacier. This suggests that the primary mechanism to explain the observed differences in the MOT results between $\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$, $\mathrm{Xe}/{\mathrm{N}}_{\mathrm{2}}$, and $\mathrm{Xe}/\mathrm{Kr}$ is a process that affects these ratios within the firn or ice, rather than the atmospheric inventories of Xe, Kr, and N₂. However, this does not rule out the existence of processes related to the latter. While the slight differences in results with different fractionation correction and between the three noble gas ratios do not affect the conclusions of the study, further investigation is necessary to gain a better understanding of these processes' influence on the MOT proxies and their associated uncertainties.

Appendix B: Ocean solubility effect on CO₂

In order to estimate the magnitude of CO₂ drawdown between 129 and 59.7 ka due to a cooling ocean, we used a simple carbon cycle model and held all model parameters constant except for ocean surface temperatures. The model communicates ocean surface temperature changes to the deep ocean boxes through circulation and mixing; thus, surface temperature changes alter the MOT. The model spin-up reached equilibrium with an atmospheric CO₂ value of 282 ppm, MOT of 4.5 ^∘C, and a mean ocean surface temperature of 18 ^∘C. Ocean surface temperatures in the six surface boxes were then allowed to vary such that the timing and relative magnitude of temperature in each box changed according to a linear scaling to the EPICA Dome C δ²H record (Jouzel et al., 2007) (Fig. 4). The δ²H record was first corrected for changes in seawater δ²H using the δ¹⁸O seawater reconstruction of Waelbroeck et al. (2002) and 8:1 scaling of δ²H:¹⁸O changes. Whole ocean salinity change was also prescribed in the model to account for the solubility effect on CO₂. Salinity was linearly scaled to the Grant et al. (2014) sea level record assuming a preindustrial salinity of 34.72 psu (practical salinity units) and a Last Glacial Maximum salinity of 35.85 psu (Adkins et al., 2002).

The absolute magnitudes of cooling in the six surface boxes between MIS 5e and MIS 4 (Table B1) were chosen such that MOT decreased by 3.1 ^∘C between MIS 5e and MIS 5a and by 0.9 ^∘C across the MIS 5a–4 transition, consistent with the ice core MOT data (Shackleton et al., 2020, and this study). Because our MOT data provide no constraints on the spatial distribution of ocean temperature change, we make relatively simple assumptions on the spatial distribution of ocean cooling. The high-latitude Southern Ocean box temperature change from MIS 5e to MIS 4 (3.5 ^∘C) is smaller than the other surface boxes (5.5 ^∘C) so that the Southern Ocean box does not approach temperatures below the freezing point of seawater. The total relative change in modeled global mean ocean surface temperature is similar to the total relative change in a stack of 136 sediment core records (Table B1, Fig. 4; Kohfeld and Chase, 2017), although the absolute value of mean ocean surface temperature is lower in the model. As discussed in the main text, the choice in the spatial distribution of the modeled ocean temperature change may influence the strength of the modeled solubility pump (Khatiwala et al., 2019). Further experimentation should consider the interplay between different spatial patterns of ocean surface cooling, ocean circulation, and disequilibrium.

Table B1 Magnitude of sea surface cooling prescribed to the carbon cycle box model between the MIS 5e maximum (129 ka) and MIS 4 minimum (67.5 ka). The magnitude of modeled global mean ocean surface cooling is also given for the period from 124 ka to 67.5 ka for comparison with published reconstructions (Kohfeld and Chase, 2017).

N/A: not applicable.