Atlantic Meridional Overturning Circulation and δ13C Variability During the Last Interglacial

The Atlantic Meridional Overturning Circulation (AMOC) is thought to be relatively vigorous and stable during Interglacial periods on multimillennial (equilibrium) timescales. However, recent proxy (δ13C benthic) reconstructions suggest that higher frequency variability in deep water circulation may have been common during some interglacial periods, including the Last Interglacial (LIG, 130–115 ka). The origin of these isotope variations and their implications for past AMOC remain poorly understood. Using iLOVECLIM, an Earth system model of intermediate complexity (EMIC) allowing the computation of δ13CDIC and direct comparison to proxy reconstructions, we perform a transient experiment of the LIG (125–115 ka) forced only by boundary conditions of greenhouse gases and orbital forcings. The model simulates large centennial‐scale variations in the interior δ13CDIC of the North Atlantic similar in timescale and amplitude to changes resolved by high‐resolution reconstructions from the LIG. In the model, these variations are caused by changes in the relative influence of North Atlantic Deep Water (NADW) and southern source water (SSW) and are closely linked to large (∼50%) changes in AMOC strength provoked by saline input and associated deep convection areas south of Greenland. We identify regions within the subpolar North Atlantic with different sensitivity and response to changes in preformed δ13CDIC of NADW and to changes in NADW versus SSW influence, which is useful for proxy record interpretation. This could explain the relatively large δ13C gradient (∼0.4%0) observed at ∼3 km water depth in the subpolar North Atlantic at the inception of the last glacial.


