Evaluating manual versus automated benthic foraminiferal δ¹⁸O alignment techniques for developing chronostratigraphies in marine sediment records

Paleoceanographic interpretations of Plio-Pleistocene climate variability over the past 5 million years rely on the evaluation of event timing of proxy changes in sparse records across multiple ocean basins. In turn, orbital-scale chronostratigraphic controls for these records are often built from stratigraphic alignment of benthic foraminiferal stable oxygen isotope (δ¹⁸O) records to a preferred dated target stack or composite. This chronostratigraphic age model approach yields age model uncertainties associated with alignment method, target selection, the assumption that the undated record and target experienced synchronous changes in benthic foraminiferal δ¹⁸O values, and the assumption that any possible stratigraphic discontinuities within the undated record have been appropriately identified. However, these age model uncertainties and their impact on paleoceanographic interpretations are seldom reported or discussed. Here, we investigate and discuss these uncertainties for conventional manual and automated tuning techniques based on benthic foraminiferal δ¹⁸O records and evaluate their impact on sedimentary age models over the past 3.5 Myr using three sedimentary benthic foraminiferal δ¹⁸O records as case studies. In one case study, we present a new benthic foraminiferal δ¹⁸O record for International Ocean Discovery Program (IODP) Site U1541 (54°13^′ S, 125°25^′ W), recently recovered from the South Pacific on IODP Expedition 383. The other two case studies examine published benthic foraminiferal δ¹⁸O records of Ocean Drilling Program (ODP) Site 1090 and the ODP Site 980/981 composite. Our analysis suggests average age uncertainties of 3 to 5 kyr associated with manually derived versus automated alignment, 1 to 3 kyr associated with automated probabilistic alignment itself, and 2 to 6 kyr associated with the choice of tuning target. Age uncertainties are higher near stratigraphic segment ends and where local benthic foraminiferal δ¹⁸O stratigraphy differs from the tuning target. We conclude with recommendations for community best practices for the development and characterization of age uncertainty of sediment core chronostratigraphies based on benthic foraminiferal δ¹⁸O records.

1 Introduction

Stratigraphic alignment, or tuning, of benthic foraminiferal stable oxygen isotope (δ¹⁸O) data to an independently dated target record, such as the LR04 benthic stack (Lisiecki and Raymo, 2005) or the Cenozoic Global Reference benthic foraminifer carbon and oxygen Isotope Dataset (CENOGRID; Westerhold et al., 2020), has been used to develop orbital-scale age models for carbonate-bearing marine sediment cores for decades. This practice arose from early work to cross-link marine sedimentary records with sparse absolute age control using stratigraphic similarities among records that are likely mechanistically linked (e.g., Imbrie et al., 1984; Martinson et al., 1987; Shackleton et al., 1990). Absolute ages for these early targets were mostly inferred from well-constrained orbital variations in Earth's radiative forcing (orbital tuning), while later targets have also been dated using astronomically forced global ice volume models (Imbrie et al., 1984; Lisiecki and Raymo, 2005; Ahn et al., 2017; Westerhold et al., 2020). When available, paleomagnetic reversals and radiometric dates have been used to provide additional absolute age control on these target records (e.g., Mankinen and Dalrymple, 1979; Stein et al., 1993; Thompson and Goldstein, 2006). Benthic-foraminiferal-δ¹⁸O-derived age models of marine sediment cores are widely relied upon in the interpretation of Plio-Pleistocene changes in climate, ocean circulation, and marine biogeochemistry (e.g., Shackleton and Hall, 1984; Ravelo et al., 2004; Elderfield et al., 2012; Hodell et al., 2023). However, benthic-foraminiferal-δ¹⁸O-derived age models are often presented without discussion of the impact of the alignment approach on assigned ages and without estimates of age model uncertainty arising from alignment method, uncertainties within the tuning target, or appropriateness of the chosen target for the specific study region.

Benthic foraminiferal δ¹⁸O stratigraphy is based on the assumption that benthic foraminifera record global variations in bottom-water δ¹⁸O caused mainly by climatic changes in temperature, salinity, and global ice volume over orbital timescales (Pisias et al., 1984; Lisiecki and Raymo, 2005; Skinner and Shackleton, 2004). Benthic foraminifera of the genus Uvigerina are considered an ideal benthic foraminiferal species for the generation of δ¹⁸O stratigraphies because they are believed to calcify in equilibrium with seawater δ¹⁸O, although they occupy a shallow infaunal habitat (Shackleton, 1974). Other widely used species for these efforts include epibenthic foraminifera of the genus Cibicidoides or Cibicides, whose stable oxygen isotope composition is generally corrected by +0.64 ‰ to match presumable equilibrium seawater δ¹⁸O values (Shackleton and Opdyke, 1973). More recent studies have found disequilibrium effects between the δ¹⁸O values of Uvigerina and Cibicidoides species that range between 0.47 ‰ (Marchitto et al., 2014) and 0.73 ‰ (Jöhnck et al., 2012), depending on local bottom-water and pore-water pH conditions. In addition, regional and temporal variations in bottom-water temperature can introduce further discrepancies between the δ¹⁸O values of benthic foraminifera and the bottom water (e.g., Marchitto et al., 2014; Elderfield et al., 2012). However, the premise of benthic δ¹⁸O stratigraphy hinges on the representation of bottom-water δ¹⁸O by the δ¹⁸O of benthic foraminiferal species of the genus Uvigerina (such as U. peregrina and U. hispida) and Cibicidoides (such as C. wuellerstorfi and C. kullenbergi), with a constant correction factor between the two (Lisiecki and Raymo, 2005). In other words, benthic δ¹⁸O stratigraphy assumes the effect of bottom-water and temperature variations in space and time to be minimal (Lisiecki and Raymo, 2005). Once a suitable benthic foraminiferal δ¹⁸O record has been generated, the undated record can be aligned to a dated benthic foraminiferal δ¹⁸O tuning target under the assumptions that the new record is stratigraphically continuous and that deep-ocean δ¹⁸O changes occur synchronously over glacial–interglacial timescales.

Each tuning target carries its own sampling biases and independent age uncertainties that impact the ages assigned to new records using benthic-foraminiferal-δ¹⁸O-based alignments. Community-preferred targets for stratigraphic benthic foraminiferal δ¹⁸O alignment have continuously evolved over the past decades with the increasing number and spatial coverage of available benthic foraminiferal δ¹⁸O records (e.g., Imbrie et al., 1984; Shackleton et al., 1995; Huybers and Wunsch, 2004; Lisiecki and Raymo, 2005; Ahn et al., 2017; Westerhold et al., 2020). Early targets, such as SPECMAP, resulting from the SPECtral Mapping Project, were built by stacking the planktonic foraminiferal δ¹⁸O records of five or fewer sediment cores and were typically limited to within the last 800 kyr of the late Pleistocene (Imbrie et al., 1984; Pisias et al., 1984; Martinson et al., 1987). These early targets were further restricted in their ability to reflect global changes in deep-seawater δ¹⁸O because they may have been biased by surface-water δ¹⁸O variations due to the inclusion of planktonic foraminiferal δ¹⁸O records. Later targets, such as LR04 (Lisiecki and Raymo, 2005) and Prob-stack (Ahn et al., 2017), improved upon these approaches by stacking benthic foraminiferal δ¹⁸O records from 57 and 180 sediment cores around the world, respectively, and extending the interval covered by the benthic foraminiferal stack to ∼5 Ma. Composite splice targets built by combining discrete high-resolution benthic foraminiferal δ¹⁸O records in series, such as the S95 composite (Shackleton, 1995), the benthic δ¹⁸O megasplice (De Vleeschouwer et al., 2017), and the CENOGRID splice (Westerhold et al., 2020), have also been generated and provide a means to examine the high-frequency variability within benthic foraminiferal δ¹⁸O data that may be smoothed over or lost in the generation of a global stack. On the other hand, such composite targets are heavily influenced by regional bottom-water δ¹⁸O and other variability at the limited number of sites included in the record.

Increasing data coverage and the recovery of more core sites across the globe have revealed regional differences in the timing of benthic foraminiferal δ¹⁸O changes across different ocean basins and water masses (Lisiecki and Raymo, 2009; Skinner and Shackleton, 2005; Stern and Lisiecki, 2014; Waelbroeck et al., 2011; Labeyrie et al., 2005), which challenge one of the basic assumptions of benthic foraminiferal δ¹⁸O stratigraphies. Early work, based on two independently dated sediment cores from the Pacific Ocean and Atlantic Ocean, identified a ∼4 kyr lag in the earliest signal of glacial termination between the deep North Atlantic and the deep equatorial Pacific (Skinner and Shackleton, 2005). Subsequent development of Pacific- and Atlantic-specific benthic foraminiferal δ¹⁸O stacks (i.e., the LR09 stacks) confirmed these results and suggested an average lag of 1.6 kyr in deglacial benthic-seawater δ¹⁸O changes in the Pacific relative to the Atlantic (Lisiecki and Raymo, 2009). More recently, the generation of seven regional stacks over the last glacial cycle from 252 sediment cores has further refined our understanding of regional asynchronicity in the timing of the benthic foraminiferal δ¹⁸O signal (Stern and Lisiecki, 2014). Stern and Lisiecki (2014) find a maximum lag of 4 kyr between the earliest signals of glacial termination detected in the intermediate South Atlantic and the latest signals detected in the deep Indian Ocean as well as shorter lags occurring between other water masses and regions. If left unconsidered, such regional differences in the timing of benthic foraminiferal δ¹⁸O changes relative to a global tuning target can introduce systematic age uncertainties into the resulting benthic foraminiferal δ¹⁸O chronostratigraphy.

