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

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 δ¹³C and Δ¹⁴C 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 ²³¹Pa and ²³⁰Th tracers into the model of intermediate complexity iLOVECLIM, which already included δ¹³C and Δ¹⁴C. 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 δ¹³C and Δ¹⁴C 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 (δ¹³C and Δ¹⁴C) 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 δ¹⁸O, Nd isotopes, Cd∕Ca or sortable silts), the sedimentary (${}^{\mathrm{231}}{\mathrm{Pa}}_{\mathrm{xs},\mathrm{0}}{/}^{\mathrm{230}}{\mathrm{Th}}_{\mathrm{xs},\mathrm{0}}$) ratio (hereafter Pa∕Th) and carbon isotopes (δ¹³C, Δ¹⁴C) are key tracers used to reconstruct and quantify past circulation patterns and water mass flow rates.

Pa∕Th is the ratio of the ²³¹Pa and ²³⁰Th 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, ²³¹Pa and ²³⁰Th 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 ²³¹Pa and 10–40 years for ²³⁰Th; Henderson and Anderson, 2003). Thus, while ²³⁰Th is rapidly transferred to the sediment, ²³¹Pa 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 ²³¹Pa and ²³⁰Th 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 δ¹³C 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 ¹²C and ¹³C (Siegenthaler and Münnich, 1981). Consequently, the surface water δ¹³C signature varies with air–sea exchange efficiency and biological activity intensity. At depth, remineralisation of organic matter releases ¹³C-depleted carbon to the water parcels, which are then mixed through large-scale ocean circulation. In the modern ocean, the global δ¹³C distribution depicts a tight relation between the apparent oxygen utilisation and the δ¹³C signature of a given water mass (Eide et al., 2017). This observation lends support to the use of the δ¹³C of benthic foraminifers as a proxy for ocean oxygen content and ventilation (Duplessy et al., 1988). However, benthic δ¹³C does not solely record deep ventilation changes. As mentioned above, the δ¹³C 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 δ¹³C 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 ¹⁴C content starts to decrease exponentially with time due to radioactive decay (half-life of ¹⁴C $=\mathrm{5730}\pm \mathrm{40}$ years; Godwin, 1962). Thus, by determining the ¹⁴C 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 ¹⁴C 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 Δ¹⁴C in ‰) in terms of past ocean ventilation and circulation changes (e.g. Adkins and Boyle, 1997).

Pa∕Th, benthic δ¹³C and Δ¹⁴C 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 ¹³C and ¹⁴C, under preindustrial or glacial conditions (Bouttes et al., 2015; Brovkin et al., 2002; Menviel et al., 2017; Tschumi et al., 2011). The ²³¹Pa and ²³⁰Th 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, δ¹³C and Δ¹⁴C.

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 ²³¹Pa and ²³⁰Th 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, δ¹³C and Δ¹⁴C. 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 1 Parameters used in the Pa∕Th module.

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 ($\mathrm{3}{}^{\circ }\times \mathrm{3}{}^{\circ }$ 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 ¹³C and ¹⁴C. Previous work has shown that the simulated oceanic δ¹³C and Δ¹⁴C distribution is in reasonable agreement with observations and with available general circulation model (GCM) simulations (Bouttes et al., 2015). We have implemented ²³¹Pa and ²³⁰Th 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, ²³¹Pa and ²³⁰Th 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 ²³¹Pa (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⁻¹. We consider that Pa and Th interact with three particles types: CaCO₃, 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 ²³¹Pa and ²³⁰Th activities are the following (Rempfer et al., 2017):

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial A_{\mathrm{d}}^j}{\partial t}}=T\left(A_{\mathrm{d}}^j\right)+{\mathit{\beta }}^j+k_{\mathrm{desorp}}^j\cdot A_{\mathrm{p}}^j-\left(K_{\mathrm{adsorp}}^j+{\mathit{\lambda }}_j\right)\cdot A_{\mathrm{d}}^j\text{(2)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial A_{\mathrm{p}}^j}{\partial t}}& =T\left(A_{\mathrm{p}}^j\right)-{\displaystyle \frac{\partial (w_{\mathrm{s}}\cdot A_{\mathrm{p}}^j)}{\partial z}}-(k_{\mathrm{desorp}}^j+{\mathit{\lambda }}_j)\cdot A_{\mathrm{p}}^j\\ & +K_{\mathrm{adsorp}}^j\cdot A_{\mathrm{d}}^j,\end{array}\end{array}$$

