Lower oceanic δ¹³C during the last interglacial period compared to the Holocene

The last time in Earth's history when high latitudes were warmer than during pre-industrial times was the last interglacial period (LIG, 129–116 ka BP). Since the LIG is the most recent and best documented interglacial, it can provide insights into climate processes in a warmer world. However, some key features of the LIG are not well constrained, notably the oceanic circulation and the global carbon cycle. Here, we use a new database of LIG benthic δ¹³C to investigate these two aspects. We find that the oceanic mean δ¹³C was ∼ 0.2 ‰ lower during the LIG (here defined as 125–120 ka BP) when compared to the Holocene (7–2 ka BP). A lower terrestrial carbon content at the LIG than during the Holocene could have led to both lower oceanic δ¹³C and atmospheric δ¹³CO₂ as observed in paleo-records. However, given the multi-millennial timescale, the lower oceanic δ¹³C most likely reflects a long-term imbalance between weathering and burial of carbon. The δ¹³C distribution in the Atlantic Ocean suggests no significant difference in the latitudinal and depth extent of North Atlantic Deep Water (NADW) between the LIG and the Holocene. Furthermore, the data suggest that the multi-millennial mean NADW transport was similar between these two time periods.

1 Introduction

The most recent and well documented warm time period is the last interglacial period (LIG), which is roughly equivalent to Marine Isotope Stage (MIS) 5e (Past Interglacials Working Group of PAGES, 2016; Shackleton, 1969). The LIG began at the end of the penultimate deglaciation and ended with the last glacial inception (∼ 129–116 ka BP; Dutton and Lambeck, 2012; Govin et al., 2015; Masson-Delmotte et al., 2013; Menviel et al., 2019). The LIG was globally warmer than the pre-industrial period (PI, ∼ 1850–1900; IPCC, 2013; Shackleton et al., 2020), with PI estimated to be ∼ 0.4 ^∘C cooler than the peak of the Holocene (10–5 ka BP) (Marcott et al., 2013). Though not an exact analogue for future warming, the LIG may still help shed light on future climates. In particular, we seek to constrain the mean LIG ocean circulation and estimate the global oceanic mean δ¹³C.

As greenhouse gas concentrations were comparable to the Holocene, the LIG was most likely relatively warm because of the high boreal summer insolation (Laskar et al., 2004). During the LIG, the atmospheric CO₂ concentration was relatively stable around ∼ 280 ppm (Bereiter et al., 2015; Lüthi et al., 2008), while during the Holocene CO₂ first decreased by about 8 ppm starting at 11.7 ka BP before increasing by ∼ 17 to 277 ppm at ∼ 2 ka BP (Fig. 1a) (Köhler et al., 2017). CH₄ and N₂O peaked at ∼ 700 and ∼ 267 ppb, respectively, during both the LIG and the Holocene (Flückiger et al., 2002; Petit et al., 1999; Spahni et al., 2005). Global sea level was 6–9 m higher at the LIG compared to PI (Dutton et al., 2015; Kopp et al., 2009), thus indicating significant ice mass loss from both Antarctica and Greenland.

Figure 1 LIG and Holocene time series of (a) CO₂ stack smoothed with a spline based on the age model AICC2012 (Köhler et al., 2017), (b) sea surface temperatures (SSTs) determined from alkenones and aligned with oxygen isotopes from the Iberian Margin (MD01-2444, blue, Martrat et al., 2007b) and the North Atlantic (GIK23414-6, green, Candy and Alonso-Garcia, 2018), (c) EPICA Dome C ice core (EDC96) deuterium measurements (orange) and estimated surface air temperature anomaly relative to the mean of the last 1 kyr (red, Bazin et al., 2013; Jouzel et al., 2007) on the AICC2012 timescale, and (d) spline of atmospheric δ¹³C from EPICA Dome C and the Talos Dome ice cores (Holocene, Eggleston et al., 2016) and Monte Carlo average of three Antarctic ice cores atmospheric δ¹³C (LIG, Schneider et al., 2013) both based on the age model AICC2012. Shading around the lines indicates 1σ. Vertical grey shading indicates the periods of analysis in this paper. Vertical dotted grey lines indicate the commencement of the LIG and Holocene.

Strong polar warming is supported by terrestrial and marine temperature reconstructions. A global analysis of sea surface temperature (SST) records suggests that the mean surface ocean was 0.5 ± 0.3 ^∘C warmer during the LIG compared to 1870–1889 (Hoffman et al., 2017), similar to another global estimate, which suggests SSTs were 0.7 ± 0.6 ^∘C higher during the LIG compared to the late Holocene (McKay et al., 2011). However, there were differences in the timing of these SST peaks in different regions compared to the 1870–1889 mean: North Atlantic SSTs peaked at +0.6 ± 0.5 ^∘C at 125 ka BP (e.g. Fig. 1b) and Southern Hemisphere extratropical SSTs peaked at +1.1 ± 0.5 ^∘C at 129 ka BP (Hoffman et al., 2017). On land, proxy records from mid-latitudes to high latitudes indicate higher temperatures during the LIG compared to PI, particularly in North America (Anderson et al., 2014; Axford et al., 2011; Montero-Serrano et al., 2011). Similarly, the EPICA DOME C record suggests that the highest Antarctic temperatures from the last 800 kyr occurred during the LIG (Masson-Delmotte et al., 2010) (Fig. 1c).

Polar warming was also associated with significant changes in vegetation. Pollen records suggest a contraction of tundra and an expansion of boreal forests across the Arctic (CAPE, 2006), in Russia (Tarasov et al., 2005), and in North America (Govin et al., 2015; Muhs et al., 2001; de Vernal and Hillaire-Marcel, 2008). The few Saharan records suggest a green Sahara period during the LIG (Drake et al., 2011; Larrasoaña et al., 2013), consistent with a stronger West African monsoon (Otto-Bliesner et al., 2021). Although these reconstructions indicate changes in vegetation distribution during the LIG, the total amount of carbon stored on land remains poorly constrained.

Recent numerical experiments of the LIG as part of the Paleomodel Intercomparison Project Phase 4 (PMIP4) simulate significant warming over Alaska and Siberia in boreal summer, with mean annual temperature anomalies of close to zero, which is in good agreement with the proxy record (Otto-Bliesner et al., 2021). Despite this and other recent data compilations and modelling efforts (including Bakker et al., 2013), there are many open questions remaining about the LIG. In particular, stronger constraints are needed on the extent of Greenland and Antarctic ice sheets; on ocean circulation and the global carbon cycle, including CaCO₃ accumulation in shallow waters; and peat and permafrost carbon storage (Brovkin et al., 2016).

