Investigation of the Natural Carbon Cycle since 6000 BC using an Intermediate Complexity Model: The Role of Southern Ocean Ventilation and Marine Ice Shelves

Abstract

The mechanisms behind the 20 ppm pre-industrial rise in atmospheric CO₂ since 6000 BC have been the focus of considerable debate in recent years. Some studies suggest that natural processes, such as a decline in global forests, calcite compensation, and warming ocean temperatures, can explain the increase in CO₂. Others have argued that, because the CO₂ increase did not occur during previous interglacial periods, it is an indication of an early human influence on the climate. In this paper, we investigate several facets of the natural carbon cycle of the past 8000 years related to ocean circulation patterns and ice shelf configuration using the University of Victoria Earth System Climate Model (v. 2.9), which includes a representation of the climate system with dynamic vegetation and an interactive carbon cycle. The fully simulated earth system for various freely evolving atmospheric carbon scenarios since 6000 BC failed to recreate the observed rise in CO₂ and consistently produced a decline in CO₂ throughout the Holocene, in keeping with projections related to previous interglacial periods. However, the extent and timing of the decline was strongly dependent on the initial state of the ocean's meridional overturning circulation and the location of marine ice shelves off the coast of Antarctica. For simulations with ice shelves, carbon-poor North Atlantic Deep Water (NADW) dominated the deep Atlantic basin, with less production of (relatively) carbon-rich deep water from the Southern Ocean. This setup led to lower net ocean carbon storage and thus elevated atmospheric CO₂ levels. These deep water distributions, which are relatively independent of orbital forcing and CO₂ concentrations but strongly dependent on Antarctic marine ice shelves, suggest that greater ice shelf extent during the Holocene relative to previous interglacial periods (Pollard & DeConto, 2009) may have contributed 5–6 ppm CO₂ to the atmosphere. The results from this study also indicate that multi-centennial scale changes in the meridional overturning circulation on the order of 2–3 Sv may lead to quasi-periodic increased atmospheric CO₂ concentrations (approximately 6 ppm) through enhanced ventilation of Pacific deep waters.

RÉSUMÉ [Traduit par la redaction] Depuis quelques années, les mécanismes qui sont intervenus dans l'augmentation du CO₂ atmosphérique de l'ordre de 20 ppm, survenue à l'ère préindustrielle font l'objet d'un vaste débat. Selon certaines études, ce phénomène serait attribuable à des processus naturels, comme le recul des forêts dans le monde, la compensation du calcite et le réchauffement de la température des océans. D'après d'autres études, en l'absence d'augmentation de la concentration de CO₂ pendant les périodes interglaciaires précédentes, il marquerait le début de l'influence des activités humaines sur le climat. Nous explorons ici plusieurs facettes du cycle naturel du carbone des 8000 dernières années dans la perspective de la configuration de la circulation océanique et de celle des plates-formes de glace flottante au moyen du modèle du système climatique de la Terre (version 2.9), qui comprend une représentation des systèmes climatiques avec une végétation dynamique et un cycle du carbone interactif. Le système terrestre entièrement simulé pour différents scénarios de carbone atmosphérique à évolution libre depuis 6000 ans avant Jésus-Christ n'a pas permis de recréer l'augmentation de la concentration de CO₂ observée; il a plutôt produit systématiquement une diminution du CO₂ pendant l'Holocène, ce qui va dans le sens des projections se rapportant aux périodes interglaciaires précédentes. Cependant, l'ampleur de la baisse et le moment où le phénomène s'est produit dépendaient fortement de l'état initial de la circulation méridienne de retournement de l'océan ainsi que de l'emplacement des plates-formes de glace flottante au large de la côte de l'Antarctique. Dans le cas des simulations avec des plates-formes de glace flottante, l'Eau profonde de l'Atlantique Nord (NADW) pauvre en carbone prédominait dans les eaux des grands fonds du bassin de l'Atlantique, la production d'eaux des grands fonds (relativement) riches en carbone en provenance de l'océan Austral étant réduite. Dans le cas de cette configuration, le stockage net de carbone dans les océans diminuait, ce qui augmentait la concentration de CO₂ dans l'atmosphère. D'après la répartition des eaux profondes, qui est relativement indépendante du forçage orbital et des concentrations de CO₂, mais qui dépend énormément des plates-formes de glace flottante dans l'Antarctique, le fait que les plates-formes de glace flottante étaient plus étendues pendant l'Holocène que durant les périodes interglaciaires précédentes (Pollard & DeConto, 2009) expliquerait peut-être un apport de 5 ou 6 ppm de CO₂ dans l'atmosphère. L’étude nous apprend également que des changements à l’échelle pluricentennale dans la circulation méridienne de retournement de l'ordre de 2 à 3 Sv peuvent se traduire par des augmentations presque systématiques des concentrations de CO₂ atmosphérique (de l'ordre d'environ 6 ppm) grâce à une meilleure ventilation des eaux profondes du Pacifique.

Keywords:

• carbon dioxide
• Holocene
• Anthropocene
• early anthropogenic hypothesis
• pre-industrial agriculture
• land-use change
• global carbon cycle

1 Introduction: a review of the Holocene carbon cycle debate

In order to understand the response of the Earth's climate to past, present, and future atmospheric carbon dioxide (CO₂) and methane (CH₄) concentrations better, earth system scientists have made a concerted effort to quantify and model the global carbon cycle. Given the complex interplay of short and long time-scale processes and feedbacks between carbon and climate dynamics, it is important to interpret major influences on the carbon cycle adequately in order to evaluate present and future trends. Fortunately, Antarctic ice core records, such as those provided by the Vostok (Petit et al., 1999) and the European Project for Ice Coring in Antarctica (EPICA) Dome Concordia (C) ice cores (Loulergue et al., 2008; Lüthi et al., 2008), have made such an inquiry possible and have been the launching point for recent paleoclimatic studies of atmospheric CO₂ and CH₄ concentrations for the past 800,000 years. The objective of this paper is to evaluate natural controls on the carbon cycle since 6000 BC using an intermediate complexity model to help interpret the CO₂ trends in the ice core data.

On shorter time scales, the Taylor Dome Antarctic ice core provides a high-resolution reconstruction of CO₂ levels for the past 11,000 years (Indermühle et al., 1999), the period encompassing the Holocene. This record demonstrates a decrease in atmospheric CO₂ by 7–8 ppmv from 9000 to 6000 BC, followed by an increase of 20–25 ppmv (henceforth ppm) from 6000 BC to the end of the pre-industrial era (Indermühle et al., 1999; the data from which are represented in the Prescribed Carbon (PC) simulation in Fig. 6a). Several hypotheses have been put forward to explain these features in the CO₂ record. Modelling and data studies (Indermühle et al., 1999; Joos, Gerber, & Prentice, 2004; Yu, 2011) tend to agree that the initial decrease in CO₂ during the early Holocene was a result of vegetation regrowth and peatland uptake over formerly ice-covered areas. However, the mechanism causing the subsequent increase in atmospheric CO₂ from approximately 260 ppm in 6000 BC to the pre-industrial level of approximately 280 ppm is still not well understood (Ruddiman 2008, Ruddiman, Kutzbach, & Vavrus, 2011).

a The Debate over Contributions from Human Land Use

Much scientific debate has focused on the contribution of human land use to the 20 ppm increase in atmospheric CO₂. Ruddiman and Thomson (2001), Ruddiman (2003), and Ruddiman et al. (2011) introduced the idea that changing human subsistence strategies during the Holocene, accompanied by notable increases in population and a transformation of the landscape, have had a significant effect on both CO₂ and CH₄ levels during the middle and late Holocene. Putting together evidence from previous interglacial periods, during which a decrease in CO₂ was observed, Ruddiman (2003) concluded that CO₂ concentrations should have also declined from the early Holocene peak of 268 ppm to a value close to 240 ppm. More recently, Ruddiman et al. (2011) confirmed that the Holocene was unique relative to the previous eight interglacial periods (using both insolation and δ¹⁸O excursion alignment methods) in terms of the substantial, continuous rise in CO₂ starting around 6000 BC.

Table 1. Summary of simulations discussed in this paper. The modern spin-up referenced in the initial conditions is a climate state for AD 1800 used in Eby et al. (2009). “Free” atmospheric CO₂ refers to a model state wherein atmospheric CO₂ responds freely to the climate state and changes in other carbon reservoirs, without being forced to follow the observed value or trend in atmospheric CO₂ from Taylor or Law Dome ice cores.

	Initial Conditions	Atmospheric CO2	Land Ice	Orbital Forcing	Winds
Equilibrium Simulations
EQ_1	Modern spin-up	Constant 6000 BC value	6000 BC configuration	6000 BC configuration	NCEP
EQ_2	Year 6200 in EQ_1	Free	6000 BC configuration	6000 BC configuration	NCEP
EQ_3	Modern spin-up	Constant 6000 BC value	6000 BC configuration	6000 BC configuration	Wind feedback
Transient Simulations
PC	Year 6200 in EQ_1	Prescribed	Prescribed	Prescribed	NCEP
BE_PC	Before event	Prescribed	Prescribed	Prescribed	NCEP
DE_PC	During event	Prescribed	Prescribed	Prescribed	NCEP
FC	Year 6200 in EQ_1	Free	Prescribed	Prescribed	NCEP
BE_FC	Before event	Free	Prescribed	Prescribed	NCEP
DE_FC	During event	Free	Prescribed	Prescribed	NCEP

Using these analogues to eliminate natural variability as an explanation for the CO₂ increase during the Holocene, Ruddiman (2003) and Ruddiman et al. (2011) suggested human activity as the key contributor to the observed CO₂ increase (the “Early Anthropogenic Hypothesis”). To support this conjecture, Ruddiman (2003) cited pollen and archaeological evidence for forest clearing in China and the Fertile Crescent/Mediterranean basin by 6000 BC because 80% of the CO₂ rise was achieved by the beginning of the Common Era (AD 1). Ruddiman (2003, 2007) also pointed to “great” deforestation by that time in southern and eastern Asia and the Mediterranean basin, as discussed in Roman chronicles and seen in pollen and erosion records from formerly forested watersheds.

However, the ability of human agricultural activities to contribute substantially to the rise in CO₂ prior to 2000 years ago seems more problematic given low global populations (5–250 million) during this period (Boyle, Gaillard, Kaplan, & Dearing, 2011). Another major challenge to the early anthropogenic hypothesis has come from the direct mechanistic modelling of human land use and the resulting net carbon release to the atmosphere. Olofsson and Hickler (2008) calculated that the atmospheric CO₂ level would have increased only a few parts per million as a result of a 26 PgC (petagrams (10¹⁵ g) of carbon) net release between 4000 BC and AD 1, not enough to account for the observed increase in CO₂ of nearly 16 ppm during the Holocene to that point. Using the HistorY Database of the global Environment (HYDE) v. 3.0 with pasture and agricultural land use for AD 1700 in the Lund-Potsdam-Jenna Dynamic Global Vegetation Model (LPJ-DGVM), Strassmann, Joos, and Fischer (2008) obtained a carbon release of 45 PgC between their land-use simulation and their spin-up run (potential natural vegetation), determining this to be the maximum pre-industrial release possible (resulting in a CO₂ increase of only a few parts per million when released gradually over the course of millennia). Further, Strassmann et al. (2008) emphasized the importance of enhanced CO₂ fertilization from forests as a negative feedback (a reduction of approximately 25%) to the land-use carbon flux to the atmosphere. Pongratz, Reick, Raddatz, and Claussen (2008) developed another land-use database that took into account increased land-use efficiency with time from AD 800 to AD 1992, the part of the Holocene that witnessed the greatest growth in human population. This dataset was applied to the European Centre Hamburg Model (ECHAM5) atmosphere-ocean general circulation model (AOGCM), and the simulations (Pongratz, Reick, Raddatz, & Claussen, 2009) produced a 5–6 ppm increase between AD 800 and AD 1850. However, an extrapolation back in time suggests that early agriculture would have only contributed about 1 ppm of the observed 16 ppm rise in CO₂ between 6000 BC and AD 1 (Pongratz et al., 2009).

