The [simple carbon project] model v1.0
- Research School of Earth Sciences, Australian National University, Canberra, Australia
Correspondence: Cameron O'Neill (email@example.com)
We construct a carbon cycle box model to process observed or inferred geochemical evidence from modern and paleo settings. The [simple carbon project] model v1.0 (SCP-M) combines a modern understanding of the ocean circulation regime with the Earth's carbon cycle. SCP-M estimates the concentrations of a range of elements within the carbon cycle by simulating ocean circulation, biological, chemical, atmospheric and terrestrial carbon cycle processes. The model is capable of reproducing both paleo and modern observations and aligns with CMIP5 model projections. SCP-M's fast run time, simplified layout and matrix structure render it a flexible and easy-to-use tool for paleo and modern carbon cycle simulations. The ease of data integration also enables model–data optimisations. Limitations of the model include the prescription of many fluxes and an ocean-basin-averaged topology, which may not be applicable to more detailed simulations.
In this paper we demonstrate SCP-M's application primarily with an analysis of the carbon cycle transition from the Last Glacial Maximum (LGM) to the Holocene and also with the modern carbon cycle under the influence of anthropogenic CO2 emissions. We conduct an atmospheric and ocean multi-proxy model–data parameter optimisation for the LGM and late Holocene periods using the growing pool of published paleo atmosphere and ocean data for CO2, δ13C, Δ14C and the carbonate ion proxy. The results provide strong evidence for an ocean-wide physical mechanism to deliver the LGM-to-Holocene carbon cycle transition. Alongside ancillary changes in ocean temperature, volume, salinity, sea-ice cover and atmospheric radiocarbon production rate, changes in global overturning circulation and, to a lesser extent, Atlantic meridional overturning circulation can drive the observed LGM and late Holocene signals in atmospheric CO2, δ13C, Δ14C, and the oceanic distribution of δ13C, Δ14C and the carbonate ion proxy. Further work is needed on the analysis and processing of ocean proxy data to improve confidence in these modelling results.
A box model divides regions of the ocean into boxes or grids based on some property of the composite water masses, such as temperature, density or chemical composition. The model equations describe the evolution of tracers in the model's boxes due to the various fluxes between each box (Fig. 1). Box models differ from more complex models such as general circulation models (GCMs), mainly due to their reduced spatial resolution (i.e. much larger grids or boxes) and with major processes and fluxes typically prescribed rather than calculated in the model. Box models range in complexity from simple, basin-averaged models (e.g. Stommel, 1961; Sarmiento and Toggweiler, 1984; Toggweiler, 1999) to more complex, multi-basin ocean models (Hain et al., 2010) and full Earth carbon cycle models (Zeebe, 2012). Box models, despite being simpler than GCMs, have been useful in illustrating key concepts in oceanography. For example, Sarmiento and Toggweiler (1984), Siegenthaler and Wenk (1984), and Knox and McElroy (1984) used simple four-box ocean–atmosphere models to show that the LGM CO2 drawdown could have resulted from increased biological productivity and/or reduced ocean overturning circulation. More recently, Hain et al. (2010) used a box model to show that a range of ocean physical and biological mechanisms were required to cause lower atmospheric CO2 levels in the LGM, and Zhao et al. (2017) used a similar model to explore ocean ventilation ages in the LGM–deglacial and Holocene periods. Despite the development of highly complex coupled atmosphere–ocean models for climate simulations, box models continue to be applied to resolve problems in the carbon cycle.
Our motivation in constructing a new box model of the full carbon cycle, the [simple carbon project] model v1.0 (SCP-M), is to contribute a simple, easy-to-use, open-access model implemented with freely available software, that is consistent with physical and biogeochemical oceanography, that caters for important features of the carbon cycle, and that has explicit avenues for data integration, optimisation and inversion. Recent advances in physical oceanography have refined our understanding of global ocean circulation and mixing fluxes. For example, Talley (2013) provided a simplified interpretation of the global ocean as a handful of large-scale processes, some of which are operating across all basins – as is the case with the global overturning circulation. De Boer and Hogg (2014) described a simple model of deep ocean mixing of water masses under the influence of sea-floor topography. These high-level concepts are easy to apply to box models. Furthermore, the growing pool of paleo proxy data across carbon isotopes and reconstructions (e.g. carbonate ion) presents an opportunity to progress model–data integrations using a number of different proxies. SCP-M caters for a range of proxies including the carbon isotopes and carbonate ion proxy, with the capacity for additional elements with minimal programming effort. The model–data experiment described in this paper provides a direct linkage between paleo data and discrete values for ocean parameters in the LGM and late Holocene periods, thus contributing to our understanding of the LGM–Holocene carbon cycle transition. Combined with the expanding dataset of paleo observations, and with advances in computing power, data-aligned models such as SCP-M have the potential to improve our understanding of past changes in climate across many other time frames. Furthermore, there are a number of features of the carbon cycle outside ocean circulation and biology which influence proxy indicators, particularly the carbon isotopes. Omitting these features could lead to erroneous modelling outcomes, particularly in the case of the terrestrial biosphere, which strongly influences atmospheric CO2 and δ13C (e.g. Francois et al., 1999). We compiled SCP-M to include a broad range of carbon cycle mechanisms, including carbonate production and dissolution, marine sediments, terrestrial biosphere, anthropogenic emissions sources, and continental weathering. Finally, SCP-M is computationally cheap and quick to run. A 10 000-year simulation takes approximately 30 s to process on a regular laptop, enabling exhaustive exploration of parameter space in optimisations that incorporate large datasets. While box models are not new, we argue that these features contribute to a new tool that is well-equipped to tackle a wide range of applications in paleoceanography, paleoclimate and the modern carbon cycle.
In this paper we describe SCP-M and illustrate its application alongside LGM and late Holocene ocean and atmosphere data, with several insights for the transition between the two periods, plus modelling of the modern and future carbon cycle under the influence of anthropogenic emissions. Emphasis is placed on the model description and configuration.
SCP-M is focussed on the ocean carbon cycle and is configured to estimate the time evolution of oceanic dissolved inorganic carbon (DIC) and its constituents: δ13C, Δ14C, plus alkalinity, phosphorus, oxygen, and atmospheric CO2, δ13C and Δ14C. It contains a simple yet realistic representation of large-scale ocean physical processes, with an overlay of ocean chemistry and biology (Fig. 1). SCP-M simulates sources and sinks of carbon across the ocean and atmosphere, marine sediments, and terrestrial biosphere. Volcanic emissions, sedimentary weathering, rivers and anthropogenic emissions are prescribed fluxes. A range of carbon cycle features are included because the concentrations of carbon in the ocean and atmosphere (and its isotopes in particular) are sensitive to many sources and sinks, and omitting them makes it difficult to compare model results with the carbon data that indelibly feature their imprint. For example, regrowth in the terrestrial biosphere imparts a clear signature on the atmosphere and ocean δ13C data profile after the LGM (e.g. Francois et al., 1999; Ciais et al., 2012; Hoogakker et al., 2016). In addition, the atmospheric radiocarbon source, marine sediments, volcanic emissions, continental weathering and now anthropogenic emissions exert important influences on carbon cycle observations.
SCP-M was designed to compare model results with data and to solve for optimal parameter combinations. Within SCP-M, realistic implementation of physical processes upon a sound biogeochemical platform enables their transmission into paleo chemical tracer signals, for which proxy data exist. Many of the key ocean physical and biological processes are prescribed in the model, allowing them to be free parameters in model–data experiments. SCP-M itself is implemented with a matrix framework which enables more boxes to be added, ocean basins to be separated, elements to be added and exploration of a range of hypotheses, all with minimal programming effort.
2.1 Model topology
The box model is mostly conceptual in nature and is designed to test high-level concepts. Therefore, excessive detail and complication are to be avoided. However, key processes that are critical to the validity of any hypothesis being tested must be represented as well as possible.
The ocean is a key part of the global carbon cycle and pre-eminent in hypotheses of glacial–interglacial carbon cycles (e.g. Kohfeld and Ridgewell, 2009; Sigman et al., 2010). Talley (2013) provided an observationally based description of ocean circulation in terms of its constituent water masses, circulation and mixing fluxes, including estimates of the present-day magnitudes of those fluxes. The Talley (2013) model builds on the models of Broecker (1991), Gordon (1991), Schmitz (1996), Lumpkin and Speer (2007), Kuhlbrodt et al. (2007), Talley (2008), and Marshall and Speer (2012). Key features of the Talley (2013) model include the following.
Atlantic thermocline water moves north, ultimately reaching the North Atlantic, driven by advection and surface buoyancy changes. High-salinity North Atlantic Deep Water (NADW) forms in the north by cooling, densification and convection and then travels south to rise up into the Southern Ocean via wind-driven upwelling and Ekman flows, forming Lower Circumpolar Deep Water (LCDW). This water comprises the upper (orange arrows) Atlantic meridional overturning circulation (AMOC) in SCP-M (Fig. 1).
A fraction of the upwelled LCDW sinks to become Antarctic Bottom Water (AABW) under the influence of cooling and brine rejection, south of the Antarctic Circumpolar Current (ACC). AABW moves northward along the ocean floor via adiabatic advection (Talley, 2013) in all basins. It upwells into deep water in the Pacific and Indian oceans and also into NADW in the Atlantic via upwelling with diapycnal diffusion (Talley, 2013). The combined LCDW–AABW–PDW–IDW–NADW global overturning circulation (GOC) is represented by the red arrows in Fig. 1.
Pacific Deep Water–Indian Deep Water (PDW–IDW) upwells at low latitudes and returns to the Southern Ocean above the NADW, forming the core of the Upper Circumpolar Deep Water (UCDW), which is identified by Talley (2013) as low-oxygen-content (old) water. A part of the upwelled PDW–IDW joins AABW formation, with the bulk of it moving north from the Southern Ocean at shallow and intermediate depths. These waters become fresher and warmer and join Antarctic Intermediate Water (AAIW) and Subantarctic Mode Water (SAMW) at the base of the subtropical thermocline (advection with surface buoyancy fluxes).
Joined thermocline waters, AAIW–SAMW, and upwelled thermocline waters from the Pacific and Indian oceans form the upper-ocean transport moving towards the North Atlantic.
A key contribution of the Talley (2013) study is that GOC is the pre-eminent process in distributing water throughout the global oceans. Talley (2013) provided a 2-D “collapsed” interpretation of a 3-D ocean layout based on the observation that similar large-scale processes (i.e. GOC) operate in all three major ocean basins, and this interpretation can directly inform a box model topology. The Talley (2013) 2-D global ocean view used in SCP-M captures the features described above in a simple ocean box model format. Talley (2013) also provided observation-based estimates of the ocean transport fluxes, which are scaled according to their ocean basin domain. For example, the GOC and AABW formation process is common to all basins and thus accounts for the largest flux of 29 Sverdrups (Sv, 106 m3 s−1). The AMOC–NADW sinking cell is confined to the Atlantic Basin and represents a smaller flux (19 Sv) of water (Talley, 2013).
The SCP-M dimensions are designed to be consistent with measured estimates of the ocean surface area, average depth, and total ocean and atmosphere volumes. The model is divided into boxes according to latitude and depth, but not by longitude. Therefore, it does not distinguish between the Atlantic, Pacific and Indian basins and does not allow for compositional variations with longitude. Each box has a surface area and depth (and therefore volume) and corresponds to a location in the global ocean with reference to latitude and average depth. It is simple to add more boxes to divide the model into ocean basins.
SCP-M contains seven ocean boxes as shown in Fig. 1. The rationale for dividing the ocean into boxes is that there are regions of the ocean that are relatively well-mixed, or at least similar in terms of their prevailing element flux behaviour. For the depth of the surface boxes, this rationale conveniently translates to the maximum wintertime mixed layer depth (MLD) (e.g. Kara et al., 2003; de Boyer Montegut et al., 2004). We choose a depth of 100 m for box 1, the low-latitude surface box, which is a reasonable approximation to the 20–150 m seasonally varying MLD for the middle and low latitudes estimated by de Boyer Montegut et al. (2004) and consistent with the depth of a similar box in the Toggweiler (1999) model. This box represents the photic zone over much of the ocean, from 40∘ S to 40∘ N. Craig (1957) estimated the depth of this layer as 75 m ± 25 m, a value used by Keeling and Bolin (1968) in their simple ocean box model. We choose 250 m of depth for the NADW box (box 2) and the subpolar surface box (box 7) as per Toggweiler (1999). These boxes are deeper than the low-latitude surface box (de Boyer Montegut et al., 2004) in order to capture the regions of deepwater upwelling (subpolar Southern Ocean) and convective downwelling (North Atlantic). The MLD in these regions can vary between 70 and >500 m of depth depending on seasonal variations (de Boyer Montegut et al., 2004). An intermediate-depth box (3) resides below the low-latitude surface box and extends from 100 m to 1000 m of depth. This box captures northward-flowing AAIW and SAMW from upwelled NADW–PDW–IDW (e.g. Talley, 2013). Box 4 is the deep ocean box, extending from 1000 m to 2500 m of depth, and incorporates the upwelling abyssal waters in all basins and downwelled NADW. This water is channelled back to the surface in the subpolar surface box and the Southern Ocean box, as per the wind-driven upwelling of Morrison and Hogg (2013) and Talley (2013). The Southern Ocean box (5) extends from 80 to 60∘ S and from the ocean surface to 2500 m of depth. This box encompasses the Southern Ocean, the ACC and deepwater formation from southward-flowing upwelled NADW–PDW–IDW (Talley, 2013). The abyssal box (6) extends the full range of the ocean, from 2500 to 4000 m of depth (our assumed average depth of the ocean). This box is the pathway for northward-flowing AABW and incorporates mixing with overlying deep water and advection–upwelling (Talley, 2013).
2.2 The model parameters, processes and equations
2.2.1 Basic features
Figure 2 shows the suite of files used to execute SCP-M. We have chosen a modular approach to reduce the complexity of each of the model files. The SCP-M routine includes data processing for the model's boxes based on geographic coordinates, model calibration to the data, model simulations, model–data optimisation, and charting and/or tabular output. SCP-M is implemented in Python 3.6, with the code and download (including user instructions) available at https://doi.org/10.5281/zenodo.1310161 (O'Neill et al., 2018).
In short, SCP-M calculates the evolution of an element or species concentration in each model box as a function of time and flux parameters (e.g. inputs and outputs to each box) or processes, such as uptake or regeneration. The model includes ocean circulation and mixing fluxes, air–sea gas exchange, chemical and biological transformations, and sources and sinks of carbon. The model equations are a set of partial differential equations, one for each element in the model. These are solved with a straightforward, first-order, Euler-forward time-stepping method with a standard time step of 1 year. We find the model to be stable and approaching steady state for most of the simulations that we undertook. However, this stability can be challenged by scenarios with strong forcing. With the Euler method, errors can propagate in proportion to the step size. This can be resolved by revising the selection of parameter inputs or starting data values or by reducing the size of the time step in each model run.
2.2.2 The ocean circulation and mixing
There are four ocean physical parameters in SCP-M. Ψ1 and Ψ2 are advection terms that represent the physical transport of water from one box to another, containing the element or species concentration of its box of origin. Ψ1 represents GOC (e.g. Sarmiento and Toggweiler, 1984; Marshall and Speer, 2012; Talley, 2013) that infiltrates all basins (Talley, 2013) and is shown by the red arrows in Fig. 1. The Ψ1 parameter allows variable allocation between transport from the deep ocean box (box 4) into the subpolar surface box (box 7) and directly into the polar box (box 5) via α. Ψ2 represents the AMOC. This is the region of northward-flowing intermediate waters in the Atlantic Ocean and the formation of NADW (Dickson and Brown, 1994; Talley, 2013), shown as orange arrows in Fig. 1. γ1 and γ2 are bidirectional mixing terms that exchange element or species concentrations between boxes without any net advection of water (blue arrows in Fig. 1). γ1 is bidirectional mixing between the deep and abyssal boxes of the form described by Lund et al. (2011) and De Boer and Hogg (2014). γ2 is a low-latitude, intermediate–shallow box “thermocline” mixing parameter, which governs the constant bidirectional exchange between these two boxes (Liu et al., 2016).
The influence of each of the ocean parameters is prescribed in box model space by matrix equations, with one matrix for each parameter. Each row and column position in the matrix corresponds to a box location. The atmosphere box is treated separately from the ocean boxes, and it does not enter the ocean parameter matrices. The volumetric circulation or mixing parameters, in Sv, are multiplied by the oceanic element concentration (mol m−3) to produce a molar flux exchanged between ocean boxes. For example, the change in concentration of carbon (as DIC) in the deep box (box 4) from ocean physical parameters is estimated by
where Ci is the concentration of DIC in each box in mol m−3 and Vi is the volume of each box in m3. There is no vertical flux between box 4 and box 3 (intermediate box) in Eq. (1). We assume that this vertical flux is small compared with the lateral transport and compared with the mixing fluxes between boxes 4 and 6 (and boxes 1 and 3 in Eq. 2 below). We assume that the boundary between boxes 3 and 4 is the divide between northward-flowing water sourced from AAIW and SAMW, which overlies southward return flow from AMOC and PDW–IDW. The fluxes out of box 4 are shown by the terms −Ψ1C4, −Ψ2C4 and −γ1C4, with the fluxes into boxes 5, 6 and 7 treated in the equations for those boxes. For the low-latitude surface box (box 1),
We assume that box 1 represents the shallow, mixed layer of the ocean (e.g. Kara et al., 2003; de Boyer Montegut et al., 2004), which is mainly influenced by surface processes. We prescribe vertical mixing between this box and the underlying intermediate box (3) via the γ2 parameter. This represents the thermocline mixing described by Liu et al. (2016). We assume that lateral transport of northward-flowing water underlies box 1, involving box 7 (subpolar Southern Ocean), box 3 and box 1 (northern ocean). This intermediate-depth water is colder and denser than the overlying mixed layer, given its provenance of AAIW, SAMW and deep-upwelled NADW–PDW–IDW (Talley, 2013). These ocean circulation and mixing operations (e.g. Eqs. 1 and 2) can be vectorised for all boxes using sparse matrices, as follows.
and T1, T2, E1 and E2 are sparse matrices defined as
Given that we have applied the global ocean interpretation of Talley (2013) to the SCP-M layout, we have also adopted the observationally based estimates for the large-scale ocean fluxes for the modern ocean from the same study: GOC (Ψ1, 29 Sv), AMOC (Ψ2, 19 Sv) and deep–abyssal mixing (γ1, 19 Sv). For thermocline mixing between boxes 1 and 3 (γ2), we have adopted the value for the corresponding flux from Toggweiler (1999) of 40 Sv.
2.2.3 Biological flux parameterisation
The biological pump (e.g. Broecker, 1982) is a descriptor of marine biological activity, whereby organisms consume nutrients in shallow waters, die, sink and then release those nutrients at depth. For example, carbon is taken up by shallow-water-dwelling phytoplankton through photosynthesis and then sequestered in deeper waters after sinking, breaking down and remineralising their nutrient load back into the water column. Volk and Hoffert (1985) made the distinction between the soft tissue pump (STP) for soft-tissued organisms and the carbonate pump (carbonate-shelled organisms). We also distinguish between the two pumps, as they have different effects on carbon and alkalinity balances and therefore pCO2 and carbonate dissolution. This section deals with the STP, and a following section deals with the carbonate pump. Most STP organic matter is remineralised in the shallow to intermediate ocean depths, leading to a decrease in the export of carbon as depth increases. According to Henson et al. (2011), only ∼15 %–25 % of organic material is exported to >100 m of depth, with most recycled in the shallower waters.
Martin et al. (1987) modelled the soft-bodied organic flux of carbon observed from sediment traps in the northeast Pacific to create a simple power rule which is easily applicable to modelling. The Martin et al. (1987) equation produces a flux of organic carbon, which is a function of depth from a base organic flux at 100 m of depth (the “Martin reference depth”). The flux at 100 m of depth was estimated by Martin et al. (1987) to be between 1.2 and 7.1 mol C m−2 yr−1 from eight station observations in the northeast Pacific Ocean. Sarmiento and Gruber (2006) estimated a range of 0.0–5.0 mol C m−2 yr−1, with some localised higher values, across the global ocean. Equation (10) shows the general form of the Martin et al. (1987) equation:
where F is a flux of carbon in mol C m−2 yr−1, F100 is an estimate of carbon flux at 100 m of depth, z is depth in metres and b is a depth scalar. In SCP-M, the Z parameter implements the Martin et al. (1987) equation. Z is an estimate of biological productivity at 100 m of depth (in mol C m−2 yr−1) and, coupled with the Martin et al. (1987) depth scalar, controls the amount of organic carbon that sinks from each model surface box to the boxes below. Each subsurface ocean box receives a flux of carbon from the box above it, at its ceiling depth (also the floor of the overlying box), and loses carbon as a function of the depth of the bottom of the box. Remineralisation in each box is accounted for as the difference between the influx and outflux of organic carbon. The biological flux out of the surface box 1 is shown by
where Z1 is the biological flux of carbon prescribed for surface box 1 in mol C m−2 yr−1, S1 is the surface area of surface box 1, d0 is the reference depth of 100 m for the Z parameter value (Martin et al., 1987), and dc and df are the ceiling and floor depths of a box, respectively. The dimensionless parameter b is the depth power function of the Martin et al. (1987) equation, which tapers biological production and export below depths of 100 m. The net biological flux for intermediate-depth box 3 is given by
The process is vectorised using sparse matrices in the following:
where Z is an array of the Zi () parameter, which varies across the surface boxes, and S is the array of surface box surface areas Si (). As with the ocean parameters, the biological flux of carbon is divided by the box volume array V to return concentrations in mol m−3. Bout and Bin are sparse matrices as follows.
The value of the parameter Z in each surface box is specified to vary as a fraction of the global value for Z in SCP-M (presently 5.0 mol C m−2 yr−1). There are higher fractions in the northern and southern oceans and smaller fractions in the low-latitude and polar oceans (e.g. Sarmiento and Gruber, 2006). During the model set-up we manually tuned the individual surface box values by multiplying the global value for Z by a scalar unique to each box. The values were tuned to align the model's output with GLODAPv2 data for DIC, phosphorous, alkalinity, the carbonate ion and carbon isotopes in each of the ocean boxes (Table 1). The resulting range for Z (1.1–5.33 mol C m−2 yr−1) compares with the observations-based range of Martin et al. (1987) of 1.2–7.1 mol C m−2 yr−1 and Sarmiento and Gruber (2006) of 0–5 mol C m−2 yr−1. We chose a value for the dimensionless b depth decay parameter of 0.75, which falls in the range of Gloege et al. (2017) of 0.68–1.13 and the error range of Berelson (2001) of 0.82±0.16. We found a global value of 0.75 to produce a better fit to the GLODAPv2 data when calibrating the model. The biological flux of other elements and species, such as phosphorous and alkalinity, are calculated from the biological carbon flux using so-called Redfield ratios (e.g. Redfield et al., 1963; Takahashi et al., 1985; Anderson and Sarmiento, 1994).
2.3 pCO2 and carbonate
Modelling air–sea gas exchange, atmospheric pCO2 and the “carbonate pump” rests on a realistic estimation of pCO2 in the ocean. For example, only a fraction of DIC in seawater can exchange with the atmosphere, and this fraction is estimated by the oceanic pCO2. DIC itself consists of three major constituents: carbonic acid, bicarbonate and carbonate. Their relative proportions depend on total DIC, alkalinity, pH, temperature and salinity (Zeebe and Wolf-Gladrow, 2001).
pCO2 is estimated by subtracting alkalinity from DIC. However, this is only accurate to ±10 % (Sarmiento and Gruber, 2006), which may cause problems for scenario analysis and sensitivity testing within such a large error band. More complex calculations can require numerous iterations and can be computationally expensive (e.g. Toggweiler and Sarmiento, 1985; Zeebe and Wolf-Gladrow, 2001; Follows et al., 2006). We apply the routine of Follows et al. (2006) in SCP-M, which is a direct solution, rather than an iterative approach to solve for pCO2 at each time step of a model run. This was demonstrated by Follows et al. (2006) to be sufficiently accurate for modelling purposes. The calculation takes inputs of DIC, alkalinity, temperature, salinity, phosphorous and silicate to estimate pCO2.
Solving for pCO2 enables the calculation of the concentrations of the three species of DIC, which further enables estimation of the dissolution and burial of carbonate in the water column and sediments. The latter is an important part of the oceanic carbon and alkalinity cycles and provides important feedbacks to atmospheric CO2 on 1000-year time frames (e.g. Farrell and Prell, 1989; Anderson et al., 2007; Mekik et al., 2012; Yu et al., 2014b).
2.3.1 The carbonate pump
According to Emerson and Hedges (2003), ∼20 %–30 % of CaCO3 formed in the ocean's surface is preserved in ocean floor sediments, with the rest dissolved in the water column. Klaas and Archer (2002) estimated that 80 % of the organic matter fluxes in the ocean below 2000 m are driven by organic matter associated with carbonate ballast. Therefore, the so-called “carbonate pump” is a relatively efficient transport of carbon and alkalinity in the ocean. According to Farrell and Prell (1989), the carbonate pump is a dynamic process, and the dissolution and burial in sediments of CaCO3 is observed to vary across (and within) glacial–interglacial cycles.
To replicate this flux of carbon and alkalinity, a term is added to the carbon cycle equation to represent the flux of calcium carbonate (shells) out of the surface boxes into the abyssal box and sediments. This is an extension of the surface organic carbon flux Z described in Eq. (13) via the “rain ratio” parameter. The rain ratio is a common term in ocean biogeochemistry (e.g. Archer and Maier-Reimer, 1994; Ridgewell, 2003) and refers to the ratio between shell-based “hard” carbon and organic “soft” carbon fluxes in the biologically driven rain of carbon from the ocean's surface. Sarmiento et al. (2002) estimated a global average value for the rain ratio of 0.06±0.03, with local maxima and minima of 0.10 and 0.02, respectively, providing a narrow range of global values. We apply the rain ratio as a parameter multiplied by the organic flux parameter Z. We chose an initial value of 0.07, which provided appropriate values for DIC, alkalinity and dissolution in the model's boxes during the model spin-up (with reference to GLODAPv2 data). The combination delivers the physical production and export of calcium carbonate at the Martin reference depth (100 m).
Once the production and export flux at the Martin reference depth is established, the distribution of calcium carbonate in the boxes below is a function of dissolution. According to Milliman et al. (1999), the theory that calcium carbonate only dissolves at great depths in carbonate-undersaturated water is “one of the oldest and most strongly held paradigms in oceanography” (e.g. Sverdrup et al., 1941). However, in nature, the alkalinity and carbonate ion concentration profiles suggest that 30 %–60 % of carbonate produced is dissolved in shallower water that is saturated (Harrison et al., 1993; Milliman et al., 1999). Theories for this outcome include the emergence of locally undersaturated waters due to remineralisation of biological carbon (Jansen et al., 2002) and dissolution by zooplankton grazing (Milliman et al., 1999). Battaglia et al. (2016) found similar skill in model results for replicating observed dissolution profiles, whether a non-saturation-dependent or saturation-dependent dissolution constant was used. Battaglia et al. (2016) recommended the use of a basic non-saturation-dependent (i.e. constant) dissolution parameter in Earth carbon system models for computing efficiency, with limited loss of accuracy. As such, we include two parts to the dissolution equation: a non-saturation-dependent dissolution constant, to reflect the “unknown” processes that likely cause the observed dissolution of calcium carbonate in waters that are saturated, and a saturation-state-dependent component using the dissolution function of Morse and Berner (1972). We include the latter to enable dynamic feedback to take place in the carbonate system after model perturbations. The saturation-dependent dissolution is a function of the average carbonate ion composition for each box relative to its temperature- and pressure-dependent saturation concentration (Morse and Berner, 1972; Millero, 1983). We choose the median depth of each box for the calculation in the ocean boxes and the floor of the abyssal box for the sediment surface dissolution. We assume 100 % of calcium carbonate takes the form of calcite. If the surface export flux of CaCO3 is greater than dissolution in the ocean boxes, then the remainder escapes to the sediments. This is a flux out of the ocean of alkalinity and carbon in the ratio of 2:1 assumed for carbonate shells (Sarmiento and Gruber, 2006). DIC and alkalinity can return to the abyssal box from the sediments via undersaturation-driven dissolution in the abyssal water overlying the sediments.
The net flux of carbonate between ocean boxes and out of the ocean and into the sediments is shown in vectorised Eq. (16):
where FCA is the rain ratio, ζ is the constant background dissolution rate, ϵ is the saturation-state-dependent dissolution function of Morse and Berner (1972) and Millero (1983), and CaCO3 is the concentration of calcium carbonate in each box. The dissolution equation of Morse and Berner (1972) operates on CaCO3, which is calculated by multiplying Ca by , where Ca is estimated from salinity in each box as per Sarmiento and Gruber (2006).
2.3.2 Air–sea gas exchange
CO2 is transported across the air–sea interface by gaseous exchange. According to Henry's law, the partial pressure of a gas [P] above a liquid in thermodynamic equilibrium will be directly proportional to the concentration of the gas in the liquid:
where KH is the solubility of a gas in mmol m−3 atm−1 and C is its concentration in the liquid. Many ocean models specify the air–sea gas exchange of CO2 as a function of the pCO2 differential between ocean and atmosphere, a CO2 solubility coefficient (e.g. Weiss, 1974), and a so-called “piston” or gas transfer velocity, which governs the rate of gas exchange, in m s−1 (e.g. Toggweiler, 1999; Zeebe, 2012; Hain et al., 2010; Watson et al., 2015). We adopt the same approach in estimating the exchange of CO2 between a surface box and the atmosphere:
where P1 is the piston velocity parameter in box 1 in m s−1. P is allowed to vary in each surface box, enabling scenario analysis such as the variation of sea-ice cover in the polar box (e.g. Stephens and Keeling, 2000; Ferrari et al., 2014). is the solubility of CO2 in mol kg−1 atm−1 (Weiss, 1974), subsequently converted into mol m−3 by multiplying by seawater density ρ. and pCO are the partial pressures of CO2 in surface ocean box 1 and the atmosphere, respectively, in ppm. The equation is vectorised as follows:
where P=Pi (), with zero values for non-surface boxes, and ().
2.4 Sea surface temperature and salinity
Ocean box temperature and salinity are forced in SCP-M, not calculated by the model. Each box has a prescribed value for temperature and salinity, which can be time dependent. During initialisation, the model takes box-averaged values for temperature and salinity from the GLODAPv2 database. The values can be varied between model experiments, for example Holocene versus LGM reconstructions. We argue that this approach is plausible given the availability of temperature and salinity proxy–data inputs for a range of paleo (e.g. Adkins et al., 2002; Kohfeld and Chase, 2017), modern (e.g. Olsen et al., 2016) and future scenarios (e.g. IPCC, 2013a). For the modern and future scenarios (Sect. 3.3), we forced temperature with time series data and projections. The temperature and salinity values feed into the calculations for ocean pCO2, which further enables the calculation of air–sea gas exchange and the species of DIC in seawater (H2CO3, HCO3− and ).
2.5 Atmosphere and terrestrial carbon cycle
SCP-M incorporates the terrestrial biosphere, continental weathering and river run-off into the ocean, plus an atmospheric radiocarbon source and volcanic and industrial emissions.
V is a constant prescribed flux of volcanic emissions of CO2 in SCP-M. Toggweiler (2008) estimated this volcanic flux of CO2 emissions at 4.98×1012 mol yr−1 using a carbon cycle model which balanced volcanic emissions with land-based weathering sinks. The weathering of carbonate and silicate rocks also creates DIC and alkalinity run-off into rivers, which finds its way into the ocean (Amiotte Suchet et al., 2003). Relative alkalinity and DIC concentrations affect ocean pCO2 and carbonate ion levels, which impacts atmospheric CO2 and the dissolution and burial of carbonates (Sarmiento and Gruber, 2006). We apply the approach of Toggweiler (2008), whereby silicate and carbonate weathering fluxes of DIC and alkalinity enter only the low-latitude surface ocean box (box 1):
where WSC is a constant silicate weathering term set at mol m−3 yr−1, WSV is a variable rate of silicate weathering per unit of atmosphere CO2 (ppm) set to 0.5 mol m−3 atm−1 CO2 yr−1 and WCV is the variable rate of carbonate weathering with respect to atmosphere CO2 set at 2 mol m−3 atm−1 CO2 yr−1 (Toggweiler, 2008).
Alkalinity is added to the ocean in the ratio of 2:1 to DIC (Toggweiler, 2008). In the case of silicate rocks, weathering is also a weak sink of CO2 (e.g. Toggweiler, 2008; Hogg, 2008). The atmospheric sink of CO2 is calculated by multiplying Eq. (20) by the volume of the low-latitude surface ocean box (box 1) and subtracting from atmospheric CO2. Equation (20) is vectorised by multiplying by a vector of boxes with only a non-zero value for box 1.
The terrestrial biosphere may act as a sink of CO2 during periods of biosphere growth (e.g. post-glacial regrowth), via carbon fertilisation, or a source of CO2 (e.g. glacial reduction) via respiration. We employ a two-part model of the terrestrial biosphere with a long-term (woody forest) and short-term (grassland) terrestrial biosphere box as described by Raupach et al. (2011) and Harman et al. (2011), with net primary productivity (NPP) and respiration parameters controlling the balance between the uptake and release of carbon. NPP is positively affected by atmospheric CO2, the so-called “carbon fertilisation” effect (e.g. Raupach et al., 2011). Respiration is assumed proportional to the carbon stock. The biosphere also preferentially partitions the lighter carbon isotope 12C, leading to a relative enrichment in δ13C in the atmosphere during net uptake of CO2. The change in atmospheric CO2 from the terrestrial biosphere in the model is given by
where Npre is NPP at a reference level (“pre”) of atmospheric CO2, and RP is the parameter to split NPP between the short-term terrestrial biosphere carbon stock (fast respiration) and the long-term stock (slow respiration) (Raupach et al., 2011). β is the parameterisation of carbon fertilisation, which causes NPP to increase (decrease) logarithmically with rising (falling) atmospheric CO2 levels, and has a typical value of 0.4–0.8 (Harman et al., 2011). Cstock1 is the short-term carbon stock and k1 is the time frame for respiration in the short-term carbon stock (in years). For the long-term terrestrial biosphere, we substitute (1−RP) in place of RP and Cstock2 and k2 for the long-term carbon stock and respiration rate, respectively. Dforest is a prescribed flux of deforestation emissions, which can be switched on or off in SCP-M. A δ13C fractionation factor is applied to the terrestrial biosphere fluxes of carbon, effecting an increase in atmospheric δ13C from biosphere growth and a decrease from respiration.
2.6 The complete carbon cycle equations
The calculation of atmospheric CO2 is
where the additional term consists of a prescribed flux of δ13C-depleted and 14C-dead CO2 to the atmosphere from human industrial emissions, which is activated by a model switch in SCP-M.
2.7 Treatment of carbon isotopes
Carbon isotopes are an important component in SCP-M because they are key sources of proxy data. Carbon isotope fluxes are treated largely the same as DIC in SCP-M, with minor modification. For example, carbon isotopes are typically reported in delta notation (δ13C and Δ14C), which is the ‰ deviation from a standard reference value in nature. SCP-M operates with a metric mol m−3 for the other ocean element concentrations and flux parameters. In order to incorporate δ13C and Δ14C into this metric for the operation of model fluxes, the method of Craig (1969) is applied to convert starting data values of δ13C and Δ14C from delta notation in ‰ into mol m−3:
where 13Ci is the 13C concentration in box i in mol m−3, δ13Ci is δ13C in ‰ in box i, R is the ratio of the standard (0.0112372 as per the Pee Dee Belemnite value) and Ci is the DIC concentration C in box i in mol m−3.
The calculation in Eq. (24) derives the fraction for the data or a model starting value, multiplies that value by the standard reference value and then by the starting model concentration for DIC (Ci) in each box. This approach rests on the assumption that the fraction is the same as . For example, there are three isotopes of carbon, each with different atomic weights. They occur in roughly the following abundances: 12C ∼98.89 %, 13C ∼1.11 % and 14C %. Therefore, the assumption that = is a valid approximation. Once converted from δ13C (‰) to 13C in mol m−3, SCP-M's ocean parameters can operate on 13C concentrations in each box, according to the same model flux equations used for DIC and CO2. The 13C model results are then converted back into δ13C notation at the end of the model run in order to compare the model output with the data, which is reported in δ13C format. The same method is applied to Δ14C. The reference standard value for is (Craig, 1969). Where fractionation of carbon isotopes takes place, such as biological or air–sea gas exchange, fractionation factors are simply added to the model flux equations.
2.7.1 Biological fractionation of carbon isotopes
Biological processes influence the carbon isotopic composition of the ocean (e.g. Fontugne and Duplessy, 1978). When photosynthetic organisms are active in shallow ocean waters, they preferentially partition 12C (the lighter carbon isotope). This activity enriches the surface ocean in 13C and relatively enriches the underlying waters in 12C when remineralisation occurs. As such, the ocean displays depletion in δ13C in the deep ocean relative to the shallow ocean (e.g. Curry and Oppo, 2005). In SCP-M, a fractionation factor, f, is simply multiplied by the biological flux in Eq. (13) to calculate marine biological fractionation of 13C and to replicate ocean distributions of δ13C:
where f is the biological fractionation factor for stable carbon (e.g. ∼0.977 in Toggweiler and Sarmiento, 1985), and Sst is the ratio of 13C to 12C in the reference standard. The typical δ13C composition of marine organisms is in the range −23 to −30‰. The same method is applied for biological fractionation of 14C, but with a different fractionation factor (Toggweiler and Sarmiento, 1985).
2.7.2 Fractionation of carbon isotopes during air–sea gas exchange
Fractionation of carbon isotopes takes place during air–sea exchange (e.g. Mook et al., 1974; Siegenthaler and Munnich, 1981; Lynch-Stieglitz et al., 1995). The lighter isotope, 12C, preferentially partitions into the atmosphere as a net effect of bidirectional gaseous exchange. This fractionation leads to the heavily depleted δ13C signature for the atmosphere relative to the ocean. The approach to capture this effect in SCP-M is per Siegenthaler and Munnich (1981):
where λ, the “kinetic fractionation effect” (Zhang et al., 1995), accounts for the slower equilibration rate of the carbon isotopes 13C and 14C across the air–sea interface compared with 12C (Zhang et al., 1995). RAt is the ratio of 13C to 12C in the atmosphere, and Ri is the ratio of 13C to 12C in surface ocean box i. pCO2At is the atmospheric pCO2, and pCO2i is the pCO2 in the surface ocean boxes. τ and π are the fractionation factors of carbon isotopes from air to sea and sea to air, respectively. These are temperature-dependent reactions and are calculated in SCP-M using the method of Mook et al. (1974).
2.7.3 Source and decay of radiocarbon
Natural radiocarbon is produced in the atmosphere from the collision of cosmic-ray-produced neutrons with nitrogen. The production rate is variable over time and can be influenced by changes in solar winds and the Earth's geomagnetic field intensity (Key, 2001). A mean production rate of 1.57 atom m−2 s−1 was estimated from the long-term record preserved in tree rings, although more recent estimates approach 2 atom m−2 s−1 (Key, 2001). For use in SCP-M, this estimate needs to be converted into mol s−1. We first convert atoms to moles by dividing by Avogadro's number (). The resultant figure is multiplied by the Earth's surface area ( cm−2) to yield a production rate of mol s−1. This source rate, divided by the molar volume of the atmosphere, is added to the equation for atmospheric radiocarbon concentration. A decay timescale for radiocarbon of 8267 years is applied to each box in the model.
The modern carbon cycle has been extensively modelled as part of efforts to understand the impact of anthropogenic emissions on climate. There are abundant data on emissions and detailed observations of the modern carbon cycle with globally coordinated ocean surveys and land-based measuring stations. In addition, numerous modelling exercises, using consensus-type emissions projection scenarios from the Intergovernmental Panel on Climate Change (IPCC), have created a body of modelling inputs and results. This provides an ideal testing ground for SCP-M. We first calibrate the model for the pre-industrial period, then simulate historical and projected human emissions under a number of scenarios.Marcott et al. (2014)Schmitt et al. (2012)Nydal and Lövseth (1996)Stuiver et al. (1998)Reimer et al. (2009)Turnbull et al. (2017)Peterson et al. (2014)Skinner and Shackleton (2004)Marchitto et al. (2007)Barker et al. (2010)Bryan et al. (2010)Skinner et al. (2010)Burke and Robinson (2012)Davies-Walczak et al. (2014)Skinner et al. (2015, 2017)Chen et al. (2015)Hines et al. (2015)Sikes et al. (2016)Ronge et al. (2016)Yu et al. (2014b, a)(Olsen et al., 2016)Broecker et al. (1980)Key (2001)Sabine et al. (2004)Eide et al. (2017)
The late Holocene is chosen as the initial model calibration due to the absence of industrial-era CO2 and bomb radiocarbon. Scripps CO2 Program data originally sourced from http://scrippsco2.ucsd.edu (last access: 10 June 2018); data currently being transitioned to http://cdiac.ess-dive.lbl.gov (last access: 10 June 2018). The Peterson et al. (2014) database incorporates ∼500 core records across the LGM and late Holocene periods.
3.1 Pre-industrial calibration
We choose the late Holocene period (6–0.2 ka) for our calibration because it has relatively good proxy data coverage (e.g. Table 2) and a relatively steady climate in the absence of perturbations such as industrial CO2 emissions, bomb radiocarbon or glacial terminations. The late Holocene is also close to the pre-industrial period (1700s) in order to act as a starting point for modern carbon cycle simulations. To calibrate the model for the late Holocene we begin with the modern-day GLODAPv2 dataset (https://www.nodc.noaa.gov/ocads/oceans/GLODAPv2/, last access: 10 June 2018), which we average into the model's boxes on depth and latitude coordinates using one of the SCP-M scripts (Fig. 2). The GLODAPv2 database incorporates data from ∼1 million seawater samples from 700 cruises over the years 1972–2013, including data from the original GLODAP dataset, plus the CARINA and PACIFICA datasets (Olsen et al., 2016). We assume an average data year of 1990 for the data accumulated over the period 1972–2013. We make adjustments to the ocean concentrations of DIC, δ13C and Δ14C for the effects of industrial emissions (the “Suess” effect) and bomb radiocarbon in the atmosphere using published estimates (Broecker et al., 1980; Key, 2001; Sabine et al., 2004; Eide et al., 2017). For example, Eide et al. (2017) establish a mathematical relationship between Suess δ13C and CFC-12 in the ocean, which we applied using GLODAPv2 CFC-12 data to correct the ocean δ13C data. We force the model with late Holocene average data for atmosphere CO2, δ13C and Δ14C (data sources in Table 2). The model's starting parameters are set from literature values (Table A1, Appendix), including the point estimates for ocean circulation and mixing fluxes from Talley (2013).
Using the Suess- and bomb-adjusted GLODAPv2 ocean dataset and late Holocene atmosphere data as the starting point, combined with the literature-determined parameter values, the model is allowed to run freely for 15 kyr in spin-up. This spin-up is ample time for model equilibration and to allow slower processes such as carbonate compensation to take effect. The resulting model equilibrium ocean and atmosphere element concentrations from the spin-up are stored and then carried forward as the starting data for subsequent late Holocene simulations. Figure 3 shows the results of the model spin-up (red stars) compared with late Holocene atmosphere data and their standard error (blue dots and error bars) across the time period. We also show the model results compared with late Holocene ocean data from various sources (Table 2), which are averaged into the box model regions for comparison.
The late Holocene calibration convincingly satisfies the atmospheric data values for CO2, δ13C and Δ14C. Model results are also in good agreement with the late Holocene ocean Δ14C, falling within error or very close for all boxes covered by data. The surface boxes (1, 2) are relatively enriched in Δ14C relative to deeper boxes, reflecting their proximity to the atmospheric source of 14C, although the spread of values across the ocean boxes is narrow. The surface boxes (1, 2 and 7) intuitively display more enriched δ13C than the intermediate (3), deep (4) and abyssal (6) boxes, mainly due to the effects of the biological pump. For most of the model's boxes, the results fall within the standard error of the late Holocene δ13C data. The Southern Ocean box (5) is an exception due to its extensive vertical coverage of 2500 m, incorporating the surface boundary with the atmosphere and the deep ocean, coupled with the sparse δ13C core data for the polar Southern Ocean (one data point, no error bars). SCP-M also exaggerates the depletion in δ13C in the deep box (4) relative to the data observation.
There is limited data coverage for the carbonate ion proxy (), although the model replicates the available data well. concentrations can be interpreted as alkalinity less DIC (Zeebe and Wolf-Gladrow, 2001; Yu et al., 2014b) for the purpose of analysing model results charts. is relatively abundant in the surface boxes (e.g. boxes 1 and 2), reflecting the higher amount of alkalinity relative to carbon due to the soft tissue biological pump, which prioritises organic carbon export over alkalinity export. is less abundant in the deep ocean (boxes 4 and 6) because there is more carbon relative to alkalinity due to the remineralisation of organic matter, a pattern that SCP-M replicates.
3.2 Sensitivity tests
We undertook parameter sensitivity tests to understand changes in atmospheric CO2, Δ14C and δ13C in SCP-M. This serves two purposes: (1) to understand the directional relationship between the parameter values and these key model outputs, and evaluate whether they make sense, and (2) to inform the LGM–Holocene model–data experiment in the following section. For example, if the GOC parameter Ψ1 displays a negative relationship with atmospheric CO2, it would make sense to probe parameter values lower than modern to replicate the lower atmospheric CO2 in the LGM. We varied parameter values around their modern-day settings in 10 kyr model runs and plotted the output against atmospheric CO2, Δ14C and δ13C (Fig. 4).
For example, Fig. 4a–d shows sensitivity variations above and below the model's modern values for ocean circulation and mixing parameters, sourced from Talley (2013) and Toggweiler (1999). Atmospheric CO2 is very sensitive to Ψ1 and Ψ2 but displays limited response to γ1 and γ2 over the ranges analysed (Fig. 4a–d). Atmospheric Δ14C and δ13C are negatively related to Ψ1 and Ψ2. The slower ocean turnover leads to a reduced rate of upwelling and surface de-gassing of Δ14C- and δ13C-depleted waters, causing higher values in the atmosphere. The effect of the mixing parameters on the atmosphere variables is muted because they have a limited impact on the upwelling regime for carbon, with any upward flux of carbon offset by a downward flux (mixing).
The soft tissue export flux parameter, Z, displays an inverse relationship with CO2 (Fig. 4e). A higher global value of Z drives the removal of carbon from the surface ocean, and the resulting CO2 flux into the ocean from the atmosphere decreases CO2. Lower Z leads to increased atmospheric CO2. δ13C is particularly sensitive to Z, moving it well away from modern (and therefore Holocene and LGM) values from a minor perturbation. The rain ratio (Fig. 4f) increases pCO2 in the surface ocean boxes, leading to de-gassing of CO2 to the atmosphere, and therefore modestly decreasing atmospheric Δ14C, as the lighter 12C is preferentially partitioned across the air–sea interface.
Increasing surface ocean box temperature (Fig. 4g) increases atmospheric CO2, an intuitive outcome given that warmer water absorbs less CO2 (Weiss, 1974), and SCP-M employs a temperature- and salinity-dependent CO2 solubility function. Air–sea fractionation of δ13C also decreases with higher temperatures, leading to higher atmospheric δ13C. According to Mook et al. (1974), air-to-sea fractionation of δ13C (making the atmosphere more depleted in δ13C) increases at a rate of approximately 0.1 ‰ ∘C−1 of cooling. SCP-M employs temperature-dependent air–sea gas δ13C fractionation factors (Mook et al., 1974). Δ14C is invariant to box temperature as the fractionation parameters employed in the model are not temperature dependent. CO2 displays a weak positive relationship with surface ocean box salinity (Fig. 4h) due to the decreasing solubility of CO2 in ocean water with increasing salinity (Weiss, 1974).
As the polar box piston velocity P slows down (Fig. 4i), atmospheric CO2 falls. At lower values of P the polar box, a region of outgassing of CO2 due to the upwelling of deep-sourced carbon-rich water in that part of the ocean, will exchange CO2 with the atmosphere at a slower rate. The reduced outgassing of δ13C-depleted carbon to the atmosphere with a lower P leads to higher δ13C values in the atmosphere. Atmospheric Δ14C increases with a slowing of P as the pathways for it to invade the ocean from its atmospheric source are slower, and there is reduced outgassing of old, low Δ14C waters.
Terrestrial biosphere NPP (Fig. 4j) is a sink of CO2 and fractionates the ratios of the isotopes of carbon, leading to higher values for δ13C and, to a lesser extent, Δ14C in the atmosphere. It is likely that NPP plays a feedback role and modulates CO2, δ13C and Δ14C (Toggweiler, 2008). Varying the ocean surface area (Fig. 4k) has modest impacts on CO2 and δ13C, but a large impact on Δ14C. Decreasing the ocean volume leads to a lower surface area for CO2 and atmospherically produced radiocarbon to enter the ocean, causing them to increase in the atmosphere. We expect that changing the ocean surface area (from sea level), and therefore volume, leads to changes in pCO2 on glacial–interglacial timescales. Increasing the fraction of deep water upwelled into the subpolar box (Fig. 4l) intuitively raises CO2, but lowers δ13C and Δ14C, by upwelling carbon-rich and isotopically depleted water to the ocean surface boxes.
3.3 Modern carbon cycle simulation
Human fossil fuel and land use change emissions contributed ∼575 Gt carbon to the atmosphere between 1751 and 2010 (Boden et al., 2017; Houghton, 2010) and up until 2014 were growing at an accelerating rate. In response, the Earth's carbon cycle continually partitions carbon between its component reservoirs, with positive and negative feedbacks. The net effect is a build-up of carbon in most reservoirs. Given the dominance of the anthropogenic emissions source in the modern global carbon cycle, a simulation model should be able to provide a plausible replication of its effects.
We modelled the effects of anthropogenic emissions and atmospheric nuclear bomb testing on the carbon reservoirs and fluxes in SCP-M. The experiment forces the late Holocene–pre-industrial SCP-M equilibrium described in Sect. 3.1, with estimates of industrial fossil fuel and land use change CO2 emissions, sea surface temperature (SST) changes and atmospheric bomb 14C fluxes from historical data dating from 1751. For the future years to 2100, we force the model with the IPCC's Representative Concentration Pathway (RCP) CO2 emissions and SST scenarios (Boden et al., 2017; Houghton, 2010; IPCC, 2013a). We compare the model results with atmospheric CO2, δ13C and Δ14C historical data and published modelling results for future years from the CMIP5 project (https://cmip.llnl.gov/cmip5/, last access: 25 June 2018).
Figure 5 shows the modern carbon cycle simulation using SCP-M compared with historical atmospheric data for CO2, δ13C and Δ14C and GLODAPv2 ocean data (estimated data year 1990). Importantly, SCP-M provides an appropriate simulation of the carbon cycle response to human emissions inputs by replicating the atmospheric patterns for CO2, δ13C and Δ14C preserved in data observations for the period 1751–2016 (a–b). The atmospheric CO2 and δ13C data are sourced from the Scripps CO2 Program (originally sourced from http://scrippsco2.ucsd.edu, last access: 10 June 2018; data currently being transitioned to http://cdiac.ess-dive.lbl.gov, last access: 10 June 2018), and Δ14C data are sourced from Nydal and Lövseth (1996), Stuiver et al. (1998), and Turnbull et al. (2017). A key feature of the historical data is the substantial increase in human emissions from circa 1950 onwards, which is accompanied by higher atmospheric CO2 and a steep drop in δ13C (Fig. 5a); this reflects δ13C-depleted anthropogenic emissions. The effect of emissions on atmospheric Δ14C (Fig. 5b) in the 20th century is largely overprinted by the influence of bomb radiocarbon. Emissions are seen as a slight downturn in the model and data Δ14C in the immediate lead-up to the release of bomb radiocarbon into the atmosphere and the downward trend from ∼2020 onwards. The spike in Δ14C during the period of bomb radiocarbon release lasts during the period 1954–1963 and then disperses as 14C is absorbed by the ocean. The simulation shows that SCP-M agrees with the GLODAPv2 ocean data by 1990 (Fig. 5c–f), with most boxes falling within the standard deviation of average data values, lending confidence to the model's simulation of the redistribution of carbon.
Figure 6 shows the emissions profile (a) and modelling results for atmospheric CO2 (b) over historical time and projected forward to 2100 for the IPCC RCPs. The SCP-M output falls below the IPCC projections for atmospheric CO2 in RCP2.6 and 4.5, but provides a close match with RCP6.0 and 8.5.
Figure 7a shows the annual uptake of CO2 by the ocean, modelled with SCP-M. The model begins the period close to a steady state between the atmosphere and surface ocean pCO2, with limited transfer across the interface. Beginning circa 1950 the ocean takes up an increased load of CO2 from the atmosphere. By 2100, SCP-M models a range of annual CO2 uptake by the ocean of 0–6 PgC annum−1 across the RCPs. This is similar to the range of values estimated by the CMIP5 models (also shown in Fig. 7a), reproduced from Jones et al. (2013). The cumulative uptake of emissions by the ocean over the period 1751–2100 (Fig. 7b) modelled by SCP-M of ∼350–750 PgC is at the upper end of the modelled range of CMIP5 models of ∼200–600 PgC over the period 1850–2100 (Jones et al., 2013). The SCP-M simulations commence in 1751 and therefore incorporate an extra 100 years of fossil fuel and land use change emissions beyond the CMIP5 model results presented in Jones et al. (2013). Wang et al. (2016) quote a range of 412–649 PgC of cumulative uptake by the ocean by 2100 from 11 CMIP5 models, a closer range to the SCP-M outcomes.
Figure 8 shows the partitioning of anthropogenic CO2 emissions into the carbon cycle reservoirs in RCP6.0 by 2100, as simulated with SCP-M and compared with modelling results presented by the IPCC for the same scenario (IPCC, 2013b). By this time, the load of human emissions is roughly 45:55 split between the atmosphere and the combined terrestrial biosphere and ocean.
By 2100, in RCP6.0 the carbon cycle is substantially changed from the pre-industrial–late Holocene state. This is the result of the accumulation of hundreds of years of human industrial CO2 emissions in the atmosphere and other reservoirs (Fig. 9). Anthropogenic CO2 emissions transfer carbon to the atmosphere, ocean and terrestrial biosphere. The fluxes between the carbon reservoirs also change. In the pre-industrial state, CO2 enters the ocean in the low latitudes and northern ocean (shown as negative fluxes in Fig. 9) and de-gasses in the Southern Ocean (positive flux) under the influence of ocean upwelling in that region. In RCP6.0, the atmospheric CO2 concentration increases to the extent that the atmosphere–ocean pCO2 gradient drives all surface ocean boxes to take carbon from the atmosphere (shown as large negative changes in the air–sea fluxes of carbon; red text in Fig. 9), despite warmer surface ocean temperatures towards the end of the projection (time series inputs). The terrestrial biosphere influx of carbon is dramatically increased by the carbon fertilisation effect, leading to a larger biomass stock, which in turn also causes more respiration – both inward and outward biosphere fluxes of CO2 are therefore greatly enhanced. The weathering of silicate rocks on the continents, a weak sink of carbon, also accelerates under the effects of burgeoning atmospheric CO2, transferring carbon from the atmosphere to the ocean via rivers. The physical fluxes of carbon within the ocean are only modestly affected, with the main exception being low-latitude thermocline mixing, which in RCP6.0 mixes a larger amount of carbon back into the surface ocean box from intermediate depths. The altered balance of DIC : alkalinity, particularly in the abyssal box, leads to a decrease in the carbonate ion concentration of abyssal waters late in the projection period, which in turn causes more dissolution of marine sediments. By 2100 this feedback brings more carbon back into the ocean, increased from 0.2 to 1.1 PgC yr−1, but also alkalinity (in a ratio of 2:1 to DIC), thereby lowering whole-ocean pCO2 – a modest negative feedback. In summary, SCP-M provides an appropriate simulation of historical atmospheric CO2, δ13C and Δ14C data when forced with anthropogenic CO2 emissions data over the same period. For the forward-looking RCP emissions projections, SCP-M falls in the range of the CMIP5 models, although the oceanic carbon uptake is exaggerated for the RCP8.5 scenario. This result suggests that a more detailed experiment, for example with non-linear representation of the piston velocity with respect to atmospheric CO2 or prescribed feedbacks from ocean circulation and biology (e.g. Meehl et al., 2007; IPCC, 2013a, b; Moore et al., 2018), might provide a closer fit to the CMIP5 models.
The LGM–Holocene transition, and glacial–interglacial variations in the carbon cycle in general, remain outstanding problems in paleoceanography (e.g. Sigman and Boyle, 2000; Kohfeld and Ridgewell, 2009; Hain et al., 2010; Ferrari et al., 2014; Kohfeld and Chase, 2017). At issue is the precise cause of 80–100 ppm variations in atmospheric CO2 across glacial and interglacial periods. These CO2 oscillations are accompanied by striking changes in ocean and atmospheric carbon isotopes, oceanic carbonate ion distributions, and other paleo chemical indicators. Of particular interest is the transition from the LGM, ∼18–24 ka (Yokoyama et al., 2000; Clark et al., 2009; Hesse et al., 2011; Hughes et al., 2013; Hughes and Gibbard, 2015), to the Holocene (11.7 ka to present) due to the growing abundance of proxy data covering that period. The causes of abrupt atmospheric CO2 rise at the termination of the LGM, and continuing up to the Holocene period, remain definitively unresolved. The ocean is likely the main driver of atmospheric CO2 on the relevant timescale due to its relative size as a carbon reservoir (e.g. Broecker, 1982; Sarmiento and Toggweiler, 1984; Kohfeld and Ridgewell, 2009; Sigman et al., 2010), alongside changes in the terrestrial biosphere stock of carbon (Francois et al., 1999; Ciais et al., 2012; Peterson et al., 2014; Hoogakker et al., 2016). Active theories within the ocean realm include changes in ocean biology (Martin, 1990; Watson et al., 2000; Martinez-Garcia et al., 2014), ocean circulation and mixing (Sarmiento and Toggweiler, 1984; Toggweiler and Sarmiento, 1985; Toggweiler, 1999; Curry and Oppo, 2005; Anderson et al., 2009; Kohfeld and Ridgewell, 2009; De Boer and Hogg, 2014; Menviel et al., 2016; Muglia et al., 2018), sea-ice cover (Stephens and Keeling, 2000), whole-ocean chemistry (Broecker, 1982; Sigman et al., 2010), and composite hypotheses (Kohfeld and Ridgewell, 2009; Ferrari et al., 2014; Kohfeld and Chase, 2017). Other mechanisms proposed include ocean temperature, the terrestrial biosphere, ocean volume and shelf carbonates (Opdyke and Walker, 1992; Trent-Staid and Prell, 2002; Ridgewell et al., 2003; Ciais et al., 2012; Annan and Hargreaves, 2013). Each hypothesis listed above is supported by site-specific tracer observations (e.g. marine carbonate cores), regional data aggregation and review, literature synthesis, or modelling. Within the spectrum of hypotheses, a simple explanation of a carbon cycle mechanism or mechanisms remains elusive. Many of the early hypotheses were presented as independent or even competing in causality for the interglacial CO2 variation (Ferrari et al., 2014).
Substantial progress has been made over the last 15 years in constraining the list of likely candidates to ocean physical and biological processes, likely in concert. The growth of paleo datasets (e.g. Oliver et al., 2010; Peterson et al., 2014; Yu et al., 2014b; Skinner et al., 2017; Zhao et al., 2017) and improvements in computing power have enabled model and model–data studies which seek to constrain the magnitude of changes in the carbon cycle across glacial–interglacial cycles (e.g. Stephens and Keeling, 2000; Toggweiler et al., 2006; Tagliabue et al., 2009; Hain et al., 2010; Bouttes et al., 2011; Hesse et al., 2011; Tschumi et al., 2011; Menviel et al., 2016; Kurahashi-Nakamura et al., 2017; Muglia et al., 2018). For example, Menviel et al. (2016) modelled slowing GOC and AMOC, with a modest increase in biological productivity in the Southern Ocean in the LGM, using δ13C data and an intermediate-complexity Earth system model. This differed from the finding of Muglia et al. (2018), who specifically examined the AMOC and Southern Ocean biological productivity. They found that a weaker AMOC and stronger biological productivity could account for the LGM and Holocene δ13C, Δ14C and 15N data. The GOC was not tested by Muglia et al. (2018). Kurahashi-Nakamura et al. (2017) contradicted both studies, diagnosing a more vigorous (but shallower) AMOC in the LGM, using a GCM with data assimilation of various proxies, notably only incorporating Atlantic data for the LGM.(Trent-Staid and Prell, 2002; Annan and Hargreaves, 2013)(Adkins et al., 2002)(Stephens and Keeling, 2000; Ferrari et al., 2014)(Adkins et al., 2002; Grant et al., 2014)(Mariotti et al., 2013)
As shown in the sensitivity tests in Fig. 4, some processes do not exert a strong influence on atmospheric CO2, but do modestly impact CO2 and strongly impact δ13C and Δ14C. Where these features are posited to vary around glacial cycles, we have incorporated them as a step change from late Holocene–modern estimates in our LGM experiment.
4.2 Model–data experiments
We illustrate SCP-M's capabilities by solving for the parameter values of best fit with late Holocene and LGM ocean and atmosphere proxy data using a comprehensive model results–data optimisation. For this illustrative example, the atmosphere and ocean data are taken from published sources (Table 2), averaged for the LGM (∼18–24 ka) and late Holocene (6.0–0.2 ka) time periods and for box coordinates in SCP-M for the ocean data (depth and latitude). The mean and variance for each box are then calculated in SCP-M. First, we probe the potential for key model parameters to drive Holocene–LGM changes in atmospheric carbon variables to focus our experiment on these parameters. It is likely that the LGM-to-Holocene carbon cycle changes were dominated by the ocean (Sigman and Boyle, 2000; Kohfeld and Ridgewell, 2009), but were also accompanied by a range of physical changes in the atmosphere and terrestrial biosphere that in aggregate could be material (e.g. Sigman and Boyle, 2000; Adkins et al., 2002; Kohfeld and Ridgewell, 2009; Ferrari et al., 2014). These changes include sea surface temperature, salinity, sea-ice cover, ocean volume and atmospheric 14C production rate. Estimates of average sea surface temperature for the LGM generally fall in the range of 3–8 ∘C cooler than the present (Trent-Staid and Prell, 2002; Annan and Hargreaves, 2013). Adkins et al. (2002) estimated that ocean salinity was 1–2 psu higher in the LGM and sea levels were ∼120 m lower (Adkins et al., 2002; Grant et al., 2014). Stephens and Keeling (2000) and Ferrari et al. (2014) highlighted the role of expanded sea-ice cover in the Southern Ocean during the LGM as a key part of the LGM CO2 drawdown. Finally, Mariotti et al. (2013) estimated that higher atmospheric radiocarbon production accounted for ‰ in atmospheric Δ14C in the LGM. Mariotti et al. (2013) simulated this variation in model experiments by increasing the radiocarbon production rate by a multiple of 1.15–1.30 (best guess 1.25) of the modern estimate in order to recreate LGM Δ14C values. Using these findings we define two background states for modelling purposes: a late Holocene state (Table A1 in the Appendix) and the LGM state (Table 3). Figure 10 shows the cumulative effect of these changes in SCP-M, within the late Holocene–LGM atmosphere 3-D CO2–δ13C–Δ14C data space.
These changes are the first stage of a model adjustment to analyse the potential for ocean circulation and biological changes to deliver the LGM atmospheric CO2, δ13C and Δ14C values and transition the model output from the red circle (late Holocene) to the black star (the LGM background settings), and then to the black circle (LGM). The decrease in ocean surface box temperatures leads to a drop in CO2 of ∼20 ppm and a lightening of δ13C by ∼0.6 ‰ owing to the increased solubility of CO2 in colder water and the increasing fractionation of δ13C with decreasing temperatures, which leaves more 12C in the atmosphere. There is limited impact on Δ14C. Increasing salinity slightly reverses these changes to CO2 and δ13C. Reducing sea surface area and volume slightly increases CO2 and increases Δ14C as the ocean's capacity to take up these elements is reduced. Slowing down the piston velocity in the polar Southern Ocean box, as a proxy for increased sea-ice cover, slightly reduces CO2 (reduced outgassing), increases Δ14C (slower rate of invasion to the ocean) and increases δ13C (reduced outgassing and sea-to-air fractionation of δ13C). Finally, increasing the rate of atmospheric radiocarbon production forces a shift in Δ14C (horizontal shift in Fig. 10) towards LGM levels (black star and circle in Fig. 10). In aggregate, these changes lead to a fall in CO2 of ∼35 ppm, a fall in δ13C of ‰ and an increase in Δ14C of ∼300 ‰.
From the black star in Fig. 10, the “LGM state”, we perform a focussed sensitivity test on key hypothesised drivers of LGM–Holocene carbon cycle changes (Kohfeld and Ridgewell, 2009; Sigman et al., 2010). These are slower GOC (Ψ1), slower AMOC (Ψ2), reduced deep–abyssal ocean mixing (γ1) and a stronger biological pump (Z). The Z global biological production parameter, varied across 5–10 mol C m−2 yr−1 (i.e. increased), can deliver the LGM CO2 changes, but steers δ13C and Δ14C away from their LGM values. γ1 drives ancillary changes in all three variables, suggesting it is not the driver of the LGM atmospheric changes but may play a modulating role. Both Ψ1 (3–29 Sv) and Ψ2 (3–19 Sv) experiments run very close to the LGM data values on their own, although neither can deliver a precise hit.
Using the literature-referenced Holocene and LGM background parameter states, and informed by the sensitivity analysis in Fig. 10, we take advantage of SCP-M's fast run time to perform thousands of multi-variant simulations over the free-floating Ψ1, Ψ2, γ1 and Z parameter spaces using the SCP-M batch module. We then perform an optimisation routine against the data for each period to solve for values for Ψ1, Ψ2, γ1 and Z. The SCP-M batch module cycles through each set of parameter combinations, with each model simulation run for 10 000 years. Table 4 shows the experiment parameter ranges for the late Holocene and LGM model–data experiments.
The parameter input ranges for the experiments were informed by the sensitivity tests shown in Figs. 4 and 10. For example, the responses of atmospheric CO2, δ13C and Δ14C to variations in Ψ1, Ψ2 and Z lead us to cater for lower values for Ψ1 and Ψ2 (weaker ocean circulation) and higher values for Z (increased biological productivity) in the LGM experiment. Where the experiments resulted in a parameter output value at the limit of the input range, the range was widened and the experiment was repeated; 16 896 and 47 616 simulations were undertaken for the late Holocene and LGM experiments, respectively.
Bold font parameter results indicate those parameters that are free-floating and determined by the model and data in the experiment. The LGM experiment shows a marked decline in the strength of global overturning circulation Ψ1 (−12 Sv) and a modest decline in Atlantic meridional overturning circulation Ψ2 to deliver the LGM atmosphere and ocean data signal. A minor increase in deep–abyssal mixing γ1 is also seen.
The SCP-M script harvests model output and performs a least-squares optimisation of the output against the LGM and late Holocene data for atmospheric CO2, atmospheric, and ocean Δ14C and δ13C, and also the oceanic carbonate ion proxy, to source the best-fit parameter values for Ψ1, Ψ2, γ1 and Z (or any parameter specified):
where Optn=1 represents the optimal value of parameters n, Ri,k is the model output for concentration of each element i in box k, Di,k is the average data concentration of each element i in box k and σi,k is the standard deviation of the data for each element i in box k. The standard deviation performs two roles. It reduces the weighting of data with high uncertainty and also normalises for the different unit scales (e.g. ppm, ‰ and µmol kg−1), which allows multiple proxies in different units to be incorporated in the optimisation. Where data are unavailable for a box, that element and box combination is automatically omitted from the optimisation routine.
The late Holocene data-optimised results for Ψ1 (30 Sv) and Ψ2 (18 Sv) show good agreement with the Talley (2013) observations for Ψ1 (29 Sv) and Ψ2 (19 Sv) from the modern ocean (Table 5). The starting global value of Z, of 5 mol C m−2 yr−1, is returned in the experiment. The experiment also successfully returns values for atmospheric CO2, δ13C and Δ14C within standard error for the late Holocene data series (Table 5).
The ocean and atmosphere SCP-M results for the LGM (bold stars) and late Holocene (transparent stars) experiments using the optimised parameter settings in Table 5 are plotted in Fig. 11 along with the corresponding data (blue dots with error bars). The experiment provides results within the error bounds of data for most of the box regions in both scenarios and an excellent fit to the change in the relative distribution of the proxies between ocean boxes and the atmosphere, which is preserved in the LGM and late Holocene data. A key feature of the ocean δ13C data is a depletion of deep ocean δ13C in the LGM, shown as a drop in δ13C values in the deep (box 4) and abyssal (box 6) boxes relative to the intermediate box (3). In the LGM δ13C data, there is a spread of 1 ‰ across these water masses, which narrows to 0.3 ‰ in the late Holocene data. The pattern is replicated in the LGM experiment, pointing to the important role of changes in abyssal–deep ocean water flows via Ψ1 in delivering the ocean δ13C data patterns. The model shift in δ13C in the deep box (box 4) of 0.6 ‰ from the LGM to late Holocene is in good agreement with a global deepwater estimate of 0.49±0.23 ‰ by Gebbie et al. (2015) and an earlier estimate of 0.46 ‰ by Curry et al. (1988). The average atmospheric δ13C value remains largely unchanged between the two periods due to the effect of the terrestrial biosphere, which causes net uptake of CO2 in the Holocene period (increases atmospheric δ13C) and net respiration of CO2 in the LGM period (decreases atmospheric δ13C).
The model results also closely replicate the reduction in deep-to-shallow ocean compositional gradient in Δ14C data moving from the LGM to Holocene period (e.g. Skinner and Shackleton, 2004; Skinner et al., 2015, 2010; Burke and Robinson, 2012; Chen et al., 2015; Hines et al., 2015; Ronge et al., 2016). The LGM data show a spread of ∼300 ‰ between abyssal (box 6) and intermediate (box 3) waters, as well as deep (box 4) versus surface (boxes 1, 2 and 7) boxes. In the late Holocene data, the spread is narrowed to ∼100 ‰. This data observation was popularly characterised as the result of increased Southern Ocean upwelling of Δ14C-depleted deep water into intermediate and shallow depths in the Holocene (e.g. Skinner et al., 2010, 2015; Burke and Robinson, 2012). A slowdown in Southern Ocean upwelling in the LGM allows Δ14C-depleted water to accumulate in the deep or abyssal ocean and a widening in the Δ14C gradient between deep and shallow waters. In SCP-M, this is simulated by lower values for Ψ1 and Ψ2. The low-latitude surface box (box 1) enrichment in Δ14C in planktonic foraminifera in the LGM is replicated by the increased atmospheric production rate of radiocarbon applied to the LGM experiment, combined with slower ocean circulation.
SCP-M results are shown for comparison with sparse carbonate proxy data (Fig. 11 bottom panel). The model results for the carbonate ion proxy mirror the limited variation in the data between the LGM and late Holocene. The changes are most pronounced in the surface boxes (boxes 1 and 2), which are under the influence of atmospheric CO2, and attenuate somewhat in the deeper boxes (boxes 4 and 6). Yu et al. (2014b) interpreted the relatively small changes in the carbonate ion in the deepest ocean (box 6) as the result of efficient buffering of deepwater pH by carbonate dissolution, most notably in the Pacific Ocean. The model result for the deep box (box 4) goes against the LGM–Holocene variation in the data, but given that there is only one data point for this part of the ocean and the variation itself is small, it is an uncertain outcome.
The LGM scenario shows important changes in the carbon redistributive behaviour of the ocean (Fig. 12). The stock of carbon increases in abyssal and deep boxes (blue text denotes the increase in PgC from late Holocene to LGM) and is reduced in the intermediate, low-latitude surface and northern surface boxes (red text denotes the decrease in PgC from late Holocene to LGM). The amount of carbon upwelled to the subpolar surface and deep boxes by GOC (Ψ1) drops by ∼5–10 PgC yr−1, with the most pronounced changes taking place at the abyssal–deep box boundary. The slower upwelling rate of carbon causes a reduced outgassing rate of CO2 from the subpolar box to the atmosphere. The weaker flux of Ψ2 also brings a reduced DIC load into the intermediate-depth ocean, the driver for lower DIC content in the intermediate and surface boxes. The optimised parameter run for the late Holocene results in a terrestrial biosphere carbon pool of 2495 PgC, which is fortuitously close to the pre-industrial estimate of Raupach et al. (2011) (2496 PgC), at the top end of acceptable values in Francois et al. (1999) and close to the “active” land carbon pool of 2370±125 estimated by Ciais et al. (2012). In the optimised LGM results, the terrestrial biosphere is reduced by 667 Pg C from the late Holocene value, which is towards the upper bound of recent estimates of this change (0–700 PgC e.g. Ciais et al., 2012; Peterson et al., 2014), but within uncertainty bounds. For example, Peterson et al. (2014) estimated a variation of 511±289 Pg C in the terrestrial biosphere carbon stock based on whole-ocean δ13C data, the same data used in this exercise. According to Francois et al. (1999), palynological and sedimentological data imply that the terrestrial biosphere carbon stock was 700–1350 PgC smaller in the LGM than the present. Ciais et al. (2012) pointed to the growth of a large inert carbon pool in steppes and tundra during the LGM, which may have modulated some of the active biosphere carbon signal (i.e. reduced terrestrial biosphere), a factor not explicitly covered in our modelling exercise. The terrestrial biosphere is clearly a key part of the LGM–late Holocene carbon cycle transition. The atmosphere-enriching fractionation of δ13C by the terrestrial biosphere during the deglacial period effectively reverses the effects of the release of δ13C-depleted carbon from the deep ocean to the atmosphere at the termination and leaves atmosphere δ13C almost unchanged from LGM values as a result (Schmitt et al., 2012). The DIC : alkalinity balances in the abyssal ocean during the LGM also drive subtle changes in the balance of carbonate outflux by sinking and influx from sediment dissolution, which build up to substantial differences in the sediment carbon stock between the LGM and Holocene simulations, mainly due to the time frames modelled in the SCP-M spin-up for each scenario (15 kyr).
5.1 Model advantages and limitations
In this paper we introduce SCP-M, a box model of the global carbon cycle. We demonstrate its application to the modern and future carbon cycle with anthropogenic emissions and reconstructing potential changes across the LGM–late Holocene carbon cycle transition. In summary, SCP-M is a simple, easy-to-use model of the carbon cycle, and its fast run time enables comprehensive scenario analysis or optimisations for scenario or hypothesis testing. It takes approximately 30 s to complete a 10 000-year simulation, making the model useful for long-term paleo reconstructions of the carbon cycle. Our LGM–late Holocene experiment (Sect. 4) includes broad variations in GOC, AMOC, deep–abyssal mixing and global biological productivity. The experiments cover up to ∼47 000 parameter combinations across the LGM and late Holocene proxy data, reducing the possibility of confirmation bias or predetermined outcomes. Furthermore, the model's simplified topology (Fig. 1), although consistent with an observationally based understanding of the ocean (e.g. Talley, 2013), makes it accessible to a wide user group and potentially useful as a teaching aid to illustrate high-level concepts in the carbon cycle. The model contains data modules that directly integrate data via box mapping and averaging processes for use in calibration and for model–data experiments (Fig. 2). It also includes a model–data optimisation routine to elicit parameter values that best fit the data, as described in Sect. 4.2.
A limitation of SCP-M v1.0 is that it does not distinguish between the Atlantic and Indo-Pacific Ocean basins, which is a large simplification. We argue that this is feasible for testing high-level hypotheses, for example involving large-scale ocean processes across the LGM–Holocene time periods, and the model is demonstrated to produce appropriate results in such an application. However, this framework may not be useful for testing localised or detailed problems. Given that SCP-M is a box model, there are other simplifications, including a rigid and perhaps somewhat arbitrary treatment of box boundaries (Fig. 1). For some experiments the box boundaries themselves may need to be a dynamic, model-determined output. In our LGM–Holocene example, we did not vary the abyssal box thickness across the time periods. However, this could be done very easily for the case of an experiment featuring shoaling or deepening water mass boundaries (e.g. Curry and Oppo, 2005). A key drawback of the model is that it can identify the cause of changes in proxy element concentrations in terms of parameterised processes, but cannot diagnose the root cause of these process changes. For example, with SCP-M we cannot directly answer the question of what causes GOC, AMOC or biological productivity to vary on glacial–interglacial cycles, but combined with data we can propose which of these varies (Sect. 4.2).
5.2 Modern carbon cycle simulations
Our simple forcing of SCP-M with historical and projected anthropogenic CO2 emissions and SST (Sect. 3.3, Figs. 6–8) shows that SCP-M can reproduce historical data and the model results of more complex CMIP5 models. The SCP-M results for atmospheric CO2 (Fig. 6), air–sea fluxes of carbon (Fig. 7) and accumulation of carbon in the various carbon reservoirs (Fig. 8) line up in the range of CMIP5 model projections. Importantly, SCP-M is shown to replicate the historical data over the period 1751–2016 for atmospheric CO2, δ13C and Δ14C (Figs. 5 and 6). The historical period is a stringent test of carbon cycle models because it incorporates the influences of anthropogenic emissions, atmospheric bomb testing and an abundance of data observations of the Earth's carbon cycle response for comparison. The radioactive decay and dispersal of bomb-produced 14C provides an excellent “time clock” for the fluxes in the carbon cycle, particularly air–sea gas exchange and ocean circulation. Our experiment incorporates forcing of atmospheric 14C during 1954–1963, and SCP-M appropriately replicates the take-up of bomb 14C by the ocean from the atmosphere in the following years (Fig. 5).
However, the SCP-M modern and future simulations are simple and fail to take account of more complex potential feedbacks in the carbon cycle. These may include a wind-shift-induced slowing of AMOC and thermocline mixing or a response of ocean biological productivity to changed pCO2, temperature and DIC in the surface ocean (e.g. Meehl et al., 2007; IPCC, 2013a, b; Moore et al., 2018). To simulate such feedbacks, the relevant parameters would need to be forced in SCP-M. A dynamical response would be expected in more complex Earth system models. The value of SCP-M is in rapidly undertaking “what-if” types of analysis to probe the effects of such changes under a variety of scenarios. For example, the high-level testing of negative emissions processes, such as alkalinity or iron seeding of the ocean, rock waste fertilisation, afforestation and/or reforestation on land, or marine fauna management, as tools for reducing atmospheric CO2 on a global scale are feasible uses of SCP-M.
5.3 LGM–late Holocene modelling
Our model–data optimisation using SCP-M and published data (Sect. 4.2) showed that variations in the strength of large-scale ocean physical processes, particularly GOC and AMOC, can account for the LGM-to-Holocene carbon cycle changes inferred in the proxy data (Table 5). Critically, the variations are accompanied by changes in SST, sea-ice cover and the terrestrial biosphere (Table 3). The result is observed on account of ocean and atmosphere data across the proxies of CO2, δ13C, Δ14C and the carbonate ion (Fig. 11). This finding is not a new one, but corroborates the model–data conclusion of Menviel et al. (2016), box modelling of Toggweiler (1999), and 14C proxy data findings of Sikes et al. (2000) and Skinner et al. (2017). The importance of GOC is at odds with Muglia et al. (2018), who found a substantially weakened AMOC and enhanced biological productivity in the Southern Ocean in the LGM, with no role examined for GOC in that study. Kurahashi-Nakamura et al. (2017) had an altogether different finding, modelling a stronger yet shallower AMOC during the LGM. Many such studies focus exclusively on the Atlantic Ocean, perhaps due to the well-understood AMOC and the more detailed proxy data coverage in that basin. For example, Curry and Oppo (2005) provided striking transects of δ13C in the LGM and late Holocene Atlantic Ocean, which provided evidence for large changes in the basin δ13C stratigraphy across the two time periods. The δ13C data compilations of Oliver et al. (2010) and Peterson et al. (2014) are also heavily skewed towards the Atlantic Basin.
Talley (2013) emphasised the importance of the Pacific and Indian Ocean overturning circulation limbs in the global ocean circulation regime, which implies that it is a critical part of the Earth's carbon cycle. This finding was corroborated by Skinner et al. (2017) in a recent review of Pacific Ocean radiocarbon data. The model–data results using SCP-M suggest that GOC was substantially reduced during the LGM (Table 5), accompanying enhanced storage of isotopically depleted carbon in the abyssal and deep ocean from atmospheric and terrestrial biosphere sources (Fig. 11). We posit that the release of volumes of carbon, that were greater than amounts stored in the deep–abyssal Atlantic alone, caused the atmospheric CO2 increase at the last glacial termination. Such a movement of carbon to and from the global abyssal ocean is likely required due to the large, opposite movement in atmospheric CO2 from the terrestrial biosphere across the same time period (Fig. 12). During the LGM, the terrestrial biosphere was reduced relative to the modern period, which was a source of CO2 to the atmosphere, and rebounded from the LGM to the Holocene, becoming a sink of CO2 during that period (Francois et al., 1999; Ciais et al., 2012; Peterson et al., 2014; Hoogakker et al., 2016). Incorporating the terrestrial biosphere into modelling experiments increases the magnitude of carbon uptake–release required from the ocean to satisfy the LGM and late Holocene atmospheric CO2 and, critically, δ13C data (even when incorporating SST, salinity, sea-ice cover proxy and ocean volume changes). The finding underscores the importance of incorporating multiple data proxies and carbon reservoirs into glacial–interglacial carbon cycle modelling.
Our model–data experiments did not find a role of changed marine biological production in the LGM–late Holocene transition (Table 5). However, this finding was the result of testing for variations in the global value of the ocean biological productivity, impacting all surface ocean boxes in SCP-M. Other studies (e.g. Menviel et al., 2016; Muglia et al., 2018) focussed specifically on the Southern Ocean biological productivity and identified its potential role in the LGM atmospheric CO2 drawdown. The Southern Ocean marine biology, in particular, is posited as a candidate for driving glacial–interglacial cycles of CO2 (e.g. Martin, 1990; Martinez-Garcia et al., 2014).
The SCP-M carbon cycle box model was constructed for the purposes of scenario or hypothesis testing (quickly and easily), model–data integration and inversion, paleo reconstructions, and analysing the distribution of anthropogenic emissions in the carbon cycle. The model contains a full ocean–atmosphere–terrestrial carbon cycle with a realistic treatment of ocean processes. Despite being relatively simple in concept and construct, SCP-M can account for a range of paleo and modern carbon cycle observations. The model applications illustrated here include integration with datasets from the present day (GLODAPv2, IPCC) and ocean paleo proxy data across the LGM and late Holocene periods. Simulations of the modern carbon cycle indicate that SCP-M provides a realistic representation of the dynamic shocks from human industrial and land use change emissions and bomb 14C. A model–data experiment using LGM and late Holocene CO2, δ13C, Δ14C and the carbonate ion proxy is able to resolve parameter values for ocean circulation, mixing and biology while reproducing model results that are very close to the proxy data for both time periods. The results indicate that the LGM-to-Holocene carbon cycle transition can be explained by variations in the strength of the global overturning circulation and Atlantic meridional overturning circulation when combined with a number of background changes, such as sea surface temperature, salinity, sea-ice cover, ocean volume and a varied atmospheric radiocarbon production rate. Further work on data quality and analysis is required to validate this finding, which is the subject of a separate paper. The results show promise in helping to further resolve the LGM-to-Holocene carbon cycle transition and point towards an ongoing application for data-constrained models such as SCP-M.
CO undertook model building, data gathering, modelling and model–data experiments. AH provided the oceanographic interpretation, supervised model design and modelling work, and designed the model–data experiments. ME provided input into model design, designed model–data experiments, and oversaw the modelling of the marine biology and isotopes. BO provided input into model design and oversaw the modelling of carbonate chemistry, marine sediments and interpretation of LGM–Holocene hypotheses. SE designed model–data experiments and oversaw the modelling of the marine biology and carbonate pump. All authors contributed to drafting and reviewing the document.
The authors declare that they have no conflict of interest.
This paper was substantially improved by input from two anonymous reviewers and the GMD topical editor. Stewart Fallon provided guidance and spreadsheet tools for the processing of radiocarbon data. Malcolm Sambridge provided input on model–data optimisation and inversions. Jimin Yu provided helpful discussions and carbonate ion proxy data.
This paper was edited by Paul Halloran and reviewed by two anonymous referees.
Amiotte Suchet, P., Probst, J., and Ludwig, W.: Worldwide distribution of continental rock lithology: Implications for the atmospheric/soil CO2 uptake by continental weathering and alkalinity river transport to the oceans, Global Biogeochem. Cy., 17, 7–1 – 7–14, 2003. a
Anderson, L. A. and Sarmiento, J. L.: Redfield ratios of remineralization determined by nutrient data analysis, Global Biogeochem. Cy., 8, 65–80, 1994. a
Anderson, R., Ali, S., Bradtmiller, L., Nielsen, S., Fleisher, M., Anderson, B., and Burckle, L.: Wind-driven upwelling in the Southern Ocean and the deglacial rise in atmospheric CO2, Science, 323, 1443–1448, 2009. a
Anderson, R. F., Fleishera, M. Q., Laoc, Y., and Winckler, G.: Modern CaCO3 preservation in equatorial Pacific sediments in the context of late-Pleistocene glacial cycles, Mar. Chem., 111, 30–46, 2007. a
Archer, D. and Maier-Reimer, E.: Effect of deep-sea sedimentary calcite preservation on atmospheric CO2 concentration, Nature, 367, 260–263, 1994. a
Barker, S., Knorr, G., Vautravers, M., Diz, P., and Skinner, L.: Extreme deepening of the Atlantic overturning circulation during deglaciation, Nat. Geosci., 3, 567–571, 2010. a
Battaglia, G., Steinacher, M., and Joos, F.: A probabilistic assessment of calcium carbonate export and dissolution in the modern ocean, Biogeosciences, 13, 2823–2848, https://doi.org/10.5194/bg-13-2823-2016, 2016. a, b
Boden, T., Marland, G., and Andres, R.: Global, Regional, and National Fossil-Fuel CO2 Emissions, Carbon Dioxide Information Analysis Center, Oak Ridge National Laboratory, US Department of Energy, 2017. a, b
Bouttes, N., Paillard, D., Roche, D. M., Brovkin, V., and Bopp, L.: Last Glacial Maximum CO2 and δ13C successfully reconciled, Geophys. Res. Lett., 38, LO2705, https://doi.org/10.1029/2010GL044499, 2011. a
Broecker, W. S.: The great ocean conveyor, Oceanography, 4, 79–89, 1991. a
Bryan, S., Marchitto, T., and Lehman, S.: The release of 14C-depleted carbon from the deep ocean during the last deglaciation: Evidence from the Arabian Sea, Earth Planet. Sc. Lett., 298, 244–254, 2010. a
Chen, T., Robinson, L., Burke, A., Southon, J., Spooner, P., Morris, P., and Ng, H.: Synchronous centennial abrupt events in the ocean and atmosphere during the last deglaciation, Science, 349, 1537–1541, 2015. a, b
Ciais, P., Tagliabue, A., Cuntz, M., Bopp, L., Scholze, M., Hoffmann, G., Lourantou, A., Harrison, S. P., Prentice, I. C., Kelley, D. I., Koven, C., and Piao, S. L.: Large inert carbon pool in the terrestrial biosphere during the Last Glacial Maximum, Nat. Geosci., 5, 74–79, 2012. a, b, c, d, e, f, g
Clark, P., Dyke, A., Shakun, J., Carlson, A., Clark, J., Wohlfarth, B., Mitrovica, J., Hostetler, S., and McCabe, A.: The Last Glacial Maximum, Science, 325, 710–714, 2009. a
Compton, J., Mallinson, D., Glenn, C., and Zanin, Y.: Variations in the global phosphorus cycle, Society of Sedimentary Geology, 66, 21–33, 2000. a
Craig, H.: The natural distribution of radiocarbon and the exchange time of carbon dioxide between atmosphere and sea, Tellus, 9, 1–17, 1957. a
Curry, W., Duplessy, J., Labeyrie, L., and Shackleton, N.: Changes in the distribution of δ13C of deep water PCO2 between the last glaciation and the Holocene, Paleoceanography, 3, 317–341, 1988. a
Curry, W. B. and Oppo, D. W.: Glacial water mass geometry and the distribution of δ13C of CO2 in the western Atlantic Ocean, Paleoceanography, 20, PA1017, https://doi.org/10.1029/2004PA001021, 2005. a, b, c, d
Davies-Walczak, M., Mix, A., Stoner, J., Southon, J., Cheseby, M., and Xuan, C.: Late Glacial to Holocene radiocarbon constraints on North Pacific Intermediate Water ventilation and deglacial atmospheric CO2 sources, Earth Planet. Sc. Lett., 397, 57–66, 2014. a
de Boyer Montegut, C., Madec, G., Fischer, A. S., Lazar, A., and Iudicone, D.: Mixed layer depth over the global ocean: An examination of profile data and a profile based climatology, J. Geophys. Res., 109, C12003, https://doi.org/10.1029/2004JC00237, 2004. a, b, c, d, e
Dickson, R. R. and Brown, J.: The production of North Atlantic Deep Water: Sources, rates, and pathways, J. Geophys. Res., 99, 12319–12341, 1994. a
Emerson, S. and Hedges, J. I.: Sediment diagenesis and benthic flux, in: Treatise on Geochemistry, edited by: Holland, H. D. and Turekian., K. K., vol. 6, chap. 9, Elsevier, Amsterdam, 2003. a
Ferrari, R., Jansen, M., Adkins, J., Burke, A., Stewart, A. L., and Thompson, A.: Antarctic sea ice control on ocean circulation in present and glacial climates, P. Natl. Acad. Sci. USA, 111, 8753–8758, 2014. a, b, c, d, e, f, g
Fontugne, M. and Duplessy, J.: Carbon Isotope Ratio of Marine Phytoplankton related to surface water masses, Earth Planet. Sc. Lett., 41, 365–371, 1978. a
Francois, L., Godderis, Y., Warnant, P., Ramstein, G., de Noblet, N., and Lorenz, S.: Carbon stocks and isotopic budgets of the terrestrial biosphere at mid-Holocene and last glacial maximum times, Chem. Geol., 159, 163–199, 1999. a, b, c, d, e, f
Gebbie, G., Peterson, C., Lisiecki, L., and Spero, H.: Global-mean marine δ13C and its uncertainty in a glacial state estimate, Quaternary Sci. Rev., 125, 144–159, 2015. a
Gloege, L., McKinley, G. A., Mouw, C. B., and Ciochetto, A. B.: Global eva- luation of particulate organic carbon flux parameterizations and implications for atmospheric pCO2, Global Biogeochem. Cy., 31, 1192–1215, 2017. a
Gordon, A.: The role of thermohaline circulation in global climate change, in: Lamont-Doherty Geological Observatory Report 1990/91, Tech. rep., Lamont-Doherty Geological Observatory of Columbia University, Palisades, New York, 1991. a
Grant, K. M., Rohling, E. J., Ramsey, C. B., Cheng, H., Edwards, R. L., Florindo, F., Heslop, D., Marra, F., Roberts, A. P., Tamisiea, M. E., and Williams, F.: Sea-level variability over five glacial cycles, Nat. Commun., 5, 1–9, 2014. a, b
Hain, M. P., Sigman, D. M., and Haug, G. H.: Carbon dioxide effects of Antarctic stratification, North Atlantic Intermediate Water formation, and subantarctic nutrient drawdown during the last ice age: Diagnosis and synthesis in a geochemical box model, Global Biogeochem. Cy., 24, GB4023, https://doi.org/10.1029/2010GB003790, 2010. a, b, c, d, e
Harman, I., Trudinger, C., and Raupach, M.: SCCM – the Simple Carbon-Climate Model: Technical Documentation, CAWCR Technical Report 047, CSIRO Centre for Australian Weather and Climate Research, CSIRO Marine and Atmospheric Research, FC Pye Laboratory, G.P.O. Box 3023, Canberra, ACT, 2601, Australia, 2011. a, b
Harrison, W., Head, E., Horne, E., Irwin, B., Li, W., Longhurst, A., Paranjape, M., and Platt, T.: The western North Atlantic bloom experiment, Deep-Sea Res. Pt. II, 40, 279–305, 1993. a
Henson, S. A., Sanders, R., Madsen, E., Morris, P. J., Moigne, F. L., and Quartly, G. D.: A reduced estimate of the strength of the ocean's biological carbon pump, Geophys. Res. Lett., 38, L04606, https://doi.org/10.1029/2011GL046735, 2011. a
Hoogakker, B. A. A., Smith, R. S., Singarayer, J. S., Marchant, R., Prentice, I. C., Allen, J. R. M., Anderson, R. S., Bhagwat, S. A., Behling, H., Borisova, O., Bush, M., Correa-Metrio, A., de Vernal, A., Finch, J. M., Fréchette, B., Lozano-Garcia, S., Gosling, W. D., Granoszewski, W., Grimm, E. C., Grüger, E., Hanselman, J., Harrison, S. P., Hill, T. R., Huntley, B., Jiménez-Moreno, G., Kershaw, P., Ledru, M.-P., Magri, D., McKenzie, M., Müller, U., Nakagawa, T., Novenko, E., Penny, D., Sadori, L., Scott, L., Stevenson, J., Valdes, P. J., Vandergoes, M., Velichko, A., Whitlock, C., and Tzedakis, C.: Terrestrial biosphere changes over the last 120 kyr, Clim. Past, 12, 51–73, https://doi.org/10.5194/cp-12-51-2016, 2016. a, b, c
Hughes, P. and Gibbard, P.: A stratigraphical basis for the Last Glacial Maximum (LGM), Quatern. Int., 383, 174–185, 2015. a
Hughes, P., Gibbard, P., and Ehlers, J.: Timing of glaciation during the last glacial cycle: evaluating the concept of a global “Last Glacial Maximum” (LGM), Earth-Sci. Rev., 125, 171–198, 2013. a
IPCC: Annex II: Climate System Scenario Tables, Climate Change 2013: The Physical Science Basis, Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, 1395–1445, Cambridge University Press, Cambridge, United Kingdom and New York, USA, 2013a. a, b, c, d
IPCC: IPCC, 2013: Climate Change 2013: The Physical Science Basis, Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, Cambridge University Press, Cambridge, United Kingdom and New York, USA, 2013b. a, b, c, d
Jansen, H., Zeebe, R., and Wolf-Gladrow, D.: Modeling the dissolution of settling CaCO3 in the ocean, Global Biogeochem. Cy., 16, 11-1–11-15, 2002. a
Jones, C., Robertson, E., Arorab, V., Friedlingstein, P., Shevliakova, E., Bopp, L., Brovkin, V., Hajima, T., Kato, E., Kawamiya, M., S., Liddicoat, Lindsay, K., Reick, C., Roelandt, C., Segschneider, J., and Tjiputra, J.: Twenty-First-Century Compatible CO2 Emissions and Airborne Fraction Simulated by CMIP5 Earth System Models under Four Representative Concentration Pathways, J. Climate, 26, 4398–4413, 2013. a, b, c, d, e
Keeling, C. and Bolin, B.: The simultaneous use of chemical tracers in oceanic studies 11. A three-reservoir model of the North and South Pacific Oceans, Tellus, 20, 17–53, 1968. a
Klaas, C. and Archer, D. E.: Association of sinking organic matter with various types of mineral ballast in the deep sea: Implications for the rain ratio, Global Biogeochem. Cy., 16, 1116, https://doi.org/10.1029/2001GB001765, 2002. a
Knox, F. and McElroy, M.: Changes in Atmospheric CO2: Influence of the Marine Biota at High Latitude, J. Geophys. Res., 89, 4269–4637, 1984. a
Kuhlbrodt, T. A., Griesel, M., Montoya, A., Levermann, M., Hofmann, M., and Rahmstorf, S.: On the driving processes of the Atlantic meridional overturning circulation, Rev. Geophys., 45, RG2001, https://doi.org/10.1029/2004RG000166, 2007. a
Lueker, T. J., Dickson, A. G., and Keeling, C. D.: Ocean pCO2 calculated from dissolved inorganic carbon, alkalinity, and equations for K-1 and K-2: validation based on laboratory measurements of CO2 in gas and seawater at equilibrium, Mar. Chem., 70, 105–119, 2000. a
Lumpkin, R. and Speer, K.: Global ocean meridional overturning, J. Phys. Oceanogr., 37, 550–562, 2007. a
Lund, D. C., Adkins, J. F., and Ferrari, R.: Abyssal Atlantic circulation during the Last Glacial Maximum: Constraining the ratio between transport and vertical mixing, Paleoceanography, 26, PA1213, https://doi.org/10.1029/2010PA001938, 2011. a
Lynch-Stieglitz, J., Stocker, T., Broecker, W., and Fairbanks, R.: The influence of air-sea exchange on the isotopic composition of oceanic carbon: Observations and modeling, Global Biogeochem. Cy., 9, 653–665, 1995. a
Marchitto, T., Lehman, S., Ortiz, J., Flückiger, J., and van Geen, A.: Marine Radiocarbon Evidence for the Mechanism of Deglacial Atmospheric CO2 Rise, Science, 316, 1456–1459, 2007. a
Marcott, S., Bauska, T., Buizert, C., Steig, E., Rosen, J., Cuffey, K., Fudge, T. J., Severinghaus, J. P., Ahn, J., Kalk, M. L., McConnell, J. R., Sowers, T., Taylor, K., White, J. W. C., and Brook, E.: Centennial-scale changes in the global carbon cycle during the last deglaciation, Nature, 514, 616–621, 2014. a
Martinez-Garcia, A., Sigman, D., Ren, H., Anderson, R., Straub, M., Hodell, D., Jaccard, S., Eglinton, T., and Haug, G.: Iron Fertilization of the Subantarctic Ocean During the Last Ice Age, Science, 343, 1347–1350, 2014. a, b
Meehl, G., Stocker, T., Collins, W., Friedlingstein, P., Gaye, A., Gregory, J., Kitoh, A., Knutti, R., Murphy, J., Noda, A., Raper, S., Watterson, I., Weaver, A., and Zhao, Z.: Global climate projections, Climate Change 2007: The Physical Science Basis, 747–845, Cambridge University Press, 2007. a, b
Mekik, F. A., Anderson, R. F., Loubere, P., François, R., and Richaud, M.: The mystery of the missing deglacial carbonate preservation maximum, Quaternary Sci. Rev., 39, 60–72, 2012. a
Menviel, L., Yu, J., Joos, F., Mouchet, A., Meissner, K. J., and England, M. H.: Poorly ventilated deep ocean at the Last Glacial Maximum inferred from carbon isotopes: A data-model comparison study, Paleoceanography, 31, 2–17, 2016. a, b, c, d, e
Milliman, J. D., Troy, P. J., Balch, W. M., Adams, A. K., Li, Y.-H., and Mackenzie, F. T.: Biologically mediated dissolution of calcium carbonate above the chemical lysocline?, Deep-Sea Res. Pt. I, 46, 1653–1669, 1999. a, b, c
Moore, J. K., Fu, W., Primeau, F., Britten, G., Lindsay, K., Long, M., Doney, S., Mahowald, N., Hoffman, F., and Randerson, J.: Sustained climate warming drives declining marine biological productivity, Science, 359, 1139–1143, https://doi.org/10.1126/science.aao6379, 2018. a, b
Morrison, A. and Hogg, A.: On the Relationship between Southern Ocean Overturning and ACC Transport, J. Phys. Oceanogr., 43, 140–148, 2013. a
Morse, J. W. and Berner, R. A.: Dissolution kinetics of calcium carbonate in sea water. II: A kinetic origin for the lysocline, Am. J. Sci., 272, 840–851, https://doi.org/10.2475/ajs.272.9.840, 1972. a, b, c, d
Mucci, A.: The solubility of calcite and aragonite in seawater at various salinities, temperatures, and one atmosphere total pressure, Am. J. Sci., 283, 780–799, 1983. a
Muglia, J., Skinner, L., and Schmittner, A.: Weak overturning circulation and high Southern Ocean nutrient utilization maximized glacial ocean carbon, Earth Planet. Sc. Lett., 496, 47–56, 2018. a, b, c, d, e, f
Nydal, R. and Lövseth, K.: Carbon-14 Measurements in Atmospheric CO2 from Northern and Southern Hemisphere Sites, 1962–1993, Tech. rep., Carbon Dioxide Information Analysis Center, Oak Ridge National Laboratory, Oak Ridge, Tennessee, 1996. a, b
Oliver, K. I. C., Hoogakker, B. A. A., Crowhurst, S., Henderson, G. M., Rickaby, R. E. M., Edwards, N. R., and Elderfield, H.: A synthesis of marine sediment core δ13C data over the last 150 000 years, Clim. Past, 6, 645–673, https://doi.org/10.5194/cp-6-645-2010, 2010. a, b
Olsen, A., Key, R. M., van Heuven, S., Lauvset, S. K., Velo, A., Lin, X., Schirnick, C., Kozyr, A., Tanhua, T., Hoppema, M., Jutterström, S., Steinfeldt, R., Jeansson, E., Ishii, M., Pérez, F. F., and Suzuki, T.: The Global Ocean Data Analysis Project version 2 (GLODAPv2) – an internally consistent data product for the world ocean, Earth Syst. Sci. Data, 8, 297–323, https://doi.org/10.5194/essd-8-297-2016, 2016. a, b, c
Opdyke, B. and Walker, J.: Return of the coral reef hypothesis: Basin to shelf partitioning of CaCO3 and its effect on atmospheric CO2, Geology, 20, 733–736, 1992. a
Raupach, M., Canadell, J., Ciais, P., Friedlingstein, P., Rayner, P., and Trudinger, C.: The relationship between peak warming and cumulative CO2 emissions, and its use to quantify vulnerabilities in the carbon-climate-human system, Tellus, 63, 145–164, 2011. a, b, c, d
Redfield, A. C., Ketchum, B. H., and Richards, F. A.: The influence of organisms on the composition of seawater, The Sea, 2, 26–77, 1963. a
Reimer, P., Baillie, M., Bard, E., Bayliss, A., Beck, J., Blackwell, P., Ramsey, C. B., Buck, C., Burr, G., Edwards, R., Friedrich, M., Grootes, P., Guilderson, T., Hajdas, I., Heaton, T., Hogg, A., Hughen, K., Kaiser, K., Kromer, B., McCormac, F., Manning, S., Reimer, R., Richards, D., Southon, J., Talamo, S., Turney, C., van der Plicht, J., and Weyhenmeyer, C.: IntCal09 and Marine09 radiocarbon age calibration curves, 0–50,000 years cal BP., Radiocarbon, 51, 1111–1150, 2009. a
Ridgewell, A., Watson, A., Maslin, M., and Kaplan, J.: Implications of coral reef buildup for the controls on atmospheric CO2 since the Last Glacial Maximum, Paleoceanography, 18, 1083, https://doi.org/10.1029/2003PA000893, 2003. a
Ronge, T., Tiedemann, R., Lamy, F., Kohler, P., Alloway, B., Pol-Holz, R. D., Pahnke, K., Southon, J., and Wacker, L.: Radiocarbon constraints on the extent and evolution of the South Pacific glacial carbon pool, Nat. Commun., 7, 11487, https://doi.org/10.1038/ncomms11487, 2016. a, b
Sabine, C., Feely, R., Gruber, N., Key, R., Lee, K., Bullister, J., Wanninkhof, R., Wong, C., Wallace, D., Tilbrook, B., Millero, F., Peng, T. H., Kozyr, A., Ono, T., and Rios, A.: The oceanic sink for anthropogenic CO2, Science, 305, 367–371, 2004. a, b, c
Sarmiento, J. L., Dunne, J., Gnanadesikan, A., Key, R. M., Matsumoto, K., and Slater, R.: A new estimate of the CaCO3 to organic carbon export ratio, Global Biogeochem. Cy., 16, 1107, https://doi.org/10.1029/2002GB001919, 2002. a, b
Schmitt, J., Schneider, R., Elsig, J., Leuenberger, D., Lourantou, A., Chappellaz, J., Köhler, P., Joos, F., Stocker, T., Leuenberger, M., and Fischer, H.: Carbon Isotope Constraints on the Deglacial CO2 Rise from Ice Cores, Science, 336, 711–714, 2012. a, b
Schmittner, A., Gruber, N., Mix, A. C., Key, R. M., Tagliabue, A., and Westberry, T. K.: Biology and air–sea gas exchange controls on the distribution of carbon isotope ratios (δ13C) in the ocean, Biogeosciences, 10, 5793–5816, https://doi.org/10.5194/bg-10-5793-2013, 2013. a
Schmitz, W. J.: On the World Ocean Circulation: Volume I. Some global features/North Atlantic circulation, Woods Hole Oceanographic Institution Technical Report, WHOI-96-03, 144, 1996. a
Siegenthaler, U. and Wenk, T.: Rapid atmospheric CO2 variations and ocean circulation, Nature, 308, 624–626, 1984. a
Sikes, E., Samson, C., Guilderson, T., and Howard, W.: Old radiocarbon ages in the southwest Pacific Ocean during the last glacial period and deglaciation, Nature, 405, 555–559, 2000. a
Sikes, E., Cook, M., and Guilderson, T.: Reduced deep ocean ventilation in the Southern Pacific Ocean during the last glaciation persisted into the deglaciation, Earth Planet. Sc. Lett., 438, 130–138, 2016. a
Skinner, L. and Shackleton, N. J.: Rapid transient changes in northeast Atlantic deep water ventilation age across Termination I, Paleoceanography, 19, PA2005, https://doi.org/10.1029/2003PA000983, 2004. a, b
Skinner, L., McCave, I., Carter, L., Fallon, S., Scrivner, A., and Primeau, F.: Reduced ventilation and enhanced magnitude of the deep Pacific carbon pool during the last glacial period, Earth Planet. Sc. Lett., 411, 45–52, 2015. a, b, c
Skinner, L., Primeau, F., Freeman, E., de la Fuente, M., Goodwin, P. A., Gottschalk, J., Huang, E., McCave, I. N., Noble, T. L., and Scrivner, A. E.: Radiocarbon constraints on the glacial ocean circulation and its impact on atmospheric CO2, Nat. Commun., 8, 16010, https://doi.org/10.1038/ncomms16010, 2017. a, b, c, d
Stommel, H.: Thermohaline convection with two stable regimes of flow, Tellus, 13, 224–230, 1961. a
Stuiver, M. and Polach, H.: Reporting of 14C data, Radiocarbon, 19, 355–363, 1977. a
Sverdrup, H., Johnson, N., and Fleming, R.: The Oceans, Prentice-Hall, Englewood Cliffs, NJ, 1941. a
Tagliabue, A., Bopp, L., Roche, D. M., Bouttes, N., Dutay, J.-C., Alkama, R., Kageyama, M., Michel, E., and Paillard, D.: Quantifying the roles of ocean circulation and biogeochemistry in governing ocean carbon-13 and atmospheric carbon dioxide at the last glacial maximum, Clim. Past, 5, 695–706, https://doi.org/10.5194/cp-5-695-2009, 2009. a
Talley, L.: Freshwater transport estimates and the global overturning circulation: Shallow, deep and throughflow components, Prog. Oceanogr., 78, 257–303, 2008. a
Talley, L.: Closure of the global overturning circulation through the Indian, Pacific, and Southern Oceans: Schematics and transports, Oceanography, 78, 257–303, 2013. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y, z, aa, ab, ac
Toggweiler, J. and Sarmiento, J.: Glacial to interglacial changes in atmospheric carbon dioxide: The critical role of ocean surface water in high latitudes, The Carbon Cycle and Atmospheric CO2: Natural Variations Archean to Present, Geophysical Monograph Series, American Geophysical Union, 32, 163–184, 1985. a, b, c, d, e, f, g
Toggweiler, J. R., Russell, J. L., and Carson, S. R.: Midlatitude westerlies, atmospheric CO2, and climate change during ice ages, Paleoceanography, 21, PA2005, https://doi.org/10.1029/2005PA001154, 2006. a
Trent-Staid, M. and Prell, W. L.: Sea surface temperature at the Last Glacial Maximum: A reconstruction using the modern analog technique, Paleoceanography, 17, 1065, https://doi.org/10.1029/2000PA000506, 2002. a, b, c
Tschumi, T., Joos, F., Gehlen, M., and Heinze, C.: Deep ocean ventilation, carbon isotopes, marine sedimentation and the deglacial CO2 rise, Clim. Past, 7, 771–800, https://doi.org/10.5194/cp-7-771-2011, 2011. a
Turnbull, J. C., Mikaloff Fletcher, S. E., Ansell, I., Brailsford, G. W., Moss, R. C., Norris, M. W., and Steinkamp, K.: Sixty years of radiocarbon dioxide measurements at Wellington, New Zealand: 1954–2014, Atmos. Chem. Phys., 17, 14771–14784, https://doi.org/10.5194/acp-17-14771-2017, 2017. a, b
Volk, T. and Hoffert, M. I.: Ocean carbon pumps: Analysis of relative strengths and efficiencies in ocean-driven atmospheric CO2 changes, in The Carbon Cycle and Atmospheric CO2: Natural Variations Archean to Present, Geophys. Monogr. Ser., 32, 99–110, 1985. a
Watson, A., Bakker, D. C. E., Ridgwell, A. J., Boyd, P. W., and Law, C.: Effect of iron supply on Southern Ocean CO2 uptake and implications for glacial atmospheric CO2, Nature, 407, 730–733, 2000. a
Watson, A., Vallis, G. K., and Nikurashin, M.: Southern Ocean Buoyancy Forcing of Ocean Ventilation and Glacial Atmospheric CO2, Nat. Geosci., 8, 861–864, 2015. a
Yokoyama, Y., Lambeck, K., Deckker, P. D., Johnston, P., and Field, K.: Timing of the Last Glacial Maximum from observed sea-level minima, Nature, 406, 713–716, 2000. a
Yu, J., Anderson, R., Jin, Z., Menviel, L., Zhang, F., Ryerson, F., and Rohling, E.: Deep South Atlantic carbonate chemistry and increased interocean deep water exchange during last deglaciation, Quaternary Sci. Rev., 90, 80–89, 2014a. a
Zhao, N., Marchal, O., Keigwin, L., Amrhein, D., and Gebbie, G.: A Synthesis of Deglacial Deep-Sea Radiocarbon Records and Their (In)Consistency With Modern Ocean Ventilation, Paleoceanography and Paleoclimatology, 33, 128–151, 2017. a, b