Early deglacial Atlantic overturning decline and its role in atmospheric CO 2 rise inferred from carbon isotopes ( δ 13 C )

The reason for the initial rise in atmospheric CO2 during the last deglaciation remains unknown. Most recent hypotheses invoke Southern Hemisphere processes such as shifts in midlatitude westerly winds. Coeval changes in the Atlantic meridional overturning circulation (AMOC) are poorly quantified, and their relation to the CO2 increase is not understood. Here we compare simulations from a global, coupled climate–biogeochemistry model that includes a detailed representation of stable carbon isotopes (δC) with a synthesis of high-resolution δC reconstructions from deepsea sediments and ice core data. In response to a prolonged AMOC shutdown initialized from a preindustrial state, modeled δC of dissolved inorganic carbon (δCDIC) decreases in most of the surface ocean and the subsurface Atlantic, with largest amplitudes (more than 1.5 ‰) in the intermediatedepth North Atlantic. It increases in the intermediate and abyssal South Atlantic, as well as in the subsurface Southern, Indian, and Pacific oceans. The modeled pattern is similar and highly correlated with the available foraminiferal δC reconstructions spanning from the late Last Glacial Maximum (LGM, ∼ 19.5–18.5 ka BP) to the late Heinrich stadial event 1 (HS1, ∼ 16.5–15.5 ka BP), but the model overestimates δCDIC reductions in the North Atlantic. Possible reasons for the model–sediment-data differences are discussed. Changes in remineralized δCDIC dominate the total δCDIC variations in the model but preformed contributions are not negligible. Simulated changes in atmospheric CO2 and its isotopic composition (δCCO2 ) agree well with ice core data. Modeled effects of AMOC-induced wind changes on the carbon and isotope cycles are small, suggesting that Southern Hemisphere westerly wind effects may have been less important for the global carbon cycle response during HS1 than previously thought. Our results indicate that during the early deglaciation the AMOC decreased for several thousand years. We propose that the observed early deglacial rise in atmospheric CO2 and the decrease in δ CCO2 may have been dominated by an AMOC-induced decline of the ocean’s biologically sequestered carbon storage.


Introduction
Earth's transition from the LGM (Last Glacial Maximum) (23-19 ka BP), into the modern warm period of the Holocene (10-0 ka BP) remains enigmatic (Denton et al., 2006).Evidence of an early warming of the Southern Hemisphere and atmospheric CO 2 increase (Petit et al., 1999;Denton et al., 2010) has prompted hypotheses of a Southern Hemisphere trigger for the deglaciation (Stott et al., 2007;Timmermann et al., 2009).But the early rise in atmospheric CO 2 , although an important forcing for deglacial global warming (Shakun et al., 2012), remains unexplained.Various mechanisms have been proposed.Prominent recent studies suggest wind (Anderson et al., 2009;Denton et al., 2010;Toggweiler et al., 2006) and/or stratification (Watson and Naveira Garrabato 2006;Schmittner et al. 2007;Sigman et al. 2007;Tschumi et al., 2011) changes in the Southern Ocean and/or changes in the North Pacific circulation (Menviel et al., 2014).
Others have suggested that the deglaciation was initiated by a collapse of the AMOC (Atlantic meridional overturning circulation) caused by the melting of Northern Hemisphere ice sheets (Clark et al., 2004;Sigman et al. 2007;Anderson et al., 2009;Denton et al., 2010;Shakun et al., 2012;He et al., 2013)

