Atlantic circulation changes across a stadial–interstadial transition

We combine consistently dated benthic carbon isotopic records distributed over the entire Atlantic Ocean with numerical simulations performed by a glacial configuration of the Norwegian Earth System Model with active ocean biogeochemistry in order to interpret the observed Cibicides δ¹³C changes at the stadial–interstadial transition corresponding to the end of Heinrich Stadial 4 (HS4) in terms of ocean circulation and remineralization changes. We show that the marked increase in Cibicides δ¹³C observed at the end of HS4 between ∼2000 and 4200 m in the Atlantic can be explained by changes in nutrient concentrations as simulated by the model in response to the halting of freshwater input in the high-latitude glacial North Atlantic. Our model results show that this Cibicides δ¹³C signal is associated with changes in the ratio of southern-sourced (SSW) versus northern-sourced (NSW) water masses at the core sites, whereby SSW is replaced by NSW as a consequence of the resumption of deep-water formation in the northern North Atlantic and Nordic Seas after the freshwater input is halted. Our results further suggest that the contribution of ocean circulation changes to this signal increases from ∼40 % at 2000 m to ∼80 % at 4000 m. Below ∼4200 m, the model shows little ocean circulation change but an increase in remineralization across the transition marking the end of HS4. The simulated lower remineralization during stadials compared to during interstadials is particularly pronounced in deep subantarctic sites, in agreement with the decrease in the export production of carbon to the deep Southern Ocean during stadials found in previous studies.

1 Introduction

During the last glacial period, surface temperatures above Greenland and in the North Atlantic region shifted between cold (stadial) and warm (interstadial) phases. Previous studies have shown that these rapid changes in surface climatic conditions were accompanied by rapid changes in ocean circulation characterized by reduced Atlantic meridional overturning circulation (AMOC) during stadials (Vidal et al., 1997; Gottschalk et al., 2015; Waelbroeck et al., 2018; Toucanne et al., 2021). However, the precise geometry and extent of Atlantic circulation changes are still under debate due to the scarcity of available observations and the large dating uncertainties of marine sediment records.

Among the proxy data that can be used to reconstruct past changes in ocean circulation, benthic foraminiferal oxygen (δ¹⁸O, expressed as an anomaly relative to a standard in ‰ versus Vienna Pee Dee Belemnite, VPDB) and carbon (δ¹³C, expressed in ‰ versus VPDB) isotopes offer the best spatial and temporal coverage. However, these proxies do not solely record ocean circulation changes and depend on several factors that complicate their interpretation. At present, the extraction of the ocean circulation signals out of benthic δ¹⁸O and δ¹³C records remains a relevant objective which has the potential to significantly improve our understanding of the ocean processes at play during the rapid climate changes of the last glacial and deglaciation periods.

Over the past 25 years, it has been shown that a decrease in benthic δ¹³C took place in the Atlantic Ocean during Heinrich stadials, whose amplitude depends on the location and water depth of the cores. This decrease in benthic δ¹³C has been classically interpreted as a reflection of a decrease in northern-sourced water (NSW) formation and a concomitant increase in the southern-sourced water (SSW) fraction (Willamowski and Zahn, 2000; Elliot et al., 2002; Skinner and Shackleton, 2006; Peck et al., 2006, 2007). However, more recent studies also suggested that this decrease in benthic δ¹³C could be due to a decrease in NSW δ¹³C preformed value (Waelbroeck et al., 2011; Crocker et al., 2016; Lund et al., 2015) or an increase in organic matter remineralization (Hoogakker et al., 2007; Lacerra et al., 2017; Voigt et al., 2017) or a combination of these factors and of an increase in the SSW fraction (Oppo et al., 2015).

Other ocean circulation proxies include sedimentary Pa/Th, sortable silt, neodymium isotopic composition, benthic foraminiferal Ba/Ca, and Cd/Ca (Lynch-Stieglitz and Marchitto, 2014). However, caveats about each of these proxies, combined with the scarcity of well-resolved records of the last glacial millennial changes, have not permitted the reconstruction of 3-D changes in ocean water masses distribution associated with the rapid climate changes of the last glacial.

In such a context, isotope-enabled models appear to be extremely useful tools for disentangling the different factors that influence the benthic δ¹³C signal and for producing 3-D pictures of the ocean circulation changes. Such models have been mainly used to simulate pre-industrial (PI) (Schmittner et al., 2013; Menviel et al., 2015) or Last Glacial Maximum (LGM; Menviel et al., 2017, 2020; Muglia et al., 2018; Muglia and Schmittner, 2021; Morée et al., 2021) conditions. A few studies with isotope-enabled models have addressed ocean circulation changes induced by hosing experiments mimicking stadial–interstadial transitions under PI boundary conditions (Schmittner and Lund, 2015; Missiaen et al., 2020). To our knowledge, there has been only one study using a transient simulation with carbon isotopes initialized under glacial conditions (Gu et al., 2021). This study focused on the mid-depth (1500–2500 m) Atlantic δ¹³C decrease observed across the transition from the LGM to Heinrich Stadial 1 and concluded that this decrease is mainly explained by increased remineralization due to AMOC slowdown, while the water mass mixture change only plays a minor role (Gu et al., 2021).

