What was the source of the atmospheric CO 2 increase during the Holocene ?

The atmospheric CO2 concentration increased by about 20 ppm from 6000 BCE to the pre-industrial period (1850 CE). Several hypotheses have been proposed to explain mechanisms of this CO2 growth based on either ocean or land carbon sources. Here, we apply the Earth system model MPI-ESM-LR for two transient simulations of climate and carbon cycle dynamics during this period. In the first simulation, atmospheric CO2 is prescribed following icecore CO2 data. In response to the growing atmospheric CO2 concentration, land carbon storage increases until 2000 BCE, stagnates afterwards, and decreases from 1 CE, while the ocean continuously takes CO2 out of the atmosphere after 4000 BCE. This leads to a missing source of 166 Pg of carbon in the ocean–land–atmosphere system by the end of the simulation. In the second experiment, we applied a CO2 nudging technique using surface alkalinity forcing to follow the reconstructed CO2 concentration while keeping the carbon cycle interactive. In that case the ocean is a source of CO2 from 6000 to 2000 BCE due to a decrease in the surface ocean alkalinity. In the prescribed CO2 simulation, surface alkalinity declines as well. However, it is not sufficient to turn the ocean into a CO2 source. The carbonate ion concentration in the deep Atlantic decreases in both the prescribed and the interactive CO2 simulations, while the magnitude of the decrease in the prescribed CO2 experiment is underestimated in comparison with available proxies. As the land serves as a carbon sink until 2000 BCE due to natural carbon cycle processes in both experiments, the missing source of carbon for land and atmosphere can only be attributed to the ocean. Within our model framework, an additional mechanism, such as surface alkalinity decrease, for example due to unaccounted for carbonate accumulation processes on shelves, is required for consistency with ice-core CO2 data. Consequently, our simulations support the hypothesis that the ocean was a source of CO2 until the late Holocene when anthropogenic CO2 sources started to affect atmospheric CO2.

Abstract.The atmospheric CO 2 concentration increased by about 20 ppm from 6000 BCE to the pre-industrial period (1850 CE).Several hypotheses have been proposed to explain mechanisms of this CO 2 growth based on either ocean or land carbon sources.Here, we apply the Earth system model MPI-ESM-LR for two transient simulations of climate and carbon cycle dynamics during this period.In the first simulation, atmospheric CO 2 is prescribed following icecore CO 2 data.In response to the growing atmospheric CO 2 concentration, land carbon storage increases until 2000 BCE, stagnates afterwards, and decreases from 1 CE, while the ocean continuously takes CO 2 out of the atmosphere after 4000 BCE.This leads to a missing source of 166 Pg of carbon in the ocean-land-atmosphere system by the end of the simulation.In the second experiment, we applied a CO 2 nudging technique using surface alkalinity forcing to follow the reconstructed CO 2 concentration while keeping the carbon cycle interactive.In that case the ocean is a source of CO 2 from 6000 to 2000 BCE due to a decrease in the surface ocean alkalinity.In the prescribed CO 2 simulation, surface alkalinity declines as well.However, it is not sufficient to turn the ocean into a CO 2 source.The carbonate ion concentration in the deep Atlantic decreases in both the prescribed and the interactive CO 2 simulations, while the magnitude of the decrease in the prescribed CO 2 experiment is underestimated in comparison with available proxies.As the land serves as a carbon sink until 2000 BCE due to natural carbon cycle processes in both experiments, the missing source of carbon for land and atmosphere can only be attributed to the ocean.Within our model framework, an additional mechanism, such as sur-face alkalinity decrease, for example due to unaccounted for carbonate accumulation processes on shelves, is required for consistency with ice-core CO 2 data.Consequently, our simulations support the hypothesis that the ocean was a source of CO 2 until the late Holocene when anthropogenic CO 2 sources started to affect atmospheric CO 2 .