and abrupt North Atlantic climate
Published by Copernicus Publications on behalf of the European Geosciences Union.
changes (Broecker et al., 1985).This idea is appealing since the AMOC is known from theory to exhibit multiple steady states with the possibility of rapid transitions between them (Stommel, 1961).Moreover, AMOC variations are consistent with the observed antiphasing of surface temperatures between the hemispheres (Crowley 1992;Blunier et al., 1998;Schmittner et al., 2003;Stocker and Johnson, 2003;EPICA community members 2006;Shakun et al., 2012) and evidence for ITCZ (intertropical convergence zone) migration (Menviel et al., 2008).However, surface temperatures and tropical rainfall patterns alone do not allow robust inferences about the AMOC (Kurahashi-Nakamura et al., 2014) and evidence from the deep ocean for circulation variations remains sparse.One widely cited record of protactinium : thorium ratios ( 231 Pa / 230 Th) from the subtropical North Atlantic has been interpreted as an AMOC collapse around 19-18 ka BP followed by a rapid resumption ∼ 15 ka BP in the warm Bølling-Allerød period (McManus et al., 2004).However, this interpretation has been questioned (Keigwin and Boyle, 2008) and a subsequent set of 231 Pa/ 320 Th records (Gherardi et al., 2009) suggested that a complete AMOC cessation during HS1 (Heinrich stadial event 1) was unlikely.Moreover, our understanding of 231 Pa / 230 Th in the modern ocean continues to evolve (Anderson and Hayes, 2013) and inferences on the basin or global scale circulation from a single site require validation with multiple proxies from a range of oceanographic locations.A quantitative deglacial AMOC reconstruction constrained by distributed interior ocean observations continues to be lacking.Here we attempt a first step towards such a reconstruction by combining model simulations with δ 13 C measurements of sediment samples.
Deep-sea reconstructions based on δ 13 C are more common than 231 Pa / 230 Th, the processes governing δ 13 C are better understood, and realistic three-dimensional models exist (e.g., Schmittner et al., 2013), providing necessary ingredients for quantitative hypothesis testing.Here we compile deep-ocean δ 13 C reconstructions at a high temporal resolution from the early deglaciation and compare them with model simulations of δ 13 C changes caused by AMOC variations in order to test the hypothesis that the AMOC was reduced during HS1.We also compare our model results to observations of atmospheric CO 2 concentrations and the δ 13 C of atmospheric CO 2 in order to assess mechanisms of the early deglacial CO 2 rise.Here we do not address the full deglaciation but restrict our investigation to its initial phase from the late LGM (∼19.5-18.5 ka BP) to the late HS1 (∼16.5-15.5 ka BP).
Various modeling studies have previously examined the effect of AMOC changes on atmospheric CO 2 , with sometimes conflicting results (Marchal et al., 1998;Marchal et al., 1999;Scholze et al., 2003;Köhler et al., 2005;Schmittner et al., 2007a;Obata, 2007;Schmittner and Galbraith, 2008;Menviel et al., 2008Menviel et al., , 2012Menviel et al., , 2014;;Bozbiyik et al., 2011).Schmittner and Galbraith (2008) found that a large AMOC reduction decreases the efficiency of the ocean's biological pump if North Atlantic Deep Water (NADW) is more depleted in preformed nutrients than water masses sourced in the south (Antarctic Bottom Water (AABW) and Antarctic Intermediate Water (AAIW)), and it thus leads to the outgassing of CO 2 into the atmosphere and gradually increasing atmospheric CO 2 over several thousand years, consistent with theory (Ito and Follows 2005;Marinov et al., 2008a, b) and ice core CO 2 reconstructions (Ahn and Brook, 2007).Some of the differences in model responses may thus be due to the simulations of preformed nutrients.Whereas Schmittner and Galbraith (2008) have demonstrated consistency of their model with modern preformed nutrient observations, such a validation is not published, to our knowledge, for other models (e.g., the LOVECLIM model used by Menviel et al., 2008Menviel et al., , 2014)).Several studies found a dependency of the results on the initial state (Köhler et al., 2005;Schmittner et al., 2007a;Menviel et al., 2008), suggesting that starting from glacial conditions may give a different answer than starting from modern conditions.However, none of these studies have validated their initial deep-ocean LGM states with reconstructions.Thus inferences from these studies regarding the sensitivity of the real ocean carbon cycle to initial conditions remain subject to considerable uncertainty.
Effects of Southern Hemisphere westerly winds on atmospheric CO 2 have also been quantified with models before (Winguth et al., 1999;Menviel et al., 2008;Tschumi et al., 2008Tschumi et al., , 2011;;d'Orgeville et al., 2010;Lee et al., 2011;Völker and Köhler, 2013).Most of these studies conclude that reasonably expected changes in the strength and/or latitude of westerly winds cannot explain a large fraction of the observed 100 ppm glacial-interglacial CO 2 amplitude (Winguth et al., 1999, Menviel et al., 2008, Tschumi et al., 2008, 2011, d'Orgeville et al., 2010, Völker and Köhler, 2013).An exception is the work by Lee et al. (2011), who find a 20-60 ppm CO 2 increase for a 25% strengthening of the westerly winds.However, their wind stress forcing was calculated from an atmosphere only model, which was forced with a very large heat flux anomaly and leads to large areas of the subtropical North Atlantic experiencing extreme cooling of more than 10 • C, much more than reconstructed (Bard et al., 2000).We will show below that more realistic coupled ocean-atmosphere model simulations of an AMOC collapse result in much smaller wind stress changes.For comparison, a complete removal of the Antarctic Ice Sheet leads only to a 50 % reduction in Southern Hemisphere westerly winds (Schmittner et al., 2011).Tschumi et al. (2011), whose simulations include δ 13 C and are forced with wind stress changes over the Southern Ocean, conclude that stratification changes there can explain the observed rise in atmospheric CO 2 and the decrease in δ 13 C CO2 during HS1.Here we propose a different mechanism, namely changes in the AMOC and its effect on the efficiency of the biological pump, as an alternative hypothesis for the ice core observations.Another difference from previous studies is that we directly and quantitatively test Anderson et al.'s (2009) hypothesis that AMOC changes affect atmospheric circulation and wind stress in the Southern Ocean to such a degree that the outgassing of CO 2 contributes importantly to the total CO 2 rise during HS1.We do this by applying realistic wind stress changes from a coupled ocean-atmosphere model simulation of an AMOC shutdown to a global carbon cycle model including isotopes.

