The [simple carbon project] model v1.0

. 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, simpliﬁed layout and matrix structure render it a ﬂexible 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 ﬂuxes and an ocean-basin-averaged topology, which may not be applicable to more detailed simulations. from the Last Glacial (LGM) the and the modern carbon inﬂuence of and ocean multi-proxy model–data parameter optimisation the late pool of published paleo atmosphere and ocean data for CO 2 , δ 13 C, (cid:49) 14 C 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 CO , δ C, (cid:49) and the oceanic distribution of δ 13 C, (cid:49) 14 C and the carbonate ion proxy. is to improve in these


Introduction
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  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 fourbox ocean-atmosphere models to show that the LGM CO 2 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 CO 2 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 mod-els 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 modeldata 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, dataaligned 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 CO 2 and δ 13 C (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 car-bon cycle under the influence of anthropogenic emissions. Emphasis is placed on the model description and configuration.

SCP-M description
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: δ 13 C, 14 C, plus alkalinity, phosphorus, oxygen, and atmospheric CO 2 , δ 13 C and 14 C. 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 δ 13 C 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.

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 Figure 1. SCP-M: configured as a seven-box ocean model plus atmosphere with marine sediments, continents and the terrestrial biosphere. The exchange of elemental concentrations, e.g. C i (i = 1, 7), occurs due to fluxes between boxes. 1 (red arrows) is global overturning circulation (GOC), 2 (orange arrows) is Atlantic meridional overturning circulation (AMOC), γ 1 (blue arrows between boxes 4 and 6) is deep-abyssal mixing, γ 2 (blue arrows between boxes 1 and 3) is low-latitude thermohaline mixing, Z (green downward arrows) is the biological pump, F CA (white downward arrows) is the carbonate pump, D CA (white squiggles) is carbonate dissolution and P (black, bidirectional arrows) is the air-sea gas exchange. Box 1: low-latitude-tropical surface ocean; box 2: northern surface ocean; box 3: intermediate ocean; box 4: deep ocean; box 5: Southern Ocean; box 6: abyssal ocean; box 7: subpolar southern surface ocean. 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 -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, 10 6 m 3 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 winddriven 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). 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).

Basic features
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.

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) andDe Boer andHogg (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 C i is the concentration of DIC in each box in mol m −3 and V i is the volume of each box in m 3 . 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 − 1 C 4 , − 2 C 4 and −γ 1 C 4 , 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. where and T 1 , T 2 , E 1 and E 2 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.

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 shallowwater-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 pCO 2 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 , F 100 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 Z 1 is the biological flux of carbon prescribed for surface box 1 in mol C m −2 yr −1 , S 1 is the surface area of surface box 1, d 0 is the reference depth of 100 m for the Z parameter value (Martin et al., 1987), and d c and d f 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 Z i (i = 1, 7) parameter, which varies across the surface boxes, and S is the array of surface box surface areas S i (i = 1, 7). 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 . B out and B in 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) (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).

pCO 2 and carbonate
Modelling air-sea gas exchange, atmospheric pCO 2 and the "carbonate pump" rests on a realistic estimation of pCO 2 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 pCO 2 . 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). pCO 2 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 pCO 2 at each time step of a model Table 1. Values for the Z biological production parameter (at 100 m of ocean depth) used in the SCP-M model calibration. A global value for Z was tuned in each surface box using scalars (column 3) to yield unique values for each surface box (column 4).

Model surface box
Global value at 100 m of ocean Scalar Model input (tuned) depth (mol C m −2 yr −1 ) (tuned) mol C m −2 yr −1 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 pCO 2 . Solving for pCO 2 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 CO 2 on 1000-year time frames (e.g. Farrell and Prell, 1989;Anderson et al., 2007;Mekik et al., 2012;Yu et al., 2014b).

The carbonate pump
According to Emerson and Hedges (2003), ∼ 20 %-30 % of CaCO 3 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 CaCO 3 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 shellbased "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-saturationdependent 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 CaCO 3 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 sed-iments 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 F CA 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 CaCO 3 is the concentration of calcium carbonate in each box. The dissolution equation of Morse and Berner (1972) operates on CaCO 3 , which is calculated by multiplying Ca by CO 2− 3 , where Ca is estimated from salinity in each box as per Sarmiento and Gruber (2006).