It is important to constrain the state of the Atlantic Meridional Overturning Circulation (AMOC) at the LIG given its significant role in modulating climate. Seven coupled climate models integrated with transient 130–115 ka BP boundary conditions simulate different AMOC trends, with some models producing a strengthening of the AMOC, while others simulate a weakening during the LIG (Bakker et al., 2013). Paleoproxy records suggest equally strong and deep North Atlantic Deep Water (NADW) during the LIG and the Holocene (e.g. Böhm et al., 2015; Lototskaya and Ganssen, 1999), with a possible southward expansion of the Arctic front related to changes in the strength of the subpolar gyre (Mokeddem et al., 2014), and AMOC weakening during a few multi centennial-scale events between 127 and 115 ka BP (e.g. Galaasen et al., 2014b; Helmens et al., 2015; Lehman et al., 2002; Mokeddem et al., 2014; Oppo et al., 2006; Rowe et al., 2019; Tzedakis et al., 2018).

Stable carbon isotopes are a powerful tool for investigating ocean circulation (e.g. Curry and Oppo, 2005; Eide et al., 2017) and the global carbon cycle (e.g. Menviel et al., 2017; Peterson et al., 2014). Since the largest carbon isotope fractionation occurs during photosynthesis, organic matter is enriched in ¹²C (low δ¹³C), while atmospheric CO₂ and surface water dissolved inorganic carbon (DIC) become enriched in ¹³C (high δ¹³C). Organic matter on land includes the terrestrial biosphere, as well as carbon stored in soils, such as in peats and permafrost. Different photosynthetic pathways (which differentiate C₃ and C₄ plants) fractionate carbon differently, producing typical signatures of about −37 ‰ to −20 ‰ for C₃ plants (Kohn, 2010) and around −13 ‰ for C₄ plants (Basu et al., 2015), though these values vary with a number of factors, including precipitation, atmospheric CO₂ concentration and δ¹³C, light, nutrient availability, and plant species (Cernusak et al., 2013; Diefendorf et al., 2010; Diefendorf and Freimuth, 2017; Farquhar, 1983; Farquhar et al., 1989; Keller et al., 2017; Leavitt, 1992; Schubert and Jahren, 2012). In the ocean, phytoplankton using the C₃ photosynthetic pathway are found to have fractionation during photosynthesis that depends on the concentration of dissolved CO₂. Thus, atmospheric δ¹³CO₂ during the LIG (Fig. 1d) is influenced by the cycling of organic carbon within the ocean, changes in the amount of carbon stored in vegetation and soils, temperature-dependent air–sea flux fractionation (Lynch-Stieglitz et al., 1995; Zhang et al., 1995), and on longer timescales by interactions with the lithosphere (Tschumi et al., 2011). The mean surface DIC is enriched by ∼ 8.5 ‰ compared to the atmosphere due to fractionation during air–sea gas exchange (Menviel et al., 2015; Schmittner et al., 2013).

NADW is characterised by low nutrients and high δ¹³C as a result of a high nutrient and carbon utilisation by marine biota and fractionation during air–sea gas exchange in the northern North Atlantic. Along its path through the Atlantic basin interior, organic matter remineralisation and mixing with southern source waters lowers δ¹³C, with δ¹³C values of ∼ 0.5 ‰ in the deep Southern Ocean.

The tight relationship between the water masses' apparent oxygen utilisation, nutrient content and δ¹³C allows δ¹³C to be used as a water mass ventilation tracer (e.g. Boyle and Keigwin, 1987; Curry and Oppo, 2005; Duplessy et al., 1988; Eide et al., 2017). The δ¹³C of benthic foraminifera shells, particularly of the species Cibicides wuellerstorfi, has been found to reliably represent the δ¹³C signature of DIC (Belanger et al., 1981; Duplessy et al., 1984; Zahn et al., 1986) and has therefore been used to better constrain the extent of different water masses. Mass balances of δ¹³C between the atmosphere, ocean, and land have been previously used to constrain changes in terrestrial carbon between the last glacial maximum (∼ 20 ka BP) and Holocene (e.g. Peterson et al., 2014). However, on longer timescales, exchanges with the lithosphere including volcanic outgassing (Hasenclever et al., 2017; Huybers and Langmuir, 2009), CaCO₃ burial in sediments and weathering, release of carbon from methane clathrates, and the net burial of organic carbon also influence the global mean δ¹³C. It has been estimated that the amount of carbon both entering and exiting the lithosphere due to weathering and burial of organic carbon fluxes could be from 0.274 to 0.344 Gt C yr⁻¹ (Schneider et al., 2013), though these vary through time (Hoogakker et al., 2006). Over timescales greater than 10 kyr, the influence of weathering and burial of carbon might dominate the δ¹³C signal (Jeltsch-Thömmes et al., 2019; Jeltsch-Thömmes and Joos, 2020), and thus a mass balance cannot be accurately applied to evaluate terrestrial carbon changes between the LIG and Holocene.

Here, we present a new compilation of benthic δ¹³C from Cibicides wuellerstorfi spanning the 130–118 ka BP time period. We use this data to compare the δ¹³C signal of the LIG with that of the Holocene and to determine the difference in average ocean δ¹³C between the two time periods. We then investigate the AMOC during the LIG with our new benthic δ¹³C database. Finally, we qualitatively explore the role of the various processes affecting the δ¹³C difference between the LIG and the Holocene.

2 Database and methods

2.1 Database

We present a new compilation of benthic δ¹³C covering the periods 130–118 and 8–2 ka BP. From these two sets of data, we select data pertaining to the LIG and compare it to data from the Holocene. Our database only includes measurements on Cibicides wuellerstorfi as no significant fractionation between the calcite shells and the surrounding DIC has been measured in this species (Belanger et al., 1981; Duplessy et al., 1984; Zahn et al., 1986).