Benthic foraminiferal δ¹⁸O records are traditionally aligned manually to the chosen target via user-defined identification and assignment of characteristic δ¹⁸O signals at specific tie points, with the assumption of a constant sedimentation rate between selected tie points. Tie points, or age markers, are typically assigned at mid-points of sharp transitions within the benthic foraminiferal δ¹⁸O record because the mid-point depth is easy to identify. However, tie points are also frequently assigned to local maxima and minima within the benthic foraminiferal δ¹⁸O record, where depth uncertainties may be higher due to limitations of sampling resolution, and no community-accepted guidelines exist for the establishment of benthic foraminiferal δ¹⁸O chronostratigraphies by manual tuning. Hence, the number of tie points identified and their exact location within each benthic foraminiferal δ¹⁸O record can vary among users, which may significantly impact the final age model of the sediment core. This user-induced variance is particularly significant in core intervals where the alignment of the benthic foraminiferal δ¹⁸O record to the tuning target is not straightforward. Further, conventional manual stratigraphic alignment efforts based on benthic foraminiferal δ¹⁸O records generally lack age uncertainty estimates, although in rare cases ad hoc uncertainties are assigned to individual tie points. Consequently, benthic-foraminiferal-δ¹⁸O-derived age models are a useful, but subjective, dating tool for establishing age models of marine sedimentary records.

While alignment software such as AnalySeries (Paillard et al., 1996) can compute quantitative correlation coefficients for user-defined alignments, such tools do not eliminate the subjectivity and related uncertainties inherent in the user-defined visual alignment approach and do not provide estimates of chronostratigraphic age uncertainties. In contrast, a probabilistic alignment algorithm like the MATLAB-based Hidden Markov Model (HMM)-Match (Lin et al., 2014) largely eliminates the subjectivity resulting from reliance on multiple user-defined tie points. In fact, such algorithms can generate user-independent age estimates for every data point within an input benthic foraminiferal δ¹⁸O record without presuming constant sedimentation rates across periods of major climatic change or other predefined intervals. HMM-Match additionally provides statistically informed estimates of alignment-derived age uncertainty based on reported benthic foraminiferal δ¹⁸O variance of the tuning target and based on differences in the benthic foraminiferal δ¹⁸O tuning target and the input δ¹⁸O record (Lin et al., 2014; Ahn et al., 2017). Consequently, HMM-Match has become an increasingly popular tool for the development of Plio-Pleistocene chronostratigraphies in marine sediment cores based on benthic foraminiferal δ¹⁸O (e.g., Hodell et al., 2023; Ford and Raymo, 2020). However, as with the manual alignment approaches, alignment algorithms like HMM-Match (Lin et al., 2014) do not currently account for the age model uncertainties associated with regional offsets between the timing of benthic foraminiferal δ¹⁸O changes at a given location compared to a globally averaged dataset.

Here, we assess the uncertainties associated with each step of the stratigraphic alignment of benthic foraminiferal δ¹⁸O records and evaluate their impact on sediment core age models over the past 3.5 Myr using three example sites from the International Ocean Discovery Program (IODP) and its predecessor, the Ocean Drilling Program (ODP; Fig. 1). As part of this analysis, we present a new benthic foraminiferal δ¹⁸O record from IODP Site U1541 in the central South Pacific Ocean, recovered on IODP Expedition 383, that is based on ∼1600 benthic foraminiferal δ¹⁸O measurements over the last 3.5 Myr. We additionally examine ODP Site 1090 in the South Atlantic and the ODP Site 980/981 composite in the North Atlantic using published benthic foraminiferal δ¹⁸O records (Venz and Hodell, 2002; Hodell et al., 2000; Oppo et al., 1998; Mcmanus et al., 1999; Flower et al., 2000; Raymo et al., 2004). We assess the impact of stratigraphic alignment approach (manual versus automated probabilistic methods) on age model output by independently aligning the benthic foraminiferal δ¹⁸O records of each example site to the LR04 stack using first the conventional manual alignment method and then the open-source HMM-Match probabilistic alignment algorithm of Lin et al. (2014). We evaluate the impact of stratigraphic target choice on age model output by aligning the benthic foraminiferal δ¹⁸O records of each example site to a suite of reference targets (LR04, LR09, Prob-stack, and the CENOGRID splice) using the HMM-Match alignment algorithm (Lin et al., 2014). We additionally assess the sensitivity of such age models to the assumptions of global synchronicity in seawater δ¹⁸O changes and stratigraphic continuity within the undated record. Finally, we present suggested best practices for the development and characterization of age uncertainty of sediment core age models based on benthic foraminiferal δ¹⁸O records.

Figure 1 Location map for sediment cores included in the LR04 (smaller black circles; Lisiecki and Raymo, 2005) and Prob-stack global benthic foraminiferal stable oxygen isotope (δ¹⁸O) tuning target stacks (larger open circles; Ahn et al., 2017). International Ocean Discovery Program (IODP) Site U1541, Ocean Drilling Program (ODP) Site 1090, and ODP Site 980/981 are denoted by stars. Color bar indicates the maximum age of the benthic foraminiferal δ¹⁸O records from each site included in the LR04 and Prob-stack. Four benthic foraminiferal δ¹⁸O records have significant data gaps and are numbered as follows: (1) MV0502-4JC data cover the intervals 2975–1575, 679–441, and 210–23 ka; (2) ODP Site 1014 data cover the last 187 kyr and 916–386 ka; (3) ODP Site 1241 data cover the intervals 5691–2443 and 2129–0 ka; (4) IODP Site U1313 data cover the intervals 3331–2414, 913–788, 700–600, and 356–549 ka.

2 Sampling and measurements

2.1 Sediment cores

IODP Site U1541 (54°13^′ S, 125°25^′ W; 3604 m water depth) was recovered on the western flank of the East Pacific Rise, ∼300 km from the ridge axis, in the sub-Antarctic Pacific Ocean on IODP Expedition 383 (Fig. 1; Winckler et al., 2021). IODP Site U1541 is currently bathed in Lower Circumpolar Deep Water (LCDW) and represents one of the first continuous sedimentary sections from the sub-Antarctic Pacific that spans the Plio-Pleistocene with orbital-scale resolution (Lamy et al., 2019; Winckler et al., 2021). As such, sediments from IODP Site U1541 can be used to investigate the dynamic atmospheric, oceanographic, and Antarctic processes impacting the large Pacific sector of the Southern Ocean across the wide range of climatic variability of the past 5 Myr and more. We present new benthic foraminiferal δ¹⁸O data from IODP Site U1541 for the past 3.5 Myr, and the age model results of our stratigraphic alignment analyses represent the first orbital-scale U1541 age models based on benthic foraminiferal δ¹⁸O tuning available for this site. Due to its location in the sparsely studied but climatically important South Pacific (Fig. 1), we use IODP Site U1541 as an end-member case study to investigate the age uncertainties associated with benthic foraminiferal δ¹⁸O stratigraphic alignment of a record from an under-sampled study region (i.e., the Pacific Ocean) to common alignment targets.

ODP Site 1090 (42°55^′ S, 8°54^′ E; 3702 m water depth) was recovered from the southern flank of the Agulhas Ridge in the sub-Antarctic Atlantic Ocean on ODP Leg 177 (Fig. 1; Gersonde et al., 1999). ODP Site 1090 is also bathed by LCDW, though it is positioned close to the lower boundary of North Atlantic Deep Water (NADW). The original stratigraphic alignments of the benthic foraminiferal δ¹⁸O records for ODP Site 1090 and its pre-site survey core, TTN057-6, were manually determined using the benthic foraminiferal δ¹⁸O record of Deep Sea Drilling Project (DSDP) Site 607 as an alignment target (Venz and Hodell, 2002; Hodell et al., 2000). Sediment chronology of the DSDP Site 607 record itself was determined via orbital tuning of its benthic foraminiferal δ¹⁸O record between paleomagnetic tie points (Raymo et al., 1989, 1990; Ruddiman et al., 1989). We use ODP Site 1090 as a case study to investigate the age uncertainties inherent in the choice of alignment technique for a site in the Atlantic Ocean that is well represented by common tuning targets and to address the impact of stratigraphic discontinuity on the resulting benthic foraminiferal δ¹⁸O chronostratigraphy.

