Articles | Volume 16, issue 3
Research article
19 May 2020
Research article |  | 19 May 2020

Carbon isotopes and Pa∕Th response to forced circulation changes: a model perspective

Lise Missiaen, Nathaelle Bouttes, Didier M. Roche, Jean-Claude Dutay, Aurélien Quiquet, Claire Waelbroeck, Sylvain Pichat, and Jean-Yves Peterschmitt

Understanding the ocean circulation changes associated with abrupt climate events is key to better assessing climate variability and understanding its different natural modes. Sedimentary Pa∕Th, benthic δ13C and Δ14C are common proxies used to reconstruct past circulation flow rate and ventilation. To overcome the limitations of each proxy taken separately, a better approach is to produce multiproxy measurements on a single sediment core. Yet, different proxies can provide conflicting information about past ocean circulation. Thus, modelling them in a consistent physical framework has become necessary to assess the geographical pattern and the timing and sequence of the multiproxy response to abrupt circulation changes.

We have implemented a representation of the 231Pa and 230Th tracers into the model of intermediate complexity iLOVECLIM, which already included δ13C and Δ14C. We have further evaluated the response of these three ocean circulation proxies to a classical abrupt circulation reduction obtained by freshwater addition in the Nordic Seas under preindustrial boundary conditions. The proxy response is shown to cluster in modes that resemble the modern Atlantic water masses. The clearest and most coherent response is obtained in the deep (> 2000 m) northwest Atlantic, where δ13C and Δ14C significantly decrease, while Pa∕Th increases. This is consistent with observational data across millennial-scale events of the last glacial. Interestingly, while in marine records, except in rare instances, the phase relationship between these proxies remains unclear due to large dating uncertainties, in the model the bottom water carbon isotope (δ13C and Δ14C) response lags behind the sedimentary Pa∕Th response by a few hundred years.

1 Introduction

Understanding rapid climate changes is key to understand the different natural modes of ocean circulation as they provide information about internal climate variability and climate sensitivity to perturbations. Indeed, during the last glacial, rapid and high-amplitude atmospheric temperature changes over Greenland of about 8 to 15 C in less than 300 years are associated with only small changes in radiative forcing (Kindler et al., 2014). Changes in the Atlantic meridional overturning circulation (AMOC) strength are thought to be one of the main drivers of these abrupt climate events (see Lynch-Stieglitz, 2017, for a review), but the underlying mechanisms remain elusive. Available AMOC records based on observations only cover a few decades, and thus identifying the processes controlling AMOC changes can only be achieved from the study of long-term variations (Smeed et al., 2014), which relies on the analysis of indirect evidence (palaeoproxies).

To date, among the numerous tracers available (e.g. benthic δ18O, Nd isotopes, Cd∕Ca or sortable silts), the sedimentary (231Paxs,0/230Thxs,0) ratio (hereafter Pa∕Th) and carbon isotopes (δ13C, Δ14C) are key tracers used to reconstruct and quantify past circulation patterns and water mass flow rates.

Pa∕Th is the ratio of the 231Pa and 230Th activities (decays per time unit) at the deposition time that are derived from the water column scavenging. The Pa∕Th ratio can be used as a kinematic circulation proxy (François, 2007; McManus et al., 2004). In short, 231Pa and 230Th are homogeneously produced in the water column at known rates with a known production ratio of 0.093 and are then transferred to the underlying sediments by particle scavenging (see François, 2007, for a review). The two isotopes have different residence times in the water column (50–200 years for 231Pa and 10–40 years for 230Th; Henderson and Anderson, 2003). Thus, while 230Th is rapidly transferred to the sediment, 231Pa can be partly transported by water mass advection along the large-scale ocean circulation. Consequently, the sedimentary Pa∕Th activity ratio can be used as a proxy of water-mass advection rate. In the Atlantic, low sedimentary Pa∕Th ratios (e.g. 0.04 for the modern North Atlantic; Yu et al., 1996) are diagnostic of an active overturning, while Pa∕Th ratios close or equal to the production ratio indicate a sluggish water mass or a marked overturning circulation slowdown (e.g. Böhm et al., 2015; François, 2007; McManus et al., 2004; Waelbroeck et al., 2018). Yet, Pa and Th scavenging to the sediment is sensitive to changes in vertical particle flux and composition; hence the sedimentary Pa∕Th circulation signal could be partly biased by a particle-related signal (e.g. Chase et al., 2002, 2003; Lippold et al., 2009). In addition, sedimentary Pa∕Th is derived from bulk sediment measurements and requires the estimation of the contributions of the detrital and authigenic fractions to the 231Pa and 230Th budgets. We have recently shown that this estimation can lead to significant uncertainties in the reconstructed patterns and amplitudes of the Pa∕Th signal, especially in locations characterised by high terrigenous inputs (Missiaen et al., 2018). These potential caveats need to be tested and could complicate the evaluation of past circulation strength changes from Pa∕Th measurements.

The carbon isotopes measured in foraminifer shells reflect the carbon isotopic content of the water mass in which they form and provide information about past water mass ventilation, in other words, the time elapsed since the tracked deep-water parcel has been isolated from the surface (e.g. Lynch-Stieglitz et al., 2014; Skinner et al., 2014).

More precisely, the water mass δ13C signature depends on biological and physical processes. Carbon is exchanged between the surface waters and the atmosphere. In the surface waters, carbon is incorporated into the organic matter through biologic productivity. Both air–sea exchange and biological activity are responsible for isotopic fractionation between 12C and 13C (Siegenthaler and Münnich, 1981). Consequently, the surface water δ13C signature varies with air–sea exchange efficiency and biological activity intensity. At depth, remineralisation of organic matter releases 13C-depleted carbon to the water parcels, which are then mixed through large-scale ocean circulation. In the modern ocean, the global δ13C distribution depicts a tight relation between the apparent oxygen utilisation and the δ13C signature of a given water mass (Eide et al., 2017). This observation lends support to the use of the δ13C of benthic foraminifers as a proxy for ocean oxygen content and ventilation (Duplessy et al., 1988). However, benthic δ13C does not solely record deep ventilation changes. As mentioned above, the δ13C signature of a deep-water mass depends on several processes: its value before the water left the surface mixed layer, the intensity of the biological activity in the mixed layer, the remineralisation intensity at depth, and finally the circulation path and strength. Thus, benthic δ13C records multiple processes, which complicates its interpretation in terms of past deep ocean ventilation and circulation changes.

Radiocarbon is produced in the upper atmosphere and enters into the ocean via air–sea exchange with surface waters. As soon as a water parcel is isolated from the surface, its 14C content starts to decrease exponentially with time due to radioactive decay (half-life of 14C =5730±40 years; Godwin, 1962). Thus, by determining the 14C age of benthic foraminifer samples of independently known calendar age, one can reconstruct past ocean ventilation, (e.g. Skinner and Shackleton, 2004; Thornalley et al., 2015). However, the interpretation of a water mass radiocarbon age is complicated by temporal variations in 14C production in the upper atmosphere and air–sea exchange efficiency. Indeed, since radiocarbon is only produced in the atmosphere and transferred to the ocean via air–sea exchange, the surface waters have an older radiocarbon age than the contemporaneous atmosphere. This radiocarbon age difference between the surface waters and the atmosphere (termed the surface reservoir age) can vary with space and time according to variations in air–sea exchanges efficiency, especially in the North Atlantic region (see Bard et al., 1994; Bondevik et al., 2006; Thornalley et al., 2011; Waelbroeck et al., 2001). Those variations are still poorly constrained and complicate the interpretation of deep water radiocarbon content (expressed as Δ14C in ‰) in terms of past ocean ventilation and circulation changes (e.g. Adkins and Boyle, 1997).

Pa∕Th, benthic δ13C and Δ14C can provide useful information about past ocean circulation flow rate, geometry and ventilation. Yet, as highlighted above, each proxy has its own caveats. To overcome the limitation of each proxy taken separately and gather more detailed information about past ocean circulation, palaeoceanographers started to conduct multiproxy studies, producing different proxy records on the same sedimentary archive. Moreover, this approach also enables one to investigate phase relationships between the different proxies (Burckel et al., 2015; Waelbroeck et al., 2018). Indeed, current dating uncertainties in marine cores are generally less than 150 years between 0 and 11 ka cal BP, increasing to ∼400 years between 11 and 30 ka cal BP and ranging from ∼600 to 1100 years between 30 and 40 ka cal BP (Waelbroeck et al., 2019). Such uncertainties usually prevent the inference of phase relationships between proxy records from different marine cores, hence limiting the benefit of the multiproxy approach.

Different proxies can bring conflicting information about past ocean circulation, and reconstructing basin-scale ocean water-mass reorganisation requires the compilation of proxy records from different locations. Despite some recent palaeoproxy compilation efforts (Lynch-Stieglitz et al., 2014; Ng et al., 2018; Zhao et al., 2018), the amount of data available remains too sparse to constrain the state and evolution of the ocean circulation in 3D across abrupt climate events. One way to better understand existing palaeo-records and overcome the scarcity of palaeodata is to use climate models that are able to simulate the evolution of different proxies in a consistent physical framework. Such work could help to explain why some events are not recorded or recorded differently at a given location. Climate models are also useful as they enable the analysis of the spatial and temporal response of a proxy to changes in ocean circulation.