Chen et al. (1995)Cortijo (2003)Members (2006)McIntyre et al. (1999)Schmiedl and Mackensen (2006)Chapman and Shackleton (1999)Holbourn et al. (2005)Hodell et al. (2008)Lyle et al. (2002)McManus et al. (1999)Oppo et al. (1998)Shackleton et al. (1990)Jansen et al. (1996)Shackleton et al. (1990)Curry and Oppo (1997)Duplessy et al. (1984)Bickert and Mackensen (2003)Lyle et al. (2002)Pisias and Mix (1997)Deaney et al. (2017)Zhang et al. (2007)Waelbroeck et al. (2001)Wang et al. (1999)Jullien et al. (2006)Labeyrie et al. (1995)Lee et al. (1999)Wei et al. (2006)Hodell and Channell (2016)Cheng et al. (2004)Ruddiman and Members (1982)Duplessy et al. (1984)Galaasen et al. (2014a)Duplessy et al. (1984)Bickert et al. (1997)Mackensen et al. (2001)Curry et al. (1995)Hodell et al. (2001)Bickert et al. (1997)McCorkle and Holder (2001)Oppo and Fairbanks (1987)Russon et al. (2009)Curry and Oppo (2005)Mix et al. (1991)Bickert and Mackensen (2003)Shackleton et al. (1995)Mix et al. (1991)Schmiedl and Mackensen (1997)Zahn et al. (1986)Sarnthein et al. (1994)Shackleton (1977)Zahn et al. (1986)Bickert and Mackensen (2003)Sarnthein et al. (1994)Shackleton (1977)Freudenthal et al. (2002)Sarnthein et al. (1994)Sarnthein et al. (1994)Duplessy (1996)Sarnthein et al. (1994)Sarnthein et al. (1994)Chapman and Shackleton (1999)Sarnthein (2003)Sarnthein et al. (1994)Sarnthein et al. (1994)Sarnthein et al. (1994)Bickert and Mackensen (2003)Curry et al. (1988)Sarnthein et al. (1994)Lynch-Stieglitz et al. (2006)Sarnthein et al. (1994)Rau et al. (2002)Ruddiman and Members (1982)Vidal et al. (1999)Raymo et al. (2004)Ruddiman and Members (1982)Oppo and Lehman (1995)Raymo et al. (1997)

Table 1 List of cores for the last interglacial period (LIG). Provided is the core name (“Core”), latitude (“Lat”, ^∘), longitude (“Long”, ^∘), depth (“Dep”, m), the region, and the reference. Regions are abbreviated as follows: NEA: northeastern Atlantic; NWA: northwestern Atlantic; SWA: southwestern Atlantic; SEA: southeastern Atlantic; SA: southern Atlantic; NP: northern Pacific; SP: southern Pacific; I: Indian. Reference abbreviations are as follows: BW96: Bickert and Wefer (1996); CL82: Curry and Lohmann (1982); dA03: de Abreu et al. (2003); KJ8994: Keigwin and Jones (1989, 1994); KS02: Keigwin and Schlegel (2002); L99: Labeyrie et al. (1999); MB99: Mackensen and Bickert (1999); OH00: Oppo and Horowitz (2000); SH84: Shackleton and Hall (1984); SS0405: Skinner and Shackleton (2004, 2005); VH02: Venz and Hodell (2002); V99: Venz et al. (1999); ZM1011: Zarriess and Mackensen (2010, 2011).

Chen et al. (1995)Sarnthein et al. (1994)Schmiedl and Mackensen (2006)Sarnthein et al. (1994)Sirocko et al. (2000)Holbourn et al. (2005)Sarnthein et al. (1994)Curry et al. (1988)Millo et al. (2006)Raymo et al. (2004)Labeyrie et al. (1996)Sowers et al. (1993)Oppo and Lehman (1995)Wang et al. (1999)Martrat et al. (2007a)Duplessy et al. (1984)Cortijo (2003)Lyle et al. (2002)Pisias and Mix (1997)McIntyre et al. (1999)Zhang et al. (2007)Sarnthein et al. (1994)Wang et al. (1999)Sarnthein et al. (1994)Wang et al. (1999)Mix and Fairbanks (1985)Lee et al. (1999)Wei et al. (2006)Lebreiro et al. (2009)Wang et al. (1999)Chapman and Shackleton (1999)Duplessy et al. (1984)Sarnthein et al. (1994)Stott et al. (2000)Sarnthein et al. (1994)Gebhardt et al. (2008)Oppo et al. (1998)McManus et al. (1999)Hodell et al. (2000)Hodell et al. (2003)Jansen et al. (1996)Hodell et al. (2001)Duplessy et al. (1984)McCorkle and Holder (2001)Curry et al. (1988)Waelbroeck et al. (2011)Curry et al. (1988)Russon et al. (2009)Curry and Oppo (1997)Mix et al. (1991)Bickert and Mackensen (2003)Shackleton et al. (1995)Keigwin et al. (1991)Mix et al. (1991)Waelbroeck et al. (2001)Samson et al. (2005)Hodell and Channell (2016)Sikes et al. (2009)Bickert et al. (1997)Pahnke and Zahn (2005)Mix and Fairbanks (1985)Curry et al. (1995)Sarnthein et al. (1994)Curry et al. (1988)Sarnthein et al. (1994)Curry et al. (1988)Sarnthein et al. (1994)Curry et al. (1988)Sarnthein et al. (1994)Curry et al. (1988)Boyle (1992)Curry et al. (1988)Sarnthein et al. (1994)Curry et al. (1988)Sarnthein et al. (1994)Bickert et al. (1997)Boyle (1992)Bickert et al. (1997)Zahn et al. (1986)Bickert et al. (1997)Sarnthein et al. (1994)Oppo and Fairbanks (1987)Zarriess et al. (2011)Hüls (1999)Zahn and Stüber (2002)Sarnthein et al. (1994)Curry and Oppo (2005)Mulitza et al. (2008)Curry et al. (1988)Sarnthein et al. (1994)Collins et al. (2011)Tjallingii et al. (2008)Sarnthein et al. (1994)Schmiedl and Mackensen (1997)Sarnthein et al. (1994)Bickert et al. (2003)Sarnthein et al. (1994)Shackleton (1977)Zahn et al. (1986)Freudenthal et al. (2002)Bickert and Mackensen (2003)Sarnthein et al. (1994)Freudenthal et al. (2002)Sarnthein et al. (1994)Sarnthein et al. (1994)Sarnthein et al. (1994)Sarnthein et al. (1994)Sarnthein et al. (1994)Curry et al. (1988)Skinner et al. (2003)Sarnthein et al. (1994)Schönfeld et al. (2003)Sarnthein et al. (1994)Boyle and Keigwin (1985)Bickert and Mackensen (2003)Sarnthein et al. (1994)Lynch-Stieglitz et al. (2006)Duplessy (1996)Rau et al. (2002)Sarnthein et al. (1994)Vidal et al. (1999)Hodell et al. (2008)Came et al. (2003)Sarnthein et al. (1994)Sarnthein et al. (1994)Sarnthein et al. (1994)Sarnthein et al. (1994)Sarnthein et al. (1994)