Methods
We have compiled 25 published deep-ocean records covering the early deglaciation at a high temporal resolution (Table 1).Mostly published age models are used, except in some cases where the radiocarbon calibration was updated as described in Lund et al. (2014).In order to be consistent with the treatment of the other cores in Lund et al. (2014), we have updated the age model of MD01-2461 by recalibrating the radiocarbon ages using INTCAL13 and reservoir ages estimated by Stern and Lisiecki (2013).The ages may have considerable (O(1 ka)) uncertainties.However, we believe that the records are of sufficient resolution and their age models are well enough constrained to evaluate multimillennial changes.The purpose of this paper is to present an initial comparison to model results focusing on model analysis.The quantification of age uncertainties and their effects on the results are beyond the scope of this paper.The sediment data compilation is available in the supplement to this paper.
We employ the Model of Ocean Biogeochemistry and Isotopes (MOBI 1.4), a coupled climate-biogeochemical system that includes δ 13 C cycling in the three-dimensional ocean, land, and atmosphere to explore the effect of AMOC variations on carbon isotopes (see the Appendix for a more detailed model description).MOBI's large-scale ocean distribution of δ 13 C DIC in dissolved inorganic carbon (DIC) is consistent with modern water column observations (Schmittner et al., 2013).It is embedded in the University of Victoria climate model of intermediate complexity version 2.9 and run to a preindustrial equilibrium with prognostic atmospheric CO 2 and δ 13 C CO 2 .Subsequently four numerical experiments, each ∼ 3500 years long, were conducted with varying amplitudes (0.05, 0.1, 0.15, and 0.2 Sv; referred to as FW0.05, FW0.1, FW0.15, and FW0.2, respectively; Sv stands for Sverdrup, i.e. 10 6 m 3 s −1 ) of a stepwise, 400year-long freshwater input to the North Atlantic between 45-65 • N and 60-0 • W (Fig. 1a).The added freshwater is not compensated for elsewhere, but it affects surface tracer concentrations though dilution.Lower salinity and increased buoyancy of surface waters causes the AMOC to slow down.Note that these are idealized experiments, designed to examine only how AMOC variations impact the global δ 13 C distribution.Realistic initial conditions for the LGM are currently not available.Thus, we do not attempt realistic deglacial simulations.However, it is well known that the δ 13 C distribution of the LGM ocean (Curry and Oppo, 2005;Gebbie, 2014)  and atmospheric CO 2 concentrations (Monnin et al., 2001;Lourantou et al., 2010;Parrenin et al., 2013;Marcott et al., 2014) were different from the preindustrial.In order to account for these differences in initial conditions, our datamodel comparison focuses on anomalies rather than absolute values.Possible sensitivity of the results to initial conditions is further discussed in the Discussion section below.Selected model data are available in the supplement to this paper.

Simulated circulation changes
The AMOC reduces in all experiments (Fig. 1b).However, in FW0.05 and FW0.1 the reduction is reversible, and, after hosing is stopped, the AMOC quickly returns to its initial state.In experiments FW0.15 and FW0.2, on the other hand, the AMOC collapses permanently (Fig. 2).Reduced salt input to the deep ocean by North Atlantic Deep Water (NADW) leads to a freshening of the deep ocean and salin-ification of the surface (not shown here; see Figs. 5 and 6 and Plate 2 in Schmittner et al., 2007a), deeper mixed layers, and decreased stratification in the Southern Ocean (Schmittner et al., 2007a) and North Pacific (Saenko et al., 2005), as illustrated in Fig. 2 by the thickening of isopycnal layers between 26.8 ≤ σ ≤ 27.63.

Simulated carbon cycle changes
The effect on atmospheric CO 2 in model simulations with partial and short AMOC reductions (FW0.05 and FW0.1) is negligible.In contrast, in the simulations with a large and prolonged AMOC decline (FW0.15 and FW0.2) CO 2 starts to rise about 500 years after the beginning of the hosing.It continues to increase gradually by ∼ 25 ppm until year 2000, after which its rate of change slows.The amplitude and rate of change of the simulated CO 2 increase agrees well with the long-term trend of measurements of HS1 air recovered from Antarctic ice (Monin et al., 2011;Parrenin et al., 2013;Marcott et al., 2014), but the model does not reproduce the rapid increase around 16 250 BP.
The simulated atmospheric CO 2 increase in FW0.15 and FW0.2 is due to a loss of ocean carbon to the atmosphere.Initially net primary production (NPP) declines within a few hundred years from 64 to 54 Pg C yr −1 , consistent with Schmittner (2005; not shown here), which reduces the production of dissolved organic carbon (DOC), whereas dissolved inorganic carbon increases initially until around year 600, after which it starts to decline.By model year 3500 the ocean has lost ∼ 120 Pg C (Fig. 3d) in FW0.15, most of which (∼ 90 Pg C) is due to DIC and less (∼ 30 Pg C) due to DOC.The ocean's DIC loss is caused by a reduced efficiency of the biological pump as indicated by the large loss of remineralized DIC (∼ 400 Pg C; Fig. 3g), most of which is due to less organic matter oxidation (DIC org ; Fig. 4), whereas it is buffered by the increase in preformed DIC due to rising surface ocean DIC and atmospheric CO 2 , consistent with previous results and theory (Schmittner and Galbraith, 2008;Ito and Follows, 2005;Marinov et al., 2008a, b).

