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

The atmospheric CO₂ 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₂ 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₂ is prescribed following ice-core CO₂ data. In response to the growing atmospheric CO₂ concentration, land carbon storage increases until 2000 BCE, stagnates afterwards, and decreases from 1 CE, while the ocean continuously takes CO₂ 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₂ nudging technique using surface alkalinity forcing to follow the reconstructed CO₂ concentration while keeping the carbon cycle interactive. In that case the ocean is a source of CO₂ from 6000 to 2000 BCE due to a decrease in the surface ocean alkalinity. In the prescribed CO₂ simulation, surface alkalinity declines as well. However, it is not sufficient to turn the ocean into a CO₂ source. The carbonate ion concentration in the deep Atlantic decreases in both the prescribed and the interactive CO₂ simulations, while the magnitude of the decrease in the prescribed CO₂ 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 CO₂ data. Consequently, our simulations support the hypothesis that the ocean was a source of CO₂ until the late Holocene when anthropogenic CO₂ sources started to affect atmospheric CO₂.

1 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₂ concentration increased by about 20 ppm between 5000 BCE and the pre-industrial period (Monnin et al., 2004; Schmitt et al., 2012; Schneider et al., 2013). Hypotheses explaining CO₂ growth in the Holocene could be roughly subdivided into ocean- and land-based. The ocean mechanisms include changes in carbonate chemistry as a result of carbonate compensation to deglaciation processes (Broecker et al., 1999, 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₂ 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₂ and δ¹³C data, Elsig et al. (2009) concluded that a significant fraction of Holocene CO₂ 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₃ and organic carbon in the ocean sediments during deglaciation and Holocene (Cartapanis et al., 2018), implicitly supporting the ocean-based mechanism of atmospheric CO₂ 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, 2017; Kaplan et al., 2011; Stocker et al., 2014, 2017). A couple of land processes (CO₂ 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 δ¹³CO₂ 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 δ¹³CO₂ 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).

2 Methods

The MPI-ESM-1.2LR used 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 ¹⁴C tree-ring data that correlate with the past changes in the solar open magnetic field (Krivova et al., 2011). CO₂, N₂O, and CH₄ 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.

Figure 1 Time series of applied forcings: (a) June–July–August insolation at 60^∘ N (W m⁻²), (b) atmospheric CO₂ concentration (ppm), and (c) N₂O and CH₄ concentrations (ppb). (d) Aerosol optical depth of volcanic eruptions and (e) global area of crop and pastures (10⁶ km²).

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

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₃), which led to a stabilization of the surface alkalinity; afterwards the alkalinity changed to keep the target pCO₂. 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₂ 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 climate–carbon cycle mode. Using equilibrium initial conditions and boundary conditions for the carbon cycle does not guarantee that the interactively simulated atmospheric CO₂ concentration will follow the trend reconstructed from ice cores. If the simulated atmospheric CO₂ 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₂ forcing, and these biases complicate the comparison of trends between the model and observations. To ensure that simulated CO₂ is close to the reconstructed time series, in the TRAFc simulation we used a CO₂ nudging technique following the approach of Gonzales and Ilyina (2016), targeting the atmospheric CO₂ record from ice-core reconstructions. If the simulated atmospheric CO₂ 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₃ sedimentation. If CO₂ 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₃, 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₃, 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₂). This explains the difference between weathering fluxes in TRAF and TRAFc experiments (Table 1).

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

The weathering flux is not accounted for in the land compartment changes (second column).

3 Results and discussion

3.1 Global response

The setup of the TRAF simulation resembles experiments performed in CMIP5 with CO₂ 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 land–atmosphere and ocean–atmosphere fluxes, respectively. Cumulative changes in the ocean and land CO₂ 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₂ 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₂ fertilization, the atmospheric CO₂ 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₂ concentration trend as reconstructed from the ice cores.

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.

Carbon cycle changes in the interactive CO₂ simulation TRAFc are shown in Fig. 2b. At the beginning of the simulation (during 6 to 5000 BCE), atmospheric CO₂ and land carbon fluctuate with an amplitude of several PgC, while the ocean becomes a small source of CO₂ 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₂ 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₂ 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₂ 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 long-duration 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₂ reveals a spike due to an excessive land source. At 3000 BCE, this spike ceases, and the simulated CO₂ 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 climate–carbon 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₃ (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₂, 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₂ concentration in the TRAFc simulation, especially during the period of lower CO₂ around 1600 CE, due to the CO₂ 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₃ 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₃, 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.

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.

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

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⁻²).

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

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₂ fertilization effects as the atmospheric CO₂ 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 δ¹³CO₂ (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 δ¹³CO₂ 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).

Figure 6 Change in land carbon density (kg C m⁻²), relative to 6000 BCE at 1 CE (a) and at 1850 CE relative to 1 CE (b) in the TRAFc simulation.