where $A_{\mathrm{d}}^j$ and $A_{\mathrm{p}}^j$ are respectively the dissolved and particle-bound activities (dpm m⁻³ yr⁻¹) of isotope j (²³¹Pa or ²³⁰Th), β^j (dpm m⁻³ yr⁻¹) is the water column production of isotope ^j by radioactive decay of its uranium parent, λ_j (yr⁻¹) is the decay constant of isotope j, w_s (m yr⁻¹) is the particle settling speed, $K_{\mathrm{adsorp}}^j$ and $k_{\mathrm{desorp}}^j$ (yr⁻¹) are the adsorption and desorption coefficient of isotope j onto particles, respectively, and T is the tracer balance evolution term (dpm m⁻³ yr⁻¹) 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 K_desorp 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):

$$\begin{array}{}\text{(3)}& K_{\mathrm{adsorp}}^j\left(\mathit{\theta },\mathrm{\Phi },z\right)={\sum }_i{\mathit{\sigma }}_{i,j}\cdot F_i(\mathit{\theta },\mathrm{\Phi },z),\end{array}$$

where θ is the latitude, Φ the longitude, z the water depth, σ_i,j the scavenging efficiency for isotope j on particle type i (m² mol⁻¹) and F_i is the particle flux (mol m⁻² yr⁻¹).

2.2 Model tuning and validation

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

$$\begin{array}{}\text{(4)}& K_{d(i,j)}={\displaystyle \frac{{\mathit{\sigma }}_{i,j}\cdot w_{\mathrm{s}}\cdot {\mathit{\rho }}_{\mathrm{sw}}}{M\left(i\right)\cdot k_{\mathrm{desorp}}}},\end{array}$$

where K_d(i,j) is the partition coefficient for isotope j for particle type i, σ_i,j is the corresponding scavenging efficiency, w_s is the settling speed, K_desorp is the desorption coefficient, M_(i) is the molar mass of particle type i (i.e. 12 g mol⁻¹ for POC, 100.08 g mol⁻¹ for CaCO₃ and 67.3 g mol⁻¹ for opal) and ρ_sw is the mean density of seawater (constant fixed to 1.03×10⁶ g m⁻³). Additionally the way each particle type fractionates the Pa and the Th is defined by the fractionation factor F(Th/Pa)_i.

$$\begin{array}{}\text{(5)}& F(\mathrm{Th}/\mathrm{Pa}{)}_i={\displaystyle \frac{K_{d_{(\mathrm{Th},i)}}}{K_{d_{(\mathrm{Pa},i)}}}}={\displaystyle \frac{{\mathit{\sigma }}_{(\mathrm{Th},i)}}{{\mathit{\sigma }}_{(\mathrm{Pa},i)}}}\end{array}$$