Because of its key role in the climate system, the carbon cycle has been modelled and heavily studied in the past decades (Friedlingstein et al., 2006; Orr et al., 2001). More recent studies simulate the evolution of carbon isotopes, and in particular 13C and 14C, under preindustrial or glacial conditions (Bouttes et al., 2015; Brovkin et al., 2002; Menviel et al., 2017; Tschumi et al., 2011). The 231Pa and 230Th tracers have also been implemented in climate models. The simplest approaches used 2D models (Luo et al., 2010; Marchal et al., 2000) or 3D models but contained oversimplifications, notably in the particles representations (Siddall et al., 2005, 2007). Latest published developments focused on an improved representation of particle fluxes and scavenging scheme (Gu and Liu, 2017; van Hulten et al., 2018; Rempfer et al., 2017). However, these recent developments either suffer from the coarse resolution of the ocean model (Rempfer et al., 2017), which contains only 36×36 grid cells (latitude–longitude), or conversely cannot simulate 1000 years in reasonable computation time (van Hulten et al., 2018). To our knowledge, there has not been any study so far that considered the geographical distribution and the temporal evolution of combined proxies such as Pa∕Th, δ13C and Δ14C.

We further investigate the spatial and temporal structure of multiproxy response to an abrupt circulation change with a climate model. For that purpose, we have implemented the production and scavenging processes of 231Pa and 230Th in the climate model of intermediate complexity iLOVECLIM. To date, the model is able to simulate the evolution of three ocean circulation proxies: Pa∕Th, δ13C and Δ14C. We also developed an analysis method that gives information about the magnitude and timing of the proxy response in the Atlantic Ocean. In this study, we address the following questions. (1) What is the response of each proxy to an imposed circulation change? (2) What is the timing and sequence of the proxy responses? How do they vary with regions and water depth in the Atlantic Ocean? (3) How can the modelled multiproxy response help to interpret the palaeoproxy records?

Table 1Parameters used in the Pa∕Th module.

Download Print Version | Download XLSX

2 Material and methods

2.1 Model description and developments

We use the Earth system model of intermediate complexity iLOVECLIM, which is a code fork of the LOVECLIM model (Goosse et al., 2010). It includes a representation of the atmosphere, ocean, sea ice, terrestrial biosphere vegetation and the carbon cycle. The ocean component (CLIO) consists of a free-surface primitive equation ocean model (3×3 horizontal grid, 20 depths layers) coupled to a dynamic-thermodynamic sea-ice model. iLOVECLIM includes a land vegetation module (VECODE) (Brovkin et al., 1997) and a marine carbon cycle model (Bouttes et al., 2015), both computing the evolution of 13C and 14C. Previous work has shown that the simulated oceanic δ13C and Δ14C distribution is in reasonable agreement with observations and with available general circulation model (GCM) simulations (Bouttes et al., 2015). We have implemented 231Pa and 230Th in the iLOVECLIM model following the approach of Rempfer et al. (2017), albeit without including an explicit parameterisation of boundary and bottom scavenging as detailed in the following. In the model, 231Pa and 230Th are homogeneously produced in the ocean by radioactive decay of their respective uranium parents. They are removed from the water column by adsorption on settling particles (reversible scavenging). Their radioactive decay is also taken into account. In its current version, iLOVECLIM does not explicitly simulate the biogeochemical cycle of biogenic silica (opal), which is thought to be an important scavenger for 231Pa (e.g. Chase et al., 2002; Kretschmer et al., 2010). Therefore, we used prescribed and constant particle concentration fields obtained running the PISCES biogeochemical model coupled to NEMO ocean model (van Hulten et al., 2018). Previous work has shown that these particle fields are reasonably consistent with present-day observations (van Hulten et al., 2018) and that a set-up with fixed particles is able to capture the major features (Gu and Liu, 2017). Like in Bern 3D (Rempfer et al., 2017; Siddall et al., 2005, 2007) and CESM (Gu and Liu, 2017), we considered one particle size with a unique sedimentation speed of 1000 m yr−1. We consider that Pa and Th interact with three particles types: CaCO3, POC and biogenic silica. The role of lithogenic particles in Pa and Th scavenging having been questioned (van Hulten et al., 2018; Siddall et al., 2005). The conservation equations for dissolved and particulate 231Pa and 230Th activities are the following (Rempfer et al., 2017):


where Adj and Apj are respectively the dissolved and particle-bound activities (dpm m−3 yr−1) of isotope j (231Pa or 230Th), βj (dpm m−3 yr−1) is the water column production of isotope j by radioactive decay of its uranium parent, λj (yr−1) is the decay constant of isotope j, ws (m yr−1) is the particle settling speed, Kadsorpj and kdesorpj (yr−1) are the adsorption and desorption coefficient of isotope j onto particles, respectively, and T is the tracer balance evolution term (dpm m−3 yr−1) resulting from the water mass advection and diffusion terms computed by the CLIO ocean model. The values used in this study for each of the above-cited parameters are compiled in Table 1.

As in Bern 3D (Rempfer et al., 2017), we chose to apply a uniform desorption coefficient denoted Kdesorp hereafter. However, the adsorption coefficient depends in our model on the particle concentration and composition of each location and is calculated with the following equation (Rempfer et al., 2017):

(3) K adsorp j θ , Φ , z = i σ i , j F i ( θ , Φ , z ) ,

where θ is the latitude, Φ the longitude, z the water depth, σi,j the scavenging efficiency for isotope j on particle type i (m2 mol−1) and Fi is the particle flux (mol m−2 yr−1).

2.2 Model tuning and validation

The scavenging efficiency σi,j is related to the partition coefficient Kd, which defines the proportion of each isotope (231Pa or 230Th) lodged in the dissolved phase or bound to particles. Therefore one Kd can be defined for each isotope and each particle type considered in the model.

(4) K d ( i , j ) = σ i , j w s ρ sw M ( i ) k desorp ,

where Kd(i,j) is the partition coefficient for isotope j for particle type i, σi,j is the corresponding scavenging efficiency, ws is the settling speed, Kdesorp is the desorption coefficient, M(i) is the molar mass of particle type i (i.e. 12 g mol−1 for POC, 100.08 g mol−1 for CaCO3 and 67.3 g mol−1 for opal) and ρsw is the mean density of seawater (constant fixed to 1.03×106 g m−3). Additionally the way each particle type fractionates the Pa and the Th is defined by the fractionation factor F(Th/Pa)i.

(5) F ( Th / Pa ) i = K d ( Th , i ) K d ( Pa , i ) = σ ( Th , i ) σ ( Pa , i )

Kd and fractionation factors, F(Th∕Pa), have been measured for both radionuclides in various areas of the modern ocean, and they show a rather large distribution (see Table S1 in the Supplement) (Chase et al., 2002; Hayes et al., 2015). Consequently, these values are currently considered as tunable parameters in modelling studies (Dutay et al., 2009; Gu and Liu, 2017; van Hulten et al., 2018; Marchal et al., 2000; Rempfer et al., 2017; Siddall et al., 2005). Considering three particle types for both radionuclides, there are thus six tunable σi,j parameters in our model.

To efficiently sample our parameter space, we used a Latin hypercube sampling (LHS) methodology (, last access: 17 April 2019). In order to only select parameters values consistent with observed F(Th∕Pa), we chose to use the couples {σTh,j, F(Th/Pa)i} as input parameters for the LHS. The value of σPa,i is then deduced from those two following equation (Eq. 5). The parameter ranges used in the LHS are given in Table S2. The LHS allowed a relatively good exploration of the parameter space with a relatively small number of model evaluations. We performed 60 tuning simulations of 1000 years each under preindustrial boundary conditions (i.e. orbital parameters, ice-sheets configuration and atmospheric gas concentrations), as described in (Goosse et al., 2010).

Table 2Best fit σi,j values and corresponding Kd values.

Download Print Version | Download XLSX

Consistent with previous modelling studies (van Hulten et al., 2018), this simulation length was sufficient for Pa and Th to reach an approximate steady state at the surface and in the deeper ocean. The model performance was evaluated by comparing outputs with present-day particulate and dissolved water column Pa and Th measurements compiled in the GEOTRACES intermediate data product 2017 (Schlitzer et al., 2018) as well as sedimentary Pa∕Th core tops data as compiled in (van Hulten et al., 2018) and references therein. We selected the ensemble member that best fits the observational constraints using five metrics corresponding to the root-mean-square error (RMSE) between observation data and the closest model grid cell average for particulate and dissolved Pa and Th as well as sedimentary Pa∕Th (see Sect. S1 in the Supplement for additional information). Table 2 presents the best-fit σi,j values. The best fit simulation is then used to investigate the multiproxy response to abrupt circulation changes.

2.3 Experimental design

