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

Abstract. Understanding the ocean circulation changes associated with last glacial abrupt climate events is key to better assess climate variability and understand 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 multi-proxy 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, the timing and sequence of the multi-proxy 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. Without a priori guess, 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 (> 2,000 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 isotopes (δ13C and Δ14C) response lags the sedimentary Pa / Th response by a few hundred years.



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 Published by Copernicus Publications on behalf of the European Geosciences Union. al., 2014), which relies on the analysis of indirect evidence (palaeoproxies).
To date, among the numerous tracers available (e.g. benthic δ 18 O, Nd isotopes, Cd/Ca or sortable silts), the sedimentary ( 231 Pa xs,0 / 230 Th xs,0 ) ratio (hereafter Pa/Th) and carbon isotopes (δ 13 C, 14 C) are key tracers used to reconstruct and quantify past circulation patterns and water mass flow rates.
Pa/Th is the ratio of the 231 Pa and 230 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, 231 Pa and 230 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 231 Pa and 10-40 years for 230 Th; Henderson and Anderson, 2003). Thus, while 230 Th is rapidly transferred to the sediment, 231 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., 2002Chase et al., , 2003Lippold 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 231 Pa and 230 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 . 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 deepwater parcel has been isolated from the surface (e.g. Lynch-Stieglitz et al., 2014;Skinner et al., 2014).
More precisely, the water mass δ 13 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 12 C and 13 C (Siegenthaler and Münnich, 1981). Consequently, the surface water δ 13 C signature varies with air-sea exchange efficiency and biological activity intensity. At depth, remineralisation of organic matter releases 13 Cdepleted carbon to the water parcels, which are then mixed through large-scale ocean circulation. In the modern ocean, the global δ 13 C distribution depicts a tight relation between the apparent oxygen utilisation and the δ 13 C signature of a given water mass (Eide et al., 2017). This observation lends support to the use of the δ 13 C of benthic foraminifers as a proxy for ocean oxygen content and ventilation (Duplessy et al., 1988). However, benthic δ 13 C does not solely record deep ventilation changes. As mentioned above, the δ 13 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 δ 13 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 14 C content starts to decrease exponentially with time due to radioactive decay (half-life of 14 C = 5730 ± 40 years; Godwin, 1962). Thus, by determining the 14 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 14 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 airsea 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 14 C in ‰) in terms of past ocean ventilation and circulation changes (e.g. Adkins and Boyle, 1997).
Pa/Th, benthic δ 13 C and 14 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 . 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 13 C and 14 C, under preindustrial or glacial conditions (Bouttes et al., 2015;Brovkin et al., 2002;Menviel et al., 2017;Tschumi et al., 2011). The 231 Pa and 230 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(Siddall et al., , 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 (latitudelongitude), 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, δ 13 C and 14 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 231 Pa and 230 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, δ 13 C and 14 C. We also developed an analysis method that gives information about the magnitude and timing of the proxy response in the At-lantic 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?

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 dynamicthermodynamic 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 13 C and 14 C. Previous work has shown that the simulated oceanic δ 13 C and 14 C distribution is in reasonable agreement with observations and with available general circulation model (GCM) simulations (Bouttes et al., 2015). We have implemented 231 Pa and 230 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, 231 Pa and 230 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 231 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., 2005Siddall et al., , 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: CaCO 3 , 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 231 Pa and 230 Th activities are the following  (Rempfer et al., 2017): where A j d and A j p are respectively the dissolved and particlebound activities (dpm m −3 yr −1 ) of isotope j ( 231 Pa or 230 Th), β 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 , w s (m yr −1 ) is the particle settling speed, K j adsorp and k j desorp (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 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): where θ is the latitude, the longitude, z the water depth, σ i,j the scavenging efficiency for isotope j on particle type i (m 2 mol −1 ) and F i is the particle flux (mol m −2 yr −1 ).

Model tuning and validation
The scavenging efficiency σ i,j is related to the partition coefficient K d , which defines the proportion of each isotope ( 231 Pa or 230 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.
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 −1 for POC, 100.08 g mol −1 for CaCO 3 and 67.3 g mol −1 for opal) and ρ sw is the mean density of seawater (constant fixed to 1.03 × 10 6 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 .
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).
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.

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 . 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.

Model-data comparison under preindustrial conditions
The main goal of our study is to assess the response of two carbon-based proxies (δ 13 C, 14 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;Mc-Manus 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. 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 3 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).

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 . 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).
We evaluate the response of Pa/Th, δ 13 C and 14 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).
The δ 13 C and 14 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 δ 13 C and 14 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 con-sistent 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 δ 13 C decreases in the first 500 m and then increases around 1000 m. The 14 C increases from the surface to 1500 m, and Pa/Th displays no clear trend. In the deep North Atlantic waters, the δ 13 C and 14 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 δ 13 C and 14 C display a dipole pattern increasing at the surface and decreasing at depth. For δ 13 C, the depth boundary is around 1500 m, while for 14 C the depth boundary is ∼ 1500 m north of 50 • S and reaches 3000 m south of 40 • S.
Looking at the proxy time of response (Fig. 4), we observe significantly different patterns for Pa/Th and the carbon isotopes. The δ 13 C and 14 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 δ 13 C has frequently more than a single response, especially at depths > 2000 m. The δ 13 C displays overall the same pattern as in the western basin except in the south deep basin (between ∼ 20 and 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. 30 • S -below 2500 m) in the case of a dual response. The 14 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 14 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.
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 δ 13 C and 14 C display marked decreases during the hosing, while the Pa/Th significantly increases. Another interesting feature is that the δ 13 C and 14 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, δ 13 C and 14 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 ∼ 30 • S), below 3000 m, δ 13 C and 14 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).

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 (δ 13 C and 14 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.

The lag of carbon isotopes, underlying physical mechanisms
One of the most interesting features of our simulations is that the δ 13 C and 14 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 δ 13 C, 14 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 trans- port 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 airsea 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 centennialto 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., 2008Menviel et al., , 2015, the freshwater addition in the Nordic Seas produces a surface cooling of about 1 • C associated with an increase in the CO 2 solubility in surface waters and a general marine productivity decrease in the North Atlantic. Subsequent to the NADW cessation, the nutrients https://doi.org/10.5194/cp-16-867-2020 Clim. Past, 16, 867-883, 2020 accumulate in the Atlantic intermediate waters, shaping a negative δ 13 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 δ 13 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).

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, dura-tion 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 δ 13 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 δ 13 C data are available (e.g. Kleiven et al., 2008), to date there are no sufficiently well-resolved 14 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 δ 13 C, 14 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 δ 13 C and 14 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(Burckel et al., , 2016Mulitza 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 δ 13 C and 14 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 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.

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 ‰ ben-thic δ 13 C decrease is observed for some Heinrich events in deep and intermediate Atlantic waters (Lynch-Stieglitz et al., 2014). Finally, several coral and foraminifer 14 C records also depict a decreased 14 C (decreased 14 C) content of deep waters and increased deep ventilation ages during HS1 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 δ 13 C and 14 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 δ 13 C during the early response in the case of a dual response, while the 14 C decreases during the single or early response. We note that the amplitude of this early δ 13 C increase is rather small compared to the late δ 13 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 δ 13 C displays a two-phase response, the actual long-term response corresponding to the late response.
In the surface and intermediate waters, the modelled δ 13 C decreases above 500 m and then increases around 1000 m, while 14 C generally increases from the surface to 1500 m (Fig. 3). This δ 13 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 δ 13 C of the photic zone decreases, and the δ 13 C of the remineralisation zone (∼ 1000 m) increases as the export from 13 C-depleted carbon is reduced during the hosing compared to the PI equilibrium state. As reported in other model studies, the increased 14 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 δ 13 C and 14 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 δ 13 C dataset on the Brazilian margin Tessin and Lund, 2013) showing increasing δ 13 C between 1000 and 1300 m and decreasing between 1600 and 2100 m water depth. These results are also fairly consistent with the modelled δ 13 C de-scribed 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 δ 13 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.

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 δ 13 C decreases by ∼ 0.5 ‰ across the Heinrich events. Similar Pa/Th and benthic δ 13 C changes, of smaller/reduced amplitude, are also observed for the DO stadials. 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 δ 13 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 δ 13 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 δ 13 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 δ 13 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 δ 13 C by about 200 years , while at the Bermuda Rise, the Pa/Th was found to lag behind δ 13 C by about 200 years (Henry et al., 2016). However, the latter result must be interpreted with caution since the crosscorrelation 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 . 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 δ 13 C and 14 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 . 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.

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