ODP Site 980 (55°29^′ N, 14°42^′ W; 2169 m water depth) and ODP Site 981 (55°29^′ N, 14°39^′ W; 2173 m water depth) were recovered from the Feni Drift in the northeastern Atlantic on ODP Leg 162 (Fig. 1; Jansen et al., 1996). ODP Sites 980 and 981 are currently bathed in lower NADW. The original stratigraphic alignments for the combined ODP Site 980/981 benthic foraminiferal δ¹⁸O records, spliced together at ∼860 ka (Raymo et al., 2004), were manually determined using the Shackleton et al. (1990) target chronology for ODP Site 677 (Oppo et al., 1998; McManus et al., 1999; Flower et al., 2000; Raymo et al., 2004). The sediment chronology of ODP Site 677 was determined via orbital tuning of its benthic foraminiferal δ¹⁸O record following the calculations of Berger and Loutre (1988) under the assumption of a constant phase relationship between astronomical radiative forcing and climatic response (Shackleton et al., 1990). Due to the high-resolution benthic foraminiferal δ¹⁸O record available for this site and its location in the densely sampled North Atlantic, we use ODP Site 980/981 as a case study with maximum influence of NADW. This site also enables investigation into the age uncertainties associated with benthic foraminiferal δ¹⁸O stratigraphies under favorable alignment conditions. In the North Atlantic, regional seawater temperature and δ¹⁸O differences, as well as associated temporal offsets in benthic foraminiferal δ¹⁸O changes, between the undated record and the global LR04 and Prob-stack targets, which are heavily weighted by North Atlantic sampling, are expected to be minimal.

2.2 Benthic foraminiferal stable oxygen isotope analyses

The benthic foraminiferal δ¹⁸O data from IODP Site U1541 were measured in the upper 85 m core composite depth below seafloor (CCSF) of the U1541 shipboard splice (with identification code: CCSF_383_U1541_ABC_20190624; Winckler et al., 2021) at a 5 cm sampling resolution. This sampling resolution yields an average temporal resolution of ∼2 kyr per sample over the 3.5 Myr record examined. These data were combined with benthic foraminiferal δ¹⁸O data spanning the last four glacial cycles from PS75/059-2 (Ullermann et al., 2016), which is the pre-site survey core of IODP Site U1541. To combine the benthic foraminiferal δ¹⁸O records from both cores, the PS75/059-2 data were mapped onto a common U1541 depth scale via a manual stratigraphic alignment of high-resolution X-ray fluorescence (XRF) Fe intensity variations in PS75/059-2 (Lamy et al., 2014) and U1541 that result in 22 tie points (Fig. S1 and Table S1 in the Supplement).

Bulk sediment samples from IODP Site U1541 were freeze-dried, sieved with deionized water over a 150 µm sieve, and dried at ∼45 °C to prepare for benthic foraminiferal δ¹⁸O measurements. The >150 µm sample fraction was then examined for benthic foraminifera under a stereomicroscope. One to five tests of benthic foraminifera of the genera Cibicidoides, Cibicides, and Uvigerina, following the morphotype description of Gottschalk et al. (2016) for Cibicidoides and Cibicides, were selected for each δ¹⁸O measurement with a wet brush. The samples were measured on a Thermo Scientific Delta V+ isotope ratio mass spectrometer with an automated Kiel IV Carbonate Device at Lamont-Doherty Earth Observatory (LDEO), Columbia University, USA. The data were corrected relative to the NBS (National Bureau of Standards)-19 calcite standard, which was analyzed every 10 samples, and are reported in δ notation with respect to the Vienna PeeDeeBee (VPDB) international standard. The long-term standard deviation for δ¹⁸O of the NBS-19 standard on this instrument is 0.06 ‰ VPDB. The mean offset between Uvigerina spp. and Cib. spp. at IODP Site U1541 is 0.49±0.19 ‰ VPDB (n=87). We therefore corrected our Cibicidoides and Cibicides δ¹⁸O record for δ¹⁸O-disequilibrium effects according to the measured mean offset. Intra-species and intra-morphotype Cibicidoides and Cibicides δ¹⁸O variability at IODP Site U1541 is similar to the analytical uncertainty of the δ¹⁸O analyses (σ<0.06 ‰ VPDB).

We use previously published benthic foraminiferal δ¹⁸O data for ODP Site 1090 (Venz and Hodell, 2002; Hodell et al., 2000). The data from ODP Site 1090 were measured on one to three specimens per sample of the species C. wuellerstorfi using a Finnigan MAT 252 mass spectrometer with an estimated 1σ analytical uncertainty of ±0.06 ‰ VPDB based on the reproducibility of repeat measurements of the NBS-19 internal standard (n=321; Venz and Hodell, 2002). The benthic foraminiferal δ¹⁸O measurements of the ODP Site 1090 survey core, TTN057-6-PC4, were made on C. wuellerstorfi and C. kullenbergi and carry an estimated 1σ analytical uncertainty of <0.10 ‰ VPDB (Hodell et al., 2000). The average sampling resolution for the combined benthic foraminiferal δ¹⁸O record of ODP Site 1090 and TTN057-6-PC4, spliced together at ∼400 ka, is ∼3 kyr over the 2.9 Myr interval examined (Venz and Hodell, 2002; Hodell et al., 2000). Where repeat analyses or measurements from both species were available, we employ the average foraminiferal δ¹⁸O value at each depth for our alignment analyses.

We use previously published benthic foraminiferal δ¹⁸O data for ODP Site 980/981. The data from ODO Site 980 were measured on one to three specimens per sample of C. wuellerstorfi or C. kullenbergi on a Finnigan MAT 252 or precision isotope ratio mass spectrometer (PRISM) with a 1σ analytical uncertainty of ±0.08 ‰ VPDB (McManus et al., 1999; Flower et al., 2000; Oppo et al., 1998). The ODP Site 981 data were measured on benthic foraminifera of the genus Cibicidoides (Raymo et al., 2004). The average sampling resolution for the combined ODP Site 980/981 record is ∼1 kyr over the past 1.8 Myr.

3 Stratigraphic alignment

3.1 Manual alignment

The manual alignment of the benthic foraminiferal δ¹⁸O record from IODP Site U1541 was generated by visual alignment of characteristic peaks and troughs in the record to the benthic foraminiferal LR04 stack using the publicly available QAnalySeries software (Kotov and Pälike, 2018), which is based on the original AnalySeries program of Paillard et al. (1996). The manual alignments of ODP Site 1090 and ODP Site 980/981 were generated by realigning the original age model tie points for each site (e.g., Venz and Hodell, 2002; Raymo et al., 2004) to the LR04 stack using the QAnalySeries software. These manual alignment approaches yielded an average temporal spacing of ∼30 kyr between tie points for IODP Site U1541 (116 tie points over the last 3.5 Myr; Fig. 2; Table S2) and ODP Site 980/981 (66 tie points over the last 1.8 Myr; Fig. 2; Table S4) and ∼45 kyr between tie points for ODP Site 1090 (65 tie points over the last 2.9 Myr; Fig. 2, Table S3). Due to the subjectivity inherent in user-defined alignments, the manual alignments presented in Tables S2–S4 represent just one of many possible manual alignment outcomes for each record.

Figure 2 Benthic foraminiferal δ¹⁸O records based on Cibicidoides and Cibicides spp. from IODP Site U1541 (purple, squares), ODP Site 1090 (black, filled circles), and ODP Sites 980/981 (green, open circles). Data plotted using age models derived from manual alignment of each record to the LR04 global stack (Lisiecki and Raymo, 2005) for (a) 0–1200 ka, (b) 1200–2400 ka, and (c) 2400–3600 ka. The benthic foraminiferal δ¹⁸O records for ODP Site 1090 and ODP Site 980/981 are plotted with a −0.5 ‰ shift and a −1.0 ‰ shift, respectively, for clarity. Stratigraphic tie points for each manual alignment are denoted on the x axis (Table S1).

3.2 Automated probabilistic alignment

HMM-Match iteratively aligns every data point of a benthic foraminiferal δ¹⁸O record along a depth scale to the benthic foraminiferal δ¹⁸O target record along an age scale using a probabilistic assessment of the implied sedimentation rates associated with each alignment. The HMM-Match algorithm accounts for natural variance in a benthic foraminiferal δ¹⁸O dataset (e.g., due to bioturbation, spatial variability, and measurement uncertainty) by assuming that the residual δ¹⁸O value between each input record and the target will fall along a normal probability distribution (Lin et al., 2014). The algorithm checks the implied relative sedimentation rate changes associated with each alignment fit against their natural likelihood based on the distribution of relative sedimentation rate changes observed in an independent compilation of 37 radiocarbon-derived sediment core chronologies over the last 40 kyr (Lin et al., 2014). The automated HMM-Match algorithm additionally generates Bayesian confidence bands for the resulting age model that provide age model uncertainties associated with the stratigraphic benthic foraminiferal δ¹⁸O alignment (Lin et al., 2014), in contrast to conventional manual stratigraphic alignment techniques.