With the best fit σi,j values (see Table 2), we ran our model for 5000 years under preindustrial (PI) conditions (as described in Goosse et al., 2010) from a simulation with an equilibrated carbon cycle (Bouttes et al., 2015). The result of this equilibrium simulation is used as a starting point to perform hosing experiments of 1200 years duration. The freshwater was added in the Nordic Seas following the approach described in (Roche et al., 2010). Each simulation contains three phases: a control phase (300 years), a hosing phase (300 years) and a recovery phase (600 years). The control phase is used to assess the natural variability of the circulation and associated proxies under the PI climate state.

3 Results

3.1 Model-data comparison under preindustrial conditions

The main goal of our study is to assess the response of two carbon-based proxies (δ13C, Δ14C) and of the Pa∕Th to an abrupt circulation change in a physically consistent framework. This work represents a first step toward a better understanding of marine multiproxy records across the last glacial abrupt events. Therefore, our model–data comparison focuses on the Atlantic Ocean, where the Pa∕Th can be used as a kinematic circulation proxy (e.g. François, 2007; McManus et al., 2004). and on the simulated bottom particulate Pa and Th activities which can be directly compared to the Pa∕Th ratio recorded in marine sediments. The water column concentration results are presented in Sect. S2 and Figs. S1 to S3 in the Supplement.

Figure 1Map showing the particulate 231Pa∕230Th activity ratio of the deepest model ocean grid cells. The simulated Pa∕Th ratio is represented in the colour background. The observations compiled in van Hulten et al. (2018) and references therein are represented as circles.

The particle-bound Pa∕Th of the deepest oceanic model grid cells is shown in Fig. 1. The observations from core-top data are superimposed as circles. Even if the observational dataset is somewhat patchy, it generally shows lower sedimentary Pa∕Th ratios (0.04 or lower) in the basin interiors compared to the coastal areas. Another interesting feature is the high sedimentary Pa∕Th values in the Southern Ocean between 50 and 75 S (opal belt), where Pa is heavily scavenged to the sediments by opal. Our model generates higher sedimentary Pa∕Th in this region compared to the deep basin interior but our modelled opal belt stands slightly northward compared to the observations (Fig. 1). In our best fit simulation, the adsorption/desorption coefficients fall in the upper range of observations (see Tables 2, S1), therefore, Pa and Th are very effectively scavenged to the sediments. Besides, while Pa is mostly scavenged by the opal particles, Th mainly reaches the sediments with CaCO3 particles (see the Supplement for detailed information). Overall, these results are comparable with GCM outputs (van Hulten et al., 2018). The modelled values display the first order characteristics observed in the modern ocean and the sediment core tops. In the following section, we test the sensitivity of the simulated sedimentary Pa∕Th to abrupt circulation slowdown (hosing experiments).

3.2 Multiproxy response to an abrupt circulation slowdown

We added a freshwater flux of 0.3 Sv (sverdrups) in the Nordic Seas over a period of 300 years, which was sufficient to cause a drastic circulation reduction (Roche et al., 2010). Under PI conditions, the maximum of the AMOC stream function is about 15 Sv in our model. During the hosing, North Atlantic water formation drops to nearly 0 Sv, and the upper overturning cell completely vanishes (see Fig. S4 in the Supplement). The AMOC recovers ∼200 years after the end of the water hosing and displays a small overshoot with the maximum Atlantic meridional stream function exceeding 20 Sv around 900 years (Fig. 5).

Figure 2Theoretical definition of the proxy response studied (dotted red horizontal line) and the time of response (dotted vertical red line) in the case of (a) a single response or (b) a dual response to the circulation change induced by the freshwater addition. In the case of dual response, we examine the early (i.e. the first) and the late (i.e. the second) proxy response and time of response. The anomaly is defined as the difference between the proxy value at the defined time of response and the average proxy value (horizontal black line) during the control period of 300 years (highlighted in grey). The horizontal dashed black lines represent the proxy variance (2σ) over the control period (highlighted in grey). The proxy variable used is purposely not specified since the interest here is to provide an illustration of the definition used in the present study with the different modelled proxy variables.


We evaluate the response of Pa∕Th, δ13C and Δ14C to the hosing in the Nordic Seas as follows. We identify, for each model grid cell and each proxy, the simulation periods exceeding 80 years during which the proxy values are outside of their natural variability range. The latter is defined as the proxy variance (2σ) under PI conditions over the first 300 years of the simulation (control phase – see Fig. 2). In most cases, zero, one or two periods of significant response are detected. In some grid cells with high proxy variability, we detect up to four periods of significant proxy response, which are difficult to relate to either the hosing or overshoot timing of our simulation. Consequently, we excluded the grid cells depicting more than two periods of significant response from the subsequent analysis. For the time series containing two or less periods of interest, we call “time of maximum response”, the simulation year for which the absolute difference between the proxy value and mean proxy value during the control phase of the simulation is maximal. The proxy value at the time of maximum response is denoted “proxy response” in the following (Fig. 2). The single or dual proxy response is compared to the control proxy value (i.e. the mean proxy value over the first 300 years of the simulation). Figure 3 represents the zonal mean proxy response in the western Atlantic basin, the delimitation between the two basins being defined by the position of the Mid-Atlantic Ridge (see Figs. S5 and S6 for eastern basin).

Figure 3Zonally averaged anomalies for Pa∕Th, δ13C and Δ14C in the western Atlantic basin in the case of a freshwater input in the Nordic Seas. The western and eastern basin are delimited by the topography increase corresponding to the Mid-Atlantic Ridge in the model grid. The anomalies are computed/defined as the proxy response minus the mean of the proxy value during the control period of 300 years under PI conditions (see text). Panels (a)(c) represent the anomalies for the three proxies in the case where exactly one proxy response has been detected. In the case of two proxy responses, (d)(f) represent the proxy anomaly value for the early (first) response, while (g)(i) represent the proxy anomaly for the late (second) proxy response. Areas left blank were not showing a single response (A) or not showing exactly two responses (B). In each subplot, the grey contours represent the ocean bottom. Single and dual responses are mutually exclusive on a per location basis. Since (A) and (B) are showing zonal averages, overlaps may arise from different locations with the same latitude but different longitudes.


The δ13C and Δ14C only display one single response in the deep western Atlantic, whereas the Pa∕Th generally displays one or two responses in this part of the basin. Generally, the single or the first response of each proxy has the same geographical pattern (Fig. 3a–f), while in the case of two distinct responses, the late response has a radically different pattern (Fig. 3g–i). For δ13C and Δ14C, the late response corresponds to a general increase in the western Atlantic basin. For Pa∕Th, the late response pattern consists of increased values in the southern Atlantic and decreased values in the North Atlantic and is the opposite of the early response pattern. For the three proxies, the late response pattern is generally consistent with the circulation overshoot observed around 900 simulated years.

Considering that the single and the first responses mostly represent the proxy response to the hosing, a consistent proxy response is observed in the following three zones of the western Atlantic basin: the surface and intermediate waters (0–1500 m), the deep North Atlantic waters (> 2000 m) and the Southern Ocean (south of 30 S). In the surface waters and intermediate waters the proxy response pattern is as follows. The δ13C decreases in the first 500 m and then increases around 1000 m. The Δ14C increases from the surface to 1500 m, and Pa∕Th displays no clear trend. In the deep North Atlantic waters, the δ13C and Δ14C decrease, while the Pa∕Th generally increases, the three proxies reflecting a reduction in the North Atlantic Deep Water (NADW) flow rate. In terms of magnitude, we observe the strongest proxy response between 60 and 40 N and between 1500 and 3000 m for the three proxies. In the Southern Ocean, the Pa∕Th generally decreases, except in the deep southern basin (between 40 and 60 S, below 3000 m). The δ13C and Δ14C display a dipole pattern increasing at the surface and decreasing at depth. For δ13C, the depth boundary is around 1500 m, while for Δ14C the depth boundary is ∼1500 m north of 50 S and reaches 3000 m south of 40 S.

Figure 4Zonally averaged times of response (years) for Pa∕Th, δ13C and Δ14C in the western Atlantic basin in the case of a freshwater input in the Nordic Seas. The times of response (years) correspond to the time of proxy maximal response after the beginning of freshwater addition (see text). (a–c) Time of response in the case where exactly one single proxy response is detected. In the case where two distinct responses are detected (B). Panels (d)(f) show the time of response of the early (first) response, and (g)(i) show the time of response corresponding to the late (second) response. Areas left blank were not showing a single response (A) or not showing exactly two responses (B). In each subplot, the grey contours represent the ocean bottom.


Looking at the proxy time of response (Fig. 4), we observe significantly different patterns for Pa∕Th and the carbon isotopes. The δ13C and Δ14C time of response increases with depth as follows: in the surface and intermediate waters the response occurs roughly 300 years after the beginning of the fresh water addition, around 3000 m the response is delayed by 150 years, and towards the ocean bottom the delay increases up to 600 years. On the contrary, the Pa∕Th displays much smaller time of response, the timing of the Pa∕Th response being generally synchronous with the minimum of the stream function 300 years after the beginning of the freshwater addition. Consequently, in most of the western Atlantic basin, the response of the carbon isotopes lags behind the Pa∕Th response by a few hundred years, especially in the deeper waters.