K_d 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 (https://CRAN.R-project.org/package=lhs, 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 2 Best fit σ_i,j values and corresponding Kd values.

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 (δ¹³C, Δ¹⁴C) 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 1 Map showing the particulate ²³¹Pa∕²³⁰Th 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 CaCO₃ 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 2 Theoretical 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, δ¹³C and Δ¹⁴C 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 3 Zonally averaged anomalies for Pa∕Th, δ¹³C and Δ¹⁴C 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 δ¹³C and Δ¹⁴C 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 δ¹³C and Δ¹⁴C, 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 δ¹³C decreases in the first 500 m and then increases around 1000 m. The Δ¹⁴C increases from the surface to 1500 m, and Pa∕Th displays no clear trend. In the deep North Atlantic waters, the δ¹³C and Δ¹⁴C 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 δ¹³C and Δ¹⁴C display a dipole pattern increasing at the surface and decreasing at depth. For δ¹³C, the depth boundary is around 1500 m, while for Δ¹⁴C the depth boundary is ∼1500 m north of 50^∘ S and reaches 3000 m south of 40^∘ S.

Figure 4 Zonally averaged times of response (years) for Pa∕Th, δ¹³C and Δ¹⁴C 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 δ¹³C and Δ¹⁴C 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 δ¹³C has frequently more than a single response, especially at depths > 2000 m. The δ¹³C 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 Δ¹⁴C 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 Δ¹⁴C 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 5 Selected 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 δ¹³C 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 δ¹³C and the blue line represents the Δ¹⁴C. 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 δ¹³C and Δ¹⁴C display marked decreases during the hosing, while the Pa∕Th significantly increases. Another interesting feature is that the δ¹³C and Δ¹⁴C 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, δ¹³C and Δ¹⁴C 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 $\sim \mathrm{30}{}^{\circ }$ S), below 3000 m, δ¹³C and Δ¹⁴C 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 (δ¹³C and Δ¹⁴C) 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 δ¹³C and Δ¹⁴C 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 δ¹³C, Δ¹⁴C 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 CO₂ 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 δ¹³C 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 δ¹³C 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 δ¹³C 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 δ¹³C data are available (e.g. Kleiven et al., 2008), to date there are no sufficiently well-resolved Δ¹⁴C 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 δ¹³C, Δ¹⁴C 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 δ¹³C and Δ¹⁴C 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 δ¹³C and Δ¹⁴C 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 δ¹³C decrease is observed for some Heinrich events in deep and intermediate Atlantic waters (Lynch-Stieglitz et al., 2014). Finally, several coral and foraminifer ¹⁴C records also depict a decreased ¹⁴C (decreased Δ¹⁴C) 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 δ¹³C and Δ¹⁴C 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 δ¹³C during the early response in the case of a dual response, while the Δ¹⁴C decreases during the single or early response. We note that the amplitude of this early δ¹³C increase is rather small compared to the late δ¹³C decrease. Moreover, in this region, the time of response is rather short (∼100–150 years). We argue here that in this particular region the δ¹³C displays a two-phase response, the actual long-term response corresponding to the late response.

In the surface and intermediate waters, the modelled δ¹³C decreases above 500 m and then increases around 1000 m, while Δ¹⁴C generally increases from the surface to 1500 m (Fig. 3). This δ¹³C 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 δ¹³C of the photic zone decreases, and the δ¹³C of the remineralisation zone (∼1000 m) increases as the export from ¹³C-depleted carbon is reduced during the hosing compared to the PI equilibrium state. As reported in other model studies, the increased Δ¹⁴C 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 δ¹³C and Δ¹⁴C 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 δ¹³C dataset on the Brazilian margin (Lund et al., 2015; Tessin and Lund, 2013) showing increasing δ¹³C between 1000 and 1300 m and decreasing between 1600 and 2100 m water depth. These results are also fairly consistent with the modelled δ¹³C 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 $\mathrm{3}{}^{\circ }\times \mathrm{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 δ¹³C, 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 δ¹³C decreases by ∼0.5 ‰ across the Heinrich events. Similar Pa∕Th and benthic δ¹³C changes, of smaller/reduced amplitude, are also observed for the DO stadials.

Figure 6 Example of the Bermuda Rise time series. (a) Pa∕Th and benthic δ¹³C (‰) 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 δ¹³C (‰).

Below, we examine our modelled time series in the Bermuda Rise basin ($\sim \mathrm{34}{}^{\circ }$ 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 δ¹³C 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 δ¹³C 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 δ¹³C 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 δ¹³C 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 δ¹³C by about 200 years (Waelbroeck et al., 2018), while at the Bermuda Rise, the Pa∕Th was found to lag behind δ¹³C 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 δ¹³C and Δ¹⁴C 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 ²³¹Pa and ²³⁰Th tracers in the climate model of intermediate complexity iLOVECLIM. The new Pa∕Th module simulates dissolved and particulate ²³¹Pa and ²³⁰Th concentrations with a quality comparable to GCMs under PI equilibrium.

To date, the model is able to simulate the evolution of δ¹³C, Δ¹⁴C 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 Δ¹⁴ show a slight increase, while the δ¹³C significantly increases. In AABW, the three proxies display consistent and synchronous responses: the Pa∕Th increases, while the δ¹³C and Δ¹⁴C decrease, corresponding to a reduced circulation and ventilation of the water mass. Finally, in NADW, the Pa∕Th increases, while the δ¹³C and Δ¹⁴C 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 δ¹³C and Δ¹⁴C decrease; (ii) distinct proxy responses for intermediate and deep waters for δ¹³C and Δ¹⁴C; 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.