Simulated carbon isotope changes
The loss of biologically sequestered, isotopically light (δ 13 C org = −20 ‰) organic carbon increases deep-ocean δ 13 C DIC by ∼ 0.06 ‰ (Fig. 3f).Accumulation of this isotopically light carbon in the surface ocean and atmosphere decreases their δ 13 C DIC (by ∼ −0.3 ‰) and δ 13 C CO 2 by ∼ −0.25 ‰, respectively (Fig. 1d).Modeled land carbon storage increases (Fig. 3a) and its average δ 13 C L decreases (Fig. 3c), implying that land cannot be the cause of the atmospheric changes.The simulated atmospheric δ 13 C CO 2 decline in models FW0.15 and FW0.2 is consistent, both in amplitude and the rate of change, with ice core measurements (Fig. 1d; Schmitt et al., 2012), but the model response may depend on boundary and/or initial conditions and this agreement may be fortuitous.
The simulated preindustrial (model year 0; Fig. 5a-c) distribution of δ 13 C DIC in the ocean is characterized by high values in the surface and deep North Atlantic and low values in the deep North Pacific, consistent with a previous model version and observations (Schmittner et al., 2013).Sinking of isotopically well-equilibrated surface waters with low nutrient and respired carbon content causes high δ 13 C DIC values in the deep North Atlantic, whereas aging and accumulation of isotopically light respired organic matter progressively decreases δ 13 C DIC as deep waters flow into the South Atlantic and further into the Indian and Pacific oceans.Thus, the modern interbasin difference in deep water δ 13 C DIC is caused by the interbasin meridional overturning circulation (MOC) (Boyle and Keigwin, 1982).Hence, as the AMOC collapses, the δ 13 C DIC difference between North Atlantic and North Pacific deep waters is reduced (Fig. 5d-f).
Differences between years 2500 and 0 (Fig. 5g-i) show the largest δ 13 C DIC decreases at intermediate depths (1-2.5 km) in the northern North Atlantic.Anomalies decrease further south, but a pronounced minimum emerges at the depth of NADW (2-3 km, Fig. 2) in the South Atlantic with positive anomalies below, at the depth of Antarctic Bottom Water, and above, at the depth of Antarctic Intermediate Water.South of 40 • S in the Atlantic as well as in the Indian and Pacific oceans, δ 13 C DIC increases everywhere below ∼ 500 m due to reduced export of 13 C-depleted carbon from the photic zone.The weakening of the biological pump causes surface ocean δ 13 C to decrease by 0.2-0.4‰ in the Indian and Pacific basins, possibly explaining the onset of planktonic δ 13 C minima on glacial terminations (Spero and Lea, 2002).The deep-ocean signal dominates the global mean δ 13 C DIC increase of 0.04 ‰ by year 2500 (Fig. 3f).In the North Pacific δ 13 C DIC shows the largest increase at a depth of around 1 km, owing to reduced stratification and intensified intermediate water formation (Saenko et al., 2005), which decreases the amount of respired carbon and increases remineralized δ 13 C (δ 13 C rem ; Fig. 6) there.δ 13 C rem increases in most of the deep Pacific, Indian and Southern oceans due to the loss of respired carbon, whereas, in the Atlantic, respired carbon accumulates, leading to a decrease in δ 13 C rem .Although changes in δ 13 C rem dominate the spatial variations of the total δ 13 C DIC changes, preformed δ 13 C (δ 13 C pre ) variations are not negligible, particularly in the Atlantic, where they decrease by more than 0.3 ‰.

Observed carbon isotope changes during HS1
Observations from the North Atlantic show large δ 13 C DIC decreases early in the deglaciation (Fig. 7a-e; Fig. 8a).
A new data set from the Brazil Margin in the South Atlantic (Figs. 7f-k, 9) (Tessin and Lund, 2013;Lund et al., 2015) shows increasing δ 13 C DIC by 0.2-0.4‰ between 1.1 and 1.3 km depth and decreasing δ 13 C DIC by ∼ 0.5 ‰ between 1.6 and 2.1 km depth, whereas deeper in the water column the data are noisier without a clear trend.Model FW0.15's initial δ 13 C DIC values at the Brazil Margin are higher than the observations mainly for two reasons (Fig. 9).First, the model does not consider the whole ocean lowering of δ 13 C DIC due to the reduction in land carbon during the LGM, and, second, it does not include the shoaling of NADW and very low δ 13 C DIC values in South Atlantic bottom waters (Curry and Oppo, 2005;Gebbie, 2014).Thus the simulated δ 13 C DIC decrease extends deeper than in the observations and shows a substantial reduction below 2.2 km.However, the reconstructed pattern of opposing δ 13 C signal between shallow-intermediate and mid-depths agrees well with the simulated changes due to large AMOC reductions (Fig. 9).The rapid initial increase at intermediate depths appears to be influenced by two factors.First, there is a reduced return flow of low δ 13 C DIC from the Indian ocean (not shown).Second, less upwelling (Schmittner et al., 2005) of low δ 13 C DIC deep water into the upper and surface Southern Ocean leads to a deepening of the high δ 13 C DIC Antarctic Intermediate and Subantarctic Mode waters, which, together with decreased stratification and deeper mixed layers (Schmittner et al., 2007a; Fig. 2), increases δ 13 C DIC by ∼ 0.3 ‰ at 1.2 km depth in all ocean basins at mid-southern latitudes (Fig. 5g-i).The simulated δ 13 C DIC increase at 1.2 km depth in the southwest Pacific (∼ 0.5 ‰) and at 1.6 km depth in the tropical Indian Ocean (∼ 0.3 ‰) agrees well with local reconstructions (Fig. 7p and o), but the simulated changes happen 1.5 ka earlier than in the sediment data.The lack of age model error estimates for the sediment data currently prevents a more detailed assessment of the simulated temporal evolution.In deep waters of the Southern and Indian oceans, the reconstructions are noisy and no clear trend can be identified (Fig. 7l-n).Core W8709A-13PC from the deep eastern North Pacific (Lund et al., 2011) shows no trend in contrast to models FW0.15 and FW0.2.More data from the deep North Pacific are needed in order to better assess model simulations there.