Here, we combine consistently dated benthic carbon isotopic records distributed over the Atlantic Ocean with numerical simulations performed by a glacial configuration of the Norwegian Earth System Model (NorESM) with active ocean biogeochemistry, and we investigate the changes observed in the Cibicides δ¹³C records across the stadial–interstadial transition corresponding to the end of Heinrich Stadial 4 (HS4). This climate transition, dated at 38.17±0.73 ka, is characterized by a rapid increase in Greenland and North Atlantic surface temperatures, leading from HS4 to Greenland Interstadial 8 (GI8) conditions (Rasmussen et al., 2014; Waelbroeck et al., 2019). This transition is particularly interesting because it is the largest and best-expressed transition in the Cibicides δ¹³C records prior to the last deglaciation. It thus offers a case study of a rapid and large climatic transition away from large changes in insolation and greenhouse gases. Therefore, one can assume that the recorded climate and ocean circulation changes across the HS4-to-GI8 transition are not driven by changes in the radiative forcing. This reduces the dynamical complexity and makes the use of a hosing experiment under constant radiative forcing adequate to interpret the observed changes in the proxy records.

Our study provides a 3-D picture of the ocean circulation changes across this stadial–interstadial transition anchored in Cibicides δ¹³C observations distributed between 2000 and 5000 m depth. We decompose the non-sea-air component of the δ¹³C change across the HS4-to-GI8 transition into a fraction due to remineralization changes and the remaining fraction due to changes in water mass origin and geometry. Our results show that the non-sea-air δ¹³C change across the HS4-to-GI8 transition due to remineralization changes decreases with increasing water depth and that, below ∼3000 m, the replacement of SSW by NSW plays a dominant role.

2 Material and methods

2.1 Observations

2.1.1 Benthic carbon isotopes

The δ¹³C isotopic ratio of the epifaunal benthic foraminifer genus Cibicides (noted as Cib. hereafter) has been shown to record the δ¹³C of bottom-water dissolved inorganic carbon (DIC), δ¹³C-DIC, with minor isotopic fractionation (Duplessy et al., 1984; Zahn et al., 1986; Schmittner et al., 2017). The initial DIC isotopic composition of a given water mass is governed by surface productivity in its formation region (i.e., the preferential consumption of ¹²C by primary productivity, thereby increasing dissolved δ¹³C), as well as by temperature-dependent air–sea exchanges (Lynch-Stieglitz et al., 1995). The δ¹³C-DIC subsequently decreases as deep water ages due to the progressive remineralization at depth of relatively ¹³C-depleted biogenic material. As a result, δ¹³C-DIC largely follows water mass structure and circulation in the modern ocean, and Cib. δ¹³C has been used to trace water masses and as a proxy of bottom-water ventilation (Duplessy et al., 1988, and numerous subsequent studies).

In short, δ¹³C-DIC is driven by the biological activity (i.e., photosynthesis and remineralization) on the one hand and by thermodynamic fractionation during air–sea exchanges on the other hand.

A recent reconstruction of the pre-industrial δ¹³C-DIC distribution and of its partitioning into biological and thermodynamic components shows that the former term prevails below 2000 m, with the air–sea δ¹³C component ranging from −0.4 ‰ VPDB in the North Atlantic to +0.2 ‰ in the Southern Ocean, which is in contrast to the biological δ¹³C component that ranges from +1.4 ‰ PDB in the North Atlantic to +0.4 ‰ in the Southern Ocean (Eide et al., 2017).

2.1.2 Cibicides δ¹³C time series

A set of 110 consistently dated foraminifer isotopic records was built during the ACCLIMATE project making use of the PARIS database structure (Lougheed et al., 2022). Among these 110 cores, 92 could be directly dated according to the approach described in Waelbroeck et al. (2019), but the age models of the remaining cores had to be established by alignment to some of the directly dated cores, thereby introducing additional dating uncertainties.

For the purposes of the present study, the dated Cib. δ¹³C records were screened in order to retain only sufficiently well-resolved (average time step of 100 to 500 years in most cases and smaller than 750 years in two cases) and well-dated Cib. δ¹³C records in the vicinity of the HS4-to-GI8 stadial–interstadial transition. Therefore, only cores with relatively precise chronological information at the end of HS4 were retained – that is, cores that presented a radiocarbon date or an alignment tie point at the end of HS4, whose calendar age has a 1σ dating uncertainty lower than 1000 years.