3.3 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₃ 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 two simulations by 1850 CE is 35 µmol kg⁻¹, similar to the difference in surface alkalinity changes shown in Fig. 7a. Accounting for 7850 years of experimental length, the alkalinity loss corresponds to 8.2 and 11.2 Tmol yr⁻¹ CaCO₃ sedimentation in TRAF and TRAFc simulations, respectively. The required excessive carbonate sedimentation in the shallow waters would be 3 Tmol yr⁻¹ in TRAFc relative to TRAF, or at the lower bound of estimates of 3.35 to 12 Tmol yr⁻¹ CaCO₃ accumulation proposed by Vecsei and Berger (2004) and Opdyke and Walker (1992). Even corresponding excessive carbonate sedimentation of 11.2 Tmol yr⁻¹ CaCO₃ 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₂ 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.

Figure 7 (a) 100-year moving average of global surface alkalinity, µmol kg⁻¹. (b) 100-year moving average of carbonate ion concentration (µmol kg⁻¹), 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 [${\mathrm{CO}}_{\mathrm{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⁻³ to kg⁻¹.

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

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 ([${\mathrm{CO}}_{\mathrm{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 [${\mathrm{CO}}_{\mathrm{3}}^{=}$], Yu et al. (2014) found that [${\mathrm{CO}}_{\mathrm{3}}^{=}$] in the deep Indian and Pacific Ocean declined by 5–15 µmol kg⁻¹ during the Holocene. Broecker et al. (1999) and Broecker and Clark (2007) suggested a similar amplitude of [${\mathrm{CO}}_{\mathrm{3}}^{=}$] changes in the deep Atlantic. Comparison of changes in [${\mathrm{CO}}_{\mathrm{3}}^{=}$] in TRAF and TRAFc simulations with [${\mathrm{CO}}_{\mathrm{3}}^{=}$] data reconstructed by Yu et al. (2013) reveals a significant difference between TRAF and TRAFc in the Atlantic (Fig. 7b). Decrease in [${\mathrm{CO}}_{\mathrm{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 [${\mathrm{CO}}_{\mathrm{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 [${\mathrm{CO}}_{\mathrm{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⁻³. Changes in the TRAFc experiment are more pronounced than in TRAF due to stronger changes in alkalinity. As expected, [${\mathrm{CO}}_{\mathrm{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₃ burial in the model (27.9–29.1, plus 13 TmolC yr⁻¹ surface removal in the TRAFc experiment) is higher than in the data (23.3 TmolC yr⁻¹); however, this is compensated by a higher modelling weathering rate (24.6–34.7 TmolC yr⁻¹) compared to 11.7 TmolC yr⁻¹ in the data. The model values are at the upper end of the data uncertainty range (11–38 TmolC yr⁻¹; see min–max range for CaCO₃ sedimentation in Table 2). Organic carbon burial in the model (10.5–11.3 Tmol yr⁻¹) is less than in the averaged data (18.3 Tmol yr⁻¹); however, the uncertainty in the burial is so high (6–58 Tmol yr⁻¹) that the model values are almost twice more than the lower end of the data range.

Table 2 Model–data comparison of average carbon fluxes from 6000 BCE to 1850 CE (Tmol yr⁻¹).

^a Including aerial volcanic CO₂ outgassing. ^b Pre-industrial fluxes according to the Fig. 1 in Cartapanis et al. (2018). ^c Subtracting volcanic outgassing.

An important question is whether the ocean on average lost carbon content over the Holocene. In equilibrium, volcanic CO₂ 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.8 TmolC yr⁻¹. Similar to the data, the model shows the loss of carbon from the water column, 13.8 and 18.7 TmolC yr⁻¹ in the TRAF and TRAFc experiments, respectively. Therefore, simulated ocean carbon losses are qualitatively (and even quantitatively for TRAFc) in line with observations.

3.4 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₂ 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 ¹³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 δ¹³CO₂. This small change is within the uncertainty bounds of δ¹³CO₂ 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 δ¹³CO₂ content. Simulated changes in biological production and export flux might have affected the atmospheric δ¹³CO₂, 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, transient deglaciation run 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₃ and POC fluxes to sediments are changing with time; this also leads to changes in the rain ratio. In the absence of factorial experiments without these changes, it is difficult to infer how these trends in nutrients and biogenic opal and carbonate fluxes affect atmospheric CO₂. These two caveats (steady-state initial conditions and fixed weathering) apply to both TRAF and TRAFc simulations. Consequently, to simulate the carbon budget correctly, models have to include interactive weathering processes.

4 Conclusions

Using the Earth system model MPI-ESM-LR, we performed two transient simulations of the climate and carbon cycle dynamics in the Holocene, one with prescribed atmospheric CO₂ and one with interactive CO₂ using nudged ocean alkalinity. In both simulations, the land is a carbon sink during the mid-Holocene (from 6000 to 2000 BCE) and a source of CO₂ after 1 CE due to land use changes. Changes in vegetation cover at 6000 BCE relative to 1 CE (enhanced vegetation cover in North Africa, northward extension of boreal forest in Asia) are in line with available pollen records.

In the prescribed CO₂ 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₂ 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 simulation with interactive CO₂ is performed in the nudging mode: we use the surface alkalinity changes as the forcing for ocean–atmosphere CO₂ 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₂ 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₂ 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 (non-anthropogenic) 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 ice-core CO₂ 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₂ until the late Holocene when anthropogenic CO₂ sources started to affect atmospheric CO₂.