Each of these modelling studies assumes that land-use efficiency before AD 1 remained roughly constant or increased conservatively with time, often projecting more recent data from the Middle Ages and the present day back to 6000 BC (land-use efficiency at 0.16–1.3 ha per person). The Boserup Hypothesis (Boserup 1965), as elaborated in Ruddiman and Ellis (2009), challenges these low estimates of the contribution of human land use to the 20 ppm increase. In particular, low land-use efficiency during the neolithic period (4 ha per capita land use) likely provided a greater carbon flux per capita to the atmosphere than modern agriculture (closer to 0.16 ha per person) as a result of the cruder technologies used in early farming.

In order to address this challenge, land-use scenarios adapted from HYDE 3.1 (Klein Goldewijk, Beusen, de Vos, & van Drecht, 2011) were incorporated mechanistically into the Bern carbon cycle (CC) model (containing the LPJ-DVGM) in Stocker, Strassmann, and Joos (2010). Their results suggest that even extreme terrestrial releases (HYDE 3.1 land use multiplied by a factor of eight) would not have been significant enough to reproduce the Holocene atmospheric CO₂ trend, leading to a contribution of only 6 ppm by AD 1. However, an alternative land use scenario of early Holocene populations summarized in Kaplan et al. (2011) places mid-Holocene land use efficiency between 5.5 and 8 ha per capita (depending on the region) with greater agriculture-related soil carbon depletion than applied in Stocker et al. (2010). This database applied to the LPJ-DVGM suggests much higher anthropogenic emissions (325–357 PgC) than previous estimates (Kaplan et al., 2011). Boyle et al. (2011) proposed a smaller neolithic land-use efficiency nearer to 1.96 ha per person (based on a combination of archaeological evidence and human dietary requirements). The authors further indicated that global population estimates during the neolithic period could have been an order of magnitude greater (near 100 million) than that proposed by McEvedy and Jones (1978) (closer to 10 million) used as the basis for most human land-use databases, including some regions in Kaplan et al. (2011). In short, given the complexity and uncertainties in reconstructing pre-industrial human land use (both in terms of global population and land-use efficiency) and the wide range of emissions estimates (from 48 to 357 PgC), the contribution of agricultural activities to the Holocene carbon cycle remains an open question (Ruddiman et al., 2011).

b Natural Explanations of the Holocene Trend

Various hypotheses involving the natural carbon cycle, however, have also been put forward to explain the pre-industrial Holocene increase in CO₂. Indermühle et al. (1999) suggested a natural vegetation release of carbon associated with the drying of the Sahara and declining solar insolation in the northern hemisphere (see Claussen (2009) for a recent review). In their comprehensive model study of the Holocene, Schurgers et al. (2006) also attributed the CO₂ increase of 10 ppm to a terrestrial release. However, Wang, Mysak, Wang, and Brovkin (2005) found that the flux of carbon to the atmosphere from North Africa was compensated for by the simulated growth of vegetation with increasing solar radiation in the southern hemisphere. Furthermore, several previous modelling studies (see Joos et al., 2004) produced estimates ranging from a 90 PgC release to a 370 PgC uptake, indicating significant uncertainty in the role of natural vegetation.

Addressing the peatland question, Wang, Roulet, Frolking, and Mysak (2009) incorporated the carbon accumulation of northern peatlands into the Green McGill Paleoclimate Model (MPM) and proposed that peatlands alone could account for a natural 160–280 PgC terrestrial sink for the atmosphere since 6050 BC. Since then, Yu, Loisel, Brosseau, Beilman, and Hunt (2010) have provided a more exhaustive compilation of peat basal dates globally and “conservatively” estimated that northern peatlands alone accumulated 547 PgC during the Holocene, with tropical wetlands and Patagonian peatlands taking up 50 PgC and 15 PgC, respectively. Yu (2011) further calculated that a net accumulation of 267 PgC should have led to a 19 ppm decrease in atmospheric CO₂ since approximately 6050 BC. Furthermore, removal of carbon from the atmosphere was likely greater (87.9 Tg yr⁻¹ from approximately 7000 to 6000 BC) during the early and mid-Holocene than during the late Holocene (16.3 TgC yr⁻¹ from approximately 1000 BC to AD 1) because of delayed decomposition of previously accumulated peat, perhaps helping to explain the accelerated decrease in CO₂ during the early Holocene as well as a gradually decreasing impact of peatlands on atmospheric CO₂ over thousands of years. This would thus allow atmospheric CO₂ to rise under the influence of other natural or anthropogenic contributors.

Another hypothesis proposes that an increase in mean oceanic sea surface temperature (SST) during the Holocene would have led to a release of CO₂ to the atmosphere. This interpretation was partly validated by Joos et al. (2004) who imposed a global SST warming of 0.6°C as a boundary condition in their model and obtained a 4–5 ppm increase in CO₂. Ruddiman et al. (2011) also attributes 9 ppm to the temperature–solubility relationship.

However, proxy data based on near-coastal alkenone (Kim & Schneider, 2004) suggest that the North Atlantic Ocean, an important region for transporting atmospheric CO₂ to the deep ocean, cooled with declining summer solar insolation in the late Holocene while other basins warmed. A modelling study using an AOGCM and an accelerated orbital forcing (Lorenz, Kim, Rimbu, Schneider, & Lohmann, 2006) yielded a similar trend in cooling (on the order of 1°–2°C) for the North Atlantic basin. Brovkin, Kim, Hofmann, and Schneider (2008) also suggested that cooling in the critical North Atlantic region may have intensified overturning in this region, leading to faster uptake of CO₂ by the deep ocean in North Atlantic deep water (NADW). This regional cooling (stimulating greater CO₂ uptake) in the model is stronger than the surface warming (associated with CO₂ release) in other basins and so dominates the global carbon cycle leading to a 1–6 ppm global net decrease (sink) of CO₂ rather than the net source postulated in Joos et al. (2004). However, Brovkin et al. (2008) does not address processes in the Southern Ocean, whereas Ruddiman (2007) and Vavrus, Ruddiman, and Kutzbach (2008) argue that feedbacks from mid-Holocene anthropogenic activities may have warmed the waters there and reduced sea ice during the late Holocene, thereby increasing air–sea exchange with the atmosphere leading to a substantial CO₂ net release of 20–24 ppm. Whether anthropogenically or naturally induced, deuterium in ice core records suggests that the Antarctic climate was unusually warm during the Holocene compared with other interglacial periods, and this might have had important feedbacks on overturning in this region (Ruddiman et al., 2011).

Another interpretation of the positive CO₂ anomaly during the Holocene relates to ocean calcite compensation and coral reef growth (Broecker et al., 1993; Broecker, Lynch-Stieglitz, Clark, Hadjas, & Bonani, 2001; Archer & Maier-Reimer, 1994; Broecker & Clark, 2003; Broecker & Stocker, 2006). Because the terrestrial biosphere (expanding forests and peatlands during deglaciation) takes carbon from the atmosphere–ocean reservoir, a deepening of the calcite lysocline should eventually lead to excess calcite deposition, reduced solubility of CO₂ in the oceans, and a corresponding recovery in atmospheric CO₂. Another way of achieving this chemical transformation of the ocean is through upslope tropical coral reef migration in response to postglacial sea level rises (Ridgwell, Watson, Maslin, & Kaplan, 2003).

A number of modelling studies also lend strong support to the ocean chemistry hypothesis (Brovkin et al., 2002; Elsig et al., 2009; Joos et al., 2004; Kleinen, Brovkin, von Bloh, Archer, & Munhoven, 2010; Ridgwell et al., 2003), albeit with a fairly large uncertainty regarding the oceanic carbon release. According to the BERN CC simulations by Joos et al. (2004), calcite compensation may have accounted for a 4–10 ppm ocean release, which explains about 25 to 50% of the rising Holocene CO₂ levels. Similarly, Elsig et al. (2009) proposed that 15 ppm CO₂ came from calcite compensation and 5 ppm from coral reef growth. On the other hand, Kleinen et al. (2010) suggested that calcite compensation has little impact on the Holocene trend, with a pre-industrial atmospheric CO₂ concentration of only 257 ppm when only calcite compensation and atmosphere–vegetation (non-peatland) carbon exchanges are considered in the CLIMate and biosphERe Model-Lund-Potsdam-Jenna (CLIMBER-LPJ) model. However, even with peatlands (a net 105 PgC uptake in their simulation), Kleinen et al. (2010) was able to approximate the Holocene carbon trend (to 278 ppm) with the inclusion of shallow water sediments (an approximation of coral reefs). However, Kleinen et al. (2010) also noted that the reduction in ocean alkalinity that produced this increase in CO₂ was much greater than that calculated by Broecker et al. (1999) and that a forced terrestrial release would be required to explain the model's divergence in δ¹³C from the observed trend after approximately 550 BC.

More recently, Goodwin, Oliver, and Lenton (2011) proposed that a reduction in the soft tissue (organic carbon) pump to the deep, possibly in conjunction with a calcite sedimentation event, could explain much of the Holocene carbon trend when taking into account δ¹³C in both ice cores (representing atmospheric carbon isotope fractionation) and marine sediments (dissolved inorganic carbon in benthic foraminifera, documenting deep-ocean carbon isotope fractionation). In particular, a synthesis of data in Goodwin et al. (2011) from mostly western Pacific sources supports the idea that deep water δ¹³C was increasing during the Holocene, which suggests that less organic carbon (strongly negative fractionations of δ¹³C) was being transported to the deep ocean.

c Arguments on the Holocene Carbon Cycle Based on δ¹³C

Many studies, including Goodwin et al. (2011), point to the small δ¹³C excursions of 0.2 ppt in Taylor Dome (since 5000–6000 BC) and an even smaller 0.05 ppt (since 4000 BC) decrease in the higher-resolution EPICA ice core (Elsig et al., 2009; Indermühle et al., 1999; Schmitt et al., 2012) as evidence that a major terrestrial release did not occur during the late Holocene. Because the biosphere generally prefers the C-12 over the C-13 carbon isotope, most vegetation carbon (and biological carbon generally) has a negative δ¹³C fractionation, and a release of terrestrial carbon to the atmosphere would accordingly lead to a measurable decrease in δ¹³C in the atmospheric CO₂ value obtained from ice cores. Based on this simple assumption, Joos et al. (2004) calculated that a decrease of 0.6 ppt in δ¹³C would be required for a 40 ppm human land-use contribution, well above the uncertainty in measured δ¹³C fractionation of CO₂ in ice cores, suggesting that the carbon flux that led to the 20 ppm increase in atmospheric CO₂ during the Holocene came from an ocean source (with approximately 0 ppt δ¹³C fractionation).

Challenges to this interpretation, however, have come to light in recent studies on high-latitude carbon cycle dynamics. In particular, Zimov et al. (2009) noted that permafrost covered large parts of Siberia, North America, and Europe during the Last Glacial Maximum (LGM). As a result of slow cooling leading up to the LGM, more and more carbon litter was frozen into the permafrost layer over the course of 100,000 years (Zimov et al., 2009). With ice sheet retreat and warming, the permafrost disappeared in many areas or retreated to greater soil depths in others, exposing large quantities of carbon to decomposition and a reduction of 36–46 kg C m⁻² in the permafrost system. This implies a release of more than 1000 PgC to the atmosphere (Zimov et al., 2009). Similarly, Ciais et al. (2012) compared fractionation of O-16 and O-18 in ice cores, which suggests global photosynthesis (GPP) approximately 57% of its pre-industrial value during the LGM; however, global carbon stocks (from δ¹³C changes in ice cores and marine sediments) imply a terrestrial carbon reservoir during the LGM of only 330 PgC (approximately 8%) less than during the pre-industrial period. The authors thus inferred that the non-photosynthizing fraction of biological carbon (approximately 700 PgC) is likely dominated by inert storage in permafrost, which would have been released during deglaciation.