Air-sea gas exchange
CO 2 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 K H 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 CO 2 as a function of the pCO 2 differential between ocean and atmosphere, a CO 2 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 CO 2 between a surface box and the atmosphere: where P 1 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). K 0 1 is the solubility of CO 2 in mol kg −1 atm −1 (Weiss, 1974), subsequently converted into mol m −3 by multiplying by seawater density ρ. pCO 2 1 and pCO 2 at are the partial pressures of CO 2 in surface ocean box 1 and the atmosphere, respectively, in ppm. The equation is vectorised as follows: where P = P i (i = 1, 7), with zero values for non-surface boxes, and K 0 = K 0 i (i = 1, 7).

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 pCO 2 , which further enables the calculation of air-sea gas exchange and the species of DIC in seawater (H 2 CO 3 , HCO 3− and CO 2− 3 ).

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 CO 2 in SCP-M. Toggweiler (2008) estimated this volcanic flux of CO 2 emissions at 4.98 × 10 12 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 pCO 2 and carbonate ion levels, which impacts atmospheric CO 2 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 W SC is a constant silicate weathering term set at 0.75 × 10 −4 mol m −3 yr −1 , W SV is a variable rate of silicate weathering per unit of atmosphere CO 2 (ppm) set to 0.5 mol m −3 atm −1 CO 2 yr −1 and W CV is the variable rate of carbonate weathering with respect to atmosphere CO 2 set at 2 mol m −3 atm −1 CO 2 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 CO 2 (e.g. Toggweiler, 2008;Hogg, 2008). The atmospheric sink of CO 2 is calculated by multiplying Eq. (20) by the volume of the low-latitude surface ocean box (box 1) and subtracting from atmospheric CO 2 .
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 CO 2 during periods of biosphere growth (e.g. post-glacial regrowth), via carbon fertilisation, or a source of CO 2 (e.g. glacial reduction) via respiration. We employ a two-part model of the terrestrial biosphere with a long-term (woody forest) and shortterm (grassland) terrestrial biosphere box as described by Raupach et al. (2011) andHarman 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 CO 2 , 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 12 C, leading to a relative enrichment in δ 13 C in the atmosphere during net uptake of CO 2 . The change in atmospheric CO 2 from the terrestrial biosphere in the model is given by where N pre is NPP at a reference level ("pre") of atmospheric CO 2 , 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) . β is the parameterisation of carbon fertilisation, which causes NPP to increase (decrease) logarithmically with rising (falling) atmospheric CO 2 levels, and has a typical value of 0.4-0.8 (Harman et al., 2011). C stock1 is the short-term carbon stock and k 1 is the time frame for respiration in the short-term carbon stock (in years). For the longterm terrestrial biosphere, we substitute (1 − RP) in place of RP and C stock2 and k 2 for the long-term carbon stock and respiration rate, respectively. D forest is a prescribed flux of deforestation emissions, which can be switched on or off in SCP-M. A δ 13 C fractionation factor is applied to the terrestrial biosphere fluxes of carbon, effecting an increase in atmospheric δ 13 C from biosphere growth and a decrease from respiration.