Table 2 List of cores for the Holocene. Provided is the core name (“Core”), latitude (“Lat”, ^∘), longitude (“Long”, ^∘), depth (“Dep”, m) the region, and the reference. Regions are abbreviated as follows: NEA: northeastern Atlantic; NWA: northwestern Atlantic; SWA: southwestern Atlantic; SEA: southeastern Atlantic; SA: southern Atlantic; NP: northern Pacific; SP: southern Pacific; I: Indian. Reference abbreviations are as follows: BW96: Bickert and Wefer (1996); CL82: Curry and Lohmann (1982); dA03: de Abreu et al. (2003); KJ8994 Keigwin and Jones (1989, 1994); KS02: Keigwin and Schlegel (2002); L99: Labeyrie et al. (1999); MB99: Mackensen and Bickert (1999); OH00: Oppo and Horowitz (2000); SH84: Shackleton and Hall (1984); SS0405: Skinner and Shackleton (2004, 2005); VH02: Venz and Hodell (2002); V99: Venz et al. (1999); ZM1011: Zarriess and Mackensen (2010, 2011).

Our compilation is predominantly based on Lisiecki and Stern (2016) (53 cores) but includes 14 cores described in Oliver et al. (2010), as well as a few other records (CH69-K09, Labeyrie et al., 2017; MD03-2664, Galaasen et al., 2014a; MD95-2042, Martrat et al., 2007a; ODP 1063, Deaney et al., 2017; and U1304, Hodell and Channell, 2016). The full core lists are provided in Tables 1 and 2 for the LIG and the Holocene, respectively.

2.2 Age models

Due to the lack of absolute age markers, such as tephra layers, the LIG age models mostly rely on alignment strategies that tie each record to a well-dated reference record. The age model tie points used in this study are taken from the original age model publications. The reference records (LS16; Lisiecki and Stern, 2016) consist of eight regional stacks (one for the intermediate and one for the deep ocean each for the North Atlantic, South Atlantic, Pacific, and Indian oceans) of benthic δ¹⁸O that were dated through alignment with other climatic archives such as ice-rafted debris records, synthetic ice core records, and speleothems. The use of regional stacks, rather than a single global stack, improved stratigraphic alignment targets and provided more robust age models. The estimated age model uncertainty (2σ) for this group of cores is 2 kyr. Please refer to Lisiecki and Stern (2016) for further details. Oliver et al. (2010) defined their age tie points assuming that sea level minima and benthic δ¹⁸O maxima are synchronous. The benthic δ¹⁸O records were aligned with each other and then tied to the Dome Fuji chronology (based on ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$) (Kawamura et al., 2007). Please refer to Shackleton et al. (2000) and Oliver et al. (2010) for an extensive method description. The age model uncertainty on this group of cores is estimated to range from 1 to 2.5 kyr.

The published age models for the additional cores were determined using similar alignment techniques: SSTs were correlated to the NGRIP Greenland ice core for CH69-K09 and MD95-2042 (Govin et al., 2012). The age model for MD03-2664 was determined by correlating MD03-2664 δ¹⁸O with previously dated MD95-2042 δ¹⁸O (Galaasen et al., 2014b). ODP 1063 and U1304 δ¹⁸O were originally aligned to the LR04 stack (Lisiecki and Raymo, 2005). In order to align all of the records, adjustments to the age models of cores from Oliver et al. (2010) and the five additional cores (CH69-K09, MD95-2042, MD03-2664, ODP 1063, and U1304) were made by aligning the δ¹⁸O minima during the LIG to the corresponding δ¹⁸O minima of the nearest LS16 stack. The δ¹⁸O data before and after the alignment is given in Fig. S1 in the Supplement.

The Holocene age models are based on planktonic foraminifera radiocarbon dates (Stern and Lisiecki, 2014; Waelbroeck et al., 2001) that have been converted into calendar ages using IntCal13 and using reservoir ages based on modern observations (Key et al., 2004), which are assumed to have remained fairly stable across the Holocene. The age uncertainty associated with these Holocene radiocarbon-based age models is generally less than 0.5 kyr. However, it is important to note that Holocene age models from Oliver et al. (2010) were derived using the same method as their LIG age models, leading to larger age uncertainties of about 1–2.5 kyr for this set of Holocene records (four cores). The tie points were used to derive a full age–depth model assuming a constant sedimentation rate between tie points (i.e. linear interpolation).

2.3 Spatial coverage

The spatial distribution of the database for the Holocene and the LIG is shown in Fig. 2, and the depth distribution in each ocean basin is shown in Fig. 3. There are more data in the Atlantic Ocean (65 LIG, 118 Holocene) than in the Pacific (15 LIG, 19 Holocene) and Indian (3 LIG, 7 Holocene) oceans. We used this database to determine (1) if there is a significant difference in the average ocean δ¹³C signal at the LIG compared to the Holocene and (2) if ocean circulation patterns were comparable. Due to the sparsity of data in the Indian and Pacific oceans, our investigation is primarily focused on the Atlantic. Additionally, the temporal uncertainties (∼ 2 kyr) do not permit an investigation of centennial-scale events, and therefore we restrict our analysis to mean LIG and Holocene conditions.

Figure 2 Global distribution of benthic foraminifera δ¹³C covering the periods studied here: the Holocene (7–2 ka BP) (a) and LIG (125–120 ka BP) (b). Symbol size indicates the number of values per core, colour indicates average δ¹³C, and the triangle direction indicates the proxy depth (upward-pointing triangle: between 1000 and 2500 m depth; downward-pointing triangle: deeper than 2500 m). Four specific regions used in Sect. 3.1 are outlined: eastern equatorial Pacific (black, grey), equatorial Atlantic (yellow, green), southeastern Atlantic (cyan, blue), and northeastern Atlantic (magenta, red). Regional boundaries used to calculate the global volume-weighted mean δ¹³C (Sect. 3.2) are indicated by dotted black lines as defined in Peterson et al. (2014).