Sensitivity to wind changes
The model results discussed above did not include the effects of wind changes.Winds enter the UVic (University of Victoria) model in three ways: 1. moisture advection velocities u q determine convergence of specific humidity and thus precipitation; 2. wind stress τ supplies momentum to the surface ocean and sea ice; 3. wind speed u modulates the air-sea exchange of heat, water, and gases (CO 2 , O 2 ).
Figure 11 shows the annual mean fields derived from the National Center for Environmental Prediction (NCEP) reanalysis (Kalnay et al., 1996) used in the above runs.
In order to test the sensitivity of our results to these variables, we performed three additional simulations, in which we use anomalies calculated from hosing experiments with the Oregon State University/University of Victorial (OSU-Vic) model (bottom panels in Fig. 11).The OSUVic model includes a fully coupled dynamical atmosphere at T42 resolution (Schmittner et al., 2011), whereas the other components are identical to the UVic model version 2.8 without dynamic vegetation.In response to an AMOC shutdown, OS-UVic simulates a large anticyclonic anomaly over the North Atlantic, a cyclonic anomaly over the North Pacific, a southward shift of the ITCZ particularly over the Atlantic, and a southward shift of Southern Hemisphere westerlies consistent with previous studies (Timmermann et al., 2007;Zhang and Delworth, 2005;Schmittner et al., 2007b).Note that the changes in Southern Hemisphere westerlies are generally less than 10 % of the absolute values of the control simulation and thus much smaller than those used by Lee et al. (2011).
The OSUVic wind anomalies are applied at model year 400 of the FW0.15 simulation (blue dashed line in panel A of Fig. 12).The wind changes have only a modest impact on simulated carbon cycle and isotope distributions (Fig. 12).The largest effect is due to changes in moisture advection velocities, which lead to a rapid decrease in vegetation and soil carbon around year 400.This causes a rapid CO 2 increase by a few parts per million and a rapid decrease of δ 13 C CO 2 by a few hundredths of a per mil.It also delays the oceanic carbon loss by a few hundred years.However, the multimillennial response and our conclusions are not impacted much by the wind changes.