Applying this screening criteria to the 110 ACCLIMATE cores, we are left with only 18 Cib. δ¹³C records (Table S1 in the Supplement). Among the selected 18 Cib. δ¹³C records, 6 are new records which were generated during the ACCLIMATE project (Sect. S1 in the Supplement). In terms of dating, six cores could not be directly dated but were aligned to one of the previously published cores (Waelbroeck et al., 2019) taken as a reference. Also, the age models of 4 of the remaining 12 cores have been improved with respect to the age model published in 2019. The reader is referred to Sect. S2 for detailed information on the dating.

2.2 Numerical experiment

In this study, we took advantage of a fast version of NorESM, NorESM1-F (Guo et al., 2019a), which features a horizontal resolution of about 2^∘ in the atmosphere and nominal 1^∘ in the ocean and sea ice. The computational efficiency of the model allows the experiments to be integrated for multiple millennia within a reasonable wall clock time. The NorESM1-F model with active ocean biogeochemistry has been used in several studies to investigate changes in the ocean circulation and carbon cycle (Kessler et al., 2018; Galaasen et al., 2020).

The NorESM glacial simulations reported in this study are based on previous work by Guo et al. (2019b) and Jansen et al. (2020). To briefly summarize, we performed an equilibrium simulation of the last glacial period forced by the boundary conditions of 38 ka, which corresponds to the beginning of GI8. The boundary conditions are different from the present day in terms of orbital forcing, greenhouse gases (CO₂, CH₄, and N₂O), land–sea mask, and the height and extent of the global ice sheets. Details of the configuration and simulation results of the physical climate are documented in Guo et al. (2019b).

The physical climate reached a stable interstadial state after 2500 model years, as shown by Guo et al. (2019b). From this point, we activated the full coupling of the physical components of NorESM1-F with its ocean biogeochemical (BGC) component (HAMOCC; Tjiputra et al., 2020) that is initialized from present-day climatology (including PO₄, NO₃, oxygen, dissolved inorganic carbon, and alkalinity; Key et al., 2004; Garcia et al., 2010a, b). Due to the long equilibration timescales of BGC tracers and the difference between glacial and present-day climate states, we continued the integration for another 2510 years (e.g., until model year 5010) to ensure that the BGC tracers are in a satisfactory quasi-equilibrium state. This ∼2500-year length of the simulation with the BGC component activated is a compromise between demanding computing time and the relatively long time period (e.g., multiple millennia) required for ocean physics and BGC fields to reach a quasi-equilibrium state. More specifically, we consider the quasi-equilibrium state of the BGC fields to be satisfactory when the residual BGC drift is substantially smaller than the signal induced by freshwater forcing (on the order of a few percent).

Next, we branched off from the simulation above and applied a freshwater-hosing method to force the interstadial climate state into a cold stadial state. The freshwater is evenly distributed between 50 and 70^∘ N, with an injection of 0.33 Sv (1 Sv = 10⁶ m³ s⁻¹) for 800 years (i.e., until model year 5810), during which time a full stadial-like state (with a weak and shallow AMOC and an extensive sea ice cover in the North Atlantic) is achieved in the early phase of the simulation with freshwater forcing. We subsequently stopped the freshwater input to study the cold-to-warm transition process of the climate system. The model simulation of recovery (i.e., without freshwater forcing) lasted for another 400 years (i.e., until model year 6210).

The freshwater simulation and transition process of the physical climate system is partly documented by Jansen et al. (2020). Note that the freshwater input lasted for only 500 years in that study; however, the dynamics of the physical climate processes during the transition are quite similar with the two forcing time lengths. In this study, we extend the forcing length to 800 years in order to (1) better match the observed length of HS4 (Wolff et al., 2010; Waelbroeck et al., 2019) and (2) achieve a better-equilibrated deep Atlantic Ocean during the stadial state.

As shown by Jansen et al. (2020), NorESM1-F provides a robust representation of the HS4-to-GI8 transition process, both in terms of the magnitude (∼10 ^∘C) and rate (nearly 1 ^∘C per decade) of temperature change when comparing to Greenland ice core records (see Fig. 3 of Jansen et al., 2020). The transition process takes about 100 years to complete, after which the climate returns to an interstadial state close to the pre-hosing period. Further dynamical processes and sequences of events in the climate system that occurred during the transition, including favorable agreement with sea surface and subsurface temperature and sea ice records in the Nordic Seas and North Atlantic, will be presented in a more detailed study (Guo, 2023). The focus of the present study is on the change in ocean circulation as revealed by the observed and simulated distribution of nutrient and BGC tracers.

2.3 Simulated tracers