The complete carbon cycle equations
Equation (22)  consists of a prescribed flux of δ 13 C-depleted and 14 C-dead CO 2 to the atmosphere from human industrial emissions, which is activated by a model switch in SCP-M.

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 (δ 13 C and 14 C), 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 δ 13 C and 14 C into this metric for the operation of model fluxes, the method of Craig (1969) is applied to convert starting data values of δ 13 C and 14 C from delta notation in ‰ into mol m −3 : where 13 C i is the 13 C concentration in box i in mol m −3 , δ 13 C i is δ 13 C in ‰ in box i, R is the 13 C 12 C ratio of the standard (0.0112372 as per the Pee Dee Belemnite value) and C i is the DIC concentration C in box i in mol m −3 .
The calculation in Eq. (24) derives the fraction 13 C 12 C 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 (C i ) in each box. This approach rests on the assumption that the fraction 13 C 12 C is the same as 13 C total carbon . For example, there are three isotopes of carbon, each with different atomic weights. They occur in roughly the following abundances: 12 C ∼ 98.89 %, 13 C ∼ 1.11 % and 14 C ∼ 1 × 10 −10 %. Therefore, the assumption that 13 C 12 C = 13 C total carbon is a valid approximation. Once converted from δ 13 C (‰) to 13 C in mol m −3 , SCP-M's ocean parameters can operate on 13 C concentrations in each box, according to the same model flux equations used for DIC and CO 2 . The 13 C model results are then converted back into δ 13 C notation at the end of the model run in order to compare the model output with the data, which is reported in δ 13 C format. The same method is applied to 14 C. The reference standard value for 14 C 12 C is 1.2×10 −12 (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.

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 12 C (the lighter carbon isotope). This activity enriches the surface ocean in 13 C and relatively enriches the underlying waters in 12 C when remineralisation occurs. As such, the ocean displays depletion in δ 13 C 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 13 C and to replicate ocean distributions of δ 13 C: where f is the biological fractionation factor for stable carbon (e.g. ∼ 0.977 in Toggweiler and Sarmiento, 1985), and S st is the ratio of 13 C to 12 C in the reference standard. The typical δ 13 C composition of marine organisms is in the range −23 to −30‰. The same method is applied for biological fractionation of 14 C, but with a different fractionation factor (Toggweiler and Sarmiento, 1985).

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, 12 C, preferentially partitions into the atmosphere as a net effect of bidirectional gaseous exchange. This fractionation leads to the heavily depleted δ 13 C 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 13 C and 14 C across the air-sea interface compared with 12 C (Zhang et al., 1995). R At is the ratio of 13 C to 12 C in the atmosphere, and R i is the ratio of 13 C to 12 C in surface ocean box i. pCO 2At is the atmospheric pCO 2 , and pCO 2i is the pCO 2 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).

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 (∼ 6.022×10 23 ). The resultant figure is multiplied by the Earth's surface area (∼ 5.1×10 18 cm −2 ) to yield a production rate of 1.3296 × 10 −5 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.