While HMM-Match can generate objective and reproducible alignments of benthic foraminiferal δ¹⁸O records, the algorithm requires user-derived depth–age tie point estimates for the start and end points for the entire record or for each alignment segment. Alignment segments for each site were chosen to maximize relative stratigraphic continuity between gaps in sediment recovery (such as coring gaps) or data availability within each record. For IODP Site U1541, six input tie points were used to align the ∼3.5 Myr record that spans two coring gaps between 31.78–32.75 and 75.67–77.12 m CCSF-A (Winckler et al., 2021; Table 1). The specific depth–age values for these tie points were initially determined visually based on peaks and troughs in the benthic foraminiferal δ¹⁸O data near the top and bottom of each segment. Subsequently, where the visual alignment of segment ends was deemed ambiguous, the visually determined segment end point ages were varied by ±40 kyr and rerun through the HMM-Match algorithm. The segment tie points of IODP Site U1541 were finalized based on the input tie points that yielded the lowest segment average age uncertainties in HMM-Match (Table 1).

Table 1 Depth–age tie points for HMM-Match-based alignments.

^∗ Tie points were excluded from the HMM-Match with Prob-stack Alignment B for ODP Site 1090 to examine sensitivity to undetected hiatuses within the record.

For ODP Sites 1090 and 980/981, the depth–age inputs for HMM-Match were selected from the first and last tie points listed for each continuous section of published chronology (Table 1; Venz and Hodell, 2002; Raymo et al., 2004). At ODP Site 1090, where previous work has identified six mid-core hiatuses (Venz and Hodell, 2002), the HMM-Match alignment was constructed from seven segments bounded by these hiatuses. At ODP Site 980/981, where no such hiatuses have been reported, the HMM-Match alignment was constructed from only two segments – one segment each for the ODP Site 980 and ODP Site 981 records (Table 1; Raymo et al., 2004). We note that the automated HMM-Match algorithm typically assigns zero age uncertainty at the depth of each tie point used to bound the segments, and thus the age uncertainty estimates reported at the start and end of each segment should be treated with caution.

Automated alignment algorithms, such as HMM-Match (Lin et al., 2014), introduce a risk of disconnecting users from nuances in the sedimentation history of a given site and possible stratigraphic discontinuities within a given benthic foraminiferal δ¹⁸O record. In order to test the HMM-Match algorithm for sensitivities to hiatuses within the input data, we use the ODP Site 1090 benthic foraminiferal δ¹⁸O data to generate an additional alignment to the Prob-stack target (ODP Site 1090 Alignment B). The ODP Site 1090 Alignment B was based on the deliberate omission of some of the identified hiatuses in the benthic foraminiferal δ¹⁸O record of ODP Site 1090, using fewer tie point constraints and resulting in only four segments (Table 1).

3.3 Tuning targets

Since 2005, the LR04 stack (Fig. 3; Lisiecki and Raymo, 2005) has been the most commonly employed reference stack for global ocean δ¹⁸O variations on glacial–interglacial timescales over the past 5.3 Myr. While the LR04 benthic foraminiferal δ¹⁸O stack consists of 57 stacked records available at the time of its development, the stack is heavily weighted towards Atlantic sediment core data due to limited core availability from the Pacific and Indian ocean basins (Fig. 1). Due to decreasing record availability with increasing distance into the past, the resolution and number of records included in LR04 decrease moving further back in time (Lisiecki and Raymo, 2005). The late Pleistocene interval of the LR04 benthic foraminiferal δ¹⁸O stack is based on the average of 30 individual sediment core records yielding a 1 kyr stack resolution, while the Pliocene interval spanning 3 to 5.3 Ma is built on the average of 15 or fewer individual records and has a 5 kyr stack resolution. Absolute ages for LR04 were determined from orbital tuning of a simple nonlinear ice volume model to the 21 June isolation curve for 65° N based on the orbital solution of Laskar et al. (1993).

Figure 3 Selected benthic foraminiferal δ¹⁸O alignment target records over the past 5 Myr. Plotted targets include Prob-stack (black; Ahn et al., 2017), LR04 (blue; Lisiecki and Raymo, 2005), the CENOGRID splice (gray; Westerhold et al., 2020), and the LR09 Atlantic (yellow) and Pacific stacks (dark red; Lisiecki and Raymo, 2009) for (a) 0–1000 ka, (b) 1000–2000 ka, (c) 2000–3000 ka, (d) 3000–4000 ka, and (e) 4000–5000 ka.

The LR09 Pacific stack is composed of 14 benthic foraminiferal δ¹⁸O records, primarily from the equatorial and North Pacific Ocean, while the LR09 Atlantic stack is composed of 20 benthic foraminiferal δ¹⁸O records from the North and South Atlantic Ocean (Lisiecki and Raymo, 2009). These stacks were generated by averaging the available Pacific and Atlantic records from LR04 and three additional sites over the past 800 kyr, with a ±1 kyr smoothing window for data points between 0 and 500 ka and a ±2 kyr smoothing window between 500 and 800 ka (Fig. 3; Lisiecki and Raymo, 2009). The LR09 regional stacks were generated to examine regional differences in the timing of bottom-water δ¹⁸O variability and only cover the last 800 kyr (Lisiecki and Raymo, 2009) but are included in our analysis to examine the impact of regional-specific targets on the HMM-Match-derived ages for our three example sites. The absolute age constraints on the LR09 Pacific and Atlantic stacks are the same as in the LR04 stack (Lisiecki and Raymo, 2009).

The Prob-stack builds upon the LR04 global benthic foraminiferal δ¹⁸O stack, covering the last 5.0 Myr with the inclusion of 123 additional records, more data coverage from the Pacific Ocean, and a realignment of each record using the HMM-Match algorithm (Figs. 1 and 3; Ahn et al., 2017). As with the LR04 benthic foraminiferal δ¹⁸O stack, the resolution and number of records included in Prob-stack decrease further back in time, from 1 kyr resolution and more than 120 records in the late Pleistocene to 5 kyr resolution and fewer than 20 records by 5 Ma (Fig. 1; Ahn et al., 2017). Although the exact structure of Prob-stack varies slightly from LR04 (Fig. 3), the Prob-stack target is based on the same absolute age constraints as the LR04 age model (Ahn et al., 2017).

In contrast to stacked targets, such as LR04 (Lisiecki and Raymo, 2005), LR09 (Lisiecki and Raymo, 2009), and Prob-stack (Ahn et al., 2017), the CENOGRID composite represents a continuous splice of the benthic foraminiferal δ¹⁸O records with the highest resolution to span the entire Cenozoic (i.e., the past 66 Myr; Westerhold et al., 2020). CENOGRID covers the past 30 Myr with at least 2 kyr resolution (Westerhold et al., 2020). The CENOGRID splice has the advantage that it resolves high-frequency variability in benthic foraminiferal δ¹⁸O and that it can be used for alignment of sedimentary records that extend beyond the Plio-Pleistocene. The drawback of using CENOGRID as an alignment target, however, is that the composite splice consists entirely of low-latitude records and its 0–5 Ma interval is heavily biased towards benthic foraminiferal δ¹⁸O records from the Ceara Rise in the equatorial Atlantic (Fig. 3; Westerhold et al., 2020). Consequently, any regional variability in benthic foraminiferal δ¹⁸O observed in an undated record from a region outside the equatorial Atlantic would complicate millennial-scale probabilistic alignment to the CENOGRID splice. The absolute age constraints for the CENOGRID composite from 0 to 5 Ma are based on astronomical tuning of sediment images data, physical property data, and benthic foraminifera δ¹⁸O data from the Ceara Rise sites (i.e., ODP Sites 925, 926, 927, 928, and 929) to the Laskar et al. (2004) orbital solution (Wilkens et al., 2017; De Vleeschouwer et al., 2017; Westerhold et al., 2020). While there is general agreement between the LR04 and CENOGRID ages over the past 3.5 Myr examined in this study, there is notable disagreement between them from 1.9 to 1.8 Ma (Wilkens et al., 2017).

The alignment uncertainties computed by HMM-Match are based on the degree of misfit between the input benthic foraminiferal δ¹⁸O record and the target, which is a function of the sampling resolution and signal-to-noise ratio of the input record as well as any uncertainties or estimated variance in the benthic foraminiferal δ¹⁸O value of the target record for each point in time (Lin et al., 2014). For alignments to Prob-stack and LR04, the published uncertainties in δ¹⁸O values for each reference target were used (Ahn et al., 2017; Lisiecki and Raymo, 2005). The 1σ uncertainty of LR04 represents the mean standard error for each point within the stack and has an average value of 0.05 ‰ to 0.06 ‰ over the past 5.3 Myr (Lisiecki and Raymo, 2005). In contrast, the reported 1σ uncertainty of Prob-stack represents the observed variance in the global records at each point within the stack, not the mean standard error, and has an average value of 0.18 ‰ over the past 5.0 Myr (Ahn et al., 2017).