Within the context of a permafrost release, the δ¹³C of dissolved inorganic carbon (DIC) in benthic foraminifera does not demonstrate a δ¹³C decrease during the late Pleistocene or early Holocene (which would indicate a flux of organic carbon to the deep ocean via the atmosphere). In fact, an overview of data from several ocean basins in Oliver et al. (2009) indicates an increase of 0.4–1.0 ppt from the LGM to 5000 BC in the Atlantic basin and an increase of 0.1–0.4 ppt in the Pacific and Indian oceans, whereas Ciais et al. (2012) calculated a global marine (benthic) δ¹³C increase of 0.34±0.13 ppt from the LGM to the pre-industrial period. Zimov et al. (2009) thus assumed that a permafrost release may have been taken up by terrestrial vegetation, with the fraction absorbed by the ocean balanced by a reduction in the soft carbon pump (transfer of sea-surface biological carbon to the deep ocean) or faster ocean ventilation during the late Pleistocene and early Holocene than during the LGM. A gradual release of permafrost carbon over this period may have transferred more of this carbon from the permafrost to northward-expanding forests and peatlands (experiencing accelerated accumulation during the early Holocene) instead of into the ocean, leading to little net biosphere extraction of CO₂ from the atmosphere–ocean carbon reservoir, thus reducing the role of calcite compensation in the Holocene CO₂ increase (Kleinen et al., 2010; Ruddiman et al., 2011).

At the beginning of the Holocene, slowly increasing atmospheric δ¹³C (approximately 0.2 ppt SST-corrected δ¹³C between 10000 BC and 4000 BC in Schmitt et al. (2012)) and marine δ¹³C have been linked to terrestrial uptake, in concert with contributions from changes in ocean circulation and oceanic biological production (Schmitt et al., 2012). Later during the middle and late Holocene, the substantial terrestrial uptake of peat (extracting C-12 from the atmosphere) documented in Yu (2011) should have led to increasingly positive δ¹³C in atmospheric CO₂ from ice core records instead of the slight decrease seen in Indermühle et al. (1999) and Elsig et al. (2009). In assuming that a virtually constant δ¹³C implies little release from the terrestrial biosphere, Joos et al. (2004) and Elsig et al. (2009) may have substantially underestimated the positive δ¹³C contribution from peatland uptake (Yu, 2011). Using the new peatland estimates established in Yu (2011), a terrestrial release on the order of 23–24 ppm (or 330 PgC) would be needed to counterbalance the peat uptake and the δ¹³C-implied atmospheric CO₂ increase of 5 ppm from terrestrial sources, which likely stems from natural (monsoon weakening and retreat of boreal forests) or anthropogenic (human land use) causes (Ruddiman et al., 2011).

d The Contribution of the UVic Model to the Holocene Carbon Cycle Debate

In this study, we use an intermediate complexity climate model to test the hypotheses involving terrestrial carbon release, sea surface temperature, and meridional overturning controls on the carbon cycle. We present simulations of the entire period of CO₂ increase (from 6000 BC to AD 1850) using the University of Victoria Earth System Climate Model (UVic ESCM) version 2.9. The model includes a fully interactive terrestrial, atmospheric, and oceanic (organic/inorganic/sediments) carbon cycle, which represents coupling and feedbacks that occur between the atmosphere, ocean, and biosphere. While the Holocene has been modelled with other intermediate-complexity models, this is the first time that this period has been extensively investigated using the UVic ESCM in continuous, fully interactive transient simulations. Furthermore, some previous studies of the Holocene climate (e.g., Renssen, Goosse, Fichefet, Masson-Delmotte, & Koç, 2005; Schurgers et al., 2006) assume unchanging pre-industrial land ice for the entire Holocene, whereas our simulations test the sensitivity of the climate system to changing ice sheets and in particular to marine ice shelves.

It is important to emphasize that these experiments only model the natural carbon cycle and potential natural vegetation changes (with no incorporated anthropogenic land use). Furthermore, the UVic ESCM does not presently include wetland or peatland processes or the upslope growth of coral reefs with increasing sea level. Without addressing these issues and their uncertainties, the UVic ESCM experiments allow a thorough investigation of carbon transport through physical ocean dynamics, a factor not investigated in studies such as Kleinen et al. (2010), and feedbacks from different ice shelf distributions. Our results suggest that the coupling between the carbon reservoirs changes substantially for different ocean thermohaline circulation states and land-ice configurations and that these would thus likely have a notable influence on how much CO₂ from human land-use emissions or other carbon sources remains in the atmosphere on millennial time scales.

The following sections are structured as follows. First, we give a description of the UVic model (Section 2), then we provide the experimental set-up (Section 3a), results for the model spin-up for 6000 BC (Sections 3a and 3b), and the set-up for transient simulations (Section 3c). In Section 4a, we evaluate the sensitivity of the carbon cycle to partial and complete ventilations of stagnant North Pacific deep water, and in Section 4b we present the results of transient simulations in which the atmospheric carbon is allowed to evolve freely from 6000 BC to the end of the pre-industrial era. In Section 4c, we compare the results of these experiments with transient runs in which atmospheric CO₂ is prescribed from observed ice core data, and in Section 4d we investigate the sensitivity of the vertical circulation of the oceans and deep ocean carbon storage to Antarctic ice shelf extent. Finally, Section 5 affords a discussion of the limitations of our model and methodology, and in Section 6 we present our major conclusions.

2 Model description

The UVic ESCM 2.9, classified as an Earth system Model of Intermediate Complexity (EMIC) (Claussen et al., 2002), has been gradually built from the original UVic model described in Weaver et al. (2001). Since its development, the UVic ESCM has been employed in a number of climate modelling efforts, including the Fourth Assessment Report of the Intergovernmental Panel on Climate Change (IPCC) (Meehl et al., 2007) and the Coupled Climate-Carbon Cycle Intercomparison Project (C⁴MIP) (Friedlingstein et al., 2006). Furthermore, the UVic model (and its incorporated Top-Down Representation of Foliage and Flora Including Dynamics (TRIFFID) dynamic global vegetation model) has also been used previously to investigate land-use impacts on climate over the past three centuries (Matthews, Weaver, Meissner, Gillett, & Eby, 2004). They showed that the total effect of land-use change for the above period was an increase in global temperature of 0.15°C. The research presented in this paper extends the application of the UVic model to include most of the Holocene epoch.

The core components of the model representing the atmosphere (a one-layer Energy Moisture Balance Model (EMBM) calculated every 1.25 model days) and the ocean (the 3D Modular Ocean Model general circulation model (GCM) v. 2.0, with 19 layers, calculated every 2.5 model days) were coupled in Fanning and Weaver (1996) and refined in Weaver et al. (2001). Coupling between the atmosphere and the ocean occurs every five model days. Moisture advection in the model relies on prescribed winds, and in this study climatological winds from the National Centers for Environmental Prediction (NCEP; Kalnay et al., 1996) were used. However, in order to test the robustness of using twentieth century winds in our simulations, all experiments were repeated with the model's wind feedback mechanism, which adjusts the NCEP wind climatology according to the real-time model-generated air temperature field (Weaver et al., 2001). Isothermal land ice and thermodynamic–dynamic sea ice were also included in the original model (Weaver et al., 2001). Land ice was prescribed in the model according to ICE-4G in Peltier (1994); it contributes no freshwater pulse to the ocean as it melts, and where it extends over the ocean in the form of ice shelves, it simply acts as an insulating lid with no latent heat exchange with waters below or beside it.

To represent the terrestrial biosphere, the TRIFFID dynamic global vegetation model (Cox 2001) was adapted to the UVic model (Matthews, Weaver, & Meissner, 2005; Meissner, Weaver, Matthews, & Cox, 2003). This module generates five plant functional types (PFTs), including broadleaf trees, needleleaf trees, C3 (mid- and high-latitude) grasses, C4 (mainly tropical) grasses, and shrubs (or small trees), which each grow under specified bioclimatic (temperature and precipitation) conditions. Where more than one PFT can grow, a dominance hierarchy (based on canopy height) prefers trees to shrubs and shrubs to grass, with grass types and tree types competing with each other (see Cox (2001) for more details). The TRIFFID model also includes a one-layer soil carbon model (depth 1 m) driven by biomass input and microbial respiration. Primary production and the resulting carbon fluxes between the land and the atmosphere are calculated by a version of the Met Office Surface Exchanges Scheme (MOSES; Cox et al., 1999; Meissner et al., 2003). The oceanic carbon cycle module was added by Schmittner, Oschlies, Giraud, Eby, and Simmons (2005); this is a classic nitrate-phytoplankon-zooplankton-detritus (NPZD) model with dissolved organic matter (but lacking iron limitation and nitrification/denitrification processes). Inorganic carbon, ocean sediment deposition and dissolution, and the evolution of the calcite lysocline occur as in the model developed in Archer (1996), which includes only the oxic respiration of sediments (Eby et al., 2009). Because there is no comprehensive weathering module in the model, the land-to-ocean weathering rate was diagnosed by the sediment burial flux for the simulations presented in this paper. This results in constant ocean alkalinity and no contribution from calcite compensation.

For all subcomponents, the model operates on a 1.8° x 3.6° (latitude-longitude) horizontal grid (Weaver et al., 2001). The model simulations presented here (see Table 1) are performed with the carbon cycle “online” with respect to the climate system, so that the two evolve interactively and the sum total of all feedbacks is considered.

3 Equilibrium simulations and transient simulation set-up

a The Model Spin-Up for 6000 BC

To spin-up the model to an equilibrium state, the UVic ESCM 2.9 was run for more than 10,000 model years, forced with the ICE-4G land-ice configuration (e.g., Fig. 1a) (Peltier, 1994), orbital forcing (Berger, 1978), and carbon dioxide (260.2 ppm) (Indermühle et al., 1999) all fixed at their 6000 BC configuration. Moisture advection and wind stress over the ocean were determined by winds from the NCEP reanalysis climatology (Kalnay et al., 1996). This simulation will henceforth be called EQ_1 (see Table 1, which describes the other equilibrium and transient simulations discussed in this study). With all forcing mechanisms being time independent, a stable equilibrium was expected to evolve. However, even with the assigned parameters, the model never reached a steady state. Instead, the meridional overturning circulation (MOC) weakened (by approximately 4–5 Sv (where 1 Sv = 10⁶ m³ s⁻¹)) and rebounded in an oscillatory fashion as shown in Fig. 2a. Many other prognostic variables also exhibited these multi-millennial oscillations (such as mean depth-integrated ocean temperature and terrestrial carbon (see Figs 2b and 2d)).

Fig. 1 Ice sheet distributions around Antarctica for (a) 6000 BC, (b) 5000 BC, (c) 4000 BC, and (d) 3000 BC according to ICE-4G (Peltier, 1994) interpolated to the resolution of the UVic ESCM. The model's land–ocean boundaries around Antarctica are provided in (e).

Fig. 2 Selected results from EQ_1 (spin-up for 6000 BC, fixed atmospheric carbon), EQ_2 (spin-up for 6000 BC, free atmospheric carbon), and EQ_3 (spin-up for 6000 BC, fixed atmospheric carbon, wind feedback). In (a) the maximum meridional overturning streamfunction (Sv) for the oceans is shown for all three simulations, and (b) illustrates mean depth-integrated ocean temperature (°C) for the same simulations; (c) shows the atmospheric carbon dioxide, and (d) illustrates the evolution of terrestrial carbon.

In order to determine the interactive exchange between the carbon cycle and the different thermohaline circulation states, another simulation (EQ_2), starting from Model Year 6200 of EQ_1, was performed with freely evolving atmospheric carbon but with land ice, winds, and orbital forcing kept fixed as in EQ_1. Despite changing atmospheric CO₂, this new simulation also produced millennial-scale oscillations in the MOC, which resulted in a 12 ppm CO₂ fluctuation in the atmosphere over about 4000 years (Fig. 2c).

In EQ_2, a similar terrestrial carbon release (approximately 7 PgC) occurred in the first hundred years of the flush, initially helping to propel atmospheric CO₂ concentrations upward. However, for the following 1000 years in EQ_2, the terrestrial carbon quickly rebounded and exceeded its original level by 15 PgC. This indicates that the elevated atmospheric CO₂ for nearly 1000 years during the millennial-scale events in the equilibrium simulation were exclusively the consequence of a slow, thorough flushing of deep ocean carbon to the atmosphere, which in turn fertilized the biosphere.