In the eastern Atlantic basin the general pattern of the proxy response is similar to that of the western basin (Figs. S4–S5). However, we note that the δ13C has frequently more than a single response, especially at depths > 2000 m. The δ13C displays overall the same pattern as in the western basin except in the south deep basin (between ∼20 and 30  S – below 2500 m) in the case of a dual response. The Δ14C shows the same pattern as in the western basin, with increased values during the hosing in surface and intermediate waters and lower values at depth. Besides, the Pa∕Th has a more complex pattern with a large increase in the arctic basin, a moderate increase in the tropical basin, and a decrease in the northern and southern basins through the entire water column. Concerning the times of response, the same pattern as in the western ocean is observed. Interestingly, we note that the Δ14C has longer times of response in the eastern basin compared to the western basin in the case of one single response (Fig. S5). This is consistent with a lower ventilation of the eastern basin associated with a lower penetration of the NADW relative to the situation in the western boundary current.

Figure 5Selected multiproxy time series for the Nordic Seas hosing experiment representing the three main Atlantic water masses. (a) Zonally averaged on the western Atlantic basin time series at 40 N, 3660 m. This time series is representative of the proxy behaviours in the NADW. (b) Zonally averaged on the eastern Atlantic basin time series at 0 N, 1225 m. This time series is representative of the intermediate waters. (c) Zonally averaged on the western Atlantic basin time series at 30 S, 4385 m. This time series is representative of the proxy response in the Antarctic Bottom Water (AABW). In the eastern basin, the δ13C generally displays two responses, first increasing (not shown) and then decreasing (as shown). In each subplot from top to bottom, the dashed black line represent the freshwater flux (FWF) applied in the Nordic Seas, the black line represents the North Atlantic maximum stream function (a–c), the grey line represents the Southern Ocean maximum stream function (c), the red line represents the particulate Pa∕Th, the green line represents the δ13C and the blue line represents the Δ14C. The thin dashed black lines represent the proxy variance (2σ) in the first 300 simulated years. The blue vertical bands indicate the timing of the proxy responses.


Overall, the objective analysis of the multiproxy response allows the identification of three regions with different reaction patterns:

  • The clearest proxy response happens in the northwestern Atlantic between 40 and 60 N, at 1000 and 5000 m. In this area, the δ13C and Δ14C display marked decreases during the hosing, while the Pa∕Th significantly increases. Another interesting feature is that the δ13C and Δ14C lag behind the Pa∕Th by about 200 years (Fig. 5a). In our model, this zone corresponds to the region of NADW formation and first portion of southward flow at depth.

  • In the tropical intermediate and surface waters, δ13C and Δ14C increase during the hosing compared to the reference value, while Pa∕Th generally displays no clear reaction or very low-amplitude increase (Fig. 5b).

  • In the Southern Ocean (south of 30 S), below 3000 m, δ13C and Δ14C display a progressive decrease, while Pa∕Th slightly increases during the hosing. No lag between the Pa∕Th and the carbon isotopes response is observed (Fig. 5c). In our model, this zone corresponds to the Antarctic Bottom Water (AABW).

4 Discussion

The need for a more direct and efficient comparison between modelled proxies, easily related to modelled physical variables (e.g. stream function and currents) and observed palaeoproxy data motivated the implementation of several tracers in climate models. In this study we have fingerprinted the response of carbon isotopes (δ13C and Δ14C) and Pa∕Th to an abrupt circulation slowdown subsequent to the freshwater addition in the Nordic Seas. In short, our simulations show that the three proxies considered do not systematically respond simultaneously and in a simple manner to a given oceanic circulation change. Instead, the proxy responses show a strong spatial variability in the Atlantic basin, tightly related to the water mass mostly bathing the considered location. Besides, in a large part of the western Atlantic, the carbon isotopes response lags behind the Pa∕Th response by about 200 years, the actual time lag displaying marked spatial variability (Figs. 4 and 5). Our simulations thus indicate that a given circulation change event is not likely to produce a synchronous response, either for a given proxy in different regions of the Atlantic basin or for different proxies taken at the same location (i.e. sediment core), potentially complicating proxy data interpretation. In this section, we discuss the potential underlying physical mechanisms and compare our results with available model results and proxy data.

4.1 The lag of carbon isotopes, underlying physical mechanisms

One of the most interesting features of our simulations is that the δ13C and Δ14C response lags behind the Pa∕Th response by a few hundred years, in particular in the northwest Atlantic, bathed by the NADW, suggesting a shorter time of response for the Pa∕Th compared to carbon isotopes.

First, let us point out that the physical mechanisms leading to the δ13C, Δ14C and Pa∕Th signals are fundamentally different. On the one hand, Pa and Th have very short residence times in the water column. In our simulations, the particles fluxes remain prescribed and constant. Therefore, any change in the particle production due to the freshwater forcing is not accounted for and does not impact the sedimentary Pa∕Th in our simulations. Instead, the sedimentary Pa∕Th in the North Atlantic purely reflects the capacity of the NADW to transport Pa further south. Thus, the NADW cessation rapidly translates in a sedimentary Pa∕Th increase, which ceases as soon as the NADW resumes and is able to transport Pa again, at the end of the hosing period (i.e. after year 600). On the other hand, the oceanic carbon isotopes depend on not only circulation changes but also the biological activity and air–sea exchanges. Carbon is further actively partitioned between different reservoirs (terrestrial biosphere, ocean and atmosphere), hence the longer time needed by carbon isotopes to adjust to any change.

In our simulations, the freshwater perturbation lasts 300 years, and the NADW formation is stopped for about 100 years (Fig. 5), which is likely too short for the isotopic systems to fully adjust to the change. However, our experimental set-up is relevant to the study of centennial- to millennial-scale events (DO or Heinrich events) across which the climate and isotopic systems likely did not have the time to fully adjust either. In line with previous hosing experiments (e.g. Kageyama et al., 2013; Mariotti et al., 2012; Menviel et al., 2008, 2015), the freshwater addition in the Nordic Seas produces a surface cooling of about 1 C associated with an increase in the CO2 solubility in surface waters and a general marine productivity decrease in the North Atlantic. Subsequent to the NADW cessation, the nutrients accumulate in the Atlantic intermediate waters, shaping a negative δ13C anomaly and its counterpart positive phosphate anomaly in the North Atlantic around 50 N (Menviel et al., 2015). In our simulations, the biological productivity or air–sea exchanges response peak around year 600, when the AMOC is the weakest (Fig. S7). Therefore, the lag between the Pa∕Th and carbon isotopes responses is not related to late biological or air–sea exchange responses themselves but rather to the time needed to transport the anomaly to the deep ocean. Indeed, in our simulations, the northern sourced δ13C anomaly progressively spreads to the deep Atlantic basin and finally reaches the bottom ocean (> 4000 m) around year 800, about 200 years after the AMOC has reached its lower value (Fig. S8).

4.2 Comparison to proxy data and previous modelling efforts

4.2.1 Context overview

From the modelling side, the response of carbon isotopes and Pa∕Th to AMOC slowdown has so far been investigated in separate climate models. Unfortunately, experimental setups for the hosing strongly vary in terms of freshwater flux, duration and location from one publication to another, complicating quantitative comparison. Nevertheless, consistently with what is described in iLOVECLIM, the models generally simulate a marked Pa∕Th increase in the northwest Atlantic (Gu and Liu, 2017; Rempfer et al., 2017; Siddall et al., 2007) and a decrease in δ13C in response to the AMOC slowdown (e.g. Menviel et al., 2015; Schmittner and Lund, 2015).

From the data side, the best analogue for our hosing experiments applied to a PI state would be the 8.2 ka cal BP event (Alley et al., 1997), attributed to the drainage of the glacial Lake Agassiz (Hoffman et al., 2012; Wiersma and Renssen, 2006). However, this event is of short duration (∼300 years), and if some benthic δ13C data are available (e.g. Kleiven et al., 2008), to date there are no sufficiently well-resolved Δ14C and Pa∕Th records covering this event.

The Heinrich events and DO-cycles are good candidates to compare with our hosing experiments because they are strongly related to freshwater fluxes in the North Atlantic (e.g. Broecker et al., 1990; Hemming, 2004) and associated with circulation changes that have been documented in δ13C, Δ14C and Pa∕Th records (Böhm et al., 2015; Henry et al., 2016; Lynch-Stieglitz, 2017; Waelbroeck et al., 2018). However, the available palaeoproxy data are relatively sparse, and the extensive time series are rare. In addition, to date there is no published combined Pa∕Th, benthic δ13C and Δ14C record available for the same sediment core. Nevertheless, remarkably comprehensive multiproxy records are available for the Iberian margin (Gherardi et al., 2005; Skinner and Shackleton, 2004), the Brazilian margin (Burckel et al., 2015, 2016; Mulitza et al., 2017; Waelbroeck et al., 2018) and the Bermuda Rise (Böhm et al., 2015; Henry et al., 2016; Lippold et al., 2009; McManus et al., 2004). A common working hypothesis is to assume that these records represent the proxy evolution of the surrounding basin. Besides, some recent compilations (e.g. Lynch-Stieglitz et al., 2014; Ng et al., 2018; Zhao et al., 2018) bring some insight about the evolution of Pa∕Th, benthic δ13C and Δ14C across the last 40 kyr in the Atlantic Ocean.