No corresponding uncertainties were listed for the LR09 Pacific and Atlantic stacks (Lisiecki and Raymo, 2009) or for the CENOGRID splice (Westerhold et al., 2020). Instead, we assign the 1σ δ¹⁸O uncertainties reported for 0–800 ka in LR04 to the LR09 Pacific and Atlantic stacks. We note that these values likely represent the upper end of the possible mean standard error of each regional stack because the LR09 Pacific stack (n=14) and the LR09 Atlantic stack (n=20) are both comprised of fewer and more regionally clustered records than the global LR04 stack (n=57). Similarly, we assign to the CENOGRID splice a 1σ δ¹⁸O uncertainty of 0.10 ‰ throughout the record, which is only slightly higher than analytical uncertainties of 0.06 ‰ to 0.08 ‰ reported for the Ceara Rise datasets used to generate the Plio-Pleistocene portion of the CENOGRID splice (ODP Sites 925, 926, 927, 928, and 929; Westerhold et al., 2020; Bickert et al., 1997; deMenocal et al., 1997; Tiedemann and Franz, 1997; Billups et al., 1998). As each stacked target record was developed with the assumption of synchronous changes in benthic foraminiferal δ¹⁸O values among the individual records that form the basis of each stack, none of the uncertainties in δ¹⁸O values reported or assigned for these stacks account for the regional variability in the timing of benthic foraminiferal δ¹⁸O changes in the global ocean. Nor do these uncertainties reflect potential offsets in the timing of benthic foraminiferal δ¹⁸O changes in an input record relative to the tuning target.

In the probabilistic HMM-Match alignment algorithm (Lin et al., 2014), tuning targets with higher designated benthic foraminiferal δ¹⁸O variance will yield stratigraphic alignment outputs with higher age model uncertainties. Given that the benthic foraminiferal δ¹⁸O uncertainties assigned to each target record are based on different assumptions and measures of variance (e.g., the 1σ δ¹⁸O values for Prob-stack reflect variance in the global dataset, while the 1σ δ¹⁸O values for CENOGRID reflect analytical uncertainty for measurements at a single site), we do not provide an in-depth comparison of the differences in HMM-Match-generated uncertainties for alignment of a single site to different tuning targets. Rather, we evaluate the impact of tuning target on resulting benthic foraminiferal δ¹⁸O chronostratigraphies by comparing the age offset between benthic foraminiferal δ¹⁸O alignments to Prob-stack with those to LR04, LR09 regional stacks, and CENOGRID. We employ Prob-stack as our benchmark tuning target because it is the most globally representative target, assuming globally synchronous changes in benthic δ¹⁸O (Ahn et al., 2017).

4 Results

For each example site (IODP Site U1541, ODP Site 1090, and ODP Site 980/981), we present the age models and age uncertainties (Tables S5–S7) generated by stratigraphic alignment of the benthic foraminiferal δ¹⁸O record via manual and automated probabilistic tuning (using the HMM-Match algorithm; Lin et al., 2014) to the LR04 stack and HMM-Match-based alignments of the benthic foraminiferal δ¹⁸O record to a selected suite of reference targets (LR04, LR09, Prob-stack, and the CENOGRID splice). At ODP Site 1090, we additionally present the age models and age uncertainties associated with the Alignment B sensitivity test for stratigraphic discontinuity.

4.1 Age models and alignment outputs for IODP Site U1541

The manual and probabilistic HMM-Match-based alignments of the benthic foraminiferal δ¹⁸O record of IODP Site U1541 to the LR04 stack yield sedimentary age models that are largely consistent on the orbital scale, except for two apparent glacial cycle offsets in alignment, generating age model offsets up to 40 and 47 kyr at ∼3.0 and ∼3.4 Ma, respectively (Fig. 4). The direction of age model offset (e.g., whether HMM-Match-based ages are younger or older than the manually derived ages) varies throughout the record. There is no systematic lead or lag associated with the age model of the HMM-Match-derived alignment relative to the age model of the manually derived alignment. Including (excluding) the two intervals of complete cycle offsets, the mean value of the absolute difference between the manually aligned and HMM-Match aligned LR04 ages for IODP Site U1541 averages 4.6 kyr (3.2 kyr) across the entire record. Periods with a high age model offset between the manually derived and HMM-Match-derived age models at IODP Site U1541 are indicative of intervals within the benthic foraminiferal δ¹⁸O record where stratigraphic similarities to the LR04 target stack are ambiguous due to regional differences in the shape of the benthic foraminiferal δ¹⁸O record, differences in temporal resolution between the input record and the target, and/or stratigraphic discontinuities (Fig. 4).

The stratigraphic alignment of the benthic foraminiferal δ¹⁸O record of IODP Site U1541 to the Prob-stack target using HMM-Match provides a sub-orbital-scale chronology spanning the past 3.5 Myr in the upper 85 m (CCSF-A) of this site (Fig. 4). The HMM-Match-generated 1σ age uncertainties associated with this alignment are on average 2 kyr across the record, although they are typically lower in the most recent 1 Myr, where benthic foraminiferal δ¹⁸O variations occur at lower frequencies and with larger amplitude (Fig. 4). The highest age uncertainties in the record, reaching a maximum unidirectional 1σ age uncertainty of 20 kyr, occur between 2190 and 2110 ka and likely result from a poorly defined Marine Isotope Stage (MIS) 82 in the benthic foraminiferal δ¹⁸O at IODP U1541 data relative to Prob-stack (Fig. 4).

Figure 4 Stratigraphic benthic foraminiferal δ¹⁸O alignments for IODP Site U1541. (a) Prob-stack benthic foraminiferal δ¹⁸O target (Ahn et al., 2017) for comparison. (b) Manual (light blue) and automated HMM-Match-based (dark blue; Lin et al., 2014) alignments of the benthic foraminiferal δ¹⁸O record from IODP Site U1541 to the LR04 benthic foraminiferal δ¹⁸O stack (Lisiecki and Raymo, 2005). (c) Age model offset between the age models for IODP Site U1541 plotted in (b). (d) Automated alignments of the benthic foraminiferal δ¹⁸O record from IODP Site U1541 to Prob-stack (black), CENOGRID (gray; Westerhold et al., 2020), and the LR09 Pacific stack (brown; Lisiecki and Raymo, 2009) using HMM-Match. (e) 1σ age model uncertainties reported by HMM-Match for the benthic foraminiferal δ¹⁸O alignment of IODP Site U1541 to Prob-stack. (f) Age offsets between the age models at IODP Site 1541 generated by the automated HMM-Match-based alignments to Prob-stack and LR04 (blue), Prob-stack and CENOGRID (gray), and Prob-stack and LR09 Pacific stack (brown). The zero line is plotted in (c), (e), and (f) for reference.

Stratigraphic alignments of the benthic foraminiferal δ¹⁸O data from IODP Site U1541 to the LR04, LR09, and CENOGRID targets yield sedimentary chronostratigraphies that are similar on the orbital scale to the Prob-stack-derived U1541 age model (Fig. 4). However, there is a clear difference between alignments between 2220 and 2110 ka, overlapping with the period of high age model uncertainty for the benthic foraminiferal δ¹⁸O alignment to Prob-stack reported by HMM-Match, when the magnitude of absolute age model offsets between the Prob-stack-based alignment and the IODP Site U1541 alignments to LR04 and CENOGRID reaches up to 60 kyr (Fig. 4). For the entire 3.5 Myr record, the average magnitude of absolute age offsets between the Prob-stack-based age model and those based on alignments to LR04 and CENOGRID are 5.2 and 5.7 kyr, respectively (Fig. 4, Table S5). Over the last 800 kyr, the period for which the LR09 benthic δ¹⁸O stack is available for the Pacific Ocean, absolute age offsets between target alignment outputs are lower. Specifically, the average absolute age offsets between Prob-stack and LR04, LR09, and CENOGRID-based alignment ages are 2.5, 2.3, and 4.7 kyr, respectively, across this interval (Table S5). The largest offsets in this interval occur between 480 and 440 ka (i.e., MIS 12 and 13), when age offsets between the manually derived and HMM-Match-derived alignments at IODP Site U1541 are also above average (Fig. 4).

The absolute values of the age offsets between the benthic foraminiferal δ¹⁸O alignment of IODP Site U1541 to Prob-stack and to the LR09 Pacific stack vary from near zero up to 16 kyr, with an average of 2.3 kyr (Fig. 4). The 16 kyr age offset occurs at 478 ka and coincides with relatively high age offsets between the benthic foraminiferal δ¹⁸O alignment of IODP Site U1541 to Prob-stack and the LR04 (16 kyr) and CENOGRID targets (14 kyr), respectively. This suggests that the age offset is unlikely to be caused by regional asynchronicity between Prob-stack and the LR09 Pacific stack over the past 800 kyr. The average age offset of 2.3 kyr between the benthic foraminiferal δ¹⁸O alignment of IODP Site U1541 to Prob-stack and the LR09 Pacific stack is consistent with temporal differences of up to 4 kyr in benthic foraminiferal δ¹⁸O records from the Pacific and the North Atlantic (Lisiecki and Raymo, 2005, 2009; Skinner and Shackleton, 2005). While tuning to the closest regional stack may mitigate age uncertainties associated with spatially heterogenous and temporally diachronous δ¹⁸O records, the closest regional stack may still not be representative of a unique local record. Limited data availability has so far inhibited the development of a regional stack suitable for LCDW in the South Pacific (e.g., Stern and Lisiecki, 2014). Further, high-resolution benthic foraminiferal δ¹⁸O records, such as that of IODP Site U1541, may include sub-orbital-scale isotopic features that are not resolved in currently available Pleistocene stacks due to the relatively lower temporal resolution of most of the records used to construct the stacks.