b A Conceptual Model of Southern Ocean Flushing Events

On running the UVic ESCM (v. 2.8), Meissner, Eby, Weaver, and Saenko (2008) discovered a similar intermillennial MOC oscillation in several equilibrium simulations under present-day forcing conditions. In their study, atmospheric CO₂ concentrations were kept constant for separate simulations, taking values between 360 ppm and 720 ppm. All simulations with fixed CO₂ levels over 400 ppm demonstrated these flushing events (total MOC weakening), with the periodicity of these deep-ocean ventilations dependent on the level of prescribed atmospheric CO₂. Furthermore, each flush was associated with deep water ventilation and a simultaneous decrease in sea-ice extent near Antarctica. This was subsequently followed by long periods of extensive sea ice and strong sea surface stratification, which served to prevent ventilation through the Southern Ocean. Changes in sea-ice extent, albeit on a shorter time scale, have been shown to be able to trigger changes in the thermohaline circulation in the North Atlantic Ocean on this order of magnitude (Mauritzen & Häkkinen, 1997).

Following the flushing events in Meissner et al. (2008) and in this study, the deep ocean warmed slowly over the course of a few thousand years through diffusion of heat, with Southern Ocean overturning by strong stratification as a result of more extensive seasonal sea ice. Eventually the deep ocean warmed to the point that it became temporarily unstable and another flushing event was initiated. Meissner et al. (2008) tested the sensitivity of these oscillations to the prescribed winds by carrying out several simulations with the model's wind feedback parameterization (see Weaver et al., 2001) turned on and off. The oscillations occurred even with the wind feedback (at a slightly different amplitude). Thus, the authors concluded that a CO₂ threshold near 400 ppm was probably responsible for the simulated millennial-scale phenomenon and might be a significant feature of a future climate with atmospheric CO₂ concentrations above 400 ppm.

In our study, a comparable flushing of the Southern Ocean and the North Pacific Ocean was produced by the model's oceanic general circulation model (OGCM). As in Meissner et al. (2008), this ventilation event occurred for a separate spin-up for 6000 BC (EQ_3) including the model's wind feedback mechanism, albeit producing flushes with a higher frequency and lower amplitude than simulations using the NCEP winds. The spatial output (not shown) from the fixed CO₂ (EQ_1 and EQ_3) spin-ups and the free carbon spin-ups (EQ_2) yielded a similar and relatively coherent picture of the ocean state, with a slow but comprehensive flushing of the deep ocean carbon completed in the one thousand years following the initiation of an event. After each flushing event, the model ocean below 1 km was remarkably “newer” (evidenced by Δ¹⁴C ages of deep-ocean water) and colder.

In summary, the flushing events produced in the three equilibrium simulations were associated with deep ventilation in the Southern Ocean and corresponded to the global replacement of very deep waters, following the same pattern as in Meissner et al. (2008). However, the presence of these cycles in this paper's spin-up for 6000 BC affords a new perspective on the model's millennial-scale variability. Our equilibrium simulations (EQ_1 and EQ_2) under mid-Holocene conditions produced these events at much lower CO₂ concentrations (260 ppm) than the 400 ppm threshold proposed by Meissner et al. (2008). Because both our results and the results of Meissner et al. (2008) suggest that the wind feedback parameterization does not reduce or eliminate the occurrence of such oscillations, the nature of mid-Holocene solar forcing and land ice play a critical role in producing the millennial-scale oscillations during the equilibrium simulation.

Another prominent feature of these flushes was a change in deep water distribution globally before (Fig. 3a) and after the flush (Fig. 3b). Figure 3 provides Δ¹⁴C profiles calculated by the model for the Atlantic basin in EQ_2, with lower Δ¹⁴C values generally indicating stagnant, older deep waters and more positive values indicating newer deep waters recently exposed and well mixed with the atmosphere. As can be seen in both Figs 3a and 3b, NADW retains characteristically higher Δ¹⁴C values, whereas Δ¹⁴C in deep waters in the Southern Ocean tends to be much more negative. The scenario just before the flush in Fig. 3a indicates that Atlantic deep waters were strongly dominated by NADW, which is relatively carbon deprived (for further discussion see Section 4d). However, after the flush (Fig. 3b), the deep waters in much of the Atlantic Basin below 2000 m depth had a much more prominent Antarctic signature and were also more enriched in carbon than the NADW that it replaced.

Fig. 3 Vertical profiles of Δ¹⁴C at grid cells in the Atlantic basin (a) immediately before the flush in EQ_2 (year 9095 in the simulation) and (b) 1400 years later at the end of the flush event in the EQ_2 (year 10495 in the simulation). The blue lines represent profiles in the North Atlantic: 54.9°N, 41.4°W, 29.7°N, 55.8°W, 0.9°N, 23.4°W; the green lines represent profiles in the South Atlantic and Southern oceans: 29.7°S, 23.4°W, 42.3°S, 30.6°W, 69.3°S, 41.4°W.

In general, DIC, the dominant contributor to total deep ocean carbon, is directly correlated with Δ¹⁴C in most regions in the model so that more negative Δ¹⁴C deep waters generally accumulate and hold more DIC. By the time deep water distributions reached those from Fig. 3b in the EQ_2 simulation, atmospheric CO₂ was decreasing. This suggests that the cold, carbon-rich, abyssal Antarctic-generated waters, henceforth broadly referenced in this paper as Antarctic Bottom Water (AABW), were being buried in the Pacific and North Atlantic basins (outside the Southern Ocean) and were no longer being mixed as frequently to the surface. As NADW slowly retook the deep Atlantic basin in the lead up to the next flush (eventually returning to the state in Fig. 3a), Antarctic-generated waters retreated to the Southern Ocean, where they were mixed to the surface more frequently along the Antarctic Circumpolar Current (ACC). As a result, much of the excess carbon buried in the deep Atlantic basin was slowly returned to the atmosphere, helping produce the gradual increase in CO₂ just prior to the flush. Then the flush itself was characterized by a ventilation of the stagnant and even more carbon-rich North Pacific Deep Water (NPDW).

It is important to note that the MOC events produced in EQ_1, EQ_2, and EQ_3 are associated with an unrealistically stable forcing (time-independent land ice, CO₂, and orbital forcing for thousands of years). However, such a millennial-scale event in the Southern Ocean has already been proposed by Ruddiman (2007) as a possible explanation for the consistent increase in atmospheric CO₂ during the late Holocene, and in nature it would represent a significant CO₂ release to the atmosphere on the order of 10 ppm. Schmittner, Brook, and Ahn (2007) and Ahn and Brook (2008) demonstrate similar, millennial-scale fluctuations in atmospheric CO₂ and Southern Ocean stratification, respectively, resulting from changes in regional overturning. Furthermore, it has also been argued that longer-term modes of meridional circulation in the Southern Ocean may be a crucial driving force in glacial-interglacial variations of atmospheric CO₂ (Stephens & Keeling, 2000). Therefore, the thermohaline circulation is an important focus of investigation in this study.

c Set-up of Transient Simulations

Considering that it would be difficult to verify 4–5 Sv variations in the strength of the MOC using proxy data, given the limited temporal and spatial resolution of sediment cores, there is a great deal of uncertainty about what constitutes the actual state of the MOC starting at 6000 BC. Given this ambiguity, any Southern Ocean ventilation state generated during the spin-up could be a valid starting point for transient simulations. The significant influence of the MOC on the evolution of atmospheric carbon in EQ_2 (Fig. 2c) thus provides the motivation behind a series of sensitivity studies described below.

The “default” starting point for several reference simulations was selected to be 6200 years following spin-up initiation in EQ_1, corresponding to maximum MOC strength (see Fig. 2a) and a deep water distribution between the two extremes in Figs 3a and 3b. Other simulations were started during a flushing event (DE), at a moment when the total MOC had stabilized at a low transport value in the spin-up and AABW constituted the dominant deep water below 3000 m depth in the oceans globally (see Figs 2a, 2b, and 3b). In addition, more experiments were started at about 200–300 years before the initiation of an event (BE) in the spin-up simulation, when warmer NADW constituted the dominant water mass in the deep Atlantic (see Fig. 3a). Altogether, these different initial conditions provided the launching points for transient simulations of the carbon cycle with prescribed and free atmospheric carbon, denoted PC and FC, respectively. The simulations summarized in Table 1 were then repeated with the model's wind feedback.

The transient forcing is provided by time dependent changes in orbital parameters (Berger, 1978) and prescribed land ice (Fig. 1), which are interpolated for the model year. Sudden transitions (rather than interpolations) from land-ice area to potential vegetative area occur in the model at 5500 and 4500 BC, which explains the instantaneous jumps evident in terrestrial carbon variables in the mid-Holocene as the last vestiges of the Laurentide ice sheet disappeared. The prescribed atmospheric carbon included in the model is based on Taylor Dome (Indermühle et al., 1999) from 6000 BC to AD 1006 and from Law Dome (Etheridge et al., 1996) from AD 1006 to the end of the simulation at AD 1850. In the prescribed carbon runs (PC, BE_PC, DE_PC), the model's atmospheric CO₂ followed the observed trend in Taylor and Law Dome cores while all other components of the carbon cycle (including vegetation/soils, organic and inorganic ocean carbon, and ocean sediments) responded directly to the forced changes in the atmospheric carbon reservoir but were free to evolve without any other constraint. By contrast, the FC simulation allowed the model's atmospheric carbon content to respond freely to the evolution of the model's interactive representation of the Earth's natural carbon cycle and climate.

In summary, the equilibrium simulations suggest that two “extreme” solutions are equally plausible for 6000 BC deep-water distributions with (1) little AABW, stagnant NPDW in the Pacific Ocean, and dominant NADW in the deep north and south Atlantic Ocean as in Fig. 3a, and with (2) newly ventilated NPDW and dominant AABW in the North Pacific and Atlantic basins below 3000 m (Fig. 3b). Furthermore, a range of solutions are possible between these two extremes. The BE initial conditions represent a starting point for case (1); the DE initial conditions represent case (2), and the default initial conditions represent a deep-water distribution between these two extremes.

4 Results

a Millennial-Scale Oscillations under Transient Forcing

As discussed in Section 3b, the presence of millennial-scale oscillations in the thermohaline circulation for the model's equilibrium simulations produced a venting of deep waters globally through the Southern Ocean and a temporary slowing of NADW formation. When atmospheric carbon was allowed to vary freely (EQ_2), this resulted in an increase in atmospheric CO₂ of 8–13 ppm, likely a result of carbon-rich deep waters from the North Pacific being ventilated at the surface (Fig. 4a) and a temporary slowing of CO₂ entrainment into the deep ocean in the North Atlantic. In Meissner et al. (2008), it was determined that a CO₂ concentration of 400 ppm was the critical threshold for such millennial-scale oscillations. However, with mid-Holocene land-ice and insolation forcing configurations, this threshold is clearly much lower (at least 260 ppm), and these oscillations continue to occur in sensitivity simulations when CO₂ (EQ_2) and winds (EQ_3) are allowed to vary. Considering their substantial impact on the carbon cycle, it is useful to identify which factor (insolation or land ice) is more important for their initiation and whether these millennial-scale oscillations can be present under transient (time-evolving) forcing conditions.

Fig. 4 Time series of (a) the maximum meridional overturning streamfunction (Sv) and (b) atmospheric CO₂ for a transient simulation including the wind feedback for the Holocene (6000 BC to AD 1850) with land ice fixed at its 6000 BC configuration.

A number of sensitivity simulations were conducted with transient orbital forcing from 6000 BC to AD 1850 but with land ice fixed at 6000 BC (see Fig. 1a) for the entire run period. The simulations, summarized in Fig. 4, start from the default initial conditions defined in Section 3c and continue to produce millennial-scale flushes (with similar characteristics to EQ_2) during the late Holocene (Fig. 4). These repetitive millennial-scale oscillations do not explicitly occur in the other transient simulations described in Table 1, which are all forced with both transiently evolving land ice and orbital forcing. This suggests that the distribution of land-ice extent (particularly Antarctic land ice prescribed by ICE-4G) in the model, and not orbital forcing, is critical to producing these millennial-scale oscillations.