Discussion
Taken together, the changes in the sedimentary deep ocean δ 13 C DIC reconstructions from the LGM to the late HS1 are most consistent with simulations of a severe and prolonged, multimillennial AMOC reduction.Model FW0.15 fits the reconstructions best, as indicated by a high correlation coefficient (r FW0.15 = 0.85; Fig. 10; Table 2), a root-mean-squared error (rms FW0.15 = 0.49) just slightly larger than for models FW0.05 and FW0.1, and a standard deviation closest to the observations (rstd FW0.15 = 1.75).However, δ 13 C DIC changes in the North Atlantic are about twice as large in the model than in the reconstructions, which causes the standard deviation in model FW0.15 to be 75 % larger than that of the observations.One reason for this discrepancy may be that AMOC changes during HS1 were smaller than those simulated here (Gherardi et al., 2009;Lund et al., 2015).A second reason could be the mismatch in initial conditions.If the LGM AMOC was weaker and shallower than in the model's preindustrial simulation, as indicated by a number of reconstructions (Lynch-Stieglitz et al., 2007;Gebbie, 2014), the model would overestimate changes in volume fluxes and perhaps carbon isotopes even if a complete AMOC collapse did occur during HS1.A third reason may be biases in the foraminifera-based δ 13 C DIC reconstructions, e.g., due to a dependency on carbonate ion concentrations (Spero et al., 1997) or dampened records of the actual δ 13 C DIC changes by smoothing and averaging due to bioturbation and/or age model errors.The former would affect particularly regions with large changes in carbonate ion concentrations such as the North Atlantic in the case of an AMOC collapse and the latter may affect particularly low-resolution sediment cores as indicated by reduced agreement with lower-resolution data from a previous study (Sarnthein et al., 1994) (r FW0.15 = 0.80; rms FW0.15 = 0.60; Fig. 10).Resolving the likelihood of these different possibilities will be an important task for future research.Due to these issues our results can only be regarded as semiquantitative.Qualitatively, they support the interpretation of McManus et al.(2004) of the 231 Pa / 230 Th record, but they cannot rule out the possibility of a reduced but not necessarily collapsed HS1 AMOC based on analyses of Atlantic carbon and oxygen isotope data (Lund et al., 2015;Oppo and Curry, 2015).More work is needed for a truly quantitative reconstruction of early deglacial AMOC changes.
Our simulations suggest that an AMOC decline during HS1 could have caused the observed rise in atmospheric CO 2 and the decrease in δ 13 C CO 2 by modulating the global efficiency of the ocean's biological pump.This is in contrast with ideas that invoke Southern Ocean (Toggweiler et al.,     2006; Anderson et al., 2009;Tschumi et al., 2011;Lee et al., 2011) and/or North Pacific (Menviel et al., 2014) mechanisms for the early deglacial CO 2 rise.However, as discussed above, one possible explanation for the overestimated North Atlantic δ 13 C DIC changes in the model is that the early deglacial AMOC changes were smaller.If this was the case, the model could possibly also overestimate the effects on atmospheric CO 2 and δ 13 C CO 2 and the agreement with the ice core observations could be fortuitous.
A critical test of our hypothesis that the AMOC reduction caused a decrease in the efficiency of the biological pump may come from additional δ 13 C DIC reconstructions from the deep Pacific and Indian oceans, which hold most of the ocean's carbon and where the model predicts δ 13 C DIC to increase but where few sedimentary data are currently available (Figs. 4 and 7).Indeed, our mechanism relies on changes in the inflow of (low preformed nutrient) Atlantic deep water into the Southern, Indian, and Pacific oceans.Currently it is not known if this inflow was weaker during the LGM.Gebbie (2014) suggests a similar AMOC to the modern one.Kwon et al. (2012) suggest an even stronger influence of North Atlantic water in the global deep ocean at the LGM.These findings may indicate that our simulations with modern initial conditions may also be applicable to the early deglacial, but a solid quantitative assessment remains to be performed.Such an assessment requires simulations with more realistic initial conditions, which will be an important task for future work.
Our wind stress experiments show much smaller impacts on the carbon cycle than those caused by the buoyancy forcing and suggest only a minor effect on the overall rise in atmospheric CO 2 during HS1, but they are subject to the same caveats as discussed above with respect to initial conditions.However, changes in tropical winds associated with ITCZ shifts impact the hydrological cycle and terrestrial carbon and cause a jump of CO 2 by a few parts per million (Fig. 12).Although this is much smaller than the 12 ppm jump recently observed around 16 250 years BP by Marcott et al. (2014), it suggests a mechanism that could explain rapid increases in atmospheric CO 2 .
We have not discussed the later parts of the deglacial period such as the Bølling-Allerød (15-13 ka BP), during which the AMOC was presumably reinvigorated (McManus et al. 2004).In our model this would lead to a decrease in atmospheric CO 2 , whereas ice core data show stable concentrations (Monin et al., 2011;Parrenin et al., 2013;Marcott et al. 2014), suggesting that an additional process counteracted the AMOC effect.We speculate that this process may be related to a deepening of the AMOC beyond LGM depths and the erosion of the deep South Atlantic reservoir of respired carbon, consistent with recent reconstructions that show that δ 13 C DIC decreases there later in the deglaciation (Lund et al., 2015).
A comparison of distributed deep-ocean δ 13 C DIC reconstructions with our model simulations suggests that, during HS1, the AMOC was substantially reduced for several thousand years.However, due to remaining model-sediment-data differences and a mismatch in initial conditions, we cannot assess the likelihood of a partial AMOC reduction versus a complete shutdown.
Agreement of simulated atmospheric CO 2 and δ 13 C CO 2 with ice core data, if not fortuitous, supports our hypothesis of an AMOC-induced reduction of the oceans' biological pump during HS1.However, this idea needs further testing with more realistic simulations in the future, improving initial conditions and transient forcing.
AMOC-induced wind changes simulated in a coupled ocean-atmosphere model only have a small impact on the carbon cycle and isotope distributions in our carbon cycle model, suggesting that wind changes were less important than previously thought in controlling atmospheric CO 2 and δ 13 C CO 2 during HS1.However, effects of wind shifts on the hydrological cycle and terrestrial carbon could explain some of the recently observed rapid CO 2 increases (Marcott et al. 2014).

Appendix A: Model description
The University of Victoria Earth System Climate Model (UVic ESCM) (Weaver et al., 2001), is used in version 2.9 (Eby et al., 2009).It consists of a coarse-resolution (1.8 • × 3.6 • , 19 vertical layers) ocean general circulation model coupled to a one-layer atmospheric energy-moisture balance model and a dynamic thermodynamic sea ice model, both at the same horizontal resolution.The model is forced with seasonally varying solar irradiance at the top of the atmosphere, cloud albedo, wind stress, wind speed, and moisture advection velocities.This seasonal forcing does not change between different years.Atmospheric CO 2 and δ 13 C are calculated in a single box assuming rapid mixing.