4.2 Age models and alignment outputs for ODP Site 1090

Over the last 1.4 Myr, the manual and automated HMM-Match-based alignments of the benthic foraminiferal δ¹⁸O record from ODP Site 1090 to the LR04 stack yield similar chronostratigraphies with an average of the absolute values of age offsets equal to 3.7 kyr between them (Fig. 5). Our manual alignment of the benthic foraminiferal δ¹⁸O record at ODP Site 1090 resulted in different user-defined tie points and estimated stratigraphic hiatus durations compared to the original tie points defined by Venz and Hodell (2002) and used for segmentation of the input record for automated HMM-Match-based alignment (Table 1). As a consequence, the age offsets between the manual and automated alignments of the benthic foraminiferal δ¹⁸O at ODP Site 1090 are larger in the older part of the record (2.9 to 1.4 Ma), where multiple stratigraphic hiatuses are thought to occur (Venz and Hodell, 2002; Fig. 5). The maximum age difference between the manual versus HMM-Match-based alignment methods of benthic foraminiferal δ¹⁸O for ODP Site 1090 across this interval is 110 kyr.

Figure 5 Stratigraphic benthic foraminiferal δ¹⁸O alignments for ODP Site 1090. (a) Prob-stack benthic foraminiferal δ¹⁸O target (Ahn et al., 2017) for comparison. (b) Manual (light blue) and automated HMM-Match-based (dark blue; Lin et al., 2014) alignments of the benthic foraminiferal δ¹⁸O record from ODP Site 1090 to the LR04 benthic foraminiferal δ¹⁸O stack (Lisiecki and Raymo, 2005). (c) Age model offset between the age models for ODP Site 1090 plotted in (b). (d) Automated alignments of the benthic foraminiferal δ¹⁸O record from ODP Site 1090 to Prob-stack (black), CENOGRID (gray; Westerhold et al., 2020), and the LR09 Atlantic stack (brown; Lisiecki and Raymo, 2009) using HMM-Match. (e) 1σ age model uncertainties reported by HMM-Match for benthic foraminiferal δ¹⁸O alignment of ODP Site 1090 to Prob-stack. (f) Age offsets between the age models at ODP Site 1090 generated by the automated HMM-Match-based alignments to Prob-stack and LR04 (blue), Prob-stack and CENOGRID (gray), and Prob-stack and LR09 Atlantic stack (brown). The zero line is plotted in (c), (e), and (f) for reference.

The HMM-Match alignment of the benthic foraminiferal δ¹⁸O record at ODP Site 1090 to the Prob-stack target generates an orbital-scale age model over the past 2.9 Myr with an average reported 1σ age uncertainty of ∼2 kyr and a maximum reported 1σ age uncertainty of 11 kyr at 2870 ka (Fig. 5). These HMM-Match estimates of alignment uncertainties at ODP Site 1090 likely underestimate the true alignment uncertainty in the bottom portion of the record (2.9 to 1.4 Ma), where the HMM-Match algorithm assigns zero age uncertainty by default to the user-defined tie points at the start and stop of each hiatus, i.e., the end points of each segment input into the alignment algorithm.

Over the past 800 kyr, the average magnitudes of absolute age offsets between the HMM-Match alignment of the benthic foraminiferal δ¹⁸O record at ODP Site 1090 to Prob-stack and similar alignments to the LR04, LR09 Atlantic, and CENOGRID tuning targets are 3.8, 3.5, and 5.7 kyr, respectively. The largest target-based age offset in the HMM-Match alignments of ODP Site 1090 is 26 kyr and occurs between the alignment to Prob-stack and the alignment to LR04 at 2870 ka (Fig. 5). This large age offset coincides with the period of highest 1σ age uncertainty in the HMM-Match-based alignment to Prob-stack for ODP Site 1090.

The HMM-Match-based Alignment B to Prob-stack, using four segments rather than seven to deliberately skip known hiatuses in the sedimentary record, provides a broader perspective on the true uncertainties associated with the benthic-foraminiferal-δ¹⁸O-derived age models at ODP Site 1090 (Fig. 6). Specifically, Alignment B was generated without user-defined constraints on three mid-section hiatuses between 50.74 and 56.08 mcd (meters composite depth; Table 1). In the absence of user-defined flags of these hiatuses into HMM-Match, Alignment B yields large age model uncertainties for ODP Site 1090 between 2.7 and 1.9 Ma, reaching up to 60 kyr between 2.4 and 2.3 Ma (Fig. 6). While the reported age model uncertainties for the standard HMM-Match alignment to Prob-stack (based on seven segments) at ODP Site 1090 are typically less than 10 kyr across this interval (Fig. 5), the age model offset between the standard alignment and Alignment B exceeds 160 kyr at 2.2 Ma (Fig. 6). This age model offset, a result of undiagnosed stratigraphic discontinuity, is the largest age model offset observed across all alignments investigated in this study.

Figure 6 Comparison of age models at ODP Site 1090 generated by automated HMM-Match-based benthic foraminiferal δ¹⁸O alignments to Prob-stack using different segmentation. (a) HMM-Match-based benthic foraminiferal δ¹⁸O alignments to Prob-stack based on seven user-defined segments (blue; i.e., the standard 1090 HMM-Match alignment as plotted in Fig. 5) versus four user-defined stratigraphic segments (orange; ODP Site 1090 Alignment B) compared to the Prob-stack target (black). (b) 1σ age model uncertainties reported by the automated HMM-Match alignment algorithm for the standard alignment and Alignment B of ODP Site 1090 to Prob-stack. (c) Age model offset between the standard alignment and Alignment B of ODP Site 1090 to Prob-stack.

4.3 Age models and alignment outputs for ODP Site 980/981

As observed at IODP Site U1541 and ODP Site 1090, the manual and automated HMM-Match-based benthic foraminiferal δ¹⁸O alignments of ODP Site 980/981 to the LR04 stack yield broadly similar chronostratigraphies on the orbital scale (Fig. 7). The absolute difference in age models between these two alignment approaches reaches a maximum value of 24 kyr at 650 ka (leading into MIS 16) and has an average value of 2.8 kyr across the 1.8 Myr record at this site. The maximum age model offset at MIS 16 coincides with an interval of high temporal resolution in the benthic foraminiferal δ¹⁸O record at ODP Site 980/981, with one sample every 0.6 kyr. Such an age offset may thus result from the inability of conventional manual alignment to resolve rapid changes in sedimentation rate across intervals with limited (i.e., smoothed) benthic foraminiferal δ¹⁸O variability in the tuning target. This observation demonstrates that age uncertainties associated with manual tuning cannot always be reduced by increasing the temporal resolution of the input data and suggests that a smoothing of a high-resolution input record prior to alignment may be beneficial for tuning.

The HMM-Match-based benthic foraminiferal δ¹⁸O alignment to the Prob-stack target at ODP Site 980/981 generates an orbital-scale age model over the past 1.8 Myr with an average and maximum reported 1σ age uncertainty of 0.8 and 5.2 kyr, respectively (Fig. 7). In comparison, the magnitudes of absolute age offsets between the age models generated by automated HMM-Match-based benthic foraminiferal δ¹⁸O alignments to the LR04, LR09 Atlantic, and CENOGRID targets at ODP Site 980/981 are larger, with a maximum value of 18 kyr between the alignments to Prob-stack and LR04 at 1.42 Ma (Fig. 7). Over the last 800 kyr, the average absolute age offsets between automated HMM-Match-based benthic foraminiferal δ¹⁸O alignments to Prob-stack and LR04, LR09 Atlantic, and CENOGRID at ODP Site 980/981 are 2.0, 1.4, and 4.7 kyr, respectively (Fig. 7).

Figure 7 Stratigraphic benthic foraminiferal δ¹⁸O alignments for ODP Site 980/981. (a) Prob-stack benthic foraminiferal δ¹⁸O target (Ahn et al., 2017) for comparison. (b) Manual (light blue) and automated HMM-Match-based (dark blue; Lin et al., 2014) alignments of the benthic foraminiferal δ¹⁸O record from ODP Site 980/981 to the LR04 benthic foraminiferal δ¹⁸O stack (Lisiecki and Raymo, 2005). (c) Age model offset between the age models for ODP Site 980/981 plotted in (b). (d) Automated alignments of the benthic foraminiferal δ¹⁸O record from ODP Site 980/981 to Prob-stack (black), CENOGRID (gray; Westerhold et al., 2020), and the LR09 Atlantic stack (brown; Lisiecki and Raymo, 2009) using HMM-Match. (e) 1σ age model uncertainties reported by HMM-Match for the benthic foraminiferal δ¹⁸O alignment of ODP Site 980/981 to Prob-stack. (f) Age offsets between the age models at ODP Site 980/981 generated by the automated HMM-Match-based alignments to Prob-stack and LR04 (blue), Prob-stack and CENOGRID (gray), and Prob-stack and LR09 Atlantic stack (brown). The zero line is plotted in (c), (e), and (f) for reference.