As explained above, δ¹³C-DIC can be decomposed into biological and air–sea exchange components. Using the approach of Broecker and Maier-Reimer (1992), we can compute the biological component (termed δ¹³C-BIO hereafter) of the δ¹³C-DIC as a linear function of the simulated phosphate concentration, PO₄:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}\text{-}\mathrm{BIO}& -{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{mean}\mathrm{o}}={\displaystyle \frac{\mathit{\alpha }\mathit{\_}\mathrm{photo}}{{\mathrm{DIC}}_{\mathrm{mean}\mathrm{o}}}}\cdot {\displaystyle \frac{\mathrm{C}}{\mathrm{P}}}\\ & \cdot \left({\mathrm{PO}}_{\mathrm{4}}-{\mathrm{PO}}_{\mathrm{4}\mathrm{mean}\mathrm{o}}\right),\end{array}\end{array}$$

where α_photo is the carbon isotopic fractionation during photosynthesis; $\mathrm{C}/\mathrm{P}$ is the ocean mean carbon-to-phosphorus ratio; and δ¹³C_mean o, DIC_mean o, and PO_4 mean o are the mean ocean δ¹³C, DIC, and PO₄, respectively.

We take α_photo = −19 ‰ (Broecker and Maier-Reimer, 1992) and $\mathrm{C}/\mathrm{P}=\mathrm{122}$ (Takahashi et al., 1985), both of which can be assumed to be constant in time. We take δ¹³C_mean o=0.14 ‰ VPDB, as estimated for the glacial ocean (Gebbie et al., 2015). We apply DIC_mean o=2188.4 µmol kg⁻¹ and PO_4 mean o=1.7483 µmol kg⁻¹, taken as the long-term mean values of the MIS3 model simulation.

The NorESM1-F model is not isotope-enabled, so in this study, we focus on water depths below 2000 m and use δ¹³C-BIO as an approximation of δ¹³C-DIC. At the pre-industrial, this approximation yields slightly overestimated δ¹³C-DIC values in the North Atlantic (by 0.0 ‰–0.4 ‰) and slightly underestimated δ¹³C-DIC values in the Southern Ocean (by 0.0 ‰–0.2 ‰; Eide et al., 2017). So, not surprisingly, we find offsets between the simulated δ¹³C-BIO and observed Cib. δ¹³C values (Fig. S1 in the Supplement). Were the model perfect, these offsets would correspond to the air–sea component of δ¹³C-DIC during the simulated interstadial and subsequent hosing experiment.

In this study, we circumvent the existence of these offsets by comparing the simulated change in δ¹³C-BIO to the observed change in Cib. δ¹³C across the HS4-to-GI8 transition rather than the computed δ¹³C-BIO and the observed Cib. δ¹³C values during the HS4 and GI8 time intervals, respectively. Furthermore, we minimize the errors resulting from the non-inclusion of the air–sea component of δ¹³C-DIC by retaining only core sites located below 2000 m.

In addition to nutrients and carbon tracers, the NorESM1-F model provides preformed tracers, such as the preformed phosphate (PO${}_{\mathrm{4}}^{\mathrm{pre}}$), DIC, or dissolved oxygen (Guo et al., 2019a; Tjiputra et al., 2020). The preformed concentration of a tracer is defined as the concentration of this tracer in a parcel of water the last time this parcel was at the ocean surface. A tracer can thus be separated into its preformed concentration and its regenerated concentration, i.e., the concentration resulting from biological remineralization in the interior of the ocean (Duteil et al., 2012; Bernardello et al., 2014). For instance, preformed phosphate reflects the amount of phosphate returning to the ocean interior by physical processes (Duteil et al., 2012).

At any location below the ocean mixed layer, the total PO₄ concentration can thus be expressed as the sum of two terms as follows:

$$\begin{array}{}\text{(2)}& {\mathrm{PO}}_{\mathrm{4}}={\mathrm{PO}}_{\mathrm{4}}^{\mathrm{pre}}+{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{rem}},\end{array}$$

where PO${}_{\mathrm{4}}^{\mathrm{rem}}$ is the amount of PO₄ locally produced by organic matter remineralization.

In the NorESM1-F model, preformed tracers are set to their respective total values in the ocean mixed layer (i.e., the upper two levels of the ocean model) at each time step. Below, the preformed tracers are advected as passive (i.e., conservative) tracers by the ocean circulation (Tjiputra et al., 2020). They can thus be used to track changes in ocean circulation and water mass geometry below the surface ocean.

In this study, we use the preformed PO₄ computed by NorESM1-F to deconvolve the GI8–HS4 δ¹³C-BIO change, Δδ¹³C-BIO, into a fraction resulting from remineralization changes and the remainder, i.e., changes in primary production (PP) in the surface ocean and changes in ocean circulation. Importantly, in what follows, we use the term ocean circulation changes to refer to all the ocean circulation changes that are reflected by water mass tracers – that is, the changes in water mass origin and geometry but not the changes in water mass flow rate. The latter are reflected in the computed ideal age of the water mass, defined as the time since the water mass last made contact with the surface (Guo et al., 2019b).

Combining Eqs. (1) and (2) and subtracting HS4 values from GI8 values, we obtain