Modelling results
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.

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 CO 2 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  (1996) data accumulated over the period 1972-2013. We make adjustments to the ocean concentrations of DIC, δ 13 C and 14 C 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 δ 13 C and CFC-12 in the ocean, which we applied using GLODAPv2 CFC-12 data to correct the ocean δ 13 C data. We force the model with late Holocene average data for atmosphere CO 2 , δ 13 C and 14 C (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 spinup. 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 CO 2 , δ 13 C and 14 C. Model results are also in good agreement with the late Holocene ocean 14 C, falling within error or very close for all boxes covered by data. The surface boxes (1, 2) are relatively enriched in 14 C relative to deeper boxes, reflecting their proximity to the atmospheric source of 14 C, although the spread of values across the ocean boxes is narrow. The surface boxes (1, 2 and 7) intuitively display more enriched δ 13 C 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 δ 13 C 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 δ 13 C core data for the polar Southern Ocean (one data point, no error bars). SCP-M also exaggerates the depletion in δ 13 C in the deep box (4) relative to the data observation.
There is limited data coverage for the carbonate ion proxy (CO 2− 3 ), although the model replicates the available data well. CO 2− 3 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. CO 2− 3 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. CO 2− 3 is less abundant in the deep ocean (boxes 4 and 6) because there is  (Table A1). (a) Model results for atmospheric δ 13 C, 14 C and CO 2 (red stars) plotted against late Holocene average data values (blue dots) with standard error bars. (b) The model results for oceanic δ 13 C, 14 C and the carbonate ion proxy (red stars) plotted against late Holocene average ocean data where available (blue dots). Data sources are shown in Table 2. more carbon relative to alkalinity due to the remineralisation of organic matter, a pattern that SCP-M replicates.

Sensitivity tests
We undertook parameter sensitivity tests to understand changes in atmospheric CO 2 , 14 C and δ 13 C 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 CO 2 , it would make sense to probe parameter values lower than modern to replicate the lower atmospheric CO 2 in the LGM. We varied parameter values around their modern-day settings in 10 kyr model runs and plotted the output against atmospheric CO 2 , 14 C and δ 13 C (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 CO 2 is very sensitive to 1 and 2 but displays limited response to γ 1 and γ 2 over the ranges analysed ( Fig. 4a-d). Atmospheric 14 C and δ 13 C are negatively related to 1 and 2 . The slower ocean turnover leads to a reduced rate of upwelling and surface de-gassing of 14 C-and δ 13 C-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 CO 2 (Fig. 4e). A higher global value of Z drives the removal of carbon from the surface ocean, and the resulting CO 2 flux into the ocean from the atmosphere decreases CO 2 . Lower Z leads to increased atmospheric CO 2 . δ 13 C 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 pCO 2 in the surface ocean boxes, leading to de-gassing of CO 2 to the atmosphere, and therefore modestly decreasing atmospheric 14 C, as the lighter 12 C is preferentially partitioned across the air-sea interface. . Univariate parameter sensitivity tests around modern-day estimated values plotted for atmospheric CO 2 , 14 C and δ 13 C versus (a) 1 , (b) 2 , (c) γ 1 , (d) γ 2 , (e) biological pump, (f) rain ratio, (g) subpolar surface box temperature, (h) subpolar surface box salinity, (i) polar box piston velocity, (j) net primary productivity, (k) ocean volume and (l) fraction of deep water upwelled into the subpolar surface box. We varied parameter input values as plotted on the x axes and show model output for atmospheric CO 2 , 14 C and δ 13 C. Atmospheric CO 2 shows the greatest sensitivity to parameters associated with ocean circulation, biology and the terrestrial biosphere. Other parameters exert less influence on atmospheric CO 2 but are important for atmospheric carbon isotope values. Modern-day estimates used in SCP-M are shown with vertical black dotted lines in each subplot (sources in the text and Appendix Table A1).
Increasing surface ocean box temperature (Fig. 4g) increases atmospheric CO 2 , an intuitive outcome given that warmer water absorbs less CO 2 (Weiss, 1974), and SCP-M employs a temperature-and salinity-dependent CO 2 solubility function. Air-sea fractionation of δ 13 C also decreases with higher temperatures, leading to higher atmospheric δ 13 C. According to Mook et al. (1974), air-to-sea fractionation of δ 13 C (making the atmosphere more depleted in δ 13 C) increases at a rate of approximately 0.1 ‰ • C −1 of cooling. SCP-M employs temperature-dependent air-sea gas δ 13 C fractionation factors (Mook et al., 1974). 14 C is invariant to box temperature as the fractionation parameters employed in the model are not temperature dependent. CO 2 displays a weak positive relationship with surface ocean box salinity (Fig. 4h) due to the decreasing solubility of CO 2 in ocean water with increasing salinity (Weiss, 1974).
As the polar box piston velocity P slows down (Fig. 4i), atmospheric CO 2 falls. At lower values of P the polar box, a region of outgassing of CO 2 due to the upwelling of deepsourced carbon-rich water in that part of the ocean, will exchange CO 2 with the atmosphere at a slower rate. The reduced outgassing of δ 13 C-depleted carbon to the atmosphere with a lower P leads to higher δ 13 C values in the atmosphere. Atmospheric 14 C 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 14 C waters.
Terrestrial biosphere NPP (Fig. 4j) is a sink of CO 2 and fractionates the ratios of the isotopes of carbon, leading to higher values for δ 13 C and, to a lesser extent, 14 C in the atmosphere. It is likely that NPP plays a feedback role and modulates CO 2 , δ 13 C and 14 C (Toggweiler, 2008). Varying the ocean surface area (Fig. 4k) has modest impacts on CO 2 and δ 13 C, but a large impact on 14 C. Decreasing the ocean volume leads to a lower surface area for CO 2 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 pCO 2 on glacial-interglacial timescales. Increasing the fraction of deep water upwelled into the subpolar box (Fig. 4l) intuitively raises CO 2 , but lowers δ 13 C and 14 C, by upwelling carbon-rich and isotopically depleted water to the ocean surface boxes.

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 CO 2 emissions, sea surface temperature (SST) changes and atmospheric bomb 14 C 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) CO 2 emissions and SST scenarios (Boden et al., 2017;Houghton, 2010;IPCC, 2013a). We compare the model results with atmospheric CO 2 , δ 13 C and 14 C 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 CO 2 , δ 13 C and 14 C 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 CO 2 , δ 13 C and 14 C preserved in data observations for the period 1751-2016 (a-b). The atmospheric CO 2 and δ 13 C data are sourced from the Scripps CO 2 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 14 C data are sourced from Nydal and Lövseth (1996), Stuiver et al. (1998), andTurnbull 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 CO 2 and a steep drop in δ 13 C (Fig. 5a); this reflects δ 13 C-depleted anthropogenic emissions. The effect of emissions on atmospheric 14 C (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 14 C in the immediate lead-up to the release of bomb radiocarbon into the atmosphere and the downward trend from ∼ 2020 onwards. The spike in 14 C during the period of bomb radiocarbon release lasts during the period 1954-1963 and then disperses as 14 C 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 CO 2 (b) over historical time and projected forward to 2100 for the IPCC RCPs. The SCP-M output falls below the IPCC projections for atmospheric CO 2 in RCP2.6 and 4.5, but provides a close match with RCP6.0 and 8.5. Figure 7a shows the annual uptake of CO 2 by the ocean, modelled with SCP-M. The model begins the period close to a steady state between the atmosphere and surface ocean pCO 2 , with limited transfer across the interface. Beginning circa 1950 the ocean takes up an increased load of CO 2 from the atmosphere. By 2100, SCP-M models a range of annual CO 2 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)   In (a-b) we plot SCP-M model results for CO 2 , δ 13 C and 14 C (lines) for the period 1751-2100 against atmospheric data for CO 2 , δ 13 C and 14 C (red dots). The SCP-M model output closely resembles the atmospheric data record. The perturbation from industrial-era, isotopically depleted (δ 13 C) and dead ( 14 C) CO 2 is clear, as is the impact of atmospheric nuclear tests on 14 C during 1954-1963. In the other rows we plot SCP-M model results (boxes as shown) versus GLODAPv2 data (dots and error bars, same colour as corresponding boxes). We assume an average data year of 1990 for the GLODAPv2 data accumulated over the period 1972-2013. For most of the SCP-M ocean boxes, the model results fall within or very close to error ranges of the GLODAPv2 data, despite large perturbations in the model and data from industrial-era emissions and bomb radiocarbon. 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 CO 2 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 CO 2 emissions in the atmosphere and other reservoirs (Fig. 9). Anthropogenic CO 2 emissions transfer carbon to the atmosphere, ocean and terrestrial biosphere. The fluxes between the carbon reservoirs also change. In the pre-industrial state, CO 2 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 CO 2 concentration increases to the extent that the atmosphere-ocean pCO 2 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 respirationboth inward and outward biosphere fluxes of CO 2 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 CO 2 , 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 pCO 2 -a modest negative feedback. In summary, SCP-M provides an appropriate simulation of historical atmospheric CO 2 , δ 13 C and 14 C data when forced with anthropogenic CO 2 emissions data over the same period. For the forwardlooking 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 CO 2 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.

Background
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 CO 2 across glacial and interglacial periods. These CO 2 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 CO 2 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 CO 2 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 . Panel (a) shows the annual uptake of CO 2 by the ocean in each of the RCPs over the period 1751-2100, modelled with SCP-M. By 2100, SCP-M estimates a range of 0-6 PgC yr −1 across the RCPs as estimated by CMIP5 models, reproduced from Jones et al. (2013). Panel (b) shows the cumulative uptake of CO 2 by the ocean over the same period modelled with SCP-M and compared with CMIP5 models (Jones et al., 2013).  . SCP-M-modelled pre-industrial carbon stocks and fluxes (in PgC in black text) compared with IPCC RCP6.0 emissions scenario by 2100 (shown as PgC changes with blue text for positive changes, red text for negative and black text for no change). Atmosphere, ocean and terrestrial biospheres take up the load of carbon from the industrial source. By 2100, carbon is fluxing into all ocean boxes, the terrestrial biosphere and continental sediment weathering-river fluxes. Pre-industrial outgassing of CO 2 in the Southern Ocean is reversed, and carbon is returned to the ocean via enhanced CaCO 3 dissolution. Box numbers on the diagram refer to ocean regions specified in Fig. 1. Negative fluxes on bidirectional air-sea exchange arrows are fluxes of CO 2 out of the atmosphere into the ocean. weiler, 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 CO 2 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 δ 13 C 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 δ 13 C, 14 C and 15 N 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.

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 14 C 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 CO 2 drawdown. Finally, Mariotti et al. (2013) estimated that higher atmospheric radiocarbon production accounted for + ∼ 200‰ in atmospheric 14 C 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 14 C 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 CO 2 -δ 13 C-14 C 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 CO 2 , δ 13 C and 14 C 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 CO 2 of ∼ 20 ppm and a lightening of δ 13 C by ∼ 0.6 ‰ owing to the increased solubility of CO 2 in colder water and the increasing fractionation of δ 13 C with decreasing temperatures, which leaves more 12 C in the atmosphere. There is limited impact on 14 C. Increasing salinity slightly reverses these changes to CO 2 and δ 13 C. Reducing sea surface area and volume slightly increases CO 2 and increases 14 C 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 CO 2 (reduced outgassing), increases 14 C (slower rate of invasion to the ocean) and increases δ 13 C (reduced outgassing and sea-to-air fractionation of δ 13 C). Finally, increasing the rate of atmospheric radiocarbon production forces a shift in 14 C (horizontal shift in Fig. 10) towards LGM levels (black star and circle in Fig. 10). In aggregate, these changes lead to a fall in CO 2 of ∼ 35 ppm, a fall in δ 13 C of ∼ −0.5 ‰ and an increase in 14 C 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 CO 2 changes, but steers δ 13 C and 14 C 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 CO 2 , δ 13 C and 14 C to variations in 1 , 2 and Z lead us to cater for lower values Table 3. Changes to ocean and atmosphere parameter settings in SCP-M to recreate the LGM background model state.

Indicator
LGM change Surface ocean box temperatures −5-6 • C (Trent-Staid and Prell, 2002;Annan and Hargreaves, 2013) Surface ocean box salinity +1.0 psu (Adkins et al., 2002) Polar ocean box piston velocity ×0.3 (Stephens and Keeling, 2000;Ferrari et al., 2014) Ocean surface area and volume −3.0 % (Adkins et al., 2002;Grant et al., 2014) Atmosphere radiocarbon production ×1.25 (Mariotti et al., 2013) As shown in the sensitivity tests in Fig. 4, some processes do not exert a strong influence on atmospheric CO 2 , but do modestly impact CO 2 and strongly impact δ 13 C and 14 C. 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. Figure 10. LGM state parameter adjustments. Using the posited LGM changes in environmental parameters in Table 3, we establish LGM foundations for exploring the impacts of varying large-scale ocean process parameters towards LGM atmospheric CO 2 -δ 13 C-14 C data space. The red circle is our starting point for the late Holocene. From the LGM state foundation (black star), variation of global overturning circulation ( 1 ), Atlantic meridional overturning circulation ( 2 ) and the soft tissue biological pump (Z) drives atmospheric CO 2 , δ 13 C and 14 C into the vicinity of their LGM data values (black circle). The biological pump Z can effect the LGM CO 2 outcome, but steers δ 13 C away from the LGM value. 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. 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.
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 CO 2 , atmospheric, and ocean 14 C and δ 13 C, 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 Opt n=1 represents the optimal value of parameters n, R i,k is the model output for concentration of each element i in box k, D i,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 CO 2 , δ 13 C and 14 C 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 δ 13 C data is a depletion of deep ocean δ 13 C in the LGM, shown as a drop in δ 13 C values in the deep (box 4) and abyssal (box 6) boxes relative to the intermediate box (3). In the LGM δ 13 C 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 abyssaldeep ocean water flows via 1 in delivering the ocean δ 13 C data patterns. The model shift in δ 13 C 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 δ 13 C value remains largely unchanged between the two periods due to the effect of the terrestrial biosphere, which causes net uptake of CO 2 in the Holocene period (increases atmospheric δ 13 C) and net respiration of CO 2 in the LGM period (decreases atmospheric δ 13 C).
The model results also closely replicate the reduction in deep-to-shallow ocean compositional gradient in 14 C data moving from the LGM to Holocene period (e.g. Skinner and Shackleton, 2004;Skinner et al., 2015Skinner et al., , 2010Burke 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 14 C-depleted deep water into intermediate and shallow depths in the Holocene (e.g. Skinner et al., 2010Skinner et al., , 2015Burke and Robinson, 2012). A slowdown in Southern Ocean upwelling in the LGM allows 14 C-depleted water to accumulate in the deep or abyssal ocean and a widening in the 14 C gradient between deep and shallow waters. In SCP-M, this is simulated by lower values for 1 and 2 . The lowlatitude surface box (box 1) enrichment in 14 C 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 CO 2 , 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 CO 2 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, at the top end 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.
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 δ 13 C 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 δ 13 C by the terrestrial biosphere during the deglacial period effectively reverses the effects of the release of δ 13 C-depleted carbon from the deep ocean to the atmosphere at the termination and leaves atmosphere δ 13 C 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).

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 largescale 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).