5 Discussion

Based on the benthic foraminiferal δ¹⁸O alignment results presented above for IODP Site U1541, ODP Site 1090, and ODP Site 980/981, we discuss and synthesize the largest sources of uncertainty in benthic-foraminiferal-δ¹⁸O-based chronostratigraphies and suggest best community practices for the application of this approach.

5.1 Manual vs. probabilistic alignment

The probabilistic stratigraphic alignment algorithm HMM-Match (Lin et al., 2014) provides three advantages over manual alignment approaches for the development of benthic foraminiferal δ¹⁸O (chrono)stratigraphies: (1) it reduces the amount of user subjectivity in assigning tie points, (2) it generates estimates of age uncertainty associated with the alignment method itself, and (3) by aligning every data point within an input benthic foraminiferal δ¹⁸O record, HMM-Match estimates high-resolution sedimentation rate changes throughout an entire record (Fig. 8). These sedimentation rate changes are more realistic than the limited average sedimentation rate estimates between user-defined tie points that result from a manual alignment of benthic foraminiferal δ¹⁸O to a tuning target. In addition, linear sedimentation rate records with high temporal resolution can provide useful information on the relationship between environmental and climatic variability in sediment accumulation, dissolution, and redistribution processes on the seafloor (e.g., Suman and Bacon, 1989; Costa et al., 2020; Francois et al., 1990).

The age offsets generated between manual and automated probabilistic (HMM-Match-based) alignments of the benthic foraminiferal δ¹⁸O record to the LR04 benthic stack at IODP Site U1541 (Fig. 4), ODP Site 1090 (Fig. 5), and ODP Site 980/981 (Fig. 7) provide a means to start to quantitatively assess the magnitude of alignment uncertainties incurred from manual alignments to user-defined tie points. These uncertainties typically average between 3 and 5 kyr across the three sites but can reach values of 40 kyr or higher near segment ends (e.g., at 2.7 Ma at ODP Site 1090) and intervals where local benthic foraminiferal δ¹⁸O changes differ from the reference target (e.g., at 3.0 and 3.4 Ma at IODP Site U1541). This range of uncertainty is comparable to the 2.5 kyr tuning uncertainty estimated by Martinson et al. (1987) based on the manual alignment of five paleoclimatological indicators from a single site (planktonic and benthic foraminiferal δ¹⁸O, radiolarian abundance, sea surface temperature, stable carbon isotopes, and carbonate content) to a stacked δ¹⁸O target. As observed at ODP Site 980/981, high temporal resolution of the benthic foraminiferal δ¹⁸O record does not guarantee a reduction in this form of alignment uncertainty (Sect. 4.3).

Figure 8 Linear sedimentation rate (LSR) estimates for IODP Site U1541, ODP Site 1090, and ODP Site 980/981. Values are based on manual alignment (blue) and automated HMM-Match alignment (gray) of benthic foraminiferal δ¹⁸O records to the LR04 benthic stack (Lin et al., 2014; Lisiecki and Raymo, 2004).

The age model uncertainties generated by HMM-Match for the benthic foraminiferal δ¹⁸O alignments at IODP Site U1541 (Fig. 4), ODP Site 1090 (Fig. 5), and ODP Site 980/981 (Fig. 7) include uncertainties associated with the reported benthic foraminiferal δ¹⁸O variance in the target reference record and the stratigraphic alignment method itself (i.e., the process of “wiggle matching”; Lin et al., 2014). Uncertainties in the latter are calculated under the assumptions of reasonable stratigraphic continuity and that user-defined ages for each segment start and end are certain (Lin et al., 2014). These values average ∼1 to 3 kyr across the three example sites and typically decrease with increased temporal resolution of the undated benthic foraminiferal δ¹⁸O record and of the tuning target (which are generally better constrained in the late Pleistocene than in the Pliocene).

The point-by-point alignment uncertainties estimated by HMM-Match are valuable for assessing changes in the confidence level of benthic foraminiferal δ¹⁸O chronostratigraphies across multi-million-year records, where long-term changes in the seawater δ¹⁸O signal or temporal resolution of the input record and tuning target may cause systematic changes in age uncertainty. Although HMM-Match assigns unrealistically low age uncertainties at user-defined tie points for the start and end of an alignment segment, this issue has been resolved in more recent evolutions of automated alignment algorithms (Lee et al., 2023). For example, the new automated alignment software package, BIGMACS, allows users to specify the uncertainties of added tie points to generate more realistic estimates of alignment uncertainty (Lee et al., 2023). However, neither HMM-Match-based nor BIGMACS-based alignment age uncertainties account for the uncertainty associated with the absolute chronology of the target reference, regional variability in bottom-water temperature, and/or the timing of local changes in seawater (and resulting benthic foraminiferal) δ¹⁸O between the input record and the selected tuning target.

While automated alignment algorithms like HMM-Match and BIGMACS provide many advantages in the generation of benthic foraminiferal δ¹⁸O chronostratigraphies, we note that depositional environments with highly variable sedimentation rate changes and benthic foraminiferal δ¹⁸O records with long data gaps or temporal variations in sampling resolution may not abide by the assumptions included in these automated algorithms. Specifically, the HMM-Match algorithm is designed to minimize large sedimentation rate changes between data points based on the probability distribution calculated from a compilation of 37 radiocarbon-dated sediment cores (Lin et al., 2014). This guiding principle may hinder the success of HMM-Match-generated alignments for benthic foraminiferal δ¹⁸O records with irregular sampling frequencies or from regions like the Antarctic Southern Ocean where very large sedimentation rate changes are expected across a deglaciation (e.g., Hasenfratz et al., 2019). Consequently, a close visual evaluation of the automated alignment outputs against their designated targets should be completed manually (i.e., by the user) to quality-check the resulting benthic foraminiferal δ¹⁸O chronostratigraphies.

5.2 Choice of tuning target and the assumption of global synchronicity

We encourage the utilization of Prob-stack as a global tuning target because Prob-stack includes the most holistic estimate to date of spatial variance in benthic foraminiferal δ¹⁸O values across the last 5 Myr, under the assumption of global synchronicity in benthic foraminiferal δ¹⁸O variability, and is therefore best suited for calculating more realistic estimates of chronostratigraphic alignment uncertainties (Ahn et al., 2017). Our analyses reveal an average age offset of 2–6 kyr between benthic foraminiferal δ¹⁸O chronostratigraphies generated for IODP Site U1541, ODP Site 1090, and ODP Site 980/981 alignments to Prob-stack, LR04, and CENOGRID using HMM-Match (Figs. 4, 5, and 7, respectively). These values suggest that the choice of tuning target, or comparing records based on alignments to different tuning targets, can introduce a level of age uncertainty comparable to the uncertainty associated with the alignment itself (Sect. 5.1). The magnitude of these age offsets is similar to that of regional differences in the timing of benthic foraminiferal δ¹⁸O changes in the global ocean across glacial–interglacial cycles (i.e., up to 4 kyr; Lisiecki and Raymo, 2009; Skinner and Shackleton, 2005; Stern and Lisiecki, 2014; Waelbroeck et al., 2011; Layberie et al., 2005).

The age model uncertainties caused by regional asynchronicity in temperature and/or seawater (i.e., benthic foraminiferal) δ¹⁸O changes across glacial cycles must also be considered when comparing event timing among benthic-foraminiferal-δ¹⁸O-tuned sediment records from different basins and water masses. The age offsets observed between the benthic foraminiferal δ¹⁸O alignments at IODP Site U1541, ODP Site 1090, and ODP Site 980/981 to Prob-stack and their respective alignments to the LR09 Pacific and Atlantic regional stacks provide a qualitative assessment of the impact of regional differences on age model uncertainties and average 1.4 to 3.5 kyr over the past 800 kyr.

Regional stacks of individual benthic foraminiferal δ¹⁸O records with high-resolution radiocarbon dates over the last glacial cycle suggest that the spatial variability in seawater δ¹⁸O and its response to changing climate may induce systematic age offsets up to 1.7 kyr between the deep South Pacific IODP Site U1541 record and ODP Sites 1090 and 980/981 in the deep South Atlantic and deep North Atlantic, respectively, during glacial terminations (Stern and Lisiecki, 2014). Significant age offsets, up to 1.0 and 3.5 kyr, should be considered when comparing the new IODP Site U1541 data to records from the intermediate Pacific and the deep Indian basin, respectively (Stern and Lisiecki, 2014). The development and refinement of a deep Pacific regional stack extending beyond the last glacial cycle would help further define the uncertainties associated with asynchronous δ¹⁸O changes in the global ocean.