$$\begin{array}{}\text{(3)}& \begin{array}{rl}\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}\text{-}\mathrm{BIO}& ={\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{photo}}}{{\mathrm{DIC}}_{\mathrm{mean}\mathrm{o}}}}\cdot {\displaystyle \frac{\mathrm{C}}{\mathrm{P}}}\cdot \left(\mathrm{\Delta }{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{pre}}+\mathrm{\Delta }{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{rem}}\right)\\ & ={\left(\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}\text{-}\mathrm{BIO}\right)}_{\mathrm{circ}+\mathrm{PP}}+{\left(\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}\text{-}\mathrm{BIO}\right)}_{\mathrm{rem}},\end{array}\end{array}$$

where (Δδ¹³C-BIO)_circ+PP is the portion of Δδ¹³C-BIO resulting from ocean circulation and surface PP changes, and (Δδ¹³C-BIO)_rem is the portion of Δδ¹³C-BIO resulting from remineralization changes in the ocean interior.

3 Results

In order to compare the simulated GI8–HS4 change in δ¹³C-BIO with the observed change in Cib. δ¹³C, we compute the mean observed and simulated values just before and just after the transition. Model outputs are averaged over a 100-year pre-transition stadial period and a 100-year post-transition interstadial period, respectively defined as model years 5700–5800 of the simulation with freshwater forcing and model years 5950–6050 of the simulation in which the freshwater input has been halted. These two time intervals correspond to relatively stable conditions, as attested by the low standard deviations of the simulated values within each interval (Table S2).

Due to the much-lower temporal resolution of the Cib. δ¹³C proxy records compared to the yearly model output, we average the Cib. δ¹³C records over 500-year periods to increase the number of records included in the study. Hence, HS4 and GI8 Cib. δ¹³C are averaged over the 38.5–39.0 and 37.5–38.0 ka intervals, respectively (Table 1). We estimate the uncertainty associated with these HS4 and GI8 Cib. δ¹³C values by combining the uncertainty resulting from the dispersion of the Cib. δ¹³C measurements within the 500-year intervals and the dating uncertainty (see Sect. S3 for more details).

Table 1 Partitioning of computed δ¹³C-BIO changes across the HS4–GI8 transition, Δδ¹³C-BIO, into (Δδ¹³C-BIO)_circ+PP and (Δδ¹³C-BIO)_rem components (Eq. 3).

^* Cores are ordered by increasing water depths and are separated into two groups: above and below 4200 m (see text). Green indicates a good trend, but the computed GI8-HS4 δ¹³C-BIO is smaller than the observed amplitude; red indicates model–data disagreement; bolding indicates significant changes.