In a simulation with constant CO₂ at 260.2 ppm and NCEP winds, a millennial-scale flush does not occur until AD 1300. When the carbon cycle was freed at AD 1200 in this simulation (not shown), a 7.5 ppm jump in atmospheric CO₂ was produced within a decade as a ventilation of NPDW was initiated. Atmospheric CO₂ then stabilized at 268 ppm for 500 years before declining slowly toward 260 ppm. In a similar simulation with land ice (Fig. 1a) and CO₂ (260.2 ppm) fixed at the 6000 BC configuration, but including the wind feedback adjustment to NCEP winds, two millennial-scale flushes occurred (amplitude of 4 Sv), one at 5500 BC and another at 1500 BC. These oscillations continued with a 3000-year period (Fig. 4a), similar to the equilibrium simulation (EQ_3, Fig. 2a) which also incorporated the wind feedback. When this experiment was repeated with atmospheric CO₂ allowed to evolve freely (Figs 4a and 4b), the earlier oscillation (5500 BC) occurred at a slightly greater amplitude (4 Sv) than the flushing event at 1500 BC (2 Sv anomaly). Both events in this free-carbon simulation led to a slow flushing of NPDW, but the first event (Fig. 5b) produced a more thorough replacement of NPDW (evidenced by Δ¹⁴C) than in the second flush at 1500 BC. Consequently, the carbon released to the atmosphere was greater during the first flush (a 10 ppm excursion comparable to the release in EQ_2) than in the second (a 6 ppm excursion).

Fig. 5 Spatial distributions of Δ¹⁴C (ppt) at 3202 m depth before and after the ventilation of NPDW in several equilibrium and transient simulations. Shown in (a) are the results from EQ_2 year 9195 (left panel) and year 10195 (right panel), (b) results from 5600 BC (left panel) and 4100 BC (right panel) in a simulation with transient orbital forcing, wind feedback, freely evolving atmospheric CO₂ and fixed 6000 BC land ice, (c) results for 5100 BC (left panel) and 4100 BC (right panel) in the default FC simulation, and (d) results for 5100 BC (left panel) and 4100 BC (right panel) in the Before Flush FC wind feedback simulation.

Fig. 6 Time series of (a) atmospheric CO₂ (ppm), (b) maximum meridional overturning streamfunction (Sv), (c) terrestrial carbon (PgC), (d) ocean carbon (PgC), and (e) sediment carbon (PgC). The results from simulations with prescribed NCEP winds for the entire simulation (as in Table 1) are given in the left panels, and the results using the wind feedback for the same initial conditions are given in the right panels. The dark blue line represents the prescribed CO₂ (PC) simulation using default initial conditions, whereas the cyan (light blue) line shows free CO₂ (FC) evolving from the same default initial conditions. The green line represents before flushing event (BE) initial conditions for free CO₂ (solid line) and prescribed CO₂ (dashed line). The red curves represent during flushing event (DE) initial conditions with free CO₂ (solid line) and prescribed CO₂ (dashed line).

The combined effect of both flushes in this simulation kept atmospheric CO₂ substantially higher (near or above 260 ppm) for a much longer interval than in the transient simulations from Table 1 (see Fig. 6a and discussion in Sections 4b and 4c). These oscillations, in turn, appear to be related to the presence of an extensive land-ice shelf prescribed around the periphery of Antarctica (Fig. 1a). The flushing events were also associated with a more frequent ventilation of very old, carbon-rich North Pacific deep water through the Southern Ocean and, as a consequence, higher atmospheric CO₂ concentrations near or above 260 ppm for much of the simulated late Holocene climate. The results thus suggest that atmospheric CO₂ is highly sensitive to land-ice extent, and this was the motivation for further sensitivity experiments concerning ice shelves in Section 4d.

Similar flushes of NPDW have been recorded by proxy evidence in nature. Galbraith et al. (2007) note that the ventilation of long-stagnant North Pacific deep waters occurred around 12650 BC and likely contributed to an abrupt 10 ppm increase in atmospheric CO₂, which was balanced, in part, by a strengthening NADW and greater CO₂ uptake through biological productivity in the surface ocean. The North Pacific Ocean thereafter became more estuarine following this event (with more frequent ventilation). The UVic model oscillations simply serve to circulate deep ocean carbon stored in the North Pacific, albeit in abrupt events that resemble the first major interglacial ventilation of NPDW at 12650 BC in the proxy data. While a flush of the same magnitude as at 12650 BC has not been documented since 6000 BC; periodic, partial ventilations of NPDW could well be a feature of the Holocene MOC.

The simulations with transiently evolving land ice elaborate on this idea of partial NPDW ventilations. None of the fully transient simulations summarized in Table 1 (see Fig. 6), with and without the model's wind feedback, produced a millennial-scale oscillation on the same scale as in the equilibrium simulations EQ_1 and EQ_2 or the experiments with land ice fixed at the 6000 BC configuration (Fig. 4). However, NPDW was still ventilated in all free carbon simulations, except with DE initial conditions, in which the NPDW had already been ventilated at the start of the transient simulation. Furthermore, both the default and BE transient simulations flushed these deep waters in an event around 5000 BC (Fig. 6c) that produced a notable loss in global ocean carbon (Fig. 6d). More dramatically, the BE IC simulation with freely varying atmospheric carbon and a wind feedback demonstrated a multi-centennial-scale weakening in the MOC near 5000 BC (Fig. 6b), which was accompanied by a more gradual and thorough (Fig. 6d) flushing of NPDW (according to Δ¹⁴C ages) and greater oceanic carbon loss (Fig. 6d) than any other simulation.

b Results from Free CO₂ Experiments

1 Atmospheric CO₂

From the transient runs summarized in Table 1, the marked contrast in the atmospheric CO₂ concentration in the FC experiment (Fig 6a, light blue line) compared to the observed atmospheric carbon from ice cores (Fig. 6a, dark blue line) may provide important implications for the early anthropogenic hypothesis. In all simulations, atmospheric CO₂ decreased to values between 245 ppm and 255 ppm by the Industrial Revolution, well below the 280 ppm observed value and in-line with the analogue prediction of decreasing Holocene CO₂ in Ruddiman (2003) and Ruddiman et al. (2011).

Initially, atmospheric CO₂ increased in the FC experiment, nearly reproducing the observed trend through natural processes (without the intervention of peatlands, coral reefs, or human land use) until about 4000 BC In the BE simulation including the wind feedback, the divergence between the model-simulated carbon cycle and the natural carbon cycle does not occur until even later (3000 BC). This suggests that the initial CO₂ increase during the mid-Holocene might be explained by natural processes represented in the model, delaying the date at which an external source of carbon would be required to recreate the Holocene carbon trend.

A later date of divergence between the Holocene's natural potential carbon cycle and the observed cycle in the proxy record may theoretically lend more weight to the revised Anthropogenic Hypothesis (Ruddiman & Ellis, 2009), because the late neolithic period was characterized by not only low land-use efficiency but also much higher human populations following several millennia of horticultural, agricultural, and pastoral subsistence. With global populations nearing 100 million by 4000–3000 BC by some estimates (Boyle et al., 2011), it can be speculated from previous studies summarized in Section 2b that the human transformation of the landscape began to have a measurable impact on global CO₂ concentrations.

On closer examination of the CO₂ trend lines in Fig. 6a anomalies (spikes) in the simulated CO₂ curve between 5000 and 4000 BC can be seen. These appear to be initiated by changes in the MOC (Fig. 6b, light blue line), which take the form of periods of deeper ventilation off the western Antarctic coast, south of Australia. However, these brief multidecadal-scale events tend to mix AABW and Antarctic Circumpolar waters locally in the Southern Ocean and do not have the same spatial and temporal scale as those observed in EQ_1 and EQ_2 (which ventilate NPDW). As a consequence, they do not appear to have the same broad influence on deep ocean temperatures and deep water carbon storage.

2 Terrestrial Carbon

The terrestrial carbon (vegetation and soil carbon) time series plotted in Fig. 6c indicate that vegetation carbon constantly increased over the course of the Holocene simulation, with a total terrestrial uptake over the period ranging from 70 to 100 PgC for the FC simulations. This led to net removal of carbon from the atmosphere–ocean system over the 8000 simulated years and is consistent with the results of Kleinen et al. (2010) and Goodwin et al. (2011), which also suggested similar net increases in terrestrial carbon of approximately 100 PgC over the course of the Holocene. However, the CLIMBER-LPJ model in Kleinen et al. (2010) simulated no net change in terrestrial carbon without peatlands: decreasing boreal vegetation during the late Holocene (carbon loss in living biomass) led simply to a transfer of carbon to litter (carbon gain to colder soils) as the tree line shifted south. Only in simulations including peatlands does the terrestrial carbon (soils and biomass) in Kleinen et al. (2010) increase by 105 PgC. The UVic ESCM 2.9 does not include an estimate of peatland accumulation but does similarly transfer shrub and tree carbon to soils during the late Holocene in the FC simulations associated with an insolation-induced southward retreat in the boreal forest line after 4000 BC (not shown). The high latitude C3 grasses that replaced trees and shrubs continued to accumulate soil carbon to produce the same net effect as the peatlands in Kleinen et al. (2010). The PC simulations, with steadily increasing temperature and greater biosphere CO₂ fertilization, produced instead a northward advance of the boreal forests during the late Holocene.

The results also indicate that the increase in terrestrial carbon is greatest from 6000 BC to 4000 BC as boreal vegetation expanded rapidly northward in Eurasia during the mid-Holocene climatic optimum and new land became available in North America as the last remnants of the Laurentian ice sheet disappeared. Because the Sahara region remained mostly bare ground in the equilibrium simulations, even in simulations incorporating the wind feedback, the model does not capture the terrestrial release in this region as proposed by Indermühle et al. (1999). However, even as boreal vegetation began to retreat southward around 4000 BC in the FC simulations, total terrestrial carbon continued to increase. This increase is primarily associated with a slowly increasing storage of carbon in high-latitude soils below sub-Arctic C3 grasses, which are not parameterized to imitate peatland dynamics.

In an approximate analogy to prairie grasses, model C3 grasses store a significant proportion (approximately 40%) of their carbon in the soil. As in areas with persistent permafrost, this storage is more pronounced at high latitudes and at high altitudes where respiration rates are low because of cold temperatures and a short summer season. Because the model shrubs do not store as much soil carbon, when shrubs advance northward and replace grass in the northern “tundra C3” region, soil respiration quickly overtakes biomass input. As a consequence, the carbon stored in the former high-latitude grass soils, which slowly accumulate carbon with time, is released to the atmosphere when overtaken by shrubs and small trees.

The tundra (not shown) nearly disappeared in pulses from 6000 BC to 4000 BC with northern hemisphere warming, then expanded southward after 4000 BC in the FC simulations as boreal forests retreated to the south in response to solar insolation changes. Clearly the accumulation in high-latitude soils more than compensated for the loss of terrestrial carbon resulting from retreating boreal vegetation, leading to a steady net increase in terrestrial carbon between 4000 BC and AD 1850.

The brief flushes in the transient simulations were associated with net terrestrial carbon releases of the order of 10 PgC (the dips in Fig. 6c), which help drive the concentration of atmospheric CO₂ upward, followed by a return of this carbon to the terrestrial biosphere a few decades later. The largest such release occurred when the NPDW waters were ventilated and ocean carbon was released around 5000 BC. Most of these brief releases occurred in high latitude regions of Eurasia and North America, where tundra carbon was being released to the atmosphere.

3 Ocean and Sediment Carbon