While age uncertainties associated with the absolute chronology of the alignment target are less of a concern for the comparison of relative event timing between sedimentary records with similarly stratigraphically determined age models, estimates of such age uncertainties are critical for comparison with terrestrial or independently dated marine records. The absolute age uncertainty associated with the LR04 stack is estimated to yield age model uncertainties of 4 kyr for 0–1 Ma, 6 kyr for 1–3 Ma, 15 kyr from 3–4 Ma, and 30 kyr or higher for 4–5.3 Ma (Lisiecki and Raymo, 2005). The absolute age uncertainty associated with the LR09 stacks and Prob-stack is expected to be similar to the values for LR04, as these stacks are based on the same time domain as LR04 (Lisiecki and Raymo, 2009; Ahn et al., 2017). Absolute age uncertainties are not specifically reported for the 0–5 Ma interval of the CENOGRID composite based on astronomical tuning of a number of high-resolution benthic foraminiferal δ¹⁸O records from Ceara Rise (i.e., ODP Sites 925, 926, 927, 928, and 929); however, a broader absolute age uncertainty estimate of 10 kyr is reported for the late Miocene to Pleistocene interval of CENOGRID (Wilkens et al., 2017; Westerhold et al., 2020). Including the absolute age uncertainties of the tuning targets in stratigraphic benthic foraminiferal δ¹⁸O alignments would significantly increase the total reported age uncertainties of marine sedimentary records.

5.3 Stratigraphic discontinuity

The largest source of age uncertainties identified by our analysis arises from stratigraphic discontinuities within the input benthic foraminiferal δ¹⁸O record, as shown by ODP Site 1090 (Figs. 5 and 6). Coring gaps, extended sampling breaks, and sedimentary hiatuses within an undated record invalidate the assumption of stratigraphic continuity on which benthic foraminiferal δ¹⁸O chronostratigraphies are built and require a sound assessment of potential hiatuses in the sedimentary record by the user. The uncertainties associated with these stratigraphic discontinuities are higher during the 40 kyr glacial cycles of the Pliocene and early Pleistocene, when there are lower-amplitude and higher-frequency benthic foraminiferal δ¹⁸O changes than during the late Pleistocene. Missing cycles are, accordingly, harder to identify within this interval (Ahn et al., 2017; Lisiecki and Raymo, 2005; Westerhold et al., 2020).

Although the HMM-Match algorithm did not independently identify the specific mid-section hiatuses reported by Venz and Hodell (2002) in the benthic foraminiferal δ¹⁸O record at ODP Site 1090 in our sensitivity test, the algorithm reported large uncertainty estimates (reaching up to ∼60 kyr) across this interval of discontinuity (Fig. 6). The increase in HMM-Match-generated alignment uncertainties thereby flags intervals of suspected stratigraphic discontinuity at ODP Site 1090 and demonstrates the utility of this algorithm to detect intervals of stratigraphic complication within a long sedimentary record (Fig. 6). Such hiatuses may otherwise be difficult to detect and could introduce tens of thousands of years of systematic error into benthic-foraminiferal-δ¹⁸O-derived sediment chronologies. We caution, however, that HMM-Match has been documented to generate alignment errors and underestimated age uncertainties in high-resolution benthic foraminiferal δ¹⁸O records with millennial-scale sampling gaps along the Iberian Margin (Lisiecki et al., 2022). Interruptions in the stratigraphic integrity of the sediment record under investigation can be assessed via core descriptions, imagery, XRF scanning data, or other high-resolution indicators of sediment stratigraphy in order to make an informed refinement of the stratigraphic alignment beyond age uncertainties reported by HMM-Match.

One may hypothesize, based on the relatively high age uncertainties for the automated HMM-Match-based benthic foraminiferal δ¹⁸O alignment to Prob-stack at IODP Site U1541 between 2190 and 2110 ka (reaching up to 20 kyr), that the benthic foraminiferal δ¹⁸O record of this MIS 82 interval may be disrupted by a mid-section hiatus (Fig. 4). However, no such stratigraphic discontinuities are observed in the high-resolution (1–2 cm) shipboard density or color records at IODP Site U1541 across this interval (Winckler et al., 2021). Rather, this section may reflect an interval where local changes in benthic foraminiferal δ¹⁸O values deviate significantly from the global stack.

5.4 Suggestions for best practice

Our HMM-Match analyses of benthic foraminiferal δ¹⁸O alignments at IODP Site U1541, ODP Site 1090, and ODP Site 980/981 to common tuning targets suggest combined age uncertainties for the HMM-Match-based alignments to Prob-stack to be ±5–6 kyr for 0 to 1 Ma, ±7–9 kyr for 1 to 3 Ma, and ±18 kyr for 3 to 3.5 Ma. These values include age uncertainties associated with the alignment process, benthic foraminiferal δ¹⁸O variance in the tuning target, and absolute age constraints on the tuning target (e.g., the absolute age uncertainties reported by Lisiecki and Raymo, 2005). The largest contribution to these uncertainty estimates is the uncertainty in absolute ages of the target record (Sect. 5.2, Lisiecki and Raymo, 2005). Continued community efforts to provide independent age control on the absolute ages of benthic foraminiferal δ¹⁸O tuning targets, especially for sediments older than 3 Ma, will be critical for the further reduction of total uncertainty associated with benthic foraminiferal δ¹⁸O chronostratigraphy.

In light of the existing caveats of benthic foraminiferal δ¹⁸O alignments, we suggest the following best practices for stratigraphically aligning benthic foraminiferal δ¹⁸O records to common targets to minimize alignment uncertainty and ensure consistent characterization of age model uncertainty during the Plio-Pleistocene.

1. Investigate possible disruptions in stratigraphic continuity using sediment core images, XRF, or other high-resolution scan data.

2. Utilize an automated and probabilistic alignment algorithm, like HMM-Match (Lin et al., 2014), to reduce user subjectivity in record alignment and to generate point-by-point estimates of age uncertainties and sedimentation rate variations resulting from the alignment.

3. Visually assess the quality of alignment outputs generated by automated algorithms for alignment mismatches or other irregularities, especially across glacial terminations.

4. Align benthic foraminiferal δ¹⁸O records to stacked records that include characterization of global or regional heterogeneity in δ¹⁸O variance (e.g., Prob-stack; Ahn et al., 2017; or an appropriate regional stack, if available) to better constrain uncertainties associated with the alignment itself.

5. Acknowledge the amount by which benthic foraminiferal δ¹⁸O changes at the study site may lead or lag the tuning target of choice. The development of updated and long-term regional stacks will help reduce these associated age uncertainties.

6. When working with multiple sediment cores with age models based on stratigraphic benthic foraminiferal δ¹⁸O alignments, ensure the age models for each record are calculated using the same approaches and tuning targets.

7. Publish sample depth data along with sample ages and measurement values so that it is transparent how published records may be updated or realigned to newer age models.

On a broader note, given the financial burden of programming platforms like MATLAB, we encourage the development or translation of automated alignment algorithms into license-free coding languages such as Python or R so that they may become more accessible to the research community.

6 Conclusions

Stratigraphic alignments based on new and published benthic foraminiferal δ¹⁸O records of IODP Site U1541, ODP Site 1090, and ODP Site 980/981 provide valuable insights into the types and range of sediment age model uncertainties associated with benthic foraminiferal δ¹⁸O chronostratigraphies. Our results suggest typical age uncertainties of 3 to 5 kyr for manual alignment efforts and 1 to 3 kyr for probabilistic HMM-Match-based alignments (Lin et al., 2014) over the past 3.5 Myr. We identify an average age offset of 2–6 kyr between stratigraphic benthic foraminiferal δ¹⁸O alignments to Prob-stack (Ahn et al., 2017), LR04 (Lisiecki and Raymo, 2005), and the CENOGRID target (Westerhold et al., 2020). The 2–6 kyr age uncertainty associated with tuning target choice is comparable to the estimated 1–4 kyr range of lag time in benthic foraminiferal δ¹⁸O changes between different ocean basins and water masses (Lisiecki and Raymo, 2009; Skinner and Shackleton, 2005; Stern and Lisiecki, 2014; Waelbroeck et al., 2011; Layberie et al., 2005). Stratigraphic discontinuities within the benthic foraminiferal δ¹⁸O input record can contribute to large age uncertainties (reaching 60 kyr and higher) in marine sediment core chronostratigraphies.

The uncertainties associated with benthic foraminiferal δ¹⁸O chronostratigraphies and their impact on paleoceanographic interpretations can be reduced using best practices for age model uncertainty characterization and alignment. These practices include investigation into the stratigraphic continuity of an undated record, utilization of probabilistic alignment algorithms (such as HMM-Match or BIGMACS; Lin et al., 2014; Lee et al., 2023), alignment targets that suitably characterize heterogeneity in δ¹⁸O variance (such as Prob-stack globally – Ahn et al., 2017, or recent regional stacks like the North Atlantic stack of Hobart et al., 2023), and acknowledgement of the amount by which benthic foraminiferal δ¹⁸O changes in the study region may lead or lag the alignment target. Continued investigation into regional asynchronicity in seawater and benthic foraminiferal δ¹⁸O changes across the Plio-Pleistocene and independent constraints on the absolute ages of tuning target intervals, particularly older than 3 Ma, should also be supported to further reduce age uncertainties based on stratigraphic benthic foraminiferal δ¹⁸O alignments of marine sedimentary records to a common tuning target.