Figure 3 Zonal distribution of benthic foraminifera δ¹³C (‰) during the LIG (125–120 ka BP; a, c, e) and the Holocene (7–2 ka BP; b, d, f) in the Atlantic Ocean (a, b), Pacific Ocean (c, d), and Indian Ocean (e, f). Symbol size indicates the number of measurements per core and colour indicates average δ¹³C.

3 Results

The δ¹³C signal varies significantly regionally and with depth. The highest average δ¹³C values are associated with NADW and are generally found at depths between ∼ 1500 and 3000 m in the North Atlantic, with organic matter remineralisation and mixing with southern source waters leading to a δ¹³C decrease along the NADW path. The lowest δ¹³C values are in the deep South Atlantic (>4000 m) because the Antarctic Bottom Water (AABW) end member is much lower than its NADW counterpart. Since the Indian and Pacific oceans are mostly ventilated from southern-sourced water masses, δ¹³C generally decreases northward in these two basins.

Since the number of cores is not consistent across the two time periods and given the high regional variability observed in δ¹³C, it is not possible to simply average all available data to determine the global mean δ¹³C. Furthermore, the spatial heterogeneity of the data density adds to the complexity of the problem. To address these points, we first analyse differences between the LIG and Holocene records for pre-defined small regions with high data density. We then calculate regional volume-weighted δ¹³C means for larger regions from which we estimate the global LIG–Holocene anomaly.

3.1 Regional reconstruction of δ¹³C

We define regions with high densities of cores to reconstruct regional mean δ¹³C (Fig. 2). These regions need to be small enough to assume reasonably small spatial variability in the δ¹³C signal and yet still have enough data to establish a reliable statistical difference between the two time periods.

Table 3 Regional summary of δ¹³C below 2500 m depth for the LIG (125–120 ka BP) and Holocene (7–2 ka BP) using a single value per core for each time slice. Shown are the non-volume-weighted means (δ¹³C, ‰), standard deviations (σ, ‰), and counts (N) for both time periods, along with the time period regional anomalies (Δδ¹³C, ‰), propagated standard deviations for the anomaly (σ, ‰), and p values from two-sample t tests between the two time periods.

Based on these requirements, four regions are selected: the northeastern Atlantic, the equatorial Atlantic, a region off the Namibian Coast (southeastern Atlantic), and a region around the Galapagos Islands (eastern equatorial Pacific). The boundaries of each region are defined in Table 3.

We then define the time periods within the LIG and the Holocene to perform our analyses. For the Holocene, as most of the available data are dated prior to 2 ka BP, we define the end of our Holocene time period as 2 ka BP. To capture as much of the Holocene data as possible, we include data back to 7 ka BP, ensuring that we do not include instability associated with the 8.2 ka BP event (Alley and Ágústsdóttir, 2005; Thomas et al., 2007). This provides a time span of 5 kyr of data that we will consider for our analysis of the Holocene.

For the LIG, we seek to avoid data associated with the end of the penultimate deglaciation, which is characterised by a benthic δ¹³C increase in the Atlantic until ∼ 128 ka BP (Govin et al., 2015; Menviel et al., 2019; Oliver et al., 2010, Fig. 4). In addition, a millennial-scale event has been identified in the North Atlantic between ∼ 127 and 126 ka BP (Galaasen et al., 2014b; Tzedakis et al., 2018). Considering the typical dating uncertainties associated with the LIG data (2 kyr), we thus decide to start our LIG time period at 125 ka BP. To ensure that the two time periods are of the same length (5 kyr), we define the LIG period for our analysis to be 125–120 ka BP. We note that our definition should also avoid data associated with the glacial inception (Govin et al., 2015; Past Interglacials Working Group of PAGES, 2016). We verify that the LIG time period has sufficient data across the selected four regions, noting that the highest density of data falls within the 125–120 ka BP time period – particularly in the equatorial Atlantic and southeastern Atlantic (Fig. 4b, c).

Figure 4 Benthic foraminifera δ¹³C (left y axis, ‰) during the LIG (left) and Holocene (right) for four defined regions: northeastern Atlantic (a), equatorial Atlantic (b), southeastern Atlantic (c), and eastern equatorial Pacific (d). Data are presented in discrete time slices spanning 1 kyr. Only cores deeper than 2500 m are shown. Circular, coloured points connected by lines show each average δ¹³C value per core per time slice. Black symbols represent δ¹³C averages per slice. Each slice has a corresponding averaged depth (right y axis, m), with 1 standard deviation on either side shown on the bars. Slices with an average depth within ±300 m of the mean core depth of all slices are represented with a square point. Slices with an average depth 300 m shallower than the mean are shown with an upward triangle, and slices that are than 300 m deeper than the mean are shown with a downward triangle. Shading shows 1 standard deviation on either side of the mean for slices where more than one point exists.

We round the data to the nearest 1 kyr, find an average per 1 kyr, and refer to this as a time slice. We consider the LIG–Holocene anomaly across these two time periods for the four regions selected, and qualitatively consider the influence of changes in the average depth in which the proxies were recorded, as indicated by the direction of the black triangles in Fig. 4.

The average δ¹³C anomaly between the LIG and Holocene periods for cores deeper than 2500 m is consistent across the different regions despite their geographic separation, suggesting a significantly lower δ¹³C during the LIG than the Holocene, with differences ranging from −0.13 ‰ in the northeastern Atlantic to −0.20 ‰ in the equatorial Pacific (Table 3). The statistical significance between the two time periods is established using a two-tailed t test on data of one mean value per core and spans all time slices (125–120 and 7–2 ka BP). The t test shows that there is a statistically significant difference everywhere except in the equatorial Atlantic, with confidence intervals varying from 0.13 in the equatorial Atlantic to 0.04 in the northeastern Atlantic. When using a single tail t test instead, the difference becomes significant in the equatorial Atlantic, giving a new p value of 0.02. Figure 4 suggests that the difficulty in determining significance in this region for cores deeper than 2500 m might be due to a singular outlier measurement in the equatorial Atlantic: a value of −0.23 ‰ from GeoB1118 at ∼ 3.5 ka BP. If this value is excluded, then an anomaly of −0.22 with a p value less than 0.005 is observed.