Among the 18 selected sites (Figs. 1 and S1, Table 1), we obtain a strict agreement between the change in simulated δ¹³C-BIO and the observed Cib. δ¹³C in 10 out of 18 sites within ±1σ (blue symbols in Fig. 1). In five out of the remaining eight sites, the simulated trend in δ¹³C-BIO is consistent with that in Cib. δ¹³C, but the model underestimates the GI8–HS4 Cib. δ¹³C change (green symbols in Fig. 1). Finally, simulated δ¹³C-BIO and observed Cib. δ¹³C changes disagree in three sites (red symbols in Fig. 1). A regression analysis of Δδ¹³C-BIO versus ΔCib. δ¹³C values for (i) all sites, (ii) only sites depicted in blue and green, and (iii) only sites depicted in blue shows that the linear regressions are all highly significant, with probabilities of no linear correlation lower than 0.007 in all three cases and correlation coefficients (Pearson's) ranging from 0.63 to 0.82 from case (i) to (iii) (Fig. S2).

Figure 1 (a) Change in δ¹³C-BIO simulated by NorESM1-F versus observed change in Cib. δ¹³C across the HS4-to-GI8 transition. Plotted error bars are ±1σ. (b, c) Position of the different core sites.

Among the three sites exhibiting model–data disagreement (red symbols), two are from the southeast Atlantic (Fig. 1). It is noteworthy that the southeast Atlantic Ocean is a problematic region for epibenthic isotopic time series for the following two reasons: (i) dating is difficult due to variable surface reservoir ages and the lack of sharp temperature changes in the Antarctic ice core records that can serve as alignment targets for sea surface temperature records; and (ii) the glacial δ¹³C value of the epibenthic genus Cibicides has been shown to be much lower for the species Cibicides kullenbergi than for the species Cibicides wuellerstorfi in certain instances (Gottschalk et al., 2016b), pointing to different habitat preferences for these two species and the presence of local phytodetritus layers on the sediment–water interface, possibly with very low δ¹³C-DIC values (Mackensen et al., 1993). The third site where observed and computed δ¹³C amplitudes disagree, CAR13-05, is located in the northeast Atlantic (Table 1). Although it passed the screening test for dating quality, its sedimentation rate is relatively low over the studied time interval (6 to 10 cm kyr⁻¹; Fig. S3). The closest core in terms of latitude, longitude, and water depth, MD03-2698 (Table 1), has a sedimentation rate of about 50 cm kyr⁻¹ (Fig. S3) and shows an excellent data–model agreement. Since marine sediment core age model precision and reliability are much improved for higher sedimentation rates, we may consider it to be that the GI8–HS4 Cib. δ¹³C change obtained for core MD03-2698 is more reliable than that obtained for core CAR13-05.

From the modeling point of view, the 1^∘ spatial resolution of the NorESM1-F model is not sufficient to correctly simulate a number of physical processes, like the coastal upwelling or deep overflows, with potential consequences for the simulated ocean circulation. Similarly, 1^∘ ocean models cannot resolve the leaking Agulhas rings, which are crucial for the heat and salinity budget in the southeast Atlantic region. The Antarctic Circumpolar Current region is also prone to model bias (Beadling et al., 2020), and spurious open-ocean convection near the Weddell Sea (Heuzé, 2021) can lead to water mass biases in the Southern Ocean and beyond. In addition to these physical uncertainties, some uncertainties could also arise from the application to the glacial period of the BGC module with the same ecosystem parameterization as for the present, whereas there is evidence that some parameters, like the remineralization rate, should be temperature dependent (Brewer and Peltzer, 2017).

However, even though two of the three sites exhibiting model–data disagreement (red symbols) are located in the southeast Atlantic, no firm conclusion can be drawn regarding possible data or model biases due to the limited number of observations in the South Atlantic. When examining the geographical distribution of the three categories of sites, we find (i) model–data agreement (blue symbols) distributed over all latitudes, longitudes, and water depths; (ii) underestimated simulated changes (green symbols) in the North Atlantic, both along the western boundary current and the mid-Atlantic ridge and on the Iberian margin, at various water depths; and (iii) model–data disagreement (red symbols) below 3700 m in the North and South Atlantic, together with several blue sites and one green site. Therefore, the different categories of model–data agreement and/or disagreement do not correspond to any geographical pattern. This lack of geographical pattern points to the absence of obvious biases in the model results, at least in the North Atlantic.

To summarize, we obtain an overall good agreement between the change in δ¹³C-BIO simulated by NorESM1-F and the observed change in Cib. δ¹³C across the HS4-to-GI8 transition below 2000 m in the North Atlantic and equatorial West Atlantic (Fig. 1, Table 1).

4 Discussion

This overall good model–data agreement suggests that the model successfully reproduces the changes in ocean nutrients across a stadial–interstadial transition below 2000 m in the Atlantic Ocean. Moreover, it warrants the use of NorESM1-F simulations to examine what caused the observed changes in Cib. δ¹³C in that portion of the world ocean.

Figure 2 NorESM1-F simulated changes in Atlantic ocean circulation and nutrients across the HS4-to-GI8 transition. Pre- and post-transition values respectively correspond to 100-year averages over model years 5700–5800 of the simulation with freshwater forcing and model years 5950–6050 of the simulation in which the freshwater input has been halted. (a, b) Zonally integrated stream function. (c–h) Section plots along the western Atlantic section depicted in the inset. Colored symbols indicate the core sites as plotted in Fig. S1 and defined in Table 1.

Prior to the stadial–interstadial transition, the simulated AMOC stream function is weak, with an upper overturning cell restrained to the upper 2000 m and a maximum overturning strength of about 10 Sv (Figs. 2a and S4). After the transition, the model results exhibit a strong overturning regime, with active deep-water formation in the northern North Atlantic and Nordic Seas leading to an upper overturning cell reaching down to ∼3000 m depth and a maximum overturning strength larger than 25 Sv (Figs. 2b and S4). This evolution of the simulated AMOC stream function is in agreement with the changes in overturning strength derived from the sedimentary Pa/Th records covering the HS4-to-GI8 transition (Henry et al., 2016; Waelbroeck et al., 2018). Moreover, even though the change in benthic δ¹⁸O across the HS4-to-GI8 transition is very small and not significant for most sites, the simulated bottom-water-density decrease resulting from a slight bottom-water warming is in agreement with the observed benthic δ¹⁸O decrease at the deep northeast Atlantic site U1308 (Sect. S4).

Figure 3 Same as Fig. 2 but for simulated PO tracer [mol O₂ m⁻³] at 2500, 3500, and 4500 m depth in the Atlantic Ocean.

We use the simulated PO tracer defined as $\mathrm{PO}=\mathrm{172}\cdot {\mathrm{PO}}_{\mathrm{4}}+{\mathrm{O}}_{\mathrm{2}}$ (Luo et al., 2018), similarly to the NO tracer first defined by Broecker (1974), to track the water mass origin. The PO tracer combines O₂ and PO₄ in such proportions that the increase in PO₄ by respiration is canceled, which makes this tracer nearly conservative. Given the difference in the preformed phosphate content of dense water produced in the northern and southern high latitudes of the Atlantic Ocean, SSWs and NSWs are characterized by high and low PO values, respectively. In our simulations of HS4 and GI8 ocean circulation, SSWs correspond to PO values above ∼0.55 mol O₂ m⁻³, and NSWs correspond to PO values below ∼0.45 mol O₂ m⁻³ (Fig. 2c and d).

Vertical sections of the PO tracer along the western basin before and after the HS4-to-GI8 transition show that NSWs were confined to the upper 2000 m before the transition, whereas active deep-water formation in the northern North Atlantic and Nordic Seas is clearly visible after the transition, with low PO values delimiting a North Atlantic Deep Water (NADW)-like water mass (Fig. 2c and d). This is further illustrated in horizontal PO maps at different water depths (Fig. 3), where the resumption of the southward-flowing deep western boundary current in the North Atlantic is reflected by low PO values down to 3500 m depth in the post-transition sections (Fig. 3e–h).

Importantly, as expected, changes in the spatial distribution of the PO tracer across the HS4-to-GI8 transition are remarkably similar to those of the preformed PO₄, both vertically (Fig. 2c–f) and horizontally (Figs. 3 and S5). This constitutes a validation of our use of PO${}_{\mathrm{4}}^{\mathrm{pre}}$ to assess the portion of the change in δ¹³C-BIO resulting from changes in ocean circulation and surface PP and further indicates that the changes in surface PP do not lead to a visible imprint on the interior ocean δ¹³C-BIO distribution. The simulated PP increase in some regions across the HS4-to-GI8 transition is indeed too weak to significantly draw down the PO${}_{\mathrm{4}}^{\mathrm{pre}}$ values and is not sufficient to prevent PO${}_{\mathrm{4}}^{\mathrm{pre}}$ from increasing in the same regions across the transition (Fig. S6). For simplicity, here, we thus simply interpret PO${}_{\mathrm{4}}^{\mathrm{pre}}$ in terms of changes in ocean circulation – that is, changes in water mass origin and geometry.

In addition, the NorESM1-F model computes the total oxygen utilization (TOU; Figs. 2g, h and 4), that is, the total oxygen consumed by organic matter remineralization, which translates into the amount of PO₄ locally produced by remineralization, PO${}_{\mathrm{4}}^{\mathrm{rem}}$ (Eq. 2). Increases in TOU or PO${}_{\mathrm{4}}^{\mathrm{rem}}$ thus correspond to decreases in δ¹³C-BIO. Although the spatial distribution of the TOU before and after the HS4-to-GI8 transition is imprinted by the low TOU values of the NSW in the North Atlantic, it is markedly different from that of ocean circulation changes shown by PO${}_{\mathrm{4}}^{\mathrm{pre}}$ or the PO tracer in the South Atlantic (Figs. 2–4).

Figure 4 Same as Fig. 3 but for simulated total oxygen utilization (mol O₂ m⁻³) at 2500, 3500, and 4500 m depth in the Atlantic Ocean.

To summarize, PO${}_{\mathrm{4}}^{\mathrm{pre}}$ (or the PO tracer) on the one hand and PO${}_{\mathrm{4}}^{\mathrm{rem}}$ (or TOU) on the other hand give us access to the partitioning of the δ¹³C-BIO change across the HS4-to-GI8 transition into the following two constituents: (i) the fraction change resulting from changes in water mass origin and geometry, and (ii) the complementary fraction caused by changes in remineralization (Eq. 3).

Based on the simulated PO${}_{\mathrm{4}}^{\mathrm{pre}}$ averaged over the previously defined HS4 and GI8 100-year intervals, we compute the portion of the GI8–HS4 δ¹³C-BIO change due to ocean circulation changes, (Δδ¹³C-BIO)_circ+PP, at the 18 selected sites. We find that the 11 sites located above 4200 m that exhibit significant GI8–HS4 δ¹³C-BIO changes all show an increase in δ¹³C-BIO across the transition. In these sites, (Δδ¹³C-BIO)_circ+PP forms 30 % to 90 % of the total GI8–HS4 δ¹³C-BIO increase (Table 1). Furthermore, we find a loose, but significant (R²=0.64), linear correlation, whereby the fraction change resulting from ocean circulation changes increases from ∼40 % at 2000 m to ∼80 % at 4000 m (Fig. S7). (Δδ¹³C-BIO)_rem thus forms 10 % to 70 % of the total GI8–HS4 δ¹³C-BIO increase in these 11 sites (Table 1), with this fraction increasing with depth from ∼20 % at 4000 m to ∼60 % at 2000 m. The model results further show that the increase in δ¹³C-BIO_rem across the transition results from a decrease in remineralization explained by the stronger overturning circulation during GI8 than during HS4, which translates into shorter water mass residence times at these 11 sites, as shown by the decrease in ideal age at these sites across the transition (Fig. S8).

This increase in δ¹³C-BIO_circ+PP at the core sites located above 4200 m corresponds to an increase in the ratio of NSW versus SSW, as illustrated by the PO maps (Fig. 3). The largest changes occur in sites where SSW is replaced by NSW as a result of the resumption of deep-water formation in the northern high latitudes and the penetration of NADW-like NSW down to 3500 m in the northwest Atlantic and toward 30^∘ S in the South Atlantic around 2500 m (Fig. 3).

Contrary to the δ¹³C-BIO increase of up to 0.80 ‰ found in sites above 4200 m, δ¹³C-BIO decreases by 0.02 ‰ to 0.12 ‰ across the HS4-to-GI8 transition in the six sites located below 4200 m (Fig. S9). At these deep sites, remineralization increases across the HS4-to-GI8 transition, leading to a 0.04 ‰ to 0.18 ‰ decrease in δ¹³C-BIO_rem (Eq. 3, Table 1). In contrast, the change in δ¹³C-BIO_circ+PP at these six deep sites is very small and barely significant (within −0.03 ‰ to +0.02 ‰), except in two southeast Atlantic sites, where it reaches +0.06 ‰ and +0.08 ‰but is counteracted by larger decreases in the δ¹³C-BIO_rem component, reaching −0.18 ‰ and −0.14 ‰, respectively (Table 1). Such a decrease in δ¹³C-BIO_rem results from an increase in respired PO₄ and hence in organic matter export from the euphotic zone to the deep ocean during GI8 with respect to HS4. In particular, in our two subantarctic sites located below 4200 m, NorESM1-F simulates higher remineralization during interstadials than during stadials, which is consistent with an increased export production from the upper ocean to the deep Southern Ocean during interstadials, as demonstrated in a number of studies (e.g., Gottschalk et al., 2016a; Martínez-García et al., 2014).

In addition to the above discussion of our model results at the selected 18 sites, the simulated maps of Δδ¹³C-BIO and of its partitioning into (Δδ¹³C-BIO)_circ+PP and (Δδ¹³C-BIO)_rem provide acceptable approximations of the change in Cib. δ¹³C and of its partitioning into a circulation and surface PP term on one hand and a remineralization term on the other hand anywhere below 2000 m in the Atlantic Ocean (Fig. S10).

It is tempting to compare our results with a recent study of the LGM-to-HS1 transition by Gu et al. (2021), who concluded that the δ¹³C-DIC decrease observed at mid-depth is mainly explained by increased remineralization due to AMOC slowdown, while the water mass mixture change plays only a minor role. However, Gu et al. (2021) used a transient simulation of the last deglaciation to examine the causes of the mid-depth Atlantic δ¹³C-DIC decrease observed across the transition from the LGM into the HS1 stadial, whereas we use a hosing experiment under constant radiative forcing to examine the transition from the HS4 stadial to the GI8 interstadial, away from the last deglaciation and its large changes in radiative forcing. The results of the two studies can thus hardly be compared. Despite these different settings, both studies show that the remineralization changes due to ventilation changes are the main factor contributing to the observed Atlantic δ¹³C-DIC change at about 2000 m water depth. But this agreement does not hold at greater depths, where our results suggest that water mass mixture change is the main factor of the Atlantic δ¹³C-DIC change between ∼2500 and 4000 m.

5 Conclusions

We show that the observed Cib. δ¹³C changes across the HS4-to-GI8 stadial–interstadial transition can be successfully approximated by the computed change in dissolved inorganic carbon δ¹³C, Δδ¹³C-BIO, assuming no change in air–sea exchanges and using the nutrient fields simulated by the NorESM1-F model across a stadial–interstadial transition.

The model results can then be used to gain insight into the causes of the observed Cib. δ¹³C change across the HS4-to-GI8 transition and to assess the fraction of the change resulting from changes in water mass origin and geometry on the one hand and that associated with changes in organic matter remineralization in the interior ocean on the other hand. Above ∼4200 m, a large fraction of the computed δ¹³C-BIO increase results from the replacement of SSW by NSW subsequent to the resumption of deep-water formation in the northern high latitudes after the freshwater input is halted. This fraction exhibits a loose linear relationship with depth, increasing from ∼40 % at 2000 m to ∼80 % at 4000 m. Conversely, the δ¹³C-BIO increase resulting from the decrease in remineralization across the HS4-to-GI8 transition in response to the faster flow rates during GI8 than during HS4 represents ∼60 % of the δ¹³C-BIO increase at 2000 m but only ∼20 % of the δ¹³C-BIO increase at 4000 m.

Below ∼4200 m, the model shows little change in δ¹³C-BIO resulting from changes in water mass origin and geometry but an increase in remineralization across the HS4-to-GI8 transition. The simulated lower remineralization during stadials than during interstadials is particularly pronounced in our deep subantarctic sites, which is in agreement with the decrease in the biological export production during stadials found in the Southern Ocean in previous studies.