Introduction
Variation in the Atlantic Meridional Overturning Circulation (AMOC) is one of the major driving force controlling climate changes as it affects the distribution of heat and carbon in the ocean, thus influencing regional climate and atmospheric CO 2 concentrations (Ganachaud & Wunsch, 2000). This has made AMOC changes a key concern in future projections, with models suggesting it will likely decrease in response to buoyancy gain in the deep-water formation regions but the extent to which remains highly uncertain (Bakker et al., 2016;Stocker et al., 2014). Meanwhile, recent observations hint that the ocean circulation may have already started changing (Caesar et al., 2018;Smeed et al., 2014).
Analyzing past variability of AMOC can inform us on its stability and tendency to change. The Interglacial periods of the late Pleistocene are particularly interesting in this regard as they shared overall similar boundary conditions with the modern climate (Berger et al., 2015). The climate of the Last Interglacial (LIG), or Marine Isotope Stage (MIS) 5e in isotopic marine series, also shared similar features with the model projections of our future climate if anthropogenic greenhouse gas emissions continue unabated: high-latitudes warming (Hoffman et al., 2017;Otto-Bliesner et al., 2013), reduction of Greenland ice-sheet, and higher sea level (Kopp et al., 2009;Otto-Bliesner et al., 2006). Hence, the high-latitude changes occurring during the warm MIS 5e could provide additional insight on AMOC variability. Data reconstructions from this period (130-115 ka) characterize this circulation with a periodically warmer and fresher than today's North Atlantic water (Hoffman et al., 2017;Otto-Bliesner et al., 2006). Reconstructions of millennial-timescale variability indicate that vigorous North Atlantic Deep Water (NADW) production persisted on millennial timescales (Adkins et al., 1997;Oppo et al., 1997), simulztaneously as southern source waters (SSWs) expanded The shorter term transient behavior of the ocean circulation during this warmer period remains, however, poorly constrained as is it difficult to robustly identify and assess such short-timescale variability in the marine record. Nevertheless, new generations of high-resolution (down to multidecadal) proxy reconstructions, more appropriate to resolve multicentennial variability, suggest the existence of centennial short-lived large perturbations of the NADW during previous interglacial periods (Galaasen et al., 2014;Hodell et al., 2009;McManus et al., 1999;Ninnemann & Charles, 2002). Data reconstructions of δ 13 C-an ocean circulation and carbon cycle tracer-from sediment cores in the North Atlantic depict abrupt centennial-scale variations on the order of 0.7%0, comparable to the deep Atlantic δ 13 C changes observed during glacial periods and associated with changes in the distribution and ventilation of NADW (Henry et al., 2016;McManus et al., 1999McManus et al., , 2004Menviel et al., 2017). This suggests possibly large and abrupt variations in ocean ventilation, which challenges the paradigm of the stability of interglacial thermohaline circulation at short timescales (Galaasen et al., 2014;Hodell et al., 2009). While the existence of such AMOC change is consistent with some observed climate variability Galaasen et al., 2015;Mokeddem et al., 2014;Tzedakis et al., 2018;Zhuravleva & Bauch, 2018), the relationship between benthic δ 13 C and ocean circulation is not straightforward and requires additional tool such as model simulations to evaluate their degree of interactions (Bakker et al., 2015).
With respect to the model-data comparison, numerous studies have analyzed the δ 13 C distribution variations under various boundary conditions using, for example, two-dimensional ocean models (Bouttes et al., 2009(Bouttes et al., , 2010(Bouttes et al., , 2012Brovkin et al., 2007) and more complex three-dimensional models (Bakker et al., 2015;Menviel et al., 2017). These studies propose linkages between the millennial-scale variations of δ 13 C and changes in a number of processes such as water mass distribution changes (which can be due to interaction with sea-ice formation), biological activity changes (e.g., due to iron fertilization), and land vegetation modification. However, to the authors' knowledge, there are no modeling studies, which simulate abrupt and large modification in the bottom water δ 13 C with the amplitude and timescale depicted by the high-resolution reconstruction data (Galaasen et al., 2014;Hodell et al., 2009), such as shown here. Such constraints are needed to assess, for example, whether it is mechanistically plausible for changes in water mass distribution to have driven abrupt centennial-scale variability in bottom water δ 13 C and, furthermore, whether such tracer field adjustments could be linked to AMOC variability. Determining the characteristic pattern or fingerprint in interior ocean δ 13 C related to AMOC changes would provide a more realistic and physically consistent framework for interpreting past circulation variability and its implication for the carbon cycling but also in understanding its recent (Caesar et al., 2018) and potential future changes (Weijer et al., 2019). Here, we span this knowledge gap by using iLOVECLIM, an Earth system model of intermediate complexity (EMIC), to perform a transient simulation of the LIG (125-115 ka) under past natural variations in natural greenhouse gases and transient orbital forcings.
The paper is organized as follows: In section 2, we describe the model, the experiment design, and the terms and metrics used to analyze the water mass geometry in the interior ocean during the transient simulation. Section 3 presents the results of the model simulations, while discussions and comparison with previous studies are presented in section 4. Finally, the study is summarized in section 5.

Paleoceanography
10.1029/2019PA003818 Maqueda (1997, 1999). The terrestrial biosphere model (VECODE) is composed of three submodels that exchange heat, stress, water, and carbon and are designed for long-term simulations (Brovkin et al., 1997). In addition to desert, the dynamic vegetation model simulates two types of plants-trees and grass-which are subdivided into four compartments (leaves, wood, litter, and soil) that exchange carbon in the biogeochemical model. The model of vegetation structure computes the fraction of plant functional type (PFT) in equilibrium with the climate where photosynthesis depends on precipitation, temperature, and atmospheric CO 2 .
The ocean carbon cycle model is divided into organic and inorganic parts that are based on the nutrientphytoplankton-zooplankton-detritus (NPZD) model described in Six and Maier-Reimer (1996). The inorganic carbon is represented by the dissolved inorganic carbon (DIC) and alkalinity (ALK), while the organic part includes phytoplankton, zooplankton, dissolved organic carbon (DOC), slow dissolved organic carbon (DOCs), particulate organic carbon (POC), and calcium carbonate (CaCO 3 ). The plankton is partially remineralized as it sinks through the water column, while all the POC and CaCO 3 are remineralized at depth. The remineralization profile follows an exponential law; however, it is adjusted to have less remineralization in the upper layers and more at depth. At the air-sea interface, the carbon flux is computed from the CO 2 partial pressure (pCO 2 ) difference between the atmosphere and the ocean at a constant gas exchange coefficient of 0.06 mol·m −2 ·yr −1 . The sea surface pCO 2 is a function of temperature, salinity, DIC, and ALK following Millero (1995). The iLOVECLIM model is using a similar carbon cycle as in CLIMBER-2 where the 13 C is computed in the ocean and terrestrial biosphere as in Brovkin et al. (2007) and is fully described in Bouttes et al. (2015).

Experiment Set Up
We have performed a transient experiment over the 125-to 115-ka period of the LIG by prescribing yearly interpolated values of greenhouse gases (CO 2 , CH 4 , and N 2 O) and orbital forcings from the third phase of the Paleoclimate Modelling Intercomparison Project (PMIP3; https://pmip3.lsce.ipsl.fr/). This experiment is branched off after two consecutive runs in order to set the initial conditions: (1) a preindustrial spin up followed by (2) an equilibrium experiment using 125-ka boundary conditions. Both experiments (1) and (2) are simulated over 5,000 model years. The preindustrial spin-up reproduces well the main ocean carbon cycle characteristics and is similar to that described in Bouttes et al. (2015). The atmospheric δ 13 C and CO 2 concentrations are −6.5%0 and 282 ppm, respectively, while the ocean and terrestrial biosphere carbon content are 38,952 and 2,094 PgC. The most important features of the different water masses are reproduced and the oceanic δ 13 C (δ 13 C DIC ) is relatively well simulated as compared to the data (Eide et al., 2017;Schmittner et al., 2013), with higher δ 13 C at the surface due to photosynthesis and lower δ 13 C in the deeper ocean where remineralization takes place ( Figure 1a). It also depicts the water mass characteristics of NADW (high δ 13 C) and Antarctic Bottom Water (AABW; low δ 13 C) and the equatorial subsurface minimum induced by strong remineralization. However, the high δ 13 C values from the NADW penetrate further south in the data reconstruction ( Figure 1b) than that simulated in our model, which is likely the result of too strong diffusion in the model (Bouttes et al., 2015). In other regions, the magnitude and the large-scale spatial pattern of the δ 13 C DIC tracer are fairly comparable to the observations.
In the 125-ka equilibrium run, that is, experiment (2), the model simulation reaches equilibrium after 3,000 model years where the ocean and land vegetation carbon content stabilize at 38,857 ± 1 PgC and 23,235 ± 4 PgC, respectively, over the proceeding 2,000 simulation years. The atmospheric CO 2 and the other greenhouse gases are kept constant (CO 2 = 276 ppmv; CH 4 = 640 ppb and N 2 O = 263 ppb). The end of the 125-ka equilibrium run is used as a starting point for the 125-to 115-ka transient simulation. Changes into the carbon reservoirs are simulated during this period corresponding into a reduction of the ocean carbon content at 125 ka as compared to 115 ka of about 360 PgC. This value is of the same order of magnitude with that suggested by previous studies analyzing the ocean carbon system during the LIG, for example, Schurgers et al. (2006) and Kessler et al. (2018) with 310 and 314 PgC, respectively. However, contrary to the data reconstruction that shows relative stable atmospheric CO 2 concentration around 270-280 ppm throughout the LIG (Lourantou et al., 2010;Schneider et al., 2013), our model simulates an increase after 123.5ka until it stabilizes again around 116 ka to 305 ppm (nearly 30 ppm above the data reconstruction). This increase is attributed to the strong reduction in land vegetation carbon content (∼420 PgC) induced by gradual cooling of the Northern Hemisphere (∼2 • C) and potentially the absence of permafrost in the model. We note that this experiment has also been used in Galaasen et al. (2020) to illustrate the similarity between simulated and reconstructed bottom water δ 13 C changes in recent interglacials. Hereafter, we provide a detailed analysis of this simulation defining the cause of the simulated δ 13 C DIC anomalies, their spatial patterns in the deep North Atlantic, and potential drivers of δ 13 C DIC variability.

13 C Tracer
The δ 13 C represents the standardized 13 C/ 12 C ratio (Zeebe & Wolf-Gladrow, 2001) and is expressed in permil units as follows: where ( 13 C/ 12 C) standard is the Pee Dee Belemnite carbon isotope standard (Craig, 1957). Hereafter, we consider the δ 13 C of DIC in the ocean, noted δ 13 C DIC . During the photosynthesis, marine and terrestrial biology preferentially use the lighter 12 C over 13 C. When organic matter is remineralized at depth, more 12 C is released as compared to 13 C, thus producing low δ 13 C DIC . On the opposite, at the surface ocean, the excess sequestration of 12 C, relative to 13 C, during photosynthesis drives an increase in δ 13 C DIC . Thus, the distribution of δ 13 C DIC is affected by air-sea gas exchange and marine biological fractionation (Eide et al., 2017;Lynch-Stieglitz et al., 1995;Schmittner et al., 2013;Zhang et al., 1995) and as other oceanic tracers, is transported by the advection-diffusion scheme in the oceanic model. The analysis of δ 13 C DIC is often used to reconstruct the variations in ocean circulation (Duplessy et al., 1988) and strength of biological pump (Crucifix, 2005;Curry & Oppo, 2005;Eide et al., 2017;Morée et al., 2018).

Phosphate Oxygen Tracer (PO)
The phosphate oxygen tracer (PO) as defined by Broecker (1974) is a measure for identifying water masses. It is computed using the simulated phosphate and oxygen fields: where r O:P is the oxygen-to-phosphate stoichiometric ratio, which is varying from 138 (surface) to 170 (bottom), in the model. For a given water mass, this tracer is presumed to be approximately constant and is based on the principle that phosphate is released while oxygen is used during remineralization, and vice versa during biological production. The distinction of water masses using PO is useful for contrasting interior water masses with very different surface PO values. Here, we mainly use PO to identify northern source water (NSW) and SSW masses in the deep ocean below 1,000 m depth characterized by low and high PO values, respectively.

Apparent Oxygen Utilization (AOU)
The apparent oxygen utilization (AOU) is a diagnostic biogeochemical tracer that provides an estimate of the amount of oxygen used since a water parcel was last at the surface and hypothetically fully saturated. It is computed according to where O sat 2 is the saturation value and O 2 is the in situ oxygen value calculated by the model. When a water parcel is advected to the interior, remineralization processes start and consume the oxygen. The older a water parcel, the more oxygen has been consumed by remineralization; hence, AOU increases.

Results
Section 3.1 describes the simulated variability in δ 13 C DIC and AMOC strength, while a composite analysis opposing strong versus weak AMOC is discussed in section 3.2. The mechanisms associated with the variability of the meridional overturning circulation are addressed in section 3.3. In order to compare and contrast our simulation with data records, we focus our discussion to the locations of the two high-resolution proxy records MD03-2664 (Galaasen et al., 2014) and Integrated Ocean Drilling Program (IODP) Site U1304 (Hodell et al., 2009), which depict centennial-scale δ 13 C perturbations during the LIG. The location of these two sites is shown, for instance, in Figure 2a.

AMOC and 13 C Transient Variability
The time series of simulated δ 13 C DIC at the two sediment core locations are shown in Figure 2a. The absolute values of δ 13 C DIC generally differ from one location to the other by about 0.4%0 but converge during peaks . Atlantic section of (a,b) overturning stream function, (c,d) PO tracer, and (d,e) δ 13 C DIC . The panels on the left (a,c,e) correspond to the average over all the years of weak AMOC (i.e., two standard deviations below the long-term mean). The right side panels (b,d,f) correspond to the strong AMOC points average (i.e., two standard deviations above the long-term mean). The black dashed-lines on the top panels represent the zero line and symbolize the separation between the upper to bottom cells. In the middle panels, the region above the black solid lines show the approximate area that is predominantly occupied by NSW. The value (∼0.36 mol m −3 ) represents the averaged surface PO over the sinking region south of Greenland.
in maximum values to a relative similar value (e.g., 116.8 and 117.8 ka). In addition, both MD03-2664 and U1304 sites show similar evolution during the transient simulation. While δ 13 C DIC is relatively stable from 125-122ka (1.28 ± 0.05%0 at MD03-2664 and 0.96 ± 0.04%0 at U1304), stronger variations are simulated during the second half of the experiment where the standard deviation is increased by a factor of three at U1304 and two at MD03-2664. These abrupt changes in bottom δ 13 C DIC are simulated at centennial timescale for both locations and show a maximum amplitude of variation of about 0.65%0 and 0.43%0 around 119 ka. This is comparable to the changes depicted by the proxy reconstruction at MD03-2664 (Galaasen et al., 2014), which also highlights similar centennial variations of about 0.7 %0 of magnitude. In addition, small decreasing trend are simulated at MD03-2664 (3.2 × 10 −5 %0·yr −1 ) and U1304 (3.5 × 10 −5 %0·yr −1 ), which is attributed to the long-term decreasing trend of the atmospheric δ 13 C and the decreasing land vegetation carbon reservoir.
The variations in bottom water δ 13 C DIC at both locations follow the same variations as depicted by AOU (Figure 2b). Site U1304 simulates in general higher AOU (∼ 40%) than that at MD03-2664, suggesting that the water mass affecting Site U1304 is more remineralized than at MD03-2664. However, during the abrupt transitions, the difference between both sites is reduced, meaning that both sites are affected by a water mass that has been similarly ventilated (i.e., similar water mass). Figure 2c shows the evolution of maximum overturning stream function at 27 • N in the Atlantic basin with an averaged value of 12.2 ± 1.8 (one standard deviation) Sv and shows similar variations to those depicted with δ 13 C DIC and AOU. High values of overturning stream function coincide with high bottom water values of δ 13 C DIC at both sediment core sites, which seem to converge toward the same δ 13 C DIC value. Conversely, a weak overturning stream function corresponds to lower values of δ 13 C DIC and lead to a divergence of δ 13 C DIC values between the two core locations. This suggests that the temporal variation of δ 13 C DIC in the North Atlantic is related mechanistically to or even directly modulated by the strength of the AMOC in our model. In order to further analyze the impact differences in AMOC strength have on the biogeochemistry and ocean δ 13 C DIC , we show a composite analysis of two contrasting states of AMOC, that is, weak versus strong AMOC. For this, we average all years where AMOC strength is two standard deviations above or below the long-term mean. The selected years with strong AMOC states are represented by the purple points ( Figure 2c; number of point n = 460), while the selected years with weak AMOC states are shown by the purple stars (n = 63).  (Figures 3c and 3d). This result is in good agreement with the deepening of the upper cell of the AMOC previously mentioned. We note that in both cases, the NADW (NSW) is sustained depicting its relative PO signature down to 2,300 (3,000) m depth under weak (strong) AMOC state. In addition, because NSW has higher δ 13 C DIC values than SSW, δ 13 C DIC -rich water masses dominate MD03-2664 and U1304 locations when the AMOC is strong (Figure 3f), which is consistent with the convergence of δ 13 C DIC values at both core sites previously mentioned (section 3.1). On the other hand, more low δ 13 C DIC water masses, resulting from the mixing of SSW and NSW, are affecting the sediment cores during periods of weaker AMOC (Figure 3e).

Composite Analysis: Weak Versus Strong AMOC
These two AMOC states are characterized by distinct differences in the location, area, and depth of convection in the North Atlantic. In the case of weak AMOC, deep convection is simulated in the Norwegian sea (>400 m depth; Figure 4a), while an area of shallower convection extends further south into the North Atlantic. Conversely, under strong AMOC periods, the Irminger Sea region including the southeast Labrador Sea also experiences deep convection down to 150 m depth (Figure 4b), while the convection depth in the Norwegian Sea slightly weakens to 300 m and extends further north under the sea-ice limit. The strongest sea-ice retreat (up to 10 • northward) is simulated in the Irminger Sea and coincides with the activation of convection in this region (Figure 4, light blue dashed lines). The changes depicted by the density field (Figures 4c and 4d) correspond to the changes in convection depth. Under simulated weak AMOC, the densest water masses are in the Norwegian Sea with values of neutral density up to 27.8 kg·m −3 where the convection occurs. When the AMOC strengthens, the surface density increases south of Greenland and in the Irminger Sea, leading to similar density values with that in the Norwegian Sea (27.8 kg·m −3 ; Figure 4d). Thus, the total area of dense water formation in these regions increases during the AMOC strong composite state. Figure 5 shows the annual mean surface temperature and salinity fields during both strong and weak AMOC states. It reveals that lower SSTs are simulated in the regions corresponding to the subpolar gyre (Irminger Sea and southeast of the Labrador Sea) and in the Norwegian Sea when the AMOC weakens (Figures 5a  and 5b). For example, the 0 • C isotherm shifts north by about 10 • between Greenland and Norway while it only moves north by about 2 • in southeastern Labrador Sea (Figure 5b). On the other hand, the warmest waters expand westward (i.e., isotherm 8 • C) when AMOC strengthens and onset/increase of convection in the Labrador Sea/Irminger Sea. This suggests strong negative buoyancy flux in these regions due to the heat loss. The SSS field also depicts strong modification in the North Atlantic, especially in the region where the convection turns off and on (Figures 5c and 5d). Under weak AMOC, the strongest SSS (>35.5 psu) are simulated around the Celtic Sea and extend northward toward the Norwegian Sea (Figure 5c), while fresher polar surface waters (<33.6 psu) are being transported along the eastern Greenland coastline following the sea-ice extent. When the AMOC strengthens, the high surface salinity spreads further north into the Norwegian Sea and west toward the Irminger Sea and the entrance of the Labrador Sea (Figure 5d). This suggests that more warm and salty surface Atlantic water is entering the Norwegian Sea under stronger ocean circulation and recirculates to the northwestern basin following the northward shift of the sea-ice edge between Greenland and Iceland (Figure 5d). These increases in SSS cause the surface density to increase, reducing the stratification of the ocean and establishing favorable conditions for deep convection to occur.
While changes in AMOC are expressed clearly in δ 13 C DIC at both locations (Figures 2a and 2c), the model predicts slightly different responses at these locations due to the position of the sites relative to major water mass pathways. MD03-2664 is located further North and West and is more influenced by NSW during both strong and weak AMOC composite states (Figures 6e and 6f). Thus, during weak composite states, bottom water δ 13 C DIC at MD03-2664 does not drop as low as at U1304, which has much greater influence of SSW than MD03-2664. Figure 6 reveals these differences by depicting the changes in δ 13 C DIC , AOU, and PO tracers at the sea floor. The pattern of these three tracers are strongly similar showing that δ 13 C DIC changes are related to water mass geometry. Under weak AMOC states, MD03-2664 depicts 105.0 μmol·kg −1 of AOU when U1304 has 146.1 μmol·kg −1 (Figure 6c), corresponding to 28% of difference. This difference is reduced to 13% when the AMOC strengthens ( Figure 6d); that is, the AOU values in MD03-2664 and U1304 depict 82.5 and 94.6 μmol·kg −1 , respectively. However, U1304 AOU values are generally higher than at MD03-2664, which suggests that the water mass affecting U1304 contains more remineralized carbon than at MD03-2664 during both strong/weak AMOC states. Consequently, the bottom water δ 13 C DIC values simulated at U1304 are lower compared to MD03-2664, and this difference is even more pronounced when the AMOC weakens ( Figure 2a). Finally, the PO tracer confirms that both sites are affected by different water masses under weak AMOC, while this difference strongly decreases when the ocean circulation strengthens. Thus, under weak AMOC, the MD03-2664 location depicts weaker PO value than at U1304 (Figure 6e; 0.35 against 0.42 mol·m −3 , respectively), while almost the same values are depicted under strong AMOC (Figure 6f; 0.32 against 0.34 mol·m −3 , respectively). Figure 7 shows the time series of the vertical density gradient for the region where the deep convection turns on and off. The oscillations of the density gradient south of Greenland covary with the AMOC strength during the centennial peak events. As the stratification decreases (density gradient weakens), the AMOC strength increases. In the model, this suggests that the stratification of the ocean south of Greenland is at least partially linked to the AMOC strength during these centennial peaks. The salinity field is the main driver for the changes in density at the surface. The amplitude of oscillations is weaker at the beginning of the simulation than toward 115 ka. The lower values remain relatively constant around 0.5 × 10 −3 kg·m −4 but the maximum values increase toward 115ka (up to 3.6 × 10 −3 kg·m −4 ), especially after 122 ka. This increasing trend in maximum of density gradient is attributed to fresher surface ocean simulated toward 115 ka in that area potentially induced by stronger variability in the growing-melting of sea-ice. Hence, the surface salinity in this region is in averaged 2 psu lower in the late LIG (116-115 ka) than at the beginning (125-124 ka). The variation of surface density are more than six times greater at the surface (27.86 ± 0.33 kg·m −3 ) than at 600 m depth (27.12 ± 0.05 kg·m −3 ) and are most likely also due to the sea-ice growth and retreat in this area.

AMOC Oscillations
In order to identify the mechanisms behind the changes in surface density, we analyze the onset and offset of the AMOC peak events. We focus on the event occurring around 116.7 ka as an example because it corresponds to the strongest variation of both AMOC and δ 13 C DIC ; however, every oscillation depicted by  Figure 2c shows similar characteristics. We note yr min as the year when the AMOC reaches its last minimum value before increasing (yr min = 116.799). Similarly, we note yr max as the year when the AMOC reaches its last maximum value before decreasing and ending the strong AMOC peak event (yr max = 116.657, 142 years after yr min ). Figure 8 depicts the salinity in latitudinal (left panels) and longitudinal (right panels) sections of the North Atlantic Ocean. The first row represents the salinity of the first 1,000 m depth in the North Atlantic Ocean at yr min . The second row represents the same sections at yr min + 11 years of simulation and the last row corresponds to the end of the onset peak at yr min + 18 years. These years represent well the changes in the water mass distribution associated to the multidecadal initiation (and offset) of the abrupt transition. Showing any other years within this interval would provide similar results. At the onset of the AMOC peak (yr min ), the AMOC is at the minimum of 9.2 Sv. A 100-to 200-m layer of fresh water extends from the high latitudes down to almost 45 • N and spreads from the Greenland coast eastward to 25 • W leading to a stratified ocean surface (Figure 8, first row panels). This fresh water is associated to sea-ice cover. The Atlantic water corresponds to the higher salinity at 45 • N below 200 m depth (Figure 8). It extents northward along the Norwegian coast and depicts salinity >35.75 psu. At yr min + 11 (panels on the second row), the volume of fresh water east of the Greenland coast reduces and starts retreating northward. Similarly, the salinity increases from the eastern side and traces of convection begin to appear around 25 • W, mixing the salinity through the water column. Hence, the AMOC slightly speeds up and increases by +1.7 Sv to reach 10.9 Sv. At the end of the onset (yr min + 18), the AMOC is near its maximum (16.7 Sv) and the salinity north of 60 • N increases. By contrast, between 50 and 55 • N, the surface water is simulated fresher and may be associated to the transport of the fresh water from the sea-ice melt. The intrusion of salty water mass seems to originate from the eastern side of the basin where the saline Atlantic waters expand from the Norwegian coast toward the Greenland coastline. As a result, saline Atlantic waters fill up the basin along the 61.5 • N latitude band creating favorable conditions for convection south of Greenland. The same ocean structure is simulated during the time of the maximum AMOC, which corresponds to a 124 years window. Figure 9 is similar to Figure 9. Similar to Figure 8, but for the offset of the upward AMOC trend. Here, yr max represents the last maximum value of the AMOC in 116.657 ka (142 years after yr min ) before it starts to decrease to finally reach its minimum at yr max + 14. Figure 8, but represents the offset of the AMOC peak, from yr max to the next minimum (yr max + 14). This decreasing AMOC phase (increasing stability) mirrors the phase of increasing AMOC (decreasing stability) and shows a gradual establishment of the freshwater layer associated with sea-ice growth. No significant differences compared to the onset could be identified.

Discussion
The ocean circulation is one of the key processes affecting the ocean carbon storage, hence regulating the atmospheric CO 2 levels on millennial timescale. The production of NADW during the LIG is vigorous and appears relatively robust at (multi)millennial timescales, that is, suggested by numerous of studies. However, new high-resolution proxy records reveal that abrupt changes in bottom water chemistry (δ 13 C) occurred on centennial timescales during the previous interglacial and are postulated to represent ocean circulation changes. Yet, it is less clear how these circulation changes relate to AMOC or, indeed, even if deep circulation changes could explain such rapid changes in deep Atlantic δ 13 C-onset and terminated in decade(s) and generally persisting for one or few centuries (Galaasen et al., 2014). In this study, we simulate a transient experiment from the last interglacial period (125-115 ka) using an Earth system model of intermediate complexity (EMIC) allowing the computation of δ 13 C DIC , an oceanic tracer for the ocean circulation, and biogeochemical processes used in paleoclimate reconstruction. Significant changes of δ 13 C DIC are simulated in the North Atlantic at short timescale (centennial), which approach in magnitude those depicted in the data reconstruction. Likewise, the transitions (increases and decreases) also occur quickly, consistent with the rapid onset and demise observed in the proxy reconstruction. The variability in bottom water δ 13 C DIC is found to be positively correlated to the AMOC strength, which oscillates between 6.3 and 19.3 Sv. These AMOC oscillations are mainly controlled by the convection south of Greenland (Labrador Sea and Irminger Sea), which turns off and on with the northward and westward expansion of highly saline Atlantic waters from the Norwegian coast to Greenland around 60 • N, associated with the northward retreat of sea-ice and fresh surface waters. The representation of the North Atlantic ventilation is coherent with that suggested by previous studies. The mid-depth North Atlantic water remains relatively well ventilated during the entire period of simulation, even during the simulated weak AMOC events (Figure 3c). This persistent production of NSW at mid-depth from 125 to 115 ka is also suggested by other studies based on proxy reconstructions (McManus et al., 2002;Mokeddem et al., 2014) and previous independent model simulations (Born et al., 2011;Wang & Mysak, 2002). Moreover, our results also support previous interpretations of proxy data concerning the deep North Atlantic water (Galaasen et al., 2014;Hodell et al., 2009), that is, that ventilation of the deep Atlantic was intermittently disturbed during the LIG and that AMOC strength and NSW-SSW shifts could plausibly explain North Atlantic deep water δ 13 C changes during interglacials. Thereby, as shown in Galaasen et al. (2020), the simulated deep Atlantic δ 13 C DIC anomalies are strongly similar in magnitude (0.5-0.7%0) and temporal characteristics (decadal scale at onset and persisting for centuries) to the variability in high-resolution δ 13 C reconstructions. Assuming that the reconstructed proxies reasonably capture the core top δ 13 C values, our model simulation suggests that the region south of Greenland would be most suitable to capture any apparent variations in the watermass properties over the LIG period while regions located further North and East (e.g., in the Nordic Seas) remains more stable. Figure 10 highlights this point by showing higher bottom water δ 13 C DIC standard deviation south of Greenland as compared to any other regions of the North Atlantic. In our simulation, when AMOC is strong, both sediment core sites are affected by the same water mass (NSW) with similar biophysical properties leading to similar δ 13 C DIC values. Conversely, when AMOC weakens, the influence of SSW at Site U1304 increases; hence, the δ 13 C DIC values of both locations diverge. This suggests that relatively large differences in tracer (δ 13 C) signals between two cores located in relatively close proximity can occur due to the changes in water mass distributions (NSW vs SSW). Using the same model for the period 132-120 ka, Bakker et al. (2015) have also shown a close correlation between North Atlantic bottom water δ 13 C DIC and AMOC strength on multimillennial timescale. Our study further supports this result on (multi)centennial timescale.
In addition to those rapid shifts, the proxy data records from MD03-2664 and Site U1304 show markedly different long-term trends in bottom water δ 13 C from 122 to 115 ka (see figure 3 in Galaasen et al., 2014). While bottom water δ 13 C at MD03-2664 increases through the late LIG, it decreases at Site U1304, resulting in a gradient of ∼0.4%0 between the western (MD03-2664) and eastern subpolar basins (Site U1304) near 115 ka and the last glacial inception. In our model simulation, the long-term transient trends deviate somewhat from the actual sites and show at both location a small decreasing trend. This discrepancy with data reconstruction can be attributed, at least in part, to the coarse spatial resolution of the model and the difficulty to properly simulate overflows in this region. However, much of our model simulation shows deep North Atlantic δ 13 C DIC gradient of ∼0.4%0 between MD03-2664 and Site U1304 (Figure 2a) similar to the proxy records toward 115 ka. In the model, this difference is caused by relatively higher influence of SSW in the eastern basin and Site U1304 (Figure 6e), which is restricted by the bathymetry from entering the western basin. Therefore, we suggest that the observed ∼0.4%0 δ 13 C gradient for the late LIG could be explained from such water mass geometry state, that is, the expansion of the SSW into the subpolar North Atlantic in the late LIG (Govin et al., 2009).
While similar in character, the forcing of the centennial-scale LIG variability likely differs between the simulated and reconstructed events. The reconstructed variability focuses around the early phase of the LIG and appears closely associated with discharge of freshwater and icebergs, including a freshwater outburst flood event associated with the final retreat of the Laurentide ice sheet (Galaasen et al., 2014;Hodell et al., 2009;Nicholl et al., 2012). In the model, the variability occurs in the mid-LIG to late LIG and is associated with changes in the distribution of sea ice and freshwater within the subpolar North Atlantic and Nordic Seas. This discrepancy from the proxy data arises mainly from the initial conditions used to begin the transient simulation-equilibrium state at 125 ka-which removes potential effects from reminiscent glacial conditions especially at the beginning of the simulation. Hence, the simulated absolute values of δ 13 C DIC may not correspond to the reality from that period. However, the climate response to orbital and greenhouse gases forcings remain sufficient to create large modifications in the interior ocean carbon distribution and water mass geometry, which comes to support the short-lived large perturbations of the climate system under warmer climate conditions than today. These rapid shifts in deep water occur as the internal ocean adjusts to switches in AMOC strength. These shifts in AMOC strength are due to changes in surface ocean convection sites and in particular involve turning on and off convection in the Nordic Seas. Similar variations in overturning geometry and strength tied to the location and strength of deep convection have been invoked previously in particular to explain abrupt, millennial-scale ocean-climate changes observed during the last glacial cycle (Rahmstorf, 2002). Our simulation provides an example of how similar shifts in convection-AMOC strength could explain centennial scale variability seen in proxy records during interglacial boundary conditions.
The multicentennial scale variability of AMOC simulated by our model is closely linked to the occurrence of deep convection events south of Greenland. Here, transport of anomalously saline subsurface water mass leads to a strong AMOC strength while fresher water leads to weaker AMOC. The importance of salinity as well as the magnitude and timescale of the AMOC changes shown here are consistent with previous studies applying ECBilt-CLIO coupled model (Jongma et al., 2007) and LOVECLIM (Friedrich et al., 2010), albeit they use slightly different boundary conditions. Furthermore, our simulation shares some important features with Friedrich et al. (2010) regarding the behavior of AMOC oscillations and the mechanisms involved in it. Among them is the AMOC preferentially operating at weak state toward 115 ka when the obliquity reduces to 22.4 and its related increase of large abrupt AMOC change occurrences during the second half of the LIG, but also changes of the Greenland-Iceland-Norwegian (GIN) sea overturning circulations (not shown here) ahead of the total AMOC changes. However, our simulation do not show any deep decoupling from subsurface temperature increase responsible for the AMOC increase, as shown in Friedrich et al. (2010), but rather atmospheric temperature anomalies above the sea level (up to 10 K, not shown here) inducing wind stress anomalies favorable to transporting anomalously saline water south of Greenland as mentioned in Schulz et al. (2007). Here, the origin of density changes and the subsequent activation of deep convection comes rather from a surface salinity change. Hence, the salinity remains to be one of the key factor in triggering the deep convection areas that switches the AMOC from weak to stronger states; the mechanisms behind their origin remain however partially unclear. Further sensitivity experiments are needed in order to fully diagnose the mechanisms responsible for the fresh versus saline water shifts south of Greenland.

Conclusion
Recent high-resolution proxy records reveal that abrupt changes in bottom water chemistry (δ 13 C) occurred on centennial timescales during the previous interglacial and are postulated to represent ocean circulation changes, thus challenging the paradigm of interglacial thermohaline circulation stability at short timescale (Galaasen et al., 2014;Hodell et al., 2009). In this study, we use the iLOVECLIM Earth system model of intermediate complexity to simulate a transient experiment of the LIG period from 125 to 115 ka to analyze possible mechanisms behind these short-lived disturbances depicted by the δ 13 C data reconstruction from that period.
Our simulation depicts large perturbations of bottom water δ 13 C DIC up to 0.43%0 and 0.65%0 in the North Atlantic at the locations corresponding to the sediment core sites MD03-2664 (Eirik Drift) and U1304 Site (Gardar Drift), respectively. These variations are strongly similar in magnitude (0.5-0.7%0) and temporal characteristics (multidecadal scale at onset and persisting for centuries) to the variability in high-resolution δ 13 C reconstruction. The simulated δ 13 C DIC at both sites generally differ by about 0.4%0, similar to the δ 13 C gradient observed in the subpolar North Atlantic during the late LIG. In the model, this gradient is attributed to the relative influence of northern source water (NADW) versus SSW to the sediment core locations. However, during the centennial scales perturbations, the simulated bottom water δ 13 C DIC values converge toward the same values, suggesting that similar water mass properties could potentially be affecting both sites, which is confirmed by our PO analysis.
The simulated variations in δ 13 C DIC at MD03-2664 and U1304 Site are correlated to the variations in AMOC strength, which ranges from 6.3 to 19.3 Sv. High values of bottom δ 13 C DIC correspond to strong AMOC and vice versa. Consequently, a simulated strong AMOC state leads to similar bottom water δ 13 C DIC values at both locations (i.e., affected by the same water mass-NADW), while when the AMOC decreases, the influence of SSW at U1304 Site increases. These large AMOC differences are linked to the activation and deactivation of convection areas south of Greenland, highlighting the importance of NADW production in this region and its role in redistributing δ 13 C DIC in the interior ocean.
In our model simulation, the mechanism associated to the activation and deactivation of the convection south of Greenland and therefore the AMOC strength is related to a weakening of the vertical density gradient mainly induced by the surface salinity in this region, which increases and decreases along with the sea-ice and associated freshwater extension. Northward intrusions of warm and saline Atlantic water into the Nordic Seas also gradually expand westward from 25 • W toward the Greenland coast. In turn, the sea-ice and associated freshwater surface layer retreats (up to ∼10 • northward) and the advected saline surface water increases the surface density. This sets favorable conditions for deep convection to occur, initiating the AMOC increase.
While the deep Atlantic δ 13 C DIC changes seen in our model simulation and the proxy data are similar in magnitude and temporal characteristics, several discrepancies and differences can be mentioned. The use as initial conditions of the equilibrium at 125 ka state sets different climate background than it may have been in reality. Consequently, (1) the absolute values of bottom water δ 13 C DIC are generally simulated higher than that depicted by the data reconstructions, and (2) the simulation does not reproduce the abrupt variations during the early LIG, potentially due to the lack effects from the final retreat of glacial ice sheet remnants (Nicholl et al., 2012). In addition, our model simulates similar decreasing long-term trends at both sediment core locations that do not correspond to what is observed in the data records depicting different trends at both locations. This is attributed to the simulated strong decline of the land vegetation carbon content (∼420 PgC) towards 115 ka and the associated decrease in atmospheric δ 13 C. Hence, terrestrial vegetation content seems to play an important role in setting the ocean δ 13 C DIC absolute values, which has also been highlighted for glacial periods (Menviel et al., 2017). Finally, the abrupt transitions occur when the model switches from a strong AMOC state (early LIG) to a weak state (late LIG) and are likely associated to the cooling (i.e., sea-ice/freshwater forcings) towards the glacial inception at 115 ka, which differs from future projected climate. Additional experiments need to be performed in order to properly address the internal mechanisms shaping the frequency and occurrence of these abrupt transitions. In the model, the typical changes in atmospheric CO 2 associated with these abrupt transitions are approximately of 7 ppm, which is in the range of the CO 2 variation observed in the data reconstruction from that period (Lourantou et al., 2010;Schneider et al., 2013), while the changes in ocean DIC are of about 15 PgC and compensated by a 30 PgC increase in the land vegetation carbon.
This study underlines that during the LIG the Labrador Sea and Irminger Sea may have been sensitive to density changes altering the strength of the AMOC over several decades by simple natural greenhouse gases and orbital forcings, hence affecting the distribution of bottom water δ 13 C DIC . These unforced oscillations may be specific to our simulation and do not exclude other potential mechanism that may be responsible for abrupt and large changes of bottom water δ 13 C depicted in some marine records. The stronger amplitude of change in δ 13 C DIC at U1304 Site as compared to MD03-2664 potentially reveals this location at a better indicator for AMOC structure and/or strength change. In addition, we highlight the need for more isotope-enabled model simulations over interglacial periods in order to better constrain and link changes observed in the marine records and potential changes in overturning circulation.