Modern carbon cycle simulations
Our simple forcing of SCP-M with historical and projected anthropogenic CO 2 emissions and SST (Sect. 3.3, shows that SCP-M can reproduce historical data and the model results of more complex CMIP5 models. The SCP-M results for atmospheric CO 2 (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 CO 2 , δ 13 C and 14 C (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 14 C 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 14 C during LGM parameter values selected from the four-parameter LGM experiment in Table 5. The LGM setting leads to a transfer of carbon from the atmosphere and terrestrial biosphere to the deep ocean. Carbon upwelled into the surface ocean falls, leading to reduced outgassing of CO 2 in the Southern Ocean boxes. Continental weathering and river fluxes of carbon are also reduced due to lower atmospheric CO 2 , leading to a change in CaCO 3 burial and dissolution in marine sediments until equilibrium is restored with river input to the oceans. Box numbers on the diagram refer to ocean regions specified in Fig. 1. [1954][1955][1956][1957][1958][1959][1960][1961][1962][1963], and SCP-M appropriately replicates the take-up of bomb 14 C 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 windshift-induced slowing of AMOC and thermocline mixing or a response of ocean biological productivity to changed pCO 2 , 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 fer-tilisation, afforestation and/or reforestation on land, or marine fauna management, as tools for reducing atmospheric CO 2 on a global scale are feasible uses of SCP-M.

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 CO 2 , δ 13 C, 14 C 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 14 C 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 wellunderstood AMOC and the more detailed proxy data coverage in that basin. For example, Curry and Oppo (2005) provided striking transects of δ 13 C in the LGM and late Holocene Atlantic Ocean, which provided evidence for large changes in the basin δ 13 C stratigraphy across the two time periods. The δ 13 C 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 CO 2 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 CO 2 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 CO 2 to the atmosphere, and rebounded from the LGM to the Holocene, becoming a sink of CO 2 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 CO 2 and, critically, δ 13 C 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 CO 2 drawdown. The Southern Ocean ma-rine biology, in particular, is posited as a candidate for driving glacial-interglacial cycles of CO 2 (e.g. Martin, 1990;Martinez-Garcia et al., 2014).

Conclusions
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 oceanatmosphere-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 14 C. A model-data experiment using LGM and late Holocene CO 2 , δ 13 C, 14 C 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. Appendix A: Parameters, data sources and dimensions  (2001) Air-sea exchange velocity (m day −1 ) 3.0 Toggweiler (1999) 13 C air-sea fractionation factors 0.9989-0.999 Mook et al. (1974) 14 C air-sea fractionation factors 0.98-0.998 Toggweiler and Sarmiento (1985) 13 C "thermodynamic" air-sea factor 0.99915 Schmittner et al. (2013) 14 C "thermodynamic" air-sea factor 0.999 Toggweiler and Sarmiento (1985) Organic δ 13 C fractionation factor ∼ 0.975 Toggweiler and Sarmiento (1985) C / P in org. Redfield ratio 130 Takahashi et al. (1985) Rain ratio (carbonate : org. in sinking particles) 0.07 Sarmiento et al. (2002) CaCO 3 dissolution rate (units day −1 ) 0.38 Hales and Emerson (1997) n order of CaCO 3 dissolution reaction rate 1.0 Hales and Emerson (1997) K sp solubility coefficient for calcite Various Mucci (1983) Carbon chemistry solubility and dissociation coefficients Various Weiss (1974), Lueker et al. (2000) Atmosphere radiocarbon production rate (atoms s −1 ) ∼ 1.6 Key (2001) Suess and bomb radiocarbon corrections Various Broecker et al. (1980), Key (2001), Sabine et al. (2004), Eide et al. (2017) Radiocarbon decay rate (yr −1 ) 1/8267 Stuiver and Polach (1977) Volcanic emissions flux CO 2 (mol C yr −1 ) 5-6 × 10 12 Modified from Toggweiler (2008) River phosphorus flux (Tg yr −1 ) 15.0 Compton et al. (2000)