We also compare the distribution of δ¹³C for cores deeper than 2500 m for three overlapping periods within the LIG (early LIG: 128–123 ka BP; LIG: 125–120 ka BP; late LIG: 123–118 ka BP). The results for the four regions are shown in Fig. 5. The statistical characteristics do not show much variation between the LIG and late LIG δ¹³C distributions. In the equatorial Pacific, the difference between the early LIG and the Holocene is smaller than between the LIG and Holocene, but this is countered with a larger difference in the equatorial Atlantic between the early LIG and Holocene. The spread in the data is generally larger during the Holocene than during the other time periods, which might be due to the greater number of measurements during the Holocene. The spread of data during the early LIG is slightly larger than during the LIG and late LIG in the equatorial and southeastern Atlantic. The equatorial Atlantic is the only region that displays significantly more points with lower δ¹³C during the early LIG. Overall, these distributions do not suggest that the LIG–Holocene anomalies that we have determined would be significantly impacted by slight variations in the selected time window. We perform an analysis of variance (ANOVA) on each region and post hoc tests on the data. We find that the Holocene data are significantly different from the three LIG periods in the northeastern Atlantic, the southeastern Atlantic, and the equatorial Pacific, while the three periods within the LIG are not significantly different from each other for any of the regions.

Figure 5 Distributions of δ¹³C for all core measurements deeper than 2500 m during the Holocene (7–2 ka BP, red), the early LIG (128–123 ka BP), the LIG (125–120 ka BP), and the late LIG (123–118 ka BP) across four regions (equatorial Pacific, equatorial Atlantic, southeastern Atlantic, northeastern Atlantic). The lower end of the box indicates quartile 1 (Q1), and the upper end indicates quartile 3 (Q3). Orange vertical lines show the median, and dotted vertical lines show the mean. The whiskers indicate the lower and upper fences of the data calculated as Q1 − 1.5 × (Q3 − Q1) and Q3 + 1.5 × (Q3 − Q1), respectively, and the clear circles are outliers.

3.2 Volume-weighted regional δ¹³C

The second approach we use to further constrain the LIG–Holocene δ¹³C anomaly is to estimate the volume-weighted regional δ¹³C. We define our regional boundaries based on the regions described in Peterson et al. (2014); however, we only include the regions where there are enough data to justify an analysis. For all the data in each of these regions, we calculate a mean value by taking the direct averages of all data. We divide the ocean basins into eight regions (Table 4, shown in Fig. 2) and calculate the volume-weighted averages δ¹³C for each of these regions. Since the Atlantic and Pacific oceans have more data than the Indian Ocean, there is greater confidence in the δ¹³C estimates for these regions. These regional averages are then used to calculate a global volume-weighted δ¹³C.

Table 4 Regional breakdown of δ¹³C data for all depths during the Holocene (7–2 ka BP) and LIG (125–120 ka BP) averaged across the 1 kyr time slices. For each region: the average number of data points (labelled as “Points”) and cores per time slice (labelled as “Cores”), the average standard deviation of δ¹³C per time slice (‰), the mean depth (m) across time slices, and the standard deviation of depth (m) between time slices (σ_depth). Regions are abbreviated as follows: NEA: northeastern Atlantic; NWA: northwestern Atlantic; SA: southern Atlantic; SEA: southeastern Atlantic; SWA: southwestern Atlantic; I: Indian, NP: northern Pacific; SP: southern Pacific.

Results for the Atlantic and Pacific oceans are given in Fig. 6 and show a mean LIG–Holocene anomaly of −0.21 ‰ and −0.27 ‰, respectively, slightly higher than the range of estimates for the four regions selected in Sect. 3.1. There is a higher offset estimated in this definition of the southwestern Atlantic (−0.45 ‰) than in Sect. 3.1; however, there are only four cores available in this region during the Holocene and two during the LIG.

Figure 6 Comparison of volume-weighted δ¹³C for the Atlantic (red) and Pacific (blue) for the LIG and Holocene, calculated using the regions from Peterson et al. (2014) from data covering all depths. Solid coloured lines indicate the mean volume-weighted δ¹³C, and the shading indicates the volume-weighted sum of square deviations from the mean. The horizontal bars indicate the mean of the stable period determined from the regional analysis as defined in Sect. 3.1 (LIG: 125–120 ka BP; Holocene: 7–2 ka BP), with the Δδ¹³C indicating the mean anomaly between these two averages and the standard deviation (σ(δ¹³C), ‰).

The estimated LIG–Holocene anomaly in the southern Pacific is relatively high at −0.39 ‰, giving a relatively large Pacific anomaly estimate of −0.27 ‰. This could be due in part to the deeper location of the LIG cores compared to the mean of the Holocene cores (439 m difference, Table 4). There is less confidence in the estimate of the Pacific volume-weighted mean since the proxy data are sparse, and the majority of cores are from the eastern equatorial Pacific as shown in Fig. 2. We also note that the average depths of cores from the Pacific Ocean (LIG: 2711 m; Holocene: 2131 m) and Indian Ocean (LIG: 2383 m; Holocene: 2303 m) are shallower than that of the Atlantic Ocean (LIG: 3531 m; Holocene: 3157 m; Fig. 3). However, as the vertical gradient below 2000 m depth in the Pacific Ocean is small (e.g. Eide et al., 2017), this might not significantly impact our results.

There is a small positive trend in the average Atlantic δ¹³C from 125 ka BP, reaching a maximum value at 118 ka BP (Fig. 6). The average core depth over the 125–120 ka BP time period does not suggest that a change in the mean depth could explain this variation. Fitting a linear regression over this period indicates an increase in δ¹³C of 0.03 ‰ kyr⁻¹ in the Atlantic, with a p value of 0.01 and an R² of 0.85 (Fig. 4a). For the Pacific, there is a ∼ 0.13 ‰ increase in δ¹³C between 7 and 5 ka BP, which could be associated with the early Holocene terrestrial regrowth (Menviel and Joos, 2012).

For the Indian Ocean, we only include four cores, as these are the only ones spanning both the LIG and Holocene. An LIG anomaly of −0.13 ‰ in the Indian Ocean compared to the Holocene is therefore associated with higher uncertainties. The whole ocean mean LIG δ¹³C anomaly is −0.25 ‰, but it is associated with higher uncertainties than each region anomaly.