However, our hosing experiments are not direct analogues of the millennial-scale climate changes in the last glacial cycle because (i) glacial millennial events occurred under glacial conditions, whereas our simulations were run under PI conditions; (ii) the Heinrich and DO events have distinct proxy patterns and cannot be entirely explained by a simple fresh water addition in the North Atlantic; and (iii) the freshwater inputs might have occurred in different locations across distinct millennial-scale events (e.g. originating from the Laurentide or Scandinavian ice sheet), while in our experiment, the freshwater was only added in the Nordic Seas. Furthermore, the sequence of mechanisms involved in Heinrich stadials is still under debate (Barker et al., 2015; Broecker, 1994), and these periods were likely subdivided in several distinct phases (Ng et al., 2018; Stanford et al., 2011).

Besides, the model set-up used in this study is rather crude. It uses fixed modelled particle fields and does not explicitly represent boundary or bottom scavenging by nepheloid layers. However, because the particle fields used show higher concentrations at the continental margins, a greater Pa and Th scavenging is achieved in the areas of observed boundary scavenging. Therefore, an additional parameterisation for boundary scavenging might not be necessary. However, the absence of representation of preferential Pa scavenging by resuspended particles at the bottom in nepheloid layers (e.g. Anderson et al., 1983; Gardner et al., 2018; Thomas et al., 2006) has not been parameterised in our model and could explain the relatively high sedimentary Pa∕Th values in the northwest Atlantic offshore Newfoundland and Florida. Given the reasons above, we will only discuss the common/divergent trends rather than precise values between the palaeodata and our model output.

4.2.2 General Atlantic features

The available palaeodata display a consistent proxy evolution in the western Atlantic and in particular in the western boundary current. Deep northwest Atlantic cores present the same Pa∕Th pattern over the last 25 kyr with a significant Pa∕Th increase during HS1 (Ng et al., 2018). A ∼0.4 ‰ benthic δ13C decrease is observed for some Heinrich events in deep and intermediate Atlantic waters (Lynch-Stieglitz et al., 2014). Finally, several coral and foraminifer 14C records also depict a decreased 14C (decreased Δ14C) content of deep waters and increased deep ventilation ages during HS1 (Chen et al., 2015; Skinner et al., 2014). While no major east–west difference has been highlighted concerning the carbon isotopes, Pa∕Th response seems less clear in the eastern Atlantic basin compared to the western boundary current. While some eastern Atlantic cores located close to the Mid-Atlantic Ridge and the Iberian margin reproduce the amplitude of Pa∕Th variations observed in the western basin, two other eastern records display no millennial-scale variability (Ng et al., 2018).

In line with the palaeoproxy data, our model results show a coherent and significant δ13C and Δ14C decrease in the deep western Atlantic. Overall the same response pattern is obtained in our model in the eastern Atlantic basin with the exception of the deep southeast Atlantic (between ∼20 and 30  S – below 2500 m). In this region, we observe a slight increase in the δ13C during the early response in the case of a dual response, while the Δ14C decreases during the single or early response. We note that the amplitude of this early δ13C increase is rather small compared to the late δ13C decrease. Moreover, in this region, the time of response is rather short (∼100–150 years). We argue here that in this particular region the δ13C displays a two-phase response, the actual long-term response corresponding to the late response.

In the surface and intermediate waters, the modelled δ13C decreases above 500 m and then increases around 1000 m, while Δ14C generally increases from the surface to 1500 m (Fig. 3). This δ13C pattern is consistent with changes in productivity reported in previous hosing experiments. Indeed, due to changes in winds and upwelling intensity and subsequent changes in nutrient availability, the productivity generally decreases in the North and equatorial Atlantic (Menviel et al., 2008). Therefore, the δ13C of the photic zone decreases, and the δ13C of the remineralisation zone (∼1000 m) increases as the export from 13C-depleted carbon is reduced during the hosing compared to the PI equilibrium state. As reported in other model studies, the increased Δ14C in the upper ocean reflects the accumulation of radiocarbon in the atmosphere and the upper ocean in the absence of carbon entrainment in deep water through North Atlantic deep-water formation (Butzin et al., 2005). The generally good ventilation state of the upper ocean indicated by both δ13C and Δ14C may also be related to an increase in the vertical mixing (and ventilation) in the intermediate waters in the absence of North Atlantic deep convection cell (Schmittner, 2007). The pattern we obtain in the intermediate waters (Figs. 3 and 5b) is consistent with a δ13C dataset on the Brazilian margin (Lund et al., 2015; Tessin and Lund, 2013) showing increasing δ13C between 1000 and 1300 m and decreasing between 1600 and 2100 m water depth. These results are also fairly consistent with the modelled δ13C described in Schmittner and Lund (2015) though our hosing experiments are notably shorter (300 vs 2000 years).

In agreement with the Pa∕Th data, the model simulates a clear Pa∕Th increase in the NADW in the western basin, while the results are more ambiguous in the eastern basin. Indeed, in the eastern basin, except at high northern latitudes, very small amplitude Pa∕Th variations are observed, with increasing values in the equatorial zone and decreasing values elsewhere. Thus, our modelled Pa∕Th response at the Iberian margin is very different from the Iberian margin record. The low amplitude of our Pa∕Th response may arise from the absence of strong Pa advection towards the southern Atlantic in relation with the east equatorial Atlantic upwellings variations across the hosing experiment (Menviel et al., 2008). Alternatively, the absence of Pa∕Th signal in the eastern basin could be due to an overestimation of the water fluxes that cross Gibraltar in the model. Indeed, the model resolution is rather coarse, and the Gibraltar Straight is represented by a full ocean grid cell of 3×3, unrealistically impacting the simulated eastern Atlantic circulation. Additionally, the modelled Pa∕Th displays interesting features in the southern Atlantic and in the surface and intermediate waters (between 40 S and 40 N). In both cases, no reliable Pa∕Th record is available because particles fluxes from the opal belt (southern Atlantic) and coastal areas (surface waters) overprint the Pa∕Th circulation signal. In our model, we observe that the single and the early responses do not match and that the response of the Pa∕Th is rather fast (∼100–150 years) in the case of a dual response. This suggests that in this area, the Pa∕Th displays a two-phase response to the freshwater addition and AMOC perturbation.

Further investigations into this two-phase response for Pa∕Th and δ13C, in particular in the Southern Ocean, would require (1) monitoring of the detailed 3D circulation pathway and carbon exchanges between the different reservoirs and (2) running experiments with increased hosing duration in order to better assess the proxy response to a sustained (∼1500 years) AMOC shutdown.

4.2.3 Bermuda Rise and Brazilian margin multiproxy records

The Brazilian margin and the Bermuda Rise are located in the western boundary current and display similar proxy patterns. The Pa∕Th increases from ∼0.06 to the production ratio (0.093) or even above while the δ13C decreases by ∼0.5 ‰ across the Heinrich events. Similar Pa∕Th and benthic δ13C changes, of smaller/reduced amplitude, are also observed for the DO stadials.

Figure 6Example of the Bermuda Rise time series. (a) Pa∕Th and benthic δ13C (‰) measured at the Bermuda Rise around HS4 (Henry et al., 2016). (b) Modelled time series corresponding to the average of the nine grid cells surrounding the model grid cell closest to the Bermuda Rise (34 N, 58 W). (a) Imposed freshwater flux, (b) Maximum North Atlantic meridional stream function (Sv), (c) simulated Pa∕Th and (d) benthic δ13C (‰).


Below, we examine our modelled time series in the Bermuda Rise basin (34 N–58 W, > 4300 m) (Fig. 6). The presented time series correspond to the average of nine model grid cells and are representative of the time series of the deep western basin (see Figs. 5 and 6). The simulated Pa∕Th significantly increases between year 350 and year 850 of the simulation, consistently with the decrease in the maximum Atlantic stream function (Fig. 6). The simulated Pa∕Th approaches the production ratio of 0.093 at year 600, and the maximum stream function is close to zero from year 400 to 600, which is consistent with a sedimentary Pa∕Th reaching the production ratio in the case of an AMOC shutdown. The simulated Pa∕Th change has a moderate amplitude of 0.015 Pa∕Th units (from ∼0.075 to 0.090) but is significant with respect to the natural variability recorded during the first 300 years of the simulation under PI conditions. Our simulated δ13C decreases from ∼0.62 ‰ to ∼0.5 ‰, reaching the minimum around year 900, i.e. 600 years after the beginning of the hosing. The simulated δ13C decrease has a moderate amplitude of ∼0.12 ‰ but is significant with respect to the natural variability recorded in the control first 300 years of the simulation. Although the trend of the simulated and observed proxy response is the same, their absolute values differ. In our simulated record, the Pa∕Th value associated with the modern circulation scheme is around 0.075, which is significantly higher than the value actually measured at the Bermuda Rise or Brazilian margin sites: ∼0.06. For both proxies, the simulated amplitude of change is much smaller than the amplitude recorded in the palaeodata across the Heinrich events; the modelled δ13C decrease is around 0.12 ‰, while it is 0.5 ‰ in the palaeodata; the modelled Pa∕Th change is ∼0.015, while it is ∼0.03 in the palaeodata. This might be a consequence of the short duration of the fresh water forcing and induced reduction in AMOC (for carbon isotopes) and of poor particle representation along the northeast American coast (for Pa∕Th).