A1 Description of land carbon isotopes (δ 13 C and δ 14 C) model
The land carbon isotopes model has not been published before.Therefore, we provide a description and evaluation here.It is based on TRIFFID, the Top-down Representation of Interactive Foliage and Flora Including Dymamics dynamic vegetation model by Cox (2001), as modified by Meissner et al. (2003) and Matthews et al. (2004), which solves prognostic equations for total vegetation carbon density C v = 12 C v + 13 C v and fractional coverage v i ∈ (0, 1) of five plant functional types (PFTs; i = 1, . .., 5): where i is net primary production (NPP) and i is litter production, which enters the soil carbon pool.Total soil carbon density is calculated according to We added prognostic equations for the heavy carbon isotopes 13 C and 14 C to both vegetation, and soil, and where fractionation during photosynthesis is indicated by factors and where is the heavy to light isotope ratio of atmospheric CO 2 .
Fractionation factors are different for C 3 and C 4 plants:  , 1988).
No fractionation occurs during litter production or respiration: 1 + β 13 , (A11) The simulated spatial distribution of average δ 13 C (Fig. 13) varies from −13 ‰ in regions dominated by C 4 grasses such as North Africa and Australia to −27 ‰ in most other regions, which are dominated by C 3 plants, due to the differences in fractionation factors for C 3 and C 4 plants used in the model.This distribution is broadly consistent with previous independent estimates (Still and Powell, 2010;Powell et al., 2012).

A2 Description of ocean carbon isotope model
We use the Model of Ocean Biogeochemistry and Isotopes (MOBI) version 1.4.The ocean carbon isotope component is described in detail in Schmittner et al. (2013).Here we only describe differences with respect to that publication.The physical UVic model version was updated to version 2.9 (Schmittner et al. (2013) used 2.8).The ocean ecosystem model has been modified by changing zooplankton grazing, using a slightly different approach to consider iron limitation of phytoplankton growth as described in detail in Keller et al. (2012).This model gives very similar results to model FeL (iron limitation), which uses a simple mask to consider iron limitation of phytoplankton growth, in Schmittner et al. (2013).
The implementation of the carbon isotope equations has been changed from the "alpha" formulation to the "beta" formulation, courtesy of Christopher Somes.In the alpha formulation, the change in the heavy (rare) isotope carbon density 13 C (in mol C m −3 ) of the product (e.g., phytoplankton) of a given process (e.g., photosynthesis), is calculated as the product of the total carbon change ∂C/∂t times the fractionation factor α for that process times the heavy to light isotope ratio of the source (e.g., sea water DIC) R 13 = 13 C/ 12 C.This formulation assumes total carbon, is equal to 12 C, which is a good approximation since 13 C is usually 2 orders of magnitude smaller than 12 C.However, the assumption (Eq.17) can be avoided by using the beta formulation, in which the heavy isotope change is calculated according to where β 13 = α 13 R 13 .
In order to convert isotope ratios to delta values, δ 13 C = (R/R std − 1), (A19) we now use the conventional standard ratio R 13 std = 0.0112372 instead of R 13 std = 1, which was used in Schmittner et al. (2013).For radiocarbon R 14 std = 1.17 × 10 −12 is used.MOBI 1.4 includes dissolved organic carbon (DOC) cycling described in Somes et al. (2014).The close agreement of the preindustrial δ 13 C DIC distributions with model FeL of Schmittner et al. (2013) suggests that none of the changes described above have a major impact on the simulated δ 13 C DIC .
The Supplement related to this article is available online at doi:10.5194/cp-11-135-2015-supplement.

Figure 1 .
Figure 1.Time series of (A) North Atlantic freshwater forcing, (B) AMOC response, (C) atmospheric CO 2 concentrations, and (D) δ 13 C of atmospheric CO 2 (solid, left axis) and global mean surface ocean δ 13 C DIC (dashed, right axis) for four model simulations (color lines).Symbols in (C) and thick black curve (error estimates are indicated by thin lines) in (D) show ice core measurements (Marcott et al., 2014, Schmitt et al., 2012) (bottom and right axes in (C)).

Figure 2 .
Figure 2. Meridional overturning stream function (black contour lines; solid and dashed lines denote negative and positive values, respectively, and clockwise and counterclockwise flow) in the Southern (A, D), Atlantic (B,E), and Indian and Pacific (C, F) oceans at year 0 (averaged from year −1000 to 0) (A, B, C) and 2500 (averaged from year 2000 to 3000) (D, E, F).Three zonally averaged potential density (σ ) isolines (27.63, 27.5, 26.8) are shown as green lines.