Both the regional analysis of our new database and our volume-weighted estimate indicate that the global mean δ¹³C was about 0.2 ‰ lower during the LIG than during the Holocene. We further test the robustness of this result in the next section.

3.3 Reconstruction of the LIG Atlantic Meridional Overturning Circulation

In this section, we analyse the spatial δ¹³C distribution in the Atlantic Ocean to assess potential changes in the penetration depth and southward expansion of NADW during the LIG, defined here as 125–120 ka BP, with respect to the Holocene. A change in NADW might influence our estimate of the mean δ¹³C, given that most of the available data are localised in the Atlantic Ocean.

We use simple statistical regression models to reconstruct NADW and AABW separately with a quadratic-with-depth and linear-with-latitude equation following the method of Bengtson et al. (2019). For consistency, the regression algorithm only includes records from cores that span both the LIG and Holocene and uses a weighted least-squares approach, where the weighting equals the number of samples per core. The modelled region is defined between 40^∘ S and 60^∘ N as this is the region where we can expect to find both the NADW and AABW δ¹³C signals.

Figure 7 Reconstructed Atlantic δ¹³C (‰) meridional section during the LIG (125–120 ka BP) and Holocene (7–2 ka BP). The circular points represent the proxy data, showing the average δ¹³C with colour and the number of points per core with size. The stars represent the proxy data that make up the end members. Background shading shows the reconstructed δ¹³C using a quadratic statistical regression of the proxy data following the method described in Bengtson et al. (2019).

The results are shown in Fig. 7. We test the robustness of our statistical model using the jackknifing technique. We systematically exclude each individual core from the database one at a time, fit the parameters using this modified database, and compare the model prediction against the core, which was excluded. This produces small variations in the average mean response of the statistical models (the standard deviations were 0.01 ‰ for both the LIG and Holocene, respectively).

We calculate end member values based on proxies located near the water mass sources. These are taken as 0.79 ‰ and 1.02 ‰ for NADW for the LIG and Holocene, respectively, and −0.09 ‰ and 0.23 ‰ for AABW for the LIG and Holocene, respectively. The end member values are calculated as the average of cores shallower than 3000 m but deeper than 1000 m and located between 50 and 70^∘ N for NADW. The NADW end member cores have an average depth of 2043 m and a standard deviation of 478 m during the LIG. For the AABW end member, the only eligible core is ODP1089, which is at ∼ 41^∘ S and 4621 m.

The mean volume-weighted δ¹³C for the Atlantic Ocean between 40^∘ S and 60^∘ N based on this interpolation is 0.53 ‰ for the LIG and 0.70 ‰ for the Holocene (Fig. 7). This suggests a 0.17 ‰ lower Atlantic δ¹³C at the LIG than the Holocene. Our statistical reconstruction points to a very similar NADW depth (∼ 2600 m) for both time periods (Fig. 7). The NADW depth is defined here as the depth of maximum δ¹³C in the North Atlantic.

Figure 8 The meridional gradient of the Atlantic Ocean benthic δ¹³C (‰). (a) Holocene (red) and LIG (blue) δ¹³C for each core (points) between 1000 and 3000 m. Dotted lines are the moving averages of the cores. Solid lines indicate the results of our statistical model at 2000 m. (b) Average δ¹³C for each record deeper than 1000 m as a proportion of the end members. A value of one indicates the NADW end member, and a value of zero indicates the AABW end member. Solid lines show the linear regressions of the records.

We also investigate the meridional gradient in δ¹³C in the Atlantic Ocean to determine whether the NADW southward penetration, transport, and remineralisation rates were significantly different during the LIG compared to the Holocene. We only consider cores that are located between depths of 1000 and 3000 m in order to stay within the main pathway of NADW (Fig. 8a). Though there is significant scatter, in accordance with our previous findings, a moving average through the Holocene and the LIG data shows that LIG δ¹³C is typically lower than the Holocene counterparts. However, the meridional δ¹³C statistical model gradients are not very different for the LIG (0.0036 ‰ per degree latitude) and the Holocene (0.0030 ‰ per degree latitude) (Fig. 8a), suggesting a similar southward penetration of NADW.

Using δ¹³C of the end members for NADW and AABW, we use a simple binary mixing model for all cores deeper than 1000 m to estimate changes in NADW penetration (Fig. 8b). The LIG and Holocene δ¹³C slopes in the Atlantic are similar, indicating similar southward penetration of NADW during both time periods. This suggests that the differences in δ¹³C between the two time periods is most likely due to change in end member values, while the mean Atlantic oceanic circulation was likely similar.

Based on our analysis, there appears to be no significant difference in the mean time-averaged AMOC between the LIG and the Holocene. Negative LIG–Holocene anomalies are found for each of the smaller regions selected (northeastern Atlantic, equatorial Atlantic, southeastern Atlantic, and eastern equatorial Pacific), with statistical significance seen in all regions except the equatorial Atlantic, where an unusual low δ¹³C value in one core is responsible for narrowing the difference between the two period means. Additionally, our volume-weighted mean δ¹³C estimates have similar anomalies in the Atlantic and Pacific oceans (−0.21 ‰ and −0.27 ‰, respectively, Fig. 6).

4 Discussion

One of the goals of our study is to assess the mean change in oceanic δ¹³C between the LIG and the Holocene. Given the uncertainties in the chronologies, avoiding data that pertains to deglaciation, and capturing the same length of time during the LIG and the Holocene, we chose the periods 125 to 120 ka BP for the LIG and 7 to 2 ka BP for the Holocene. Using a similar geographical distribution of data points for both periods, we find that the oceanic δ¹³C was ∼ 0.2 ‰ lower during the LIG than the Holocene.

Our analysis of the δ¹³C signal suggests consistent LIG–Holocene δ¹³C anomalies in different regions of the Atlantic basins, as well as in the Pacific and Indian oceans, even if there are significant uncertainties with the later due to fewer available records. The δ¹³C distribution in the Atlantic Ocean suggests that there was no significant mean change in the southward penetration or depth of NADW during the LIG (125–120 ka BP) compared to the Holocene (7–2 ka BP). A statistical reconstruction of the early LIG (128–123 ka BP) δ¹³C compared to our 125–120 ka BP reconstruction does not reveal a significant difference in either the NADW core depth or NADW extent, as indicated by the meridional δ¹³C gradients (Fig. S2). The volume-weighted average δ¹³C during the early LIG is 0.06 ‰ lighter than during the LIG period considered here (125–120 ka BP). Since both time slices (128–123 and 125–120 ka BP) are 5 kyr averages and include dating uncertainties of ∼ 2 kyr, it is not possible to resolve potential centennial-scale oceanic circulation changes (e.g. Galaasen et al., 2014b; Tzedakis et al., 2018).