The lead–lag relationship between Pa∕Th and benthic δ13C was previously examined on data both at the Brazilian margin and the Bermuda Rise with opposite conclusions: at the Brazilian margin the Pa∕Th was found to lead δ13C by about 200 years (Waelbroeck et al., 2018), while at the Bermuda Rise, the Pa∕Th was found to lag behind δ13C by about 200 years (Henry et al., 2016). However, the latter result must be interpreted with caution since the cross-correlation method used to evaluate the lead or lag relationships in (Henry et al., 2016) has been designed for signals that are stationary in time and is thus is not suitable to analyse nonstationary climatic signals (Waelbroeck et al., 2018). Therefore, the palaeodata seem to suggest a lead of the Pa∕Th response with respect to the carbon isotopes response in the western boundary current. However, this observation may simply reflect the impact of bioturbation on sediment archives. Indeed, even if the two proxies are recorded in a single sediment core, it is important to note that both proxies are not hosted by the same sediment fraction; the Pa and Th being preferentially adsorbed on the fine grain particles (< 100 µm, Chase et al., 2002; Kretschmer et al., 2010; Thomson et al., 1993), while δ13C and Δ14C are measured on foraminifer shells that correspond on larger particle sizes (> 150 µm). It has been shown that bioturbation could affect different particle sizes differently (e.g. Wheatcroft, 1992). Therefore Pa∕Th and carbon isotopes could be affected by bioturbation in a different way, and the 200-year Pa∕Th response lead on carbon isotopes observed in Brazilian margin sediments could then solely be explained by sediment bioturbation, as suggested in (Waelbroeck et al., 2018). In contrast, in our model we observe a systematic lead of the Pa∕Th response with respect to the carbon isotopes in the NADW. Therefore, our model suggests that this Pa∕Th lead may be a feature of the proxy response to millennial-scale variability and is not necessarily an artefact due to the marine core bioturbation.

5 Conclusions and perspectives

We have implemented the 231Pa and 230Th tracers in the climate model of intermediate complexity iLOVECLIM. The new Pa∕Th module simulates dissolved and particulate 231Pa and 230Th concentrations with a quality comparable to GCMs under PI equilibrium.

To date, the model is able to simulate the evolution of δ13C, Δ14C and Pa∕Th over thousands of years in a consistent physical framework and in a reasonable computation time (∼800 years per 24 h). We tested and fingerprinted the response of these three proxies to an imposed and abrupt circulation change by adding freshwater in the Nordic Seas. The proxy response displays high spatial and temporal variability, but three patterns of proxy response, tightly related to the main Atlantic water masses, can be identified. In intermediate waters, the Pa∕Th and Δ14 show a slight increase, while the δ13C significantly increases. In AABW, the three proxies display consistent and synchronous responses: the Pa∕Th increases, while the δ13C and Δ14C decrease, corresponding to a reduced circulation and ventilation of the water mass. Finally, in NADW, the Pa∕Th increases, while the δ13C and Δ14C decrease, with a lag of about 200 years between the Pa∕Th and carbon isotopes responses. Our simulations show that this lag is induced by the different mechanisms producing the proxy records. Indeed, since our model set-up uses prescribed particles fields, the sedimentary Pa∕Th in the North Atlantic basin is purely driven by the capacity of the ocean circulation to transport Pa further south. On the contrary, the carbon isotopes- are impacted by multiple processes, such as biological productivity, remineralisation and air–sea exchanges, and the anomaly formed in the intermediate North Atlantic takes time to be exported toward the deep basin.

Despite the crudeness of the model set-up and incomplete representation of the processes governing Pa and Th scavenging in the water column, the main features of our modelled Pa∕Th and carbon isotopes responses are consistent with the available palaeoproxy record. We observe (i) coherent proxy response along the Atlantic western boundary current for the three proxies consisting in a significant Pa∕Th increase and δ13C and Δ14C decrease; (ii) distinct proxy responses for intermediate and deep waters for δ13C and Δ14C; and (iii) a constitutive lag between the Pa∕Th and carbon isotopes, consistent with a multiproxy record from the Brazilian margin. Our study therefore suggests that a given circulation change event does not likely produce a synchronous response, either for a given proxy in different regions of the Atlantic basin or for different proxies taken at the same location (i.e. sediment core), potentially complicating proxy data interpretation. Future work is required to evaluate the multiproxy response in more realistic numerical experiments, using glacial boundary conditions and coupled/interactive particle fields. In addition, a more complete dataset containing multiproxy records is needed to achieve a more complete model–data comparison.

Code and data availability