The ocean carbon reservoir (Fig. 6d) also played a critical role in the simulated Holocene carbon cycle. Between 5000 BC and 3500 BC, sediments gradually absorbed about 6 PgC (Fig. 6e, light blue line) before levelling off. However, this does not account for the 35 PgC decline in total ocean carbon (NPDW ventilation), most of which occurred between 5500 BC and 4500 BC (Fig. 6d, light blue line). Following this pronounced loss early in the simulation, ocean carbon remained virtually constant or in slow decline thereafter (with no substantial aging of NPDW during the late Holocene). The relatively stable, declining ocean carbon in the FC simulations is in contrast to the PC simulations, which showed an increase in ocean carbon in response to the forced increase in atmospheric carbon (Fig. 6d). The contrast between the FC and PC simulations clearly indicates that the ocean carbon in the prescribed runs primarily responded to atmospheric carbon and was not forced by the internal ocean chemistry.

Furthermore, a comparison of Fig. 6a with Fig. 6d suggests that the substantial increase in atmospheric CO₂ around 6000 BC, on the order of approximately 8 ppm (excluding the additive effect from terrestrial carbon), is associated first and foremost with a release of carbon from NPDW to the atmosphere. However, different ocean carbon losses between the simulations suggest differing degrees of “thoroughness” and rapidity (as discussed in Section 4a) in the flushing of NPDW. It is further interesting that each simulation with default and BE initial conditions produced some kind of NPDW ventilation around 5000 BC, which is surprising considering that one would expect an earlier ventilation from the BE initial conditions and a much more delayed event from the default (FC) initial conditions. This suggests that the orientation and transient evolution of marine ice shelves (for further discussion see Section 4d) is more important for the initiation of these NPDW flushing events than initial temperature and salinity of the deep water.

However, initial conditions and winds can clearly influence the thoroughness of the flushes of NPDW (and thus the quantity and duration of the carbon release to the atmosphere). The BE initial conditions with the wind feedback produced the most thorough ventilation of NPDW and the longest period of elevated CO₂ in the atmosphere. Default initial conditions resulted in a less thorough ventilation of NPDW (as evidenced by model Δ¹⁴C ages) and a shorter period of elevated atmospheric CO₂. The abrupt rise in simulated CO₂ to 270 ppm by 5000 BC seems exaggerated compared to the smooth, gradual increase in atmospheric CO₂ suggested by ice core records. However, if considering the source data for these proxies, ice in Antarctica must be given time to compact under newer layers of slowly accumulating snow above it to form isolated ice bubbles from which paleo-CO₂ concentrations can be extracted. Thus, the proxy data do not represent one single year (as do the model data) but rather decades if not centuries of mixed air bubbles, depending on the accumulation rate at the core site. A comparable running average of yearly data of the events produced in our model (not shown) produced lower-amplitude “humps” than is presented in Fig. 6a. This is more in keeping with the observed CO₂ trend from ice sheets, as humps (periods of accelerated CO₂ increases over a thousand years of approximately 5 ppm, followed by decreases on the order of 2–3 ppm) are visible in EPICA Dome C between approximately 3571 and 2787 BC and between approximately 968 BC and AD 198. Millennial-scale flushes of NPDW of a smaller magnitude than EQ_2 (periodic, incomplete ventilations on millennial time scales) thus may be embedded in the observed ice core data with other processes leading to the net increase in atmospheric CO₂ over the course of the Holocene.

The ocean carbon reservoir also interacts with centennial-scale changes in ventilation. During the periods when the Southern Ocean ventilation was locally active (characterized by dips in the maximum MOC streamfunction in Fig. 6b), relatively carbon-rich deep waters in the Southern Ocean were circulated to the surface and NADW uptake slowed. The brief loss of carbon from the ocean to the atmosphere during these events was minimal under NCEP winds but pronounced for the wind feedback simulations (see the zig-zag pattern in Fig. 6d). These centennial events occurred with relative frequency and are associated with elevated CO₂, higher Antarctic temperatures, and much greater Antarctic precipitation rates.

After 4000 BC, the observed and simulated natural atmospheric carbon diverged (Fig. 6a), with atmospheric carbon dropping rapidly between 4000 BC and 2000 BC. The most significant decline in atmospheric CO₂ concentration in the FC simulation occurred during and 1000 years after a permanent 2–3 Sv reduction in the maximum overturning streamfunction (Fig. 6b). Over this period, the cycles in Southern Ocean ventilation and Antarctic sea ice essentially stopped and mean-annual sea ice remained virtually unchanged for thousands of years.

The restabilization of the model MOC appears to be related to the disappearance of a major fraction of the Antarctic Ice Shelf at 3500 BC (interpolation between Figs 1c and 1d), which corresponds to the single greatest loss of southern hemisphere land-ice area during the 19000 BC to AD 2000 period in the interpolated ICE-4G dataset (it should be reinforced that the model does not experience this event as a freshwater pulse). At first, a large portion of the Weddell Sea was free of ice, which initially stimulated intense ventilation between 3500 BC and 3400 BC and a corresponding spike in atmospheric CO₂. Following this event, however, annual Southern Ocean sea ice expanded and reached a steady equilibrium, with areal coverage of sea ice (and thus winter production rates of sea ice in the Southern Ocean) changing little during the following millennia. Accordingly, virtually no periods of intense Southern Ocean ventilation were experienced for the rest of the Holocene simulation. Furthermore, with no additional abrupt NPDW ventilations, atmospheric CO₂ fell quickly to adjust to the lowered ocean carbon content (from 2000 years previously) as the MOC approached a new equilibrium in the transient simulation.

Following the stabilization of the MOC between 3000 BC and 2000 BC, simulated atmospheric CO₂ continued to drop slowly, approaching a limit of 252 ppm at the end of the pre-industrial era. The BE_FC simulation levelled out at a similar value (251 ppm) and the DE_FC simulation stabilized near 245 ppm. These values are significantly (approximately 25–35 ppm) lower than the ice core record of atmospheric CO₂ concentrations but are above the 240 ppm value projected by Ruddiman (2003) and used in the GCM studies of Vavrus et al. (2008) and Kutzbach, Ruddiman, Vavrus, and Philippon (2009). The range of 245–251 ppm by AD 1800 are, however, close to the average decrease (just below 250 ppm) of the six previous interglacial periods aligned by their (δ¹⁸O-determined) start date (Ruddiman et al., 2011).

The decline of atmospheric carbon appears to respond to the continuous uptake of carbon in the terrestrial biosphere. This increase is dominated by soil uptake as the high-latitude C3 grasses begin to expand southward, in part, as a result of declining northern hemisphere solar insolation after 4000 BC. However, the terrestrial uptake of approximately 30 PgC only accounts for the relatively small net 8 PgC (3 ppm) decline in atmospheric CO₂ from 2000 BC to AD 1850, approximately equivalent to the contribution of terrestrial release to the atmosphere from human land use proposed in Pongratz et al. (2009). Figure 6 reveals that much of the carbon gains in the biosphere during this stable MOC period after 2000 BC appear to be more or less balanced by a corresponding slow loss in ocean carbon.

c Comparison of Free Carbon Results with Prescribed Carbon Runs

As described above, a series of experiments, labelled PC, DE_PC, and BE_PC scenarios, were performed with prescribed atmospheric carbon, and the results are also summarized in Figs 6a to 6e. We have noted that the PC and FC carbon experiments yielded similar results for the first 1500–2000 years of both simulations; however, after 4000 BC a pronounced divergence is evident. Furthermore, we established in Sections 4a and 4b that the MOC state plays a particularly crucial role in the evolution of atmospheric carbon. On the surface, the evolution of the MOC (Fig. 2b) seems quite similar in the FC and PC runs, although the greatest variations between the two appear during times of greater ice shelf extent between 6000 and 3000 BC.

This suggests that emissions of carbon to the atmosphere might feed back to the MOC at times when Southern Ocean ventilation is more variable and net downwelling of Antarctic-generated waters is reduced. Small differences in atmospheric CO₂ between the individual PC and FC simulations has the largest impact on the MOC strength during this period of greater Antarctic ice shelf extent. When Antarctic ice shelves retreated after 3400–3300 BC, the MOC states between the PC and the FC simulations approximately converged. We thus infer that some other effect (such as the 3500 BC Antarctic ice shelf disappearance) and not CO₂ was driving MOC variability and deep water distributions. It is also clear from a comparison of the FC and PC simulations that the increasingly large divergence of atmospheric CO₂ (and the corresponding radiative forcing) between the two experiments had virtually no effect on the stabilized MOC after 3500 BC. In summary, our results suggest that the MOC is only notably altered by the CO₂ difference between the PC and FC simulations during periods when ice shelves are more extensive.

Our model results therefore indicate that the greatest impact resulting from small perturbations in CO₂ (such as those associated with terrestrial releases as proposed by Ruddiman (2007)) would have occurred before 3300 BC, the date of final MOC stabilization in all simulations. Although this date is very early in human agricultural history, it also may have been a time of low-efficiency (greater per-capita) land use related to the cruder technologies of the neolithic period (Ruddiman & Ellis, 2009). Furthermore, this period (before the disappearance of Antarctic marine ice shelves) is characterized by low-efficiency ocean uptake of carbon resulting from faster local circulation of Antarctic-generated deep waters to the surface and thus reduced long-term carbon storage in these waters.

However, it is not clear whether the multi-millennial steadiness obtained for the model MOC after 3300 BC is accurate (see Section 4d). In particular, proxy evidence of δ¹⁵N off the coast of Chile (Robinson, Mix, & Martinez, 2007), taken as a measure of the stratification of Southern Ocean surface waters (e.g., Ahn & Brook, 2008), shows variations during the Holocene that are nearly on the same order as those during Heinrich events earlier in the proxy record. Also, abundances of penguin rookeries described in Baroni and Orombelli (1994) and Emslie, Berkman, Ainley, Coats, and Polito (2003) suggest that 2000–1000 BC and AD 700–1400 were periods of low Antarctic sea-ice extent, interspersed with more extensive sea ice (for example, during the period from 1000 BC to AD 1). Furthermore, our modelling results suggest that Antarctic warming and cooling are closely coupled to Southern Ocean ventilation, and abrupt warmings and coolings during the Holocene have been documented in Antarctic proxies (for example, Kulbe, Melles, Verkulich, & Pushina, 2001; Shevenell, Ingalls, Domack, & Kelly, 2011).

The proxy evidence for significant variability in sea ice (and probably Southern Ocean overturning) implies that the production of sea ice off the coast of Antarctica following 3300 BC was not as steady as indicated by our results. Furthermore, the MOC naturally undergoes oscillations on the same order of magnitude and temporal scale as those seen in the FC and PC simulations between 5300 BC and 3300 BC, making circumspect the near constancy of the MOC for thousands of years after 3300 BC. Therefore, it is quite plausible that Southern Ocean ventilation continued to be locally active in the Southern Ocean during the late Holocene, and correspondingly the increasing extent of human land use with time might have mildly influenced the net terrestrial carbon release to the atmosphere observed for these events. The potential magnitude of these contributions should be investigated by future UVic simulations for different land-use scenarios using adaptations of HYDE 3.1 (Klein Goldewijk, 2001; Klein Goldewijk et al., 2011) and land-use scenarios from Kaplan et al. (2011).

The most compelling evidence for some kind of external carbon forcing not represented in the model, and potentially provided by human land use, calcite compensation, and/or coral reef regrowth, is demonstrated in Fig. 7. As seen in all the FC simulations, total earth system carbon is generally conserved by the UVic ESCM 2.9 when the carbon cycle is allowed to evolve freely. It is also approximately conserved for the PC and BE_PC simulations for the first 2000 years.

Fig. 7 The total carbon stored in the atmosphere, ocean, sediments, and terrestrial biosphere/soils for all simulations in Table 1.

However, after 4000 BC, total system carbon began to increase in these two PC experiments. A comparison with Fig. 6d reveals that much of this carbon is stored in the ocean so the atmosphere is clearly forcing the ocean carbon to increase in the prescribed runs. Furthermore, Fig. 6c reveals that terrestrial carbon also increased as did atmospheric CO₂ during the PC simulations as a result of greater CO₂ fertilization and more large (needleleaf) trees in subpolar regions replacing small trees (shrubs). The combined warming and fertilization effect from higher CO₂ concentrations enhanced carbon storage in the terrestrial biosphere much more than the greater soil carbon storage in high-latitude C3 grasses from the FC simulations (Figs 6c and 8c), which conversely responded to both decreasing CO₂ levels (cooling) and decreasing insolation after 4000 BC.

Therefore, some external forcing to the model's climate system is clearly necessary in order to reproduce the observed increasing atmospheric carbon trend shown in Indermühle et al. (1999) and Lüthi et al. (2008). This unaccounted-for input into the climate system, which increases with time, may be determined by the model's sediment and ocean chemistry initial conditions (see the limitations described in Section 5) or unaccounted-for volcanic activity. Peatlands (a sink for atmospheric CO₂ estimated to be as great as 35 ppm by Frolking and Roulet (2007) and slightly less by Yu (2011)) and coral reef-related calcite compensation (a source of atmospheric CO₂ estimated to be as much as 40 ppm by Ridgwell et al. (2003) and approximately 28 ppm by Kleinen et al. (2010)) are also not represented in the model. Although these two carbon pools might have balanced each other during previous interglacial periods, leading to little net influence on the atmospheric carbon concentrations, a slight lag or imbalance between the two during the Holocene could also be responsible for some of the increase in total system carbon. It is also worth noting that this increase in earth system carbon in the PC simulation exceeds 400 PgC, which is considerably larger than many estimates of the potential carbon flux to the atmosphere associated with early agriculture (Stocker et al., 2010) but not far from the emissions estimates (as high as 327–357 PgC) associated with crop rotations in Kaplan et al. (2011).

d Sensitivity of the Simulated Holocene Climate to Land-Ice Configuration

As established in Section 4a, the model's ability to ventilate long-stagnant, carbon-rich NPDW is strongly dependent on the land-ice configuration between 6000 BC and 5000 BC. Although this purging of deep ocean waters in the Pacific Ocean led to higher atmospheric CO₂ (largely between 258 and 270 ppm in transient simulations), it failed to explain the permanent increase in CO₂ to 280 ppm during the Holocene. However, these millennial-scale flushes kept atmospheric CO₂ near or above 260 ppm for the entire late Holocene, as opposed to 245–252 ppm in model simulations without comparable flushes. Therefore, we might speculate that a greater transfer of NPDW carbon to the atmosphere, associated with MOC sensitivity to land-ice extent, might explain part of the trend during the Holocene.

Furthermore, the results (Sections 4b and 4c) from transient simulations in Fig. 6 suggest that the land-ice transition between 4000 BC and 3000 BC led to a sudden stabilization of the MOC for the remainder of the simulated late Holocene, regardless of the initial conditions of the simulation or adjustments to the model winds. Comparing Fig. 1c with Fig. 1d indicates that this transition coincided with the complete disappearance of large ice shelves off the coast of Antarctica in the ICE-4G dataset; this land ice had only been melting slowly between previous time slices.

Because an abrupt and permanent stabilization of the model's MOC and climate state may not represent the reality of the Holocene; the transient simulations in Table 1 (with and without the wind feedback) were repeated by eliminating this transition. This was done very simply by fixing the land-ice configuration after 4000 BC to the configuration in Fig 1c. The motivation for these simulations is that the interpolation of ICE-4G to the model grid cell eliminates nearly all marine ice shelves (Fig. 1d compared with Fig. 1e), in particular, the prominent Ross and Ronne-Filchner ice shelves in the Ross and Weddell Seas, respectively, as well as extensive marine ice shelves off Graham Land (Antarctic Peninsula) and Queen Maud Land (east of Weddell sea) that still exist today (Vaughan 2008; their Fig. 1). Because sea ice, temperature, and ventilation depth are particularly variable in these regions during the modelled millennial- and centennial-scale flushes in the Southern Ocean (Sections 4a and 4b), these new experiments test the sensitivity of the global carbon cycle to the presence of Antarctic marine ice shelves during the late Holocene.

Prescribed ice shelves in the UVic ESCM 2.9 essentially act as a lid over the ocean grid cells below, cutting off fluxes (heat, moisture, carbon, and tracer) between the ocean and the atmosphere. As evident in a comparison of Figs 1a to 1d with Fig. 1e, grid cells are either open ocean or marine ice shelf (with no partial ice shelf coverage). With these ice sheets prescribed above the sea surface, grounding (marine ice shelves extending to the seafloor along continental shelves) is not represented. The open ocean grid cells adjacent to ice shelf margins are influenced by the albedo of the nearby ice sheet, which affects regional surface temperature and, correspondingly, precipitation and runoff to the ocean. In addition, the presence of the ice sheet raises the elevation of the grid cell (which also feeds back to regional temperature, precipitation, and runoff) and cuts off the ocean from overlying winds (which will in turn influence overturning). However, there is no exchange of latent heat between the ice sheet and the nearby or underlying ocean waters (which would serve to intensify convection along ice shelf margins), and the retreat of ice shelves does not contribute to a freshwater pulse to the ocean. Additionally, because ice shelves are only represented as a lid, ocean convection and circulation can continue beneath them.

Taking these considerations into account, Fig. 8 indicates there is a non-negligible difference in the carbon cycle dynamics between the simulations without ice shelves after 3000 BC (as for all the simulations represented in Fig. 6) and those with ice shelves (unchanging land ice after 4000 BC). Atmospheric CO₂ was notably higher for simulations with the ice shelves (Fig. 8a), yielding a pre-industrial value approaching 258 ppm with BF initial conditions, approximately 5–6 ppm higher than the same simulations without ice shelves. This increase was smaller (2–3 ppm) with the wind feedback because atmospheric CO₂ was already 1–2 ppm higher with the wind feedback in the no-ice-shelf experiments (Fig. 6a). Furthermore, Fig. 8b shows a divergence in MOC strength between the default PC and default FC simulations with ice shelves. This indicates that changing CO₂ forcing may influence the MOC, but only when marine ice shelves are present in the model (refer to Section 4a3).

Fig. 8 As in Fig. 6, but with results for simulations (dashed lines) when Antarctic ice shelves are kept in the model (unchanging land-ice configuration after 4000 BC).

This excess carbon in the atmosphere cannot be associated with a terrestrial release, because the biosphere and soils actually held more carbon (approximately 10 PgC) in the simulations with ice shelves (Fig. 8c). Furthermore, North Atlantic overturning was slightly stronger (Fig. 8b) for the simulations with ice shelves, which would normally suggest a greater CO₂ uptake efficiency and reduced atmospheric CO₂. However, a more substantial net loss of carbon in the ocean (Figs 8d and 8e) explained the increased atmospheric CO₂ concentrations for the simulations with ice shelves.

The model's spatial output further revealed that this loss was most pronounced in the deep ocean, particularly in the Atlantic basin between 2500 and 5500 m depth. Figure 9 compares temperature and total ocean carbon at 3202 m depth for the BE simulation (NCEP winds) initial conditions with and without ice shelves, and it indicates that deep waters in the Atlantic were significantly (1.5°C) colder and held substantially more carbon in the simulation without marine ice shelves. However, the simulation with ice shelves had much more stagnant NPDW and thus stored more carbon in the North Pacific Ocean (Fig. 9a). As a result, in the simulation with ice shelves, the net loss of ocean carbon from a zonally averaged perspective (not shown) was in deep waters south of the equator, concentrated in the tropical and subtropical South Atlantic Ocean.

Fig. 9 The difference (BE_FC without Antarctic ice shelves - BE_FC with Antarctic ice shelves) at AD 295 and 3202 m depth of the spatial fields for (a) total ocean carbon concentration and (b) temperature. Similar difference fields appear in (c) and (d) but with results from the default initial conditions PC simulations.

The PC simulation (with default initial conditions) was repeated with constant ice shelves after 4000 BC, and a similar profile was produced (Figs 9c and 9d). This suggests that the difference in the distribution of water masses between the two simulations was dependent on ice sheet configuration (not CO₂ concentration) for the range of CO₂ values (245–280 ppm) considered in this study. Furthermore, the MOC states (as implied in Fig. 8b) and characteristics of water masses, both with and without ice shelves, were relatively stable over time, and thus this pattern changed little with increasing southern hemisphere insolation during the late Holocene.

Figure 10 further clarifies the change in water mass characteristics at depth between the two simulations. As discussed in Section 3b, the NADW is characterized by more positive Δ¹⁴C values and is generally more oxygenated, whereas deep waters in the Southern Ocean are less oxygenated and have more negative Δ¹⁴C values. However, Fig. 10 indicates that waters with NADW characteristics extended much deeper (3000–5000 m) in the North and South Atlantic Ocean basins for the simulation with ice shelves. However, without ice shelves, the tongue of NADW was much narrower in the Atlantic Ocean and remained mostly above 3000 m depth, with more extensive colder, carbon-rich, oxygen-deprived AABW at depth.

Fig. 10 The AD 295 vertical profiles of Δ¹⁴C at grid cells in the Atlantic basin (a) in the transient BE_FC simulation without ice shelves after 3000 BC and (b) in the BE_FC with ice shelves (unchanging land-ice configuration after 4000 BC). The blue lines represent profiles in the North Atlantic Ocean: 54.9°N, 41.4°W, 29.7°N, 55.8°W, 0.9°N, 23.4°W; the green lines represent profiles in the South Atlantic and Southern oceans: 29.7°S, 23.4°W, 42.3°S, 30.6°W, 69.3°S, 41.4°W. In (c), the field difference (BE_FC without Antarctic ice shelves - BE_FC with Antarctic ice shelves) is represented at AD 295 for latitudinal averages of oxygen (O₂) concentration with depth.

Taken together, this evidence suggests that AABW is much more prominent in both the North and South Atlantic Ocean basins when ice shelves are not present in the model. More AABW formation in the Weddell Sea would also help explain the slightly faster replacement of NPDW in the simulations without ice shelves (Fig. 9a). When Antarctic ice shelves were kept in the simulation, AABW was limited mostly to the southern hemisphere and deep trenches in the North Atlantic Ocean, and warmer, newer, carbon-poor NADW appeared dominant between 2500 and 3500 m depth. The characteristics of deep water masses in the model simulations with ice shelves thus closely resemble today's deep water distributions.

Greater AABW formation and burial for the no-shelf experiments can be linked to greater sea-ice production and greater sea surface salinity, because shifting the open ocean closer to the Antarctic continent allows a greater influence of downsloping Antarctic winds. When ice shelves were more extensive in the model simulations, overlying winds (and associated heat and moisture fluxes) were weaker, less sea ice was produced (and is generally more variable), and sinking was initiated at lower latitudes. This resulted in a less extensive AABW, a more uniform Antarctic circumpolar water mass, and dominant NADW at depth in much of the Atlantic basin. Because NADW holds less carbon than other water masses, its prominence during the simulation with ice shelves lowers ocean-wide DIC, which may further lead to a long-term increase in pCO₂ in the global ocean and thus decreased oceanic CO₂ solubility, following Galbraith et al. (2007).

Another interesting feature in the simulations with ice shelves is the more stagnant NPDW and greater build-up of carbon in the North Pacific Ocean. This deep water was never flushed over the course of these simulations and has the potential to increase atmospheric CO₂ substantially if it were to be ventilated, as in EQ_2 and transient simulations with land ice fixed at the 6000 BC configuration (Fig. 4). The land ice at this earlier date (Fig. 1a) was slightly more extensive than land ice at 4000 BC (Fig. 1c), further restricting AABW formation to even lower latitudes with less influence from continental downsloping Antarctic winds. Therefore, with very little AABW being generated, the NPDW remained even more stagnant and isolated and warmed slowly through diffusion until it could be ventilated (through the Southern Ocean and equatorial upwelling) in a millennial-scale flushing event. With waters in both the Atlantic and North Pacific basins now newer and carbon-deprived, ocean DIC was even lower in this experiment than the simulations with ice shelves fixed at 4000 BC, and oceanic CO₂ solubility should be reduced even more on longer time scales.

This change in deep water distribution for varying Antarctic ice shelf extents may be significant for carbon cycle comparisons between different interglacial periods. Huybrechts (2002) suggested that the edge of the Eemian (approximately 118000 BC) Ronne-Filchner ice sheet was 50–150 km further south than today. Taking into account sea level changes, insolation, and δ¹⁸O proxies, Pollard and DeConto (2009) used a combined ice sheet–ice shelf model to reconstruct changes in Antarctic marine ice sheets during the past five million years. Their results indicate that there were much less extensive marine ice sheets during the last three interglacial periods, and the Ross and Ronne-Filchner ice sheets virtually disappeared (their reconstruction also correctly reproduced the present more extensive Ross and Ronne-Filchner ice sheets during the Holocene interglacial period).

Our results on deep water distribution, which are independent of orbital forcing changes and CO₂ concentrations but strongly dependent on prescribed ice shelf extent, suggest that this is one potential contributing factor to the higher atmospheric CO₂ values during the late Holocene. Indeed, the model estimate of 5–6 ppm through the ice shelf mechanism by itself would almost entirely account for the 8 ppm contribution that Ruddiman et al. (2011) attribute to the Southern Ocean, and increasing evidence suggests that marine ice shelves during the Holocene were, and remain, unusually extensive relative to many previous interglacial periods.

5 Limitations and future work

It is necessary to underscore the limitations of our model and the methodology that we employed. Our modelling experiments lack representations of coral reefs and peatlands, both of which changed substantially and likely played critical roles in the natural carbon cycle over the course of the Holocene. As for the appropriate spin-up state to use, this has no unique answer because the Earth's climate system was not likely in steady state with orbital, CO₂, and land-ice forcing fixed at 6000 BC, especially given the climate system's abrupt reaction to the approximately 6200 BC (8.2 kyr before present) cooling event/thermohaline weakening prior to this date. This is also especially true for ocean sediments (Fig. 6e), which at 6000 BC had likely been responding to decreasing atmospheric CO₂ since the early Holocene CO₂ maximum (268 ppm) (Broecker et al., 1993). Thus, in order to evaluate the full response of the ocean sediments to calcite compensation, a starting point from the post-glacial CO₂ maximum (or earlier) would be the more appropriate choice.

However, regardless of changes in ocean chemistry, the results presented here do suggest that changes in the ocean carbon reservoir depend strongly on vertical overturning in the ocean. Regarding the MOC itself, the inability of low-resolution OGCMs (in this study 1.8° latitude × 3.6° longitude) to capture downsloping currents may provide a further limitation to the interpretation of sea-ice induced instabilities in the thermohaline circulation (Winton, Hallberg, & Gnanadesikan, 1998).

In all of the simulations described above, the wind forcing (important for ocean surface wind stress and moisture advection) followed the seasonal cycle of the modern climatology. This is a significant caveat in the interpretation of our results, because the wind forcing from 1958 to 1998 is clearly not representative of the actual winds for the entire Holocene. For this reason, the model was unable to generate (for example) a green Sahara or a green Arabian Peninsula during the mid-Holocene. Although GCM time-slice winds have been produced for different periods of the Holocene, because these periods have climate states that do not correspond exactly to those in our model, we tested our results with NCEP winds using the wind parameterization feedback available in the UVic ESCM (Weaver et al., 2001). Sensitivity studies (Figs 6 and 8, right panels) suggest that the wind feedback (for all initial conditions) only leads to a net 1–2 ppm increase in atmospheric CO₂ compared to the NCEP climatology simulations, mostly after AD 1 as a result of increasing overturning in the Southern Ocean.

6 Conclusions

We used the UVic ESCM v. 2.9 to investigate various aspects of the natural carbon cycle from 6000 BC to the pre-industrial period. Because the spin-up (EQ_1) did not reach a steady climate state but instead produced millennial-scale oscillations in the thermohaline circulation, we used different meridional overturning conditions (default, before a flushing event (BE), and during a flushing event (DE)) as initial conditions for transient simulations. These different initial conditions explore the possibility of different deep water mass distributions (Section 3c), with BE initial conditions representing a dominant NADW in the deep Atlantic basin while DE initial conditions start from a state with extensive, newly produced AABW dominating the deep Atlantic basin.

We performed several experiments with atmospheric CO₂ prescribed to observed values (PC simulations), whereas in others the carbon cycle was allowed to freely evolve (FC). For all transient simulations, terrestrial carbon storage increased continuously during the Holocene, placing this model's output in a league with studies such as Kaplan, Prentice, Knorr, and Valdes (2002) and Kleinen et al. (2010), which showed an increase in natural terrestrial carbon storage during the Holocene. In general, the values we obtained (approximately 100 PgC) are in the positive part of the −90 to +370 PgC range of terrestrial carbon uptake summarized in Joos et al. (2004).

Under a steady MOC with the Southern Ocean capped by sea ice and no Antarctic ice shelves, terrestrial carbon sequestration continued during the late Holocene with a corresponding decrease in ocean carbon storage, leading to little net change (on the order of a few parts per million) in atmospheric CO₂. This relatively constant MOC and atmospheric CO₂ is a feature of every UVic free carbon simulation for the late Holocene with and without ice shelves after 4000 BC. This long-term stabilization occurs largely as an indirect consequence of an unchanging land-ice configuration and, consequently, a steady pattern of deep-water distribution.

As discussed in Section 4d, the presence of ice shelves off the coast of Antarctica led to less Antarctic deep water formation, slower circulation of deep waters in the North Pacific, and the dominance of newer, carbon-poor NADW in the Atlantic basin. Significant AABW formation occurred mostly during the centennial-scale ventilation events (the dips in Fig. 6b associated with brief southward retreats in the sea-ice margin), and AABW remained largely constrained to the Southern Ocean. At the same time, atmospheric CO₂ increased during these periods of short-term ventilation in response to ocean carbon losses (mostly from the deep Southern Ocean) and also to climatically induced releases of terrestrial carbon at high latitudes associated with these events. These periods of greater ventilation were accompanied by a weakened NADW (and thus slower uptake), causing the ocean to become a weakened sink or even a brief source of carbon with respect to the atmosphere. Thus, during periods of more localized deep circulation in the Southern Ocean, atmospheric CO₂ tended to increase while oceanic carbon decreased slightly.

The net effect is that the relatively carbon-rich AABW was regionally constrained and better locally ventilated in the transient simulations described in Table 1 before the disappearance of marine ice shelves in the model between 4000 BC and 3000 BC Correspondingly, with less high-carbon AABW in the deep Atlantic, atmospheric CO₂ stabilized at a higher value (257–259 ppm). This suggests that, for more extensive ice shelves and a distribution of deep water masses comparable to their known orientations today, the natural carbon cycle as simulated by the UVic ESCM produced no net change in atmospheric CO₂ after 4000 BC.

On the other hand, the transient simulations without ice shelves after 3500 BC produced more carbon-rich Antarctic deep water that spread globally through the deep ocean. This led to a slightly better-circulated North Pacific (with less deep-water carbon storage in this region) but significantly greater carbon storage in the deep Atlantic basin. Because more carbon remained in the deep ocean globally as a result of the greater prevalence of deep Antarctic-generated waters, this caused CO₂ to stabilize at lower values near 251–252 ppm for the FC and BE_FC simulations. This deep water orientation and carbon storage pattern stayed relatively unchanged during the late Holocene, regardless of the evolution of insolation. Furthermore, the prescribed carbon (PC) simulation without ice shelves after 3500 BC had significantly more AABW formation than a PC simulation with ice shelves (Figs 9c to 9d), suggesting that the deep water distribution is not sensitive to the range of CO₂ values during the Holocene. With the “ice shelf effect” on the global thermohaline circulation being relatively independent of both CO₂ and orbital forcing, we speculate that the Holocene's more extensive Antarctic ice shelves compared to the three previous interglacial periods (Pollard & DeConto, 2009) may be a contributing factor to the higher observed CO₂ concentrations.

Furthermore, the Holocene transient simulation with fixed ice after 6000 BC (in the case of slightly more extensive Antarctic ice shelves from Fig. 1a) produced AABW that was even more local in scale. In this case, NPDW waters stagnated, accumulated carbon and warmed diffusively, while even the mid-latitude South Atlantic Ocean retained a deep and prominent NADW signature. This led to a diffusive warming of the NPDW, which was eventually flushed out through the equatorial Pacific and Southern Oceans when it became vertically unstable in millennial-scale flushing events. Thus, the more localized centennial-scale flushes were replaced with a deep flush of carbon to the atmosphere on 1000-year time scales, which raised atmospheric CO₂ values for more than a millennium.

Although not sensitive to insolation, EQ_3, EQ_2, and Section 4a demonstrate that the periodicity of these millennial-scale events, as well as the “thoroughness” of the flushing of the North Pacific's deep carbon reservoir to the atmosphere, is sensitive to the winds used. In particular, the use of NCEP winds led to longer periods with less frequent flushes whereas the wind feedback adjustment to the climate state suggests more frequent, moderate flushes of NPDW carbon during the late and mid-Holocene. For the transient simulations with transiently evolving ice shelves, the thoroughness of NPDW ventilation around 5000 BC appeared to be dependent both on winds (with the wind feedback leading to more thorough flushes) as well as to some extent on initial conditions (with BE initial conditions producing the most complete ventilation of NPDW). In turn, the much lower deep-ocean DIC in both the North Pacific and the Atlantic Ocean basins for more thorough NPDW ventilations would help reduce ocean CO₂ solubility globally on long time-scales, potentially contributing to an increase in atmospheric CO₂. Therefore, the physical dynamics of the MOC and the source region characteristics of deep water masses highlighted in this study likely feed back to other processes controlling CO₂ solubility in the oceans (such as calcite compensation, shallow water sediments, and temperature changes).

Previous studies have also noted that the efficiency of calcite compensation, SST changes, and human terrestrial release mechanisms might have depended on the thermohaline circulation (Broecker & Ellis, 2007; Ruddiman, 2007). The results presented here indicate that the millennial-scale fluctuations in the MOC, which may only have an amplitude of 2–3 Sv globally but are associated with significantly different distributions of water masses, control atmospheric CO₂ concentrations during the Holocene to a non-negligible degree. Such patterns influence both the release of oceanic and terrestrial carbon to the atmosphere as well as the long-term oceanic absorption of atmospheric CO₂. Schmittner and Galbraith (2008) also highlighted the importance of similar Southern Ocean ventilation patterns in controlling global atmospheric CO₂ levels, although on a much longer time scale than that considered here.

This strong sensitivity to ocean circulation and paleo-ice conditions is a factor not considered in previous simulations using present-day ice, such as in Schurgers et al. (2006), and thus requires further exploration using the UVic ESCM and other models. Our results also underscore the fact that millennial-scale variability and data assimilation of proxy evidence (land ice in this case), two major sources of uncertainty in modelling studies (Crucifix, Loutre, & Berger, 2005), may have played a critical role during the Holocene carbon cycle.

However, with all the simulations considered, the model only produced no net change (approximately 260 ppm) in atmospheric CO₂ since 6000 BC. Therefore, regardless of the initial conditions used or Antarctic ice-sheet extent, no simulations produced a permanent increase in CO₂ over the course of the Holocene as in the observed record. Figure 7 suggests that an external release of 400 PgC into the earth system would be required to reproduce the increasing CO₂ levels (to 280 ppm) during the Holocene. Thus, the inability of the model to conserve carbon for the prescribed CO₂ simulations suggests that some kind of external forcing not represented in the model's carbon cycle after 4000–3000 BC, such as human land use and coral reefs, may be responsible for much of the Holocene's observed CO₂ increase.

Acknowledgements

The support of this work by Natural Sciences and Engineering Research Council (NSERC) Discovery Grants awarded to L.A.M and H.D.M. is gratefully acknowledged. H.D.M. is also thankful for the support of a Canadian Foundation for Climate and Atmospheric Sciences (CFCAS) Project Grant for this research. In addition, we readily thank Eric Galbraith, McGill University, for his helpful comments on a draft of this paper and Alex Matveev, Université du Québec à Montréal, for contributing to the discussion on soil carbon. The anonymous reviewers are also thanked for their constructive input used to improve this paper.