Figure 3 .
Figure 3. Simulated changes in global land (left) and ocean (right) carbon inventories (in Pg C; 1 ppm = 2.1 Pg C) and their averaged δ 13 C (in per mil) in model FW0.15.Changes in (A) land carbon C L = C V + C S (black) are due to vegetation C V (green) and soil C S (red) changes.Vegetation is composed of C 3 and C 4 plants (C V = C 3 + C 4 ), but C 4 plants contribute only a small fraction (B) to the total.Changes in C 4 plant biomass (blue line in A) are negligible compared to those in C 3 plants (difference between green and blue lines).Panel (C) shows biomass weighted mean δ 13 C of the land δ 13 C L = ( i C v,i × δ 13 C i + C s × δ 13 C s )/( i C v,i + C s ) (black), vegetation δ 13 C v = ( i C v,i × δ 13 C i )/( i C v,i ) (green), and soil δ 13 C s (red).Ocean carbon changes C O = DIC + DOC + POC Figure 3. Simulated changes in global land (left) and ocean (right) carbon inventories (in Pg C; 1 ppm = 2.1 Pg C) and their averaged δ 13 C (in per mil) in model FW0.15.Changes in (A) land carbon C L = C V + C S (black) are due to vegetation C V (green) and soil C S (red) changes.Vegetation is composed of C 3 and C 4 plants (C V = C 3 + C 4 ), but C 4 plants contribute only a small fraction (B) to the total.Changes in C 4 plant biomass (blue line in A) are negligible compared to those in C 3 plants (difference between green and blue lines).Panel (C) shows biomass weighted mean δ 13 C of the land δ 13 C L = ( i C v,i × δ 13 C i + C s × δ 13 C s )/( i C v,i + C s ) (black), vegetation δ 13 C v = ( i C v,i × δ 13 C i )/( i C v,i ) (green), and soil δ 13 C s (red).Ocean carbon changes C O = DIC + DOC + POC∼0

Figure 4 .
Figure 4. Vertical profiles of globally horizontally averaged ocean DIC (top left) and δ 13 C (bottom left) at years 0 (black) and 2500 (red) of experiment FW0.15.Right panels show changes (year 2500 minus year 0) in DIC = DIC pre + DIC rem and δ 13 C = δ 13 C pre + δ 13 C rem (green) as well as preformed DIC pre and remineralized DIC rem = DIC org + DIC CaCO 3 .See Schmittner et al. (2013) for the calculation of the individual terms.The differences between the blue and green lines are due to changes in preformed DIC and δ 13 C.

Figure 5 .
Figure 5. Zonally averaged distributions of δ 13 C DIC (color shading and black isolines) as a function of latitude and depth simulated by model FW0.15 in the Atlantic (left), Indian (center) and Pacific (right) ocean basins at model years 0 (A-C) and 2500 (D-F; A-F use top color scale), and the difference (G-I; bottom color scale).Symbols in bottom panels denote locations of observations shown in Fig. 7.

Figure 6 .
Figure 6.Impact of AMOC collapse on δ 13 C pre (top) and δ 13 C rem (bottom).Zonally averaged changes between year 2500 and year 0 of model run FW0.15 in the Atlantic (left), Indian Ocean (center), and Pacific (right).Note the different color scales and isoline differences used.For absolute values of δ 13 C pre and δ 13 C rem , see Figs. 6 and 12 in Schmittner et al. (2013).

Figure 7 .
Figure 7.Comparison of simulated (lines as in Fig.1; left and top axes) and observed (symbols as in Fig.5; right and bottom axes) δ 13 C DIC time series in the North Atlantic (A-E), South Atlantic (F-L), Indian (M-O), and Pacific (P) oceans.If no numbers are given on the right axis, the scale is identical to the left axis.If numbers are given on the right axis, the scale is different but the range (max-min) is identical to that of the left axis.Note that different ranges of the vertical axis are used for the different columns, whereas, within each column they are the same.

Figure 8 .
Figure 8. Heinrich Stadial 1 (16.5-15.5 ka BP) minus LGM (19.5-18.5 ka BP) difference in δ 13 C DIC in the Atlantic (left) and Indian and Pacific (right) basins from our high-resolution synthesis of reconstructions averaged on the model grid (top) compared to model FW0.15 results (bottom; averages of model years 2000 to 3000 minus averages of model years −1000 to 0.).

Figure 9 .
Figure 9. Simulated (solid; model FW0.15) and observed (dashed) vertical profiles of δ 13 C DIC at the Brazil Margin in the South Atlantic before (black) and after (red) the AMOC collapse.Observations show 1 ka averages of smoothed (2 ka) data.Results for model FW0.2 are very similar to FW0.15 (not shown).However, models FW0.05 and FW0.1 show almost no changes from their initial (year −500) distribution (not shown).

Figure 10 .
Figure 10.HS1 minus LGM change in δ 13 C from Sarnthein et al. (1994; blue) our high-resolution compilation (red) vs. changes between years −500 and 2500 (1000-year centered averages) predicted by model experiment FW0.15 at the same locations.The diagonal 1 : 1 line corresponds to a perfect match.

Figure 12 .
Figure12.Sensitivity to changes in winds.Experiment FW0.15 (red) is repeated with changes in moisture advection velocities u q (light blue), u q plus wind stress τ (green), and u q + τ plus wind speed u (dark blue) calculated from the OSUVic model.See Fig.11and text for more details.

Figure 13 .
Figure 13.Simulated average preindustrial land δ 13 C distribution (model year 0).Each pool's (five vegetation plant functional types, PFTs, and one soil, S, carbon compartment) δ 13 C value was weighted by its mass in calculating the average as explained in the caption for Fig. 2. Desert regions with negligible vegetation carbon (< 10 g m −2 ) are shown in white.

Table 1 .
Sediment cores used in this study.Note that cores 1 and 2 have been averaged for the comparison with the model simulations shown in Figs.7-10 and Table2.