The iLOVECLIM source code is based on the LOVECLIM model version 1.2 whose code is accessible at (last access: 17 April 2019, UCLouvain, 2019). The developments in the iLOVECLIM source code are hosted at (last access: 17 April 2020, LUDUS, 2020) but are not publicly available due to copyright restrictions. Access can be granted on demand by request to Didier M. Roche ( The model output related to this article is available on Pangaea (, Missiaen, 2019).


The supplement related to this article is available online at:

Author contributions

LM, CW and DMR designed the research. LM, DMR, NB, JCD and AQ developed the iLOVECLIM model. LM performed the simulations. LM and JYP developed the postprocessing algorithm. JCD and SP contributed with expert knowledge of Pa∕Th. LM wrote the paper with the inputs from all the coauthors.

Competing interests

The authors declare that they have no conflict of interest.


This is a contribution to ERC project ACCLIMATE; the research leading to these results has received funding from the European Research Council under the European Union Seventh Framework Programme (FP7/2007-2013), ERC grant agreement 339108. Lise Missiaen acknowledges funding from the Australian Research Council grant DP180100048. We thank Santiago Moreira and Fanny Lhardy for their help with Python. Elisabeth Michel is thanked for expert knowledge discussion of 14C. This is a LSCE contribution 6949.

Financial support

This research has been supported by the European Research Council (ACCLIMATE (grant no. 339108)) and the Australian Research Council (grant no. DP180100048).

Review statement

This paper was edited by Ed Brook and reviewed by two anonymous referees.


Adkins, J. F. and Boyle, E. A.: Changing atmospheric Δ14C and the record of deep water paleoventilation ages, Paleoceanography, 12, 337–344,, 1997. 

Alley, R. B., Sowers, T., Mayewski, P. A., Stuiver, M., Taylor, K. C., and Clark, P. U.: Holocene climatic instability: A prominent, widespread event 8200 yr ago, Geology, 25, 483–486,<0483:HCIAPW>2.3.CO;2, 1997. 

Anderson, R. F., Bacon, M. P., and Brewer, P. G.: Removal of 230Th and 231Pa at ocean margins, Earth Planet. Sc. Lett., 66, 73–90,, 1983. 

Bard, E., Arnold, M., Mangerud, J., Paterne, M., Labeyrie, L., Duprat, J., Mélières, M.-A., Sønstegaard, E., and Duplessy, J.-C.: The North Atlantic atmosphere-sea surface 14C gradient during the Younger Dryas climatic event, Earth Planet. Sc. Lett., 126, 275–287,, 1994. 

Barker, S., Chen, J., Gong, X., Jonkers, L., Knorr, G., and Thornalley, D.: Icebergs not the trigger for North Atlantic cold events, Nature, 520, 333–336, 2015. 

Böhm, E., Lippold, J., Gutjahr, M., Frank, M., Blaser, P., Antz, B., Fohlmeister, J., Frank, N., Andersen, M. B., and Deininger, M.: Strong and deep Atlantic meridional overturning circulation during the last glacial cycle, Nature, 517, 73–76, 2015. 

Bondevik, S., Mangerud, J., Birks, H. H., Gulliksen, S., and Reimer, P.: Changes in North Atlantic Radiocarbon Reservoir Ages During the Allerød and Younger Dryas, Science, 312, 1514,, 2006. 

Bouttes, N., Roche, D. M., Mariotti, V., and Bopp, L.: Including an ocean carbon cycle model into iLOVECLIM (v1.0), Geosci. Model Dev., 8, 1563–1576,, 2015. 

Broecker, W. S.: Massive iceberg discharges as triggers for global climate change, Nature, 372, 421–424,, 1994. 

Broecker, W. S., Bond, G., Klas, M., Bonani, G., and Wolfli, W.: A salt oscillator in the glacial Atlantic? 1. The concept, Paleoceanography, 5, 469–477,, 1990. 

Brovkin, V., Ganopolski, A., and Svirezhev, Y.: A continuous climate-vegetation classification for use in climate-biosphere studies, Ecol. Model., 101, 251–261,, 1997. 

Brovkin, V., Bendtsen, J., Claussen, M., Ganopolski, A., Kubatzki, C., Petoukhov, V., and Andreev, A.: Carbon cycle, vegetation, and climate dynamics in the Holocene: Experiments with the CLIMBER-2 model, Global Biogeochem. Cy., 16, 1139,, 2002. 

Burckel, P., Waelbroeck, C., Gherardi, J. M., Pichat, S., Arz, H., Lippold, J., Dokken, T., and Thil, F.: Atlantic Ocean circulation changes preceded millennial tropical South America rainfall events during the last glacial, Geophys. Res. Lett., 42, 411–418,, 2015. 

Burckel, P., Waelbroeck, C., Luo, Y., Roche, D. M., Pichat, S., Jaccard, S. L., Gherardi, J., Govin, A., Lippold, J., and Thil, F.: Changes in the geometry and strength of the Atlantic meridional overturning circulation during the last glacial (20–50 ka), Clim. Past, 12, 2061–2075,, 2016. 

Butzin, M., Prange, M., and Lohmann, G.: Radiocarbon simulations for the glacial ocean: The effects of wind stress, Southern Ocean sea ice and Heinrich events, Earth Planet. Sc. Lett., 235, 45–61,, 2005. 

Chase, Z., Anderson, R. F., Fleisher, M. Q., and Kubik, P. W.: The influence of particle composition and particle flux on scavenging of Th, Pa and Be in the ocean, Earth Planet. Sc. Lett., 204, 215–229,, 2002. 

Chase, Z., Anderson, R. F., Fleisher, M. Q., and Kubik, P. W.: Scavenging of 230Th, 231Pa and 10Be in the Southern Ocean (SW Pacific sector): the importance of particle flux, particle composition and advection, Deep-Sea Res. Pt. II, 50, 739–768,, 2003. 

Chen, T., Robinson, L. F., Burke, A., Southon, J., Spooner, P., Morris, P. J., and Ng, H. C.: Synchronous centennial abrupt events in the ocean and atmosphere during the last deglaciation, Science, 349, 1537,, 2015. 

Duplessy, J. C., Shackleton, N. J., Fairbanks, R. G., Labeyrie, L., Oppo, D., and Kallel, N.: Deepwater source variations during the last climatic cycle and their impact on the global deepwater circulation, Paleoceanography, 3, 343–360,, 1988. 

Dutay, J.-C., Lacan, F., Roy-Barman, M., and Bopp, L.: Influence of particle size and type on 231Pa and 230Th simulation with a global coupled biogeochemical-ocean general circulation model: A first approach, Geochem. Geophy. Geosy., 10, Q01011,, 2009. 

Eide, M., Olsen, A., Ninnemann, U. S., and Johannessen, T.: A global ocean climatology of preindustrial and modern ocean δ13C, Global Biogeochem. Cy., 31, 515–534,, 2017. 

François, R.: Paleoflux and paleocirculation from sediment 230Th and 231Pa∕230Th. Proxies in Late Cenozoic Paleoceanography, Elsevier, Amsterdam, 2007. 

Friedlingstein, P., Cox, P., Betts, R., Bopp, L., von Bloh, W., Brovkin, V., Cadule, P., Doney, S., Eby, M., and Fung, I.: Climate–carbon cycle feedback analysis: results from the C4MIP model intercomparison, J. Climate, 19, 3337–3353, 2006. 

Gardner, W. D., Richardson, M. J., and Mishonov, A. V.: Global assessment of benthic nepheloid layers and linkage with upper ocean dynamics, Earth Planet. Sc. Lett., 482, 126–134, 2018. 

Gherardi, J.-M., Labeyrie, L., McManus, J. F., Francois, R., Skinner, L. C., and Cortijo, E.: Evidence from the Northeastern Atlantic basin for variability in the rate of the meridional overturning circulation through the last deglaciation, Earth Planet. Sc. Lett., 240, 710–723,, 2005. 

Godwin, H.: Half-life of radiocarbon, Nature, 195, 984–984, 1962. 

Goosse, H., Brovkin, V., Fichefet, T., Haarsma, R., Huybrechts, P., Jongma, J., Mouchet, A., Selten, F., Barriat, P.-Y., Campin, J.-M., Deleersnijder, E., Driesschaert, E., Goelzer, H., Janssens, I., Loutre, M.-F., Morales Maqueda, M. A., Opsteegh, T., Mathieu, P.-P., Munhoven, G., Pettersson, E. J., Renssen, H., Roche, D. M., Schaeffer, M., Tartinville, B., Timmermann, A., and Weber, S. L.: Description of the Earth system model of intermediate complexity LOVECLIM version 1.2, Geosci. Model Dev., 3, 603–633,, 2010. 

Gu, S. and Liu, Z.: 231Pa and 230Th in the ocean model of the Community Earth System Model (CESM1.3), Geosci. Model Dev., 10, 4723–4742,, 2017. 

Hayes, C. T., Anderson, R. F., Fleisher, M. Q., Vivancos, S. M., Lam, P. J., Ohnemus, D. C., Huang, K.-F., Robinson, L. F., Lu, Y., Cheng, H., Edwards, R. L., and Moran, S. B.: Intensity of Th and Pa scavenging partitioned by particle chemistry in the North Atlantic Ocean, Mar. Chem., 170, 49–60,, 2015. 

Hemming, S. R.: Heinrich events: Massive late Pleistocene detritus layers of the North Atlantic and their global climate imprint, Rev. Geophys., 42, RG1005,, 2004. 

Henderson, G. M. and Anderson, R. F.: The U-series toolbox for paleoceanography, Rev. Mineral. Geochem., 52, 493–531, 2003. 

Henry, L. G., McManus, J. F., Curry, W. B., Roberts, N. L., Piotrowski, A. M., and Keigwin, L. D.: North Atlantic ocean circulation and abrupt climate change during the last glaciation, Science, 353, 470–474, 2016. 

Hoffman, J. S., Carlson, A. E., Winsor, K., Klinkhammer, G. P., LeGrande, A. N., Andrews, J. T., and Strasser, J. C.: Linking the 8.2 ka event and its freshwater forcing in the Labrador Sea, Geophys. Res. Lett., 39, L18703,, 2012. 

Kageyama, M., Merkel, U., Otto-Bliesner, B., Prange, M., Abe-Ouchi, A., Lohmann, G., Ohgaito, R., Roche, D. M., Singarayer, J., Swingedouw, D., and X Zhang: Climatic impacts of fresh water hosing under Last Glacial Maximum conditions: a multi-model study, Clim. Past, 9, 935–953,, 2013. 

Kindler, P., Guillevic, M., Baumgartner, M., Schwander, J., Landais, A., and Leuenberger, M.: Temperature reconstruction from 10 to 120 kyr b2k from the NGRIP ice core, Clim. Past, 10, 887–902,, 2014. 

Kleiven, H. F., Kissel, C., Laj, C., Ninnemann, U. S., Richter, T. O., and Cortijo, E.: Reduced North Atlantic Deep Water Coeval with the Glacial Lake Agassiz Freshwater Outburst, Science, 319, 60–64,, 2008. 

Kretschmer, S., Geibert, W., Rutgers van der Loeff, M. M., and Mollenhauer, G.: Grain size effects on 230Thxs inventories in opal-rich and carbonate-rich marine sediments, Earth Planet. Sc. Lett., 294, 131–142,, 2010. 

Lippold, J., Grützner, J., Winter, D., Lahaye, Y., Mangini, A., and Christl, M.: Does sedimentary 231Pa/230Th from the Bermuda Rise monitor past Atlantic Meridional Overturning Circulation?, Geophys. Res. Lett., 36, L12601,, 2009. 

LUDUS:, last access: 17 April 2020. 

Lund, D. C., Tessin, A. C., Hoffman, J. L., and Schmittner, A.: Southwest Atlantic water mass evolution during the last deglaciation, Paleoceanography, 30, 477–494,, 2015. 

Luo, Y., Francois, R., and Allen, S. E.: Sediment 231Pa∕230Th as a recorder of the rate of the Atlantic meridional overturning circulation: insights from a 2-D model, Ocean Sci., 6, 381–400,, 2010. 

Lynch-Stieglitz, J.: The Atlantic Meridional Overturning Circulation and Abrupt Climate Change, Annu. Rev. Mar. Sci., 9, 83–104,, 2017. 

Lynch-Stieglitz, J., Schmidt, M. W., Gene Henry, L., Curry, W. B., Skinner, L. C., Mulitza, S., Zhang, R., and Chang, P.: Muted change in Atlantic overturning circulation over some glacial-aged Heinrich events, Nat. Geosci., 7, 144–150, 2014. 

Marchal, O., François, R., Stocker, T. F., and Joos, F.: Ocean thermohaline circulation and sedimentary 231Pa∕230Th ratio, Paleoceanography, 15, 625–641,, 2000. 

Mariotti, V., Bopp, L., Tagliabue, A., Kageyama, M., and Swingedouw, D.: Marine productivity response to Heinrich events: a model-data comparison, Clim. Past, 8, 1581–1598,, 2012. 

McManus, J. F., Francois, R., Gherardi, J.-M., Keigwin, L. D., and Brown-Leger, S.: Collapse and rapid resumption of Atlantic meridional circulation linked to deglacial climate changes, Nature, 428, 834–837,, 2004. 

Menviel, L., Timmermann, A., Mouchet, A., and Timm, O.: Meridional reorganizations of marine and terrestrial productivity during Heinrich events, Paleoceanography, 23, PA1203,, 2008. 

Menviel, L., Mouchet, A., Meissner, K. J., Joos, F., and England, M. H.: Impact of oceanic circulation changes on atmospheric δ13CO2, Global Biogeochem. Cy., 29, 1944–1961,, 2015. 

Menviel, L., Yu, J., Joos, F., Mouchet, A., Meissner, K. J., and England, M. H.: Poorly ventilated deep ocean at the Last Glacial Maximum inferred from carbon isotopes: A data-model comparison study, Paleoceanography, 32, 2–17,, 2017. 

Missiaen, L., Pichat, S., Waelbroeck, C., Douville, E., Bordier, L., Dapoigny, A., Thil, F., Foliot, L., and Wacker, L.: Downcore Variations of Sedimentary Detrital (238U∕232Th) Ratio: Implications on the Use of 230Thxs and 231Paxs to Reconstruct Sediment Flux and Ocean Circulation, Geochem. Geophy. Geosy., 19, 2560–2573,, 2018. 

Missiaen, L.: Output from iLOVECLIM Pa∕Th module: Pre-industrial equilibrium and Nordic Seas hosing experiment, PANGAEA,, 2019. 

Mulitza, S., Chiessi, C. M., Schefuß, E., Lippold, J., Wichmann, D., Antz, B., Mackensen, A., Paul, A., Prange, M., Rehfeld, K., Werner, M., Bickert, T., Frank, N., Kuhnert, H., Lynch-Stieglitz, J., Portilho-Ramos, R. C., Sawakuchi, A. O., Schulz, M., Schwenk, T., Tiedemann, R., Vahlenkamp, M., and Zhang, Y.: Synchronous and proportional deglacial changes in Atlantic meridional overturning and northeast Brazilian precipitation, Paleoceanography, 32, 622–633,, 2017. 

Ng, H. C., Robinson, L. F., McManus, J. F., Mohamed, K. J., Jacobel, A. W., Ivanovic, R. F., Gregoire, L. J., and Chen, T.: Coherent deglacial changes in western Atlantic Ocean circulation, Nat. Commun., 9, 2947,, 2018. 

Orr, J. C., Maier-Reimer, E., Mikolajewicz, U., Monfray, P., Sarmiento, J. L., Toggweiler, J. R., Taylor, N. K., Palmer, J., Gruber, N., Sabine, C. L., Le Quéré, C., Key, R. M., and Boutin, J.: Estimates of anthropogenic carbon uptake from four three-dimensional global ocean models, Global Biogeochem. Cy., 15, 43–60,, 2001. 

Rempfer, J., Stocker, T. F., Joos, F., Lippold, J., and Jaccard, S. L.: New insights into cycling of 231Pa and 230Th in the Atlantic Ocean, Earth Planet. Sc. Lett., 468, 27–37,, 2017. 

Roche, D. M., Wiersma, A. P., and Renssen, H.: A systematic study of the impact of freshwater pulses with respect to different geographical locations, Clim. Dynam., 34, 997–1013,, 2010. 

Schlitzer, R., Anderson, R. F., Dodas, E. M., Lohan, M., Geibert, W., Tagliabue, A., Bowie, A., Jeandel, C., Maldonado, M. T., and Landing, W. M.: The GEOTRACES intermediate data product 2017, Chem. Geol., 493, 210–223, 2018. 

Schmittner, A. and Lund, D. C.: Early deglacial Atlantic overturning decline and its role in atmospheric CO2 rise inferred from carbon isotopes (δ13C), Clim. Past, 11, 135–152,, 2015. 

Schmittner, A. B.: Impact of the Ocean's Overturning Circulation on Atmospheric CO2, American Geophysical Union, 2007. 

Siddall, M., Henderson, G. M., Edwards, N. R., Frank, M., Müller, S. A., Stocker, T. F., and Joos, F.: 231Pa∕230Th fractionation by ocean transport, biogenic particle flux and particle type, Earth Planet. Sc. Lett., 237, 135–155,, 2005. 

Siddall, M., Stocker, T. F., Henderson, G. M., Joos, F., Frank, M., Edwards, N. R., Ritz, S. P., and Müller, S. A.: Modeling the relationship between 231Pa∕230Th distribution in North Atlantic sediment and Atlantic meridional overturning circulation, Paleoceanography, 22, PA2214,, 2007. 

Siegenthaler, U. and Münnich, K. O.: C∕12C fractionation during CO2 transfer from air to sea, Carbon Cycle Modelling. Bolin B(ed) Wiley, New York, 249–257, 1981. 

Skinner, L. C. and Shackleton, N. J.: Rapid transient changes in northeast Atlantic deep water ventilation age across Termination I, Paleoceanography, 19, PA2005,, 2004. 

Skinner, L. C., Waelbroeck, C., Scrivner, A. E., and Fallon, S. J.: Radiocarbon evidence for alternating northern and southern sources of ventilation of the deep Atlantic carbon pool during the last deglaciation, P. Natl. Acad. Sci. USA, 111, 5480,, 2014. 

Smeed, D. A., McCarthy, G. D., Cunningham, S. A., Frajka-Williams, E., Rayner, D., Johns, W. E., Meinen, C. S., Baringer, M. O., Moat, B. I., Duchez, A., and Bryden, H. L.: Observed decline of the Atlantic meridional overturning circulation 2004–2012, Ocean Sci., 10, 29–38,, 2014. 

Stanford, J. D., Rohling, E. J., Bacon, S., Roberts, A. P., Grousset, F. E., and Bolshaw, M.: A new concept for the paleoceanographic evolution of Heinrich event 1 in the North Atlantic, Quaternary Sci. Rev., 30, 1047–1066,, 2011. 

Tessin, A. C. and Lund, D. C.: Isotopically depleted carbon in the mid-depth South Atlantic during the last deglaciation, Paleoceanography, 28, 296–306,, 2013. 

Thomas, A. L., Henderson, G. M., and Robinson, L. F.: Interpretation of the 231Pa/230Th paleocirculation proxy: New water-column measurements from the southwest Indian Ocean, Earth Planet. Sc. Lett., 241, 493–504,, 2006. 

Thomson, J., Colley, S., Anderson, R., Cook, G. T., MacKenzie, A. B., and Harkness, D. D.: Holocene sediment fluxes in the northeast Atlantic from 230Thexcess and radiocarbon measurements, Paleoceanography, 8, 631–650,, 1993. 

Thornalley, D. J. R., Barker, S., Broecker, W. S., Elderfield, H., and McCave, I. N.: The Deglacial Evolution of North Atlantic Deep Convection, Science, 331, 202–205,, 2011. 

Thornalley, D. J. R., Bauch, H. A., Gebbie, G., Guo, W., Ziegler, M., Bernasconi, S. M., Barker, S., Skinner, L. C., and Yu, J.: A warm and poorly ventilated deep Arctic Mediterranean during the last glacial period, Science, 349, 706–710,, 2015. 

Tschumi, T., Joos, F., Gehlen, M., and Heinze, C.: Deep ocean ventilation, carbon isotopes, marine sedimentation and the deglacial CO2 rise, Clim. Past, 7, 771–800,, 2011. 

UCLouvain: LOVECLIM coupled model, version 1.2, available at:, last access: 17 April 2019. 

van Hulten, M., Dutay, J.-C., and Roy-Barman, M.: A global scavenging and circulation ocean model of thorium-230 and protactinium-231 with improved particle dynamics (NEMO–ProThorP 0.1), Geosci. Model Dev., 11, 3537–3556,, 2018.  

Waelbroeck, C., Duplessy, J.-C., Michel, E., Labeyrie, L., Paillard, D., and Duprat, J.: The timing of the last deglaciation in North Atlantic climate records, Nature, 412, 724–727,, 2001. 

Waelbroeck, C., Pichat, S., Böhm, E., Lougheed, B. C., Faranda, D., Vrac, M., Missiaen, L., Vazquez Riveiros, N., Burckel, P., Lippold, J., Arz, H. W., Dokken, T., Thil, F., and Dapoigny, A.: Relative timing of precipitation and ocean circulation changes in the western equatorial Atlantic over the last 45 kyr, Clim. Past, 14, 1315–1330,, 2018. 

Waelbroeck, C., Lougheed, B. C., Riveiros, N. V., et al.: Consistently dated Atlantic sediment cores over the last 40 thousand years, Sci. Data, 6, 165, 2019. 

Wheatcroft, R. A.: Experimental tests for particle size-dependent bioturbation in the deep ocean, Limnol. Oceanogr., 37, 90–104,, 1992. 

Wiersma, A. P. and Renssen, H.: Model–data comparison for the 8.2 ka BP event: confirmation of a forcing mechanism by catastrophic drainage of Laurentide Lakes, Quaternary Sci. Rev., 25, 63–88,, 2006. 

Yu, E.-F., Francois, R. and Bacon, M. P.: Similar rates of modern and last-glacial ocean thermohaline circulation inferred from radiochemical data, Nature, 379, 689–694,, 1996. 

Zhao, N., Marchal, O., Keigwin, L., Amrhein, D., and Gebbie, G.: A Synthesis of Deglacial Deep-Sea Radiocarbon Records and Their (In)Consistency With Modern Ocean Ventilation, Paleoceanogr. Paleocl., 33, 128–151,, 2018.