Explanations for the 0.2 ‰ lower δ¹³C anomaly in the ocean may include a redistribution between the ocean–atmosphere system. Such a redistribution can result from a change in end member values (Fig. 8). As fractionation during air–sea gas exchange is temperature dependent, globally higher SSTs at the LIG could lead to a lower oceanic δ¹³C. However, the effect of this is likely small (Brovkin et al., 2002) and would also lead to a higher atmospheric δ¹³CO₂ at the LIG, which is inconsistent with Antarctic ice core measurements that suggest an anomaly of −0.3 ‰ (Schneider et al., 2013). Lower nutrient utilisation in the North Atlantic would decrease surface ocean δ¹³C and thus the δ¹³C end members. However, this would also imply that less organic carbon would be remineralised at depth. Therefore, it is unlikely that the lower average oceanic mean δ¹³C results from a change in end members through lower surface ocean nutrient utilisation. Currently, there is still a lack of constraints on nutrient utilisation in these end member regions during the LIG compared to the Holocene. Therefore, the lower δ¹³C in the ocean–atmosphere system cannot be explained by a simple redistribution of δ¹³C between the atmosphere and the ocean.

An alternative explanation for the anomaly is a change in the terrestrial carbon storage, which has a typical signature of approximately −37 ‰ to −20 ‰ for C₃-derived plant material (Kohn, 2010) and −13 ‰ for C₄-derived plant material (Basu et al., 2015). The total land carbon content at the LIG is poorly constrained. Proxies generally suggest extensive vegetation during the LIG compared to the Holocene (CAPE, 2006; Govin et al., 2015; Larrasoaña et al., 2013; Muhs et al., 2001; Tarasov et al., 2005; de Vernal and Hillaire-Marcel, 2008), which would imply a greater land carbon store. However, other terrestrial carbon stores including peatlands and permafrost may also have differed during the LIG compared to the Holocene. With an estimated ∼ 550 Gt C stored in peats today (mean δ¹³C ∼ −28 ‰; Dioumaeva et al., 2002; Novák et al., 1999) and ∼ 1000 Gt C in the active layer in permafrost, which may have been partially thawed during the LIG (Reyes et al., 2010; Schuur et al., 2015; Stapel et al., 2018), less carbon stored in peat and permafrost at the LIG could have led to a lower total land carbon store compared to the Holocene. However, it is not possible to infer this total land carbon change from the oceanic and atmospheric δ¹³C anomalies because it cannot be assumed that the mass of carbon and ¹³C is preserved within the ocean–atmosphere–land biosphere system on glacial–interglacial timescales. There is indeed continuous exchange of carbon and ¹³C between the lithosphere and the coupled ocean, atmosphere, and land biosphere carbon reservoirs. Isotopic perturbations associated with changes in the terrestrial biosphere are communicated to the burial fluxes of organic carbon and CaCO₃ and are therefore attenuated on multi-millennial time scales (Jeltsch-Thömmes et al., 2019; Jeltsch-Thömmes and Joos, 2020). Nevertheless, when hypothetically neglecting any exchange with the lithosphere, we find that the change in terrestrial carbon needed to explain the difference in δ¹³C would be on the order of 310 ± 44 Gt C less during the LIG than the Holocene (Text S1).

In addition, due to the warmer conditions at the LIG than during the Holocene, there could have been a release of methane clathrates which would have added isotopically light carbon (δ¹³C: ∼ −47 ‰) to the ocean-atmosphere system. However, available evidence suggests that geological CH₄ sources are rather small (Bock et al., 2017; Dyonisius et al., 2020; Hmiel et al., 2020; Saunois et al., 2020) making this explanation unlikely, although we cannot completely exclude the possibility that the geological CH₄ source was larger at the LIG than the Holocene. Similarly, since δ¹³CO₂ from volcanic outgassing has a similar value to atmospheric δ¹³CO₂ (Brovkin et al., 2016) and modelling suggests volcanic outgassing likely only had a minor impact on δ¹³CO₂ (Roth and Joos, 2012), it is unlikely that volcanic outgassing of CO₂ played a significant role in influencing the mean oceanic δ¹³C.

While changes in terrestrial carbon could have impacted the oceanic δ¹³C at the LIG, the LIG–Holocene differences in the isotopic signal of both the atmosphere and ocean were most likely due to a long-term imbalance between the isotopic fluxes to and from the lithosphere, including the net burial (or redissolution) of organic carbon and CaCO₃ in deep-sea sediments, and changes in shallow water sedimentation and coral reef formation (Jeltsch-Thömmes and Joos, 2020).

5 Conclusions

We present a new compilation of benthic δ¹³C from 130 to 115 ka BP covering the LIG. Over this time period, benthic δ¹³C generally display a maximum value at ∼ 121 ka BP (±2 kyr), in phase with the maximum atmospheric δ¹³CO₂ (LIG value of −6.5 ‰ at ∼ 120 ka BP). As there are significant chronological uncertainties associated with LIG records, we analyse data between 125 and 120 ka BP to avoid data associated with millennial-scale events and the deglaciation. We compare this LIG benthic δ¹³C data to a similar database covering the Holocene (7–2 ka BP). We find that during these specific time periods, LIG oceanic δ¹³C was about 0.2 ‰ lower than during the Holocene. This anomaly is consistent across different regions in the Atlantic Ocean. Even though there are fewer records available, benthic δ¹³C data from the Pacific Ocean also support an anomaly of about 0.2 ‰.

An analysis of δ¹³C gradients across the Atlantic Ocean suggests that there were no significant changes in mean long-term ocean circulation between the two intervals. While reduced high northern latitude peat and permafrost caused by higher temperatures at the LIG than during the Holocene (Otto-Bliesner et al., 2021) could have lead to a lower atmospheric and oceanic δ¹³C, the most likely explanation for the lower LIG oceanic δ¹³C is a long-term imbalance in the weathering and burial of carbon. Additional studies are required to further constrain the LIG carbon balance.