Introduction
The recent interglacial period, the Holocene, began in about 9700 BCE (Before Common Era) and is characterized by a relatively stable climate.In geological archives, the Holocene is the best recorded period, making it possible to reconstruct changes in climate and vegetation in remarkable detail (e.g.Wanner et al., 2008).Proxy-based reconstructions suggest a decrease in sea surface temperatures in the North Atlantic (Marcott et al., 2013;Kim et al., 2004) simultaneous with an increase in land temperature in western Eurasia (Baker et al., 2017), so that net changes in the global temperature are small.From Antarctic ice-core records, we know that the atmospheric CO 2 concentration increased by about 20 ppm between 5000 BCE and the preindustrial period (Monnin et al., 2004;Schmitt et al., 2012;Schneider et al., 2013).Hypotheses explaining CO 2 growth in the Holocene could be roughly subdivided into oceanand land-based.The ocean mechanisms include changes in carbonate chemistry as a result of carbonate compensation to deglaciation processes (Broecker et al., 1999(Broecker et al., , 2001;;Joos et al., 2004), redistribution of carbonate sedimentation from the deep ocean to shelves, mostly due to coral reef regrowth (Ridgwell et al., 2003;Kleinen et al., 2016), CO 2 degassing due to an increase in sea surface temperatures, predominantly in the tropics (Indermühle et al., 1999;Brovkin et al., 2008), and a decrease in the marine soft tissue pump in response to circulation changes (Goodwin et al., 2011).Using a deconvolution approach based on ice-core CO 2 and δ 13 C data, Elsig et al. (2009) concluded that a significant fraction of Holocene CO 2 changes were attributed to carbonate compensation effects during deglaciation.Recent synthesis of carbon burial in the ocean during the last glacial cycle suggests excessive accumulation of CaCO 3 and organic carbon in the ocean sediments during deglaciation and Holocene (Cartapanis et al., 2018), implicitly supporting the ocean-based mechanism of atmospheric CO 2 growth in response to a decrease in the ocean alkalinity.
The land-based explanations suggest a reduction in natural vegetation cover, such as decreased boreal forests' area and increased desert in North Africa (e.g.Foley, 1994;Brovkin et al., 2002) and anthropogenic land cover changes, mainly deforestation (Ruddiman, 2003(Ruddiman, , 2017;;Kaplan et al., 2011;Stocker et al., 2014Stocker et al., , 2017)).A couple of land processes (CO 2 fertilization, peat accumulation) should have led to terrestrial carbon increase in the Holocene (Yu, 2012;Stocker et al., 2017), complicating the land source explanation.The land-based hypotheses, as well as the ocean soft tissue pump explanation, also have difficulty conforming with the atmospheric δ 13 CO 2 reconstructions that show no substantial changes in the Holocene (Schmitt et al., 2012;Schneider et al., 2013), contrary to the expectation of a significant atmospheric δ 13 CO 2 decrease due to the release of isotopically light biological carbon.The deconvolution approach by Elsig et al. (2009) suggested a land carbon uptake of 290 GtC between 9000 and 3000 BCE.For a detailed overview of process-based hypotheses, see, for example, Brovkin et al. (2016).
While many numerical experiments have been done to test the above-mentioned hypotheses with intermediate complexity models, this class of models is limited in spatial and temporal resolution and consequently does not resolve well climate patterns and variability.Here, we apply the full-scale Earth system model MPI-ESM-LR for simulations of the coupled climate and carbon cycle during the period from 6000 BCE to 1850 CE.The focus of this paper is on the carbon cycle dynamics for terrestrial and marine components, while changes in climate are considered in the companion paper by Bader et al. (2019).

Methods
The MPI-ESM-1.2LRused in the Holocene simulations consists of the coupled general circulation models for the atmosphere and the ocean, ECHAM6 (Stevens et al., 2013) and MPIOM (Jungclaus et al., 2013), respectively, the land surface model JSBACH (Raddatz et al., 2007;Reick et al., 2013), and the marine biogeochemistry model HAMOCC5 (Ilyina et al., 2013).In comparison to the MPI-ESM-LR model used in the Climate Model Intercomparison Project Phase 5 (CMIP5) simulations (Giorgetta et al., 2013), the model has been updated with several new components for the land hydrology and carbon cycle.In addition to the previously developed dynamic vegetation model (Brovkin et al., 2009), the new JSBACH component includes the soil carbon model YASSO (Goll et al., 2015), a five-layer hydrology scheme (Hagemann and Stacke, 2015), and an interactive albedo scheme (Vamborg et al., 2011).HAMOCC was updated with prognostic nitrogen fixers (Paulsen et al., 2017) and a parameterization for a depth-dependent detritus settling velocity and remineralization of organic material (Martin et al., 1987).
We used a combination of several forcings as boundary conditions for transient simulations (Fig. 1).The orbital forcing (Fig. 1a) follows a reconstruction by Berger (1978).We use a solar irradiance forcing that was reconstructed from 14 C tree-ring data that correlate with the past changes in the solar open magnetic field (Krivova et al., 2011).CO 2 , N 2 O, and CH 4 forcings stem from ice-core reconstructions (Fortunat Joos, personal communication, 2016, Fig. 1b,c); see also comment by Köhler (2019).We applied a new reconstruction of volcanic forcing (Bader et al., 2019) based on the GISP2 ice-core volcanic sulfate record (Zielinski et al., 1996) (Fig. 1d).To maintain consistency with the CMIP5 simulations, we included the land use changes based on the LUH1.0 product by Hurrt et al. (2011) and Pongratz et al. (2011), provided for the period after 850 CE.To avoid an abrupt change in the land use forcing at 850 CE, we interpolated the land use map from 850 CE linearly backwards to no land use conditions at 150 BCE (Fig. 1e).The land use areas in 1750 CE are in line with the most recent update of the HYDE dataset (Klein Goldewijk et al., 2017), although our interpolation underestimates crop and pasture areas earlier than 1750 CE.
The spin-up simulation 8KAF started from initial conditions for the pre-industrial climate and continued with boundary conditions for 6000 BCE for more than 1000 years in order to establish an equilibrium of the climate and carbon cycle with the boundary conditions.During the spin-up period the atmospheric CO 2 concentration was kept at a constant level of 260 ppm.The weathering flux was not changed compared to pre-industrial conditions.The change in the boundary conditions to 6000 BCE led to slightly different sedimentation fluxes, which resulted in a slow decline of alkalinity in 8KAF.Afterwards, the model was run with an interactive carbon cycle to ensure a dynamic equilibrium between land, ocean, and atmospheric carbon cycle components (simulation 8KAFc).In the 8KAFc simulation, the equilibration procedure in HAMOCC followed the CMIP5 spin-up procedure (Ilyina et al., 2013): Throughout the equilibration process, weathering fluxes and CaCO 3 content in sediments have been changed, which led to changes in total alkalinity (TA).This would have occurred naturally, without leading to excess TA, had the biogeochemistry model been given a long enough spin-up time to equilibrate its sediments.Along with change in TA, also DIC changed (with the molar ratio 2 : 1) in order to maintain the correct pCO 2 .
In the interactive carbon spin-up, the model firstly used the same weathering rates as in 8KAF for ∼ 300 years, then it stabilized the system by increasing all the weathering fluxes (Si, organic matter, CaCO 3 ), which led to a stabilization of the surface alkalinity; afterwards the alkalinity changed to keep the target pCO 2 .For the last few hundred years, weathering was adjusted, which led to alkalinity stabilization.In total, the 8KAFc spin-up took more than 1000 model years.
We performed two transient simulations from 6000 BCE to 1850 CE, a commonly defined onset of the industrial period.The first transient simulation, TRAF, was initiated from the end of the spin-up simulation 8KAF.In the TRAF simulation, the atmospheric CO 2 concentration is prescribed from ice-core reconstructions.Consequently, land and ocean carbon uptakes do not interfere with each other, and the total mass of carbon in the land-ocean-atmosphere system is not conserved.
The second transient simulation, TRAFc, started from the end of the 8KAFc simulation in the interactive climatecarbon cycle mode.Using equilibrium initial conditions and boundary conditions for the carbon cycle does not guarantee that the interactively simulated atmospheric CO 2 concentration will follow the trend reconstructed from ice cores.If the simulated atmospheric CO 2 concentration differs substantially from reconstructed data, both climate and carbon cycle components deviate from results which would be obtained if the model were driven by the reconstructed CO 2 forcing, and these biases complicate the comparison of trends between the model and observations.To ensure that simulated CO 2 is close to the reconstructed time series, in the TRAFc simulation we used a CO 2 nudging technique following the approach of Gonzales and Ilyina (2016), targeting the atmospheric CO 2 record from ice-core reconstructions.If the simulated atmospheric CO 2 dropped below the target, the surface ocean total alkalinity and dissolved inorganic carbon (DIC) concentrations were decreased in a 2 : 1 ratio to mimic the process of CaCO 3 sedimentation.If CO 2 was higher than the target, the alkalinity and DIC were not changed.This alkalinity-forced approach is supported by evidence of decreasing deep ocean carbonate ion concentration in the course of the Holocene (Yu et al., 2014) and excessive coral reef buildup on shelves during deglaciation (Opdyke and Walker, 1992;Vecsei and Berger, 2004) as well as recent synthesis by Cartapanis et al. (2018).The alkalinity decline resembles excessive shallow-water carbonate sedimentation such as coral reefs, which are not included in HAMOCC.
HAMOCC includes a module of sediment processes (Heinze et al., 1999).Interactive simulation of the sediment porewater chemistry and accumulation of solid sediment components, such as CaCO 3 , particulate organic carbon (POC), opal, and clay, is a necessary condition to calculate changes in the ocean biogeochemistry on millennial timescales.To compensate for the POC, CaCO 3 , and opal losses due to sedimentation, the fluxes into the sediment over the last 300 years of the spin-up runs were analysed, and a globally uniform weathering input for silicate, alkalinity, nutrients in the form of organic matter, and dissolved inorganic carbon was prescribed.Note that the weathering flux calculated using this approach is sensitive to changes in the model setups (prescribed versus interactive CO 2 ).This explains the difference between weathering fluxes in TRAF and TRAFc experiments (Table 1).The weathering flux is not accounted for in the land compartment changes (second column).
3 Results and discussion

Global response
The setup of the TRAF simulation resembles experiments performed in CMIP5 with CO 2 concentrations prescribed from Representative Concentration Pathway (RCP) scenarios.Similar to the RCP simulations, changes in carbon pools on land and in the ocean could be estimated from landatmosphere and ocean-atmosphere fluxes, respectively.Cumulative changes in the ocean and land CO 2 fluxes reveal that, by the end of the simulation, the ocean is a sink of 152 PgC, while the land is a source of 39 PgC (Fig. 2a, Table 1).Accounting for an increase in the atmospheric carbon pool by 53 PgC, the total carbon budget has a deficit of 166 PgC by 1850 BCE (Fig. 3).Assuming that the CO 2 airborne fraction on a millennial timescale is about one-sixth (Maier-Reimer and Hasselmann, 1987), or even less if we account for the land response due to CO 2 fertilization, the atmospheric CO 2 concentration by the end of the simulation would be 11-13 ppm less than observed (286 ppm).Therefore, carbon budget changes in the TRAF experiment imply that other boundary conditions are necessary to obtain the amplitude of the simulated CO 2 concentration trend as reconstructed from the ice cores.Carbon cycle changes in the interactive CO 2 simulation TRAFc are shown in Fig. 2b.At the beginning of the simulation (during 6 to 5000 BCE), atmospheric CO 2 and land carbon fluctuate with an amplitude of several PgC, while the ocean becomes a small source of CO 2 to the atmosphere.Between 5000 and 2000 BCE, atmospheric carbon storage increases by about 30 PgC, and the land takes about 60 PgC due to CO 2 fertilization, while the ocean releases ca.90 PgC.Between 2000 BCE and 1 CE (start of the Common Era; note that the year 0 does not exist in the CE system), land is a source of 10 PgC to the atmosphere.After 1 CE, land carbon losses accelerate due to land use changes, and by 1850 CE land carbon decreases by an additional 60 PgC.In 1850 CE, the land and ocean are sources of 17 and 39 PgC, respectively, while the atmosphere gains 56 PgC.The atmospheric minimum in CO 2 around 1600 CE apparent in the reconstruction is not reproduced by the model, confirming that the abrupt uptake of carbon by land or ocean is difficult to attribute to internal variability in the coupled climate-carbon system (Pongratz et al., 2011), and external forcing scenarios, such as an abrupt reforestation of tropical America (Kaplan et al., 2011), might need to be accounted for.
Decadal-scale excursions in ocean and land carbon storages in Fig. 2b are mainly explained by responses to surface cooling resulting from volcanic eruptions.The most visible example of this CO 2 response is during the period around 3200 BCE, when reconstructed aerosol optical depth shows an enhancement which is moderate in magnitude but is of long duration (Fig. 1d), potentially resulting from a longduration high-latitude eruption or from contamination of the volcanic record by biogenic sulfate (Zielinski et al., 1994).In response to this applied forcing, the land takes up carbon due to decreased respiration (see, for example, Brovkin et al., 2010;Segschneider et al., 2013), while the alkalinity adjustment in the ocean counteracts the land carbon uptake, leading to carbon release from the ocean.After a few decades, the land turns into a source of carbon due to reduced productivity, ocean carbon uptake restores, and atmospheric CO 2 reveals a spike due to an excessive land source.At 3000 BCE, this spike ceases, and the simulated CO 2 continues to fluctuate around the ice-core time series.Although the volcanic forcing included in the simulations at 3200 BCE is likely an overestimate, this case illustrates the response of the climatecarbon system to an extreme volcanic aerosol forcing, which leads to pronounced cooling of the land and ocean surfaces.
Changes in carbon density on land and in the ocean in the course of both TRAF (not shown) and TRAFc simulations reveal complex patterns (Fig. 4).In the ocean, the vertically integrated DIC is decreasing everywhere, causing negative change patterns dominating in the Southern Ocean and North Pacific.Carbon sedimentation is high in upwelling zones, mainly in coastal areas and the tropical Pacific, and that causes strong accumulation patterns.These sedimentation patterns are typical for the HAMOCC model with interactive sediments (Heinze et al., 1999); they are generally well comparable with observed sedimentation patterns for organic carbon and CaCO 3 (Seiter et al., 2004;Archer, 1996).The land has a mixed pattern of increased carbon density, mostly in South America and in central North America, with decreased densities in Africa and East Asia.This is caused by an interplay between climate, CO 2 , and land use effects on soil and biomass storages.
Changes in the carbon budget components over the experimental period are provided in Table 1 and Fig. 3.For the atmosphere, the difference between the TRAF and TRAFc simulations is minor (3 PgC).For the land, a difference of 22 PgC is caused mainly by the relatively higher CO 2 concentration in the TRAFc simulation, especially during the period of lower CO 2 around 1600 CE, due to the CO 2 fertilization effect on the plant productivity.The ocean-atmosphere cumulative fluxes (−152 and 39 PgC for TRAF and TRAFc, respectively) are minor in comparison with the ocean carbon budget components, and the difference of 191 PgC is explained by the applied surface alkalinity removal in the TRAFc simulation.The carbon inventory of the water column that predominantly includes dissolved inorganic carbon (DIC) loses 1324 and 1799 PgC in the TRAF and TRAFc runs, respectively.Sediments accumulate more than 3500 PgC in the form of CaCO 3 and organic carbon, mainly compensated by the weathering flux from land.In the TRAFc experiment, 1224 PgC was removed from the ocean surface in the form of CaCO 3 , effectively reducing the weathering flux (3270 PgC) to a scale below the TRAF experiment (2137 PgC).In total, despite large changes in the cumulative fluxes of weathering and sedimentation, the net cumulative ocean-atmosphere flux is minor.

Land carbon and vegetation
Natural changes in vegetation and tree cover are most pronounced for the period before 1 CE, before the start of substantial land use forcing.Comparing with 6000 BCE, vegetation cover becomes much less dense in Africa, mainly due to decreased rainfall in response to the decreasing summer radiation in the Northern Hemisphere (Fig. 5a).Boreal forests moved southward in both North America and Eurasia (Fig. 5b).The southward shift of vegetation in North Africa, and of the treeline in Eurasia from the mid-Holocene to the pre-industrial period, as well as the increase in vegetation and tree cover in central North America, is in line with pollen evidence (Prentice et al., 2000).The southward retreat of the boreal forest in North America is much less pronounced than in Eurasia (Fig. 5b).This is also in line with reconstructions, as there is no evidence for a significant shift of the treeline in North America (Bigelow et al., 2003), likely due to the cooling effects of the remains of the Laurentide ice sheet, which is not accounted for as a forcing in our simulations.
Simulated changes in vegetation cover are reflected in the carbon density changes (Fig. 6a).From 6000 BCE to 1 CE, the carbon density decreases in North Africa, East Asia, northern South America and above 60 • N slightly in Eurasia.In most of the rest of the land ecosystems, the carbon density increases, mostly due to CO 2 fertilization effects as the atmospheric CO 2 concentration increases by about 15 ppm by 1 CE.A strong increase in the Southern Hemisphere and central North America is also due to increased vegetation density.After 1 CE, land carbon declines due to land use changes, predominantly deforestation (Fig. 6b).Patterns of carbon decrease after 1 CE reflect land use patterns except in South America, southern Africa, and central North America.The simulated increase in land carbon storage before 2000 CE and decrease afterwards is consistent with the changes in atmospheric δ 13 CO 2 (Elsig et al., 2009;Schmitt et al., 2012).The deconvolution approach (Fig. 3 in Elsig et al., 2009) resulted in the land uptake of ca.140 PgC from 6000 to 3000 BCE, divided rather equally into ca.70 PgC from 6000 to 5000 BCE and 70 PgC between 5000 and 3000 BCE.The 70 PgC uptake from 6000 to 5000 BCE deduced from the increase in atmospheric δ 13 CO 2 is not reproduced in our experiments, likely because it is a non-equilibrium land response which can be captured only in transient simulations during the last deglaciation.In our TRAF and TRAFc experiments, land accumulates about 60 PgC between 5000 and 2000 BCE, comparable with land uptake of 70 PgC between 5000 and 3000 BCE inferred by Elsig et al. (2009).We can conclude that after 5000 BCE, the land carbon dynamics in MPI-ESM (uptake of 60 PgC by 2000 BCE, release of 80-100 PgC by 1850 CE, predominantly due to land use) is similar to the land carbon changes estimated by Elsig et al. (2009).

Ocean carbon
Simulated physical ocean fields, including sea surface temperatures and the Atlantic meridional overturning, do not change substantially in the Holocene.The main reason for the declining carbon storage in the water column (Fig. 4, Table 1) is a decrease in ocean alkalinity (Figs. 7a;8a).This is explained by the applied surface ocean alkalinity forcing and also by the response of the ocean carbonate chemistry to changes in carbonate production.The global CaCO 3 export from the surface to the aphotic layer increases by about 5 % between 6000 and 2000 BCE in both TRAF and TRAFc simulations and returns to the 6000 BCE level by the end of the simulation.Comparing TRAF and TRAFc simulations, the difference in the globally averaged ocean alkalinity in these  quired excessive carbonate sedimentation in the shallow waters would be 3 Tmol yr −1 in TRAFc relative to TRAF, or at the lower bound of estimates of 3.35 to 12 Tmol yr −1 CaCO 3 accumulation proposed by Vecsei and Berger (2004) and Opdyke and Walker (1992).Even corresponding excessive carbonate sedimentation of 11.2 Tmol yr −1 CaCO 3 in the TRAFc simulation would fall into this observational range, although at the higher bound.Let us note that in the 8KAF and TRAF experiments the weathering was not adjusted to changes in boundary conditions, and this likely caused surface alkalinity decrease in the transient run (Fig. 7c).In particular, this alkalinity drift in TRAF explains the initial decrease in the ocean carbon storage until ca.4500 BCE (Fig. 2a), despite an increase in atmospheric CO 2 concentration.If TRAF had started from an equilibrated system as TRAFc did, the beginning of TRAF would have been more similar to TRAFc, and the ocean carbon uptake would have started earlier.
Besides the estimate of applied carbonate accumulation forcing, another way to address the plausibility of simulated alkalinity trends is to compare changes in the carbonate ion concentration ([CO   2013) reveals a significant difference between TRAF and TRAFc in the Atlantic (Fig. 7b).Decrease in [CO =  3 ] in the TRAFc simulation is more significant than in the TRAF experiment, presumably due to a stronger decrease in alkalinity in the former simulation.Interestingly, changes in [CO = 3 ] at the Pacific site are not significant in both simulations, while the data propose a slight decrease in carbonate ion concentration.The difference between the Atlantic and Pacific responses is visible in Fig. 8.In both experiments, simulated changes in [CO =  3 ] in the Atlantic and Southern Ocean are stronger than in the Indo-Pacific.At a depth of 4 km, comparable with the depth of the cores by Broecker and Clark (2007), changes in the tropical oceans in both simulations are in the range of 0-15 mol m −3 .Changes in the TRAFc experiment are more pronounced than in TRAF due to stronger changes in alkalinity.As expected, [CO = 3 ] changes are more pronounced for depths of 3 km than for 4 km (Fig. 8).Comparison of simulated ocean carbon budget with recent carbon data synthesis (Cartapanis et al., 2018) is shown in Table 2.In comparison with mean data values, CaCO 3 burial in the model (27.9-29.1,plus 13 TmolC yr −1 surface removal in the TRAFc experiment) is higher than in the data (23.3TmolC yr −1 ); however, this is compensated by a higher modelling weathering rate (24.6-34.7 TmolC yr −1 ) compared to 11.7 TmolC yr −1 in the data.The model values are at the upper end of the data uncertainty range (11-38 TmolC yr −1 ; see min-max range for CaCO 3 sedimentation in Table 2).Organic carbon burial in the model (10.5-11.3Tmol yr −1 ) is less than in the averaged data (18.3Tmol yr −1 ); however, the uncertainty in the burial is so (6-58 Tmol yr −1 ) that the model values are almost twice more than the lower end of the data range.
An important question is whether the ocean on average lost carbon content over the Holocene.In equilibrium, volcanic CO 2 outgassing (not accounted for explicitly in the model), both aerial and submarine, compensates for weathering; therefore for proper comparison it needs to be included in the table as the ocean-atmosphere budget.In that case, averaged ocean water column losses in the data are 20.8TmolC yr −1 .Similar to the data, the model shows the loss of carbon from the water column, 13.8 and 18.7 TmolC yr −1 in the TRAF and TRAFc experiments, respectively.Therefore, simulated ocean carbon losses are qualitatively (and even quantitatively for TRAFc) in line with observations.

Limitation of the model setup
There are certain limitations of the carbon cycle models used in the study.Firstly, the applied version of JSBACH does not include wetland and peatland processes.If the Holocene peat accumulation of several hundred gigatonnes of carbon (GtC; Yu, 2012) were accounted for, the land would be a stronger sink of carbon during 6000 to 2000 BCE.This might require an even stronger ocean source.On the other hand, we neglect other sources of atmospheric CO 2 which might at least partly compensate for the peatland growth, for example, emissions due to ongoing thermokarst formation and erosion of permafrost soils, especially close to the Arctic coast (Lindgren et al., 2018).Secondly, HAMOCC does not include coral reefs as a process-based component.This is one of the reasons why the surface alkalinity was forced directly in the TRAFc simulation.Thirdly, the applied versions of JSBACH and HAMOCC do not simulate carbon isotope changes, in particular 13 C changes.For land carbon, the increase in the carbon storage on land by 50-60 PgC by 2000 BCE and its decrease by 80 PgC by 1850 CE would be translated into a ca.0.05 ‰ decrease in atmospheric δ 13 CO 2 .This small change is within the uncertainty bounds of δ 13 CO 2 reconstructed from ice cores (Elsig et al., 2009;Schmitt et al., 2012).For ocean carbon, carbonate changes would not significantly modify the ocean and atmospheric δ 13 CO 2 content.Simulated changes in biological production and export flux might have affected the atmospheric δ 13 CO 2 , but the scale will likely be small.
In addition, two limitations are intrinsic to the setups of spin-up and transient simulations.Firstly, assuming that all carbon cycle components are initially in equilibrium with boundary conditions is a simplification.Changes in climate due to slowly changing boundary conditions, such as orbital or greenhouse gas forcing, are occurring on timescales similar to long-term processes in the carbon system (soil buildup on land, carbonate compensation in the ocean).Therefore, the carbon cycle is never in full equilibrium, and memory in the carbon cycle processes, due for example to carbonate compensation in the ocean during deglaciation, affects the carbon dynamics afterwards.A proper way to account for the memory effect is to set spin-up simulations millennia before the Holocene, e.g. at the last glacial maximum (19 000 BCE), and perform a transient deglaciation simulation.Simulations with intermediate complexity models suggested that the impact of the memory effect from deglaciation on Holocene carbon dynamics, in particular due to carbonate compensation, is significant (e.g.Menviel and Joos, 2012).However, In the prescribed CO 2 experiment TRAF, the ocean is a sink of carbon.This strengthens the argument that neither changes in circulation nor in sea surface temperatures are capable of explaining CO 2 growth in the Holocene.In the coupled land-ocean-atmosphere system, there is a total deficit of 166 PgC by the end of the experiment.The TRAFc sim-ulation with interactive CO 2 is performed in the nudging mode: we use the surface alkalinity changes as the forcing for ocean-atmosphere CO 2 flux.In response to this forcing, the ocean serves as a source of carbon over the Holocene.The alkalinity decline is within the bounds of proposed changes in the carbonate sedimentation in shallow waters and consistent with available proxies for carbonate ion decrease in the deep sea.
There are several limitations of our simulations related to initial conditions and forcings.We cannot simply overcome them by repeating runs in a different setup or by doing additional sensitivity experiments due to the high computational costs of full-scale ESMs.Despite these limitations, we can make several conclusions on the potential source of CO 2 to the atmosphere during the last 8000 years.Regarding the land source, experiments demonstrate that natural carbon dynamics led to an increase in the land carbon storage during the first half of the simulation (until 2000 BCE).This is in line with previous simulations performed with intermediate complexity models (e.g.Kaplan et al., 2002;Kleinen et al., 2016) and with ice-core deconvolution studies (Elsig et al., 2009;Schmitt et al., 2012).During 6000 to 2000 BCE, the atmospheric CO 2 increase is about two-thirds of the estimated 20 ppm increase.Although the TRAF and TRAFc simulations do not account for land use changes during this period, assuming that land use was a source of carbon to the atmosphere requires about 100 PgC to compensate for natural land, ocean, and atmospheric carbon content increase during this time.If we account for the peat carbon accumulation (neglected in TRAFc), emissions from land use would need to be higher (about 200 PgC over the period 6000 to 2000 BCE).This is not absolutely impossible (Kaplan et al., 2011), but such a high-end land use emission scenario for the end of the Neolithic period, when agriculture was not yet widespread in Europe and America, is rather unlikely.
Regarding the ocean source, both TRAF and TRAFc simulations show a decrease in ocean alkalinity.Even if this decrease is a result of a drift in the carbonate system due to imperfect initialization of the balance between sedimentation and weathering, in both simulations the model is capable of producing a decrease in the carbonate ion concentrations in the Atlantic, which is in the direction proposed by proxy data (Fig. 8b).The magnitude of the decrease in the TRAF experiment is underestimated compared to the proxy data, while it is in line with the data for the TRAFc experiment.As land serves as a carbon sink until 2000 BCE due to natural (nonanthropogenic) carbon cycle processes in both experiments, the missing source of carbon for land and atmosphere could only be attributed to ocean.Within our model framework, an additional mechanism is required for consistency with icecore CO 2 data, such as surface alkalinity decrease, due, for example, to unaccounted for carbonate accumulation processes on shelves supported by observational evidence.Finally, our simulations support the hypothesis that the ocean was a source of CO 2 until the late Holocene when anthropogenic CO 2 sources started to affect atmospheric CO 2 .

Figure 2 .
Figure 2. Changes in the cumulative fluxes for major carbon cycle components (land, ocean, and atmosphere) (PgC), from 6000 BCE to 1850 CE in the TRAF (a) and TRAFc (b) simulations.The colour legend is as follows: land is denoted in cyan, ocean in blue, atmosphere in magenta, and ice-core reconstruction in orange.

Figure 3 .
Figure 3. Changes in the carbon cycle compartments (land, ocean, and atmosphere) (PgC), and cumulative fluxes between them (PgC), from 6000 BCE to 1850 CE in the TRAF (a) and TRAFc (b) simulations.Flux from the lithosphere include weathering (TRAF) and weathering minus shallow-water sedimentation (TRAFc).Carbon budget has a negative 166 Pg disbalance in the TRAF simulation, while it is closed in the TRAFc run.

Figure 4 .
Figure 4. Combined map of changes in the ocean carbon storage (vertically integrated ocean water column plus sediments minus weathering) and land (soil plus vegetation) at the end of the TRAFc simulation (1850 CE) relative to 6000 BCE (in kgC m −2 ).

Figure 5 .
Figure 5. Change in vegetation fraction (a) and tree cover fraction (b) at 1 CE relative to 6000 BCE in the TRAFc simulation.

Figure 6 .
Figure 6.Change in land carbon density (kg C m −2 ), relative to 6000 BCE at 1 CE (a) and at 1850 CE relative to 1 CE (b) in the TRAFc simulation.
= 3 ]) in the deep Atlantic and Pacific Ocean with available reconstructions of carbonate ion concentrations.Using the benthic foraminiferal B/Ca proxy for deep water [CO = 3 ], Yu et al. (2014) found that [CO = 3 ] in the deep Indian and Pacific Ocean declined by 5-15 µmol kg −1 during the Holocene.Broecker et al. (1999) and Broecker

Figure 7 .
Figure 7. (a) 100-year moving average of global surface alkalinity, µmol kg −1 .(b) 100-year moving average of carbonate ion concentration (µmol kg −1 ), averaged for nine neighbouring model grid cells centred in the deep Atlantic (12 • N, 60 • W; 3400 m) and deep Pacific (1 • S, 160 • W; 3100 m).The data (circles) are [CO = 3 ] for sites VM28-122 and GGC48 reconstructed by Yu et al. (2013), which appropriately correspond to the ocean grid cells accounting for the model bathymetry mask.Data uncertainties (1σ ) reported for the Atlantic by Yu et al. (2013) are indicated by whiskers.35 ‰ salinity is used for model unit conversion from m −3 to kg −1 .

Figure 8 .
Figure 8. Differences in carbonate ion concentration (mmol m −3 ), between TRAFc and TRAF simulations in 1850 CE at the depth of 3 km (a) and 4 km (b).

Table 1 .
Changes in compartments and cumulative fluxes at 1850 CE relative to 6000 BCE (PgC).

Table 2 .
Ilyina et al., 2013;Heinze et al., 2016)fluxes from 6000 BCE to 1850 CE (Tmol yr −1 ).Including aerial volcanic CO 2 outgassing.bPre-industrialfluxesaccording to the Fig.1inCartapanis et al. (2018).cSubtractingvolcanic outgassing.transientdeglaciationrun is presently too challenging for full-scale ESMs due to high computational costs.Secondly, the weathering fluxes are assumed to be constant during the transient simulation.While this is a common practice for ocean biogeochemistry simulations (e.g.Ilyina et al., 2013;Heinze et al., 2016), it results in a mismatch between weathering and sedimentation under changing boundary conditions in transient simulations.As land climate evolves, weathering fluxes are changing due to their dependence on runoff and temperature.This causes a shift in land-ocean fluxes of carbon, alkalinity, and nutrients, leading to inventory changes and a possible drift in ocean-atmosphere fluxes.In particular, as POC fluxes to the sediment are not properly compensated by the fixed weathering, this leads to changes in nutrient inventory in transient simulations.Both CaCO 3 and POC fluxes to sediments are changing with time; this also leads to changes in the rain ratio.In the absence of factorial experi- a