A Computer Model

The use of models for the study and prediction of MRA extends back to Craig (Reference Craig1957), Arnold (Reference Arnold1957), Revelle and Suess (Reference Revelle and Suess1957) and was updated later by Oeschger et al. (Reference Oeschger, Siegenthaler, Schotterer and Gugelmann1975) with the introduction of a diffusion coefficient to simulate mixing in the deep ocean box (eddy diffusivity). This so-called box-diffusion model has been fundamental to the construction of the marine radiocarbon age calibration curve since 1986 (Stuiver et al. Reference Stuiver, Pearson and Braziunas1986). Over that time, the complexity of these models has increased as we have gained more knowledge of both the Earth systems and past climate constraints. Recent approaches include more advanced carbon cycle box models which allow transient forward simulations such as BICYCLE (Köhler et al. Reference Köhler, Muscheler and Fischer2006) through to full three-dimensional ocean general circulation models such as the LSG OGCM (Butzin et al. Reference Butzin, Prange and Lohmann2005, Reference Butzin, Prange and Lohmann2012, Reference Butzin, Köhler and Lohmann2017, Reference Butzin, Köhler, Heaton and Lohmann2020 in this issue). For all of these models, the basic principle is the same: given a specific user-defined set of inputs/parameters, they provide an estimate/prediction of the variable of interest, in our case the non-polar global-average marine Δ¹⁴C, where this average in the following is restricted to the area between 40°S and either 50°N (Atlantic) or 40°N (Pacific) to avoid the potential complications of including regions covered with sea ice (e.g. Butzin et al. Reference Butzin, Köhler and Lohmann2017).

The most useful computer model for a particular situation will depend upon the specific question one wishes to answer. For Marine20, we wish to not only obtain a reliable estimate of the global-average marine ¹⁴C concentration over time incorporating known changes in both the carbon cycle and atmospheric ¹⁴C production rate, but also to be able to quantify our uncertainty on that estimate. Such quantification of uncertainty is fundamental to radiocarbon calibration. When calibrating a radiocarbon date, we need not only a best “point estimate” for its calendar age but also a plausible range, i.e. to understand the uncertainty on our calibrated calendar age. These two competing aims of reliability and uncertainty quantification typically require a trade-off in model selection. We require a computer model that can incorporate desired physical processes and changes to the carbon cycle that proxy data tell us occurred; and yet also allows us to understand and propagate model uncertainty, specifically in the values of the input parameters, through to the final marine Δ¹⁴C estimate. In the case of our chosen BICYCLE computer model, the various inputs are illustrated in Figure 2.

Figure 2 Propagating uncertainty through the BICYCLE model. For each simulation, we generate different inputs for BICYCLE by drawing from the variables on which we consider uncertainty (shown in the yellow box). AMOC denotes the Atlantic meridional overturning circulation and piston velocity the rate of air-sea gas exchange, see Section 4 for further details. These random inputs are combined with the inputs for which no uncertainties are incorporated (shown in green) and entered into the BICYCLE model. Each simulation of BICYCLE therefore provides a different, deterministic historic global-average estimate of marine surface Δ¹⁴C from 0–55 cal kBP. From the resulting ensemble we generate a Monte Carlo estimate of the mean and variance at any calendar age which is then used for the Marine20 radiocarbon age calibration curve.

A Brief Introduction to Quantifying Model Uncertainty

The majority of computer models simulating changes in the carbon cycle are deterministic. This means that they do not internally introduce randomness—if one drives the model with the same inputs one will always get identical outputs. Let us assume that such a computer model takes as inputs a range/vector of parameters φ which we specify. In our context, these model parameters are any input needed by the model—they could be a forcing, a user chosen setting or a tuned variable. In the case of marine ¹⁴C reservoir age simulations, these input parameters typically include the historic levels of atmospheric Δ ¹⁴C and CO₂ concentration; ocean circulation; and air-sea gas exchange amongst others. Having specified values for these inputs, the model prediction at any calendar age θ is then represented by the deterministic function f(θ; φ).

In the case of a deterministic computer model, output uncertainty arises from two sources (see Section 2 of Kennedy and O’Hagan Reference Kennedy and O’Hagan2001). The first is model parameter uncertainty, i.e. that, typically, we do not know the true value of the model’s input parameters. We specify our uncertainty on our computer model inputs in terms of a distribution $\phi \sim \pi \left( \phi \right)$ that encapsulates our current understanding about their values. The second source of uncertainty is model discrepancy, i.e. that any computer model is not an entirely accurate representation of reality and will not capture the complete carbon cycle of the Earth. This type of uncertainty is much harder to evaluate but one might expect a more complex and detailed model to typically have a smaller model discrepancy. We do not attempt to formally incorporate model discrepancy into Marine20 but we are able to incorporate model parameter uncertainty.

For calibration purposes, we wish to estimate, for any calendar age θ both the expectation ${\mu _{Mar}}\left( \theta \right)$ and the variance $\sigma _{Mar}^2\left( \theta \right)$ in the global-average level of marine $\Delta _{MarineAv}^{14}C\left( \theta \right)$. To propagate our model parameter uncertainty, we therefore need to evaluate both:

$${\mu _{Mar}}\left( \theta \right) = E\left[ {\Delta _{MarineAv}^{14}C\left( \theta \right)} \right] = \int f\left( {\theta ;\phi } \right)\pi \left( \phi \right)d\phi ,$$

and

$$\sigma _{Mar}^2\left( \theta \right) = Var\left[ {\Delta _{MarineAv}^{14}C\left( \theta \right)} \right] = {\rm{ }}\left\{ {\int {{f^2}} \left( {\theta ;\,\phi } \right){\rm{ }}\pi \left( \phi \right)d\phi } \right\}{\rm{ }} - \mu _{Mar}^2\left( \theta \right),$$

Where $\pi \left( \phi \right)$ is the density of the model’s input parameters that specifies our prior beliefs/uncertainties about their true values, and f(θ; φ) is the model output at calendar year θ cal kBP when it has been run with specific inputs φ.

Monte-Carlo Estimation

Due the complexity of the computer model, it is not possible to calculate either of these integrals explicitly, but we can estimate them via a Monte-Carlo approach. To do so, we sample a range of N plausible inputs φ ₁, …, φ_N from our density $\pi \left( \phi \right)$. The model is then run with each of these inputs giving us a large ensemble of N outputs $f\left( {\theta ;{\phi _1}} \right), \ldots ,f\left( {\theta ;{\phi _N}} \right)$. The sample mean and variance of this ensemble then provide estimates for ${\mu _{Mar}}\left( \theta \right)$ and $\sigma _{Mar}^2\left( \theta \right)$ respectively, i.e. we estimate

$${\hat \mu _{Mar}}\left( \theta \right) = {1 \over N}\mathop \sum \limits_{i=1}^N f\left( {\theta ;{\phi _i}} \right)\;{\rm{and}}$$

$$\hat \sigma _{Mar}^2\left( \theta \right) = {\rm{}}{1 \over {N - 1}}\mathop \sum \limits_{i = 1}^N {\left( \,{f\left( {\theta ;{\phi _i}} \right) - {{\hat \mu }_{Mar}}\left( \theta \right)} \right)^2},\;{\rm{where}}\;{\rm{each}}\,{{\phi}_i}\,\,\sim\,\,\,\,\pi (\phi ).$$

Importantly, as most computer models are highly non-linear, we stress that only running them at extreme “end-member” inputs is not sufficient to reliably represent output uncertainty. Equally, running a computer model with the value of an input parameter set at its mean is not sufficient to estimate the mean of the output when input parameter uncertainty is properly incorporated. Both these points are illustrated with a simple example in the next subsection.

To perform accurate Monte-Carlo estimation, it is necessary to generate a sufficiently large ensemble N of outputs in order that they provide reliable estimates. Here, we selected N=500 based on a smaller pilot study which indicated such a sample size would provide Monte-Carlo estimates that were within ±1% of the true model mean, and ±5% of the true model variance. This requires repeated model simulations and so it is necessary to have a model that is sufficiently fast to permit this. BICYCLE provides an appropriate compromise between reducing model discrepancy while still running with sufficient computational speed to enable creation of a large ensemble of outputs and hence accurately propagate model parameter uncertainty. As we show later in Figure 4 and Section 6, there are only small differences between the global-average MRA obtained from individual BICYCLE simulations and from the more complex 3D LSG OGCM alternative, yet each BICYCLE simulation can be performed in seconds rather than the weeks required for the LSG OGCM.

Finally, we note that, as described in Heaton et al. (Reference Heaton, Blaauw, Blackwell, Bronk Ramsey, Reimer and Scott2020 in this issue) laying out the statistical methodology for the atmospheric calibration curve, one could calibrate directly against the ensemble of N model outputs if covariance information on the marine calibration curve is desired. However, since current calibration tools use pointwise calibration curve means and variances, we do not discuss this further here.

A Toy Monte-Carlo Example

To illustrate the need for the Monte-Carlo approach, if we wish to accurately estimate the mean output of a computer model and capture the model uncertainty, we consider a simple example. While somewhat pathological it does illustrate the key points that running a computer model at the extreme values of its inputs does not ensure we have captured the full uncertainty in the model output; and also that setting the inputs at their mean values does not result in the mean model output.

Specifically, let us suppose we have a computer model that takes a single input x and returns the value $g\left( x \right) = {x^2}$. Further, assume that we do not know x but can quantify our prior belief that it lies uniformly between −1 and 1, i.e. $X \sim Unif\left[ { - 1,1} \right]$. With such a simple model we can work out exactly what the distribution of the outputs should be. This is shown by the solid blue line in Figure 3. All values in the range [0,1] are clearly possible outputs. Further, we can calculate the true mean output $E\left[ {g\left( X \right)} \right] = {1 \over 3}$ (dot-dashed blue line).

Figure 3 Toy example to illustrate the need for the Monte Carlo approach to accurately estimate both the mean and variability in the output of a computer model.

However, if we run the model at the extreme ends of its possible inputs (either x=1 or −1), in both cases we observe the output g(x) = 1 (green dashed line). If we were to assume that these end-member outputs, generated by the extreme inputs, informed us about the full range of potential outputs then we would falsely assume that the output of the model had no uncertainty. Similarly, to find the true mean of the model output we cannot simply run the model with the input x set to its mean value, i.e. E(X] = 0. Doing so would provide the output g(0) = 0 (red dashed line) rather than the correct mean of ${1 \over 3}$.

On the other hand, if we create a Monte Carlo sample of 500 independent ${X_i} \sim Unif\left[ { - 1,1} \right]$ values and propagate each of these through the model we obtain the ensemble of outputs shown by the histogram. Comparing this against the true output density in solid blue, we can see that our random ensemble is able to capture the correct uncertainty in the model’s output. Further, the mean of these 500 outputs (orange dashed line) shows that, by Monte-Carlo, we are able to accurately estimate the true mean output.

Why Do We Need Monte-Carlo for Marine20?

To be of use for calibration, it is essential to accurately represent the full range of uncertainty in the level of marine radiocarbon. When calibrating a radiocarbon determination, we need to discover all calendar ages which are potentially consistent with that level of ¹⁴C. As such we are required to understand the variability in the output of our, complex and non-linear, carbon cycle computer model. This is somewhat different from much modeling which aims to test specific scenarios. As described in Section 2, in order to capture the full variability in the computer model output we must explore the full range of input scenarios. Understanding this uncertainty in computer model output requires us to consider not just the extremes but also those intermediate scenarios, i.e. a Monte Carlo approach.

Limitations of a Constant ΔR Approach

It is expected that the most significant temporal changes in MRA will be shared and hence occur on a global scale rather than in individual regions, i.e. changes in values over time will predominantly be seen in ${R^{GlobalAv}}\left( \theta \right)$ as opposed to the regional ΔR(θ) adjustments (Stuiver et al. Reference Stuiver, Pearson and Braziunas1986). These global-scale MRA changes are incorporated, by construction, in our global-average Marine20 curve and hence automatically taken into account for all calibration users wishing to use the Marine20 curve.

However, in assuming a constant ΔR we do not permit for potential smaller region-specific temporal changes. This assumption of temporally invariant ΔR values should be adopted with caution in the context of a changing climate and a changing marine carbon cycle. Any global marine radiocarbon calibration product that makes use of modern ΔR values should therefore be used with similar caution.

How Were Marine Data Incorporated into IntCal20?

For time periods where sufficient tree-ring ¹⁴C determinations exist, i.e. continuously from 0–13,913 cal BP, no marine data were used in the creation of the IntCal20 curve. Further back in time however, marine observations did contribute alongside speleothems, floating tree rings and macrofossils from Lake Suigetsu (see Figure 1 and Reimer et al. Reference Reimer, Austin, Bard, Bayliss, Blackwell, Bronk Ramsey, Butzin, Cheng, Edwards, Friedrich, Grootes, Guilderson, Hajdas, Heaton, Hogg, Hughen, Kromer, Manning, Muscheler, Palmer, Pearson, van der Plicht, Reimer, Richards, Scott, Southon, Turney, Wacker, Adophi, Büntgen, Capano, Fahrni, Fogtmann-Schulz, Friedrich, Kudsk, Miyake, Olsen, Reinig, Sakamoto, Sookdeo and Talamo2020 in this issue). To incorporate this marine data into the atmospheric IntCal20, estimates of the MRA for each constituent datum were required. With the exception of Cariaco Basin (Hughen and Heaton Reference Hughen and Heaton2020 in this issue), for which a different approach was taken due to its unusual topography, to obtain these MRA estimates the following approach was taken.

To begin, an interim atmospheric ¹⁴C curve based upon only the Hulu Cave speleothems (Southon et al. Reference Southon, Noronha, Cheng, Edwards and Wang2012; Cheng et al. Reference Cheng, Edwards, Southon, Matsumoto, Feinberg, Sinha, Zhou, Li, Li and Xu2018) was constructed using the same Bayesian spline with errors-in-variables methodology as the final IntCal20 curve (Heaton et al. Reference Heaton, Blaauw, Blackwell, Bronk Ramsey, Reimer and Scott2020 in this issue). This interim Hulu-based ¹⁴C history was then entered as a forcing to the LSG OGCM to get the coarse regional MRA estimates under its GS scenario (Butzin et al. Reference Butzin, Köhler, Heaton and Lohmann2020 in this issue) shown in Figure 3 of Reimer et al. (Reference Reimer, Austin, Bard, Bayliss, Blackwell, Bronk Ramsey, Butzin, Cheng, Edwards, Friedrich, Grootes, Guilderson, Hajdas, Heaton, Hogg, Hughen, Kromer, Manning, Muscheler, Palmer, Pearson, van der Plicht, Reimer, Richards, Scott, Southon, Turney, Wacker, Adophi, Büntgen, Capano, Fahrni, Fogtmann-Schulz, Friedrich, Kudsk, Miyake, Olsen, Reinig, Sakamoto, Sookdeo and Talamo2020 in this issue).

Since LSG OGCM estimates are intended for the open-ocean as opposed to the more coastal locations of the marine IntCal20 data, a further constant coastal radiocarbon reservoir age adjustment/shift was however required. For each marine dataset, an independent prior was placed on this required coastal shift based upon the observed offset between the ¹⁴C ages of the more recent observations from that specific marine dataset and IntCal20’s atmospheric dendrodated trees (which extend back to ca. 14,190 cal BPFootnote ²).

Intuitively, this coastal shift approach used the relative changes (i.e. the shape) in the regional LSG OGCM open-ocean MRA estimates but not their absolute values. The internally estimated coastal shift is equivalent to estimating a temporally invariant pseudo-ΔR for the specific coastal location relative to the nearby open-ocean. These coastal pseudo-ΔR shifts from the GS scenario of LSG OGCM ranged from approximately 12 ¹⁴C yr in Vanuatu up to over 280 ¹⁴C yr in Kiritimati indicating the difficulty in using open-ocean OGCM estimates directly for coastal data.

These preliminary Hulu-Cave based MRA estimates are expected to be overly smooth (since they are based upon the relatively smooth Hulu Cave record) however they should provide robust first-order approximations that enable marine data to contribute to IntCal20. The region-specific prior MRA estimates were then used as inputs for the Bayesian approach to IntCal20’s curve construction allowing for some further fine-scale MRA variability not captured by the coarse preliminary Hulu-based estimates. In this curve construction process, the coastal shift priors were further updated to resolve differences between the diverse pre-Holocene data sets, see Heaton et al. (Reference Heaton, Blaauw, Blackwell, Bronk Ramsey, Reimer and Scott2020 in this issue) for complete details.

Importantly, as shown in Figure 3 of Reimer et al. (Reference Reimer, Austin, Bard, Bayliss, Blackwell, Bronk Ramsey, Butzin, Cheng, Edwards, Friedrich, Grootes, Guilderson, Hajdas, Heaton, Hogg, Hughen, Kromer, Manning, Muscheler, Palmer, Pearson, van der Plicht, Reimer, Richards, Scott, Southon, Turney, Wacker, Adophi, Büntgen, Capano, Fahrni, Fogtmann-Schulz, Friedrich, Kudsk, Miyake, Olsen, Reinig, Sakamoto, Sookdeo and Talamo2020 in this issue), for locations of the marine data used for making IntCal20, the relative differences in the shapes of the Hulu-based regional open-ocean MRA estimates provided by the LSG OGCM were small (in the order of ±5% of the absolute value of MRA). Furthermore, these regional open-ocean MRA estimates had similar shapes to the non-polar global-average MRA. Hence, in terms of the final MRAs used for the marine data contributing to IntCal20, there would be little difference in having applied individual temporally invariant coastal shifts to the non-polar global-average MRA (i.e. the approach we recommend for Marine20 to restrict temporal changes in MRA to the global-average) of the LSG OGCM rather than the region-specific output. Our proposed global-average Marine20 approach to calibration is therefore consistent with that taken to integrate the marine data into IntCal20.

Why Not Use OGCM Output Directly to Provide Regional Calibration Curves?

Much of the marine data which users wish to calibrate come from coastal regions. With our current state of knowledge, most global OGCMs do not have the resolution to simulate these coastal regions accurately and their results may be flawed at such sites. This includes the LSG OGCM. We are therefore still recommending, for the current update, the use of a global-average radiocarbon age calibration curve with the adjustment of a constant ΔR estimated independently based upon the offset between local observations and the Marine20 curve.

In Figure 4, we plot the pre-Holocene LSG OGCM regional open-ocean MRA estimates obtained when forced by the pointwise mean of the final IntCal20 curve (Reimer et al. Reference Reimer, Austin, Bard, Bayliss, Blackwell, Bronk Ramsey, Butzin, Cheng, Edwards, Friedrich, Grootes, Guilderson, Hajdas, Heaton, Hogg, Hughen, Kromer, Manning, Muscheler, Palmer, Pearson, van der Plicht, Reimer, Richards, Scott, Southon, Turney, Wacker, Adophi, Büntgen, Capano, Fahrni, Fogtmann-Schulz, Friedrich, Kudsk, Miyake, Olsen, Reinig, Sakamoto, Sookdeo and Talamo2020 in this issue) under the GS scenario. These estimates are the closest LSG OGCM analogue to our Marine20 approach but do not incorporate any uncertainty propagation as Monte-Carlo is infeasible due to the LSG OGCM’s long run time. Shown are the LSG OGCM’s MRA estimates for the open-ocean sites lying closest to each of the marine locations used in IntCal20. Also shown are the LSG OGCM (50°S–50°N) globally averaged estimate with the same IntCal20 pointwise mean forcing and GS scenario, and also the Monte-Carlo mean of our BICYCLE/Marine20 globally averaged MRA discussed in detail later.

Figure 4 Regional open-ocean MRA estimates generated by the LSG OGCM when forced by the mean of the IntCal20 (Reimer et al. Reference Reimer, Austin, Bard, Bayliss, Blackwell, Bronk Ramsey, Butzin, Cheng, Edwards, Friedrich, Grootes, Guilderson, Hajdas, Heaton, Hogg, Hughen, Kromer, Manning, Muscheler, Palmer, Pearson, van der Plicht, Reimer, Richards, Scott, Southon, Turney, Wacker, Adophi, Büntgen, Capano, Fahrni, Fogtmann-Schulz, Friedrich, Kudsk, Miyake, Olsen, Reinig, Sakamoto, Sookdeo and Talamo2020 in this issue) reconstruction of NH atmospheric Δ¹⁴C. These estimates are obtained under LSG’s GS scenario (Butzin et al. Reference Butzin, Köhler, Heaton and Lohmann2020 in this issue) and correspond to the open-ocean sites nearest to the locations of the IntCal20 marine data. Also shown is the LSG (50°S–50°N) global-average MRA estimate in the GS scenario; and also the mean of BICYCLE’s global-average MRA for comparison.

A potential approach to adapt this open-ocean LSG OGCM output for calibration in wider and/or coastal locations would be that taken for inclusion of marine data in IntCal20. One would initially identify the regional open-ocean LSG OGCM site that lies closest to the particular location of interest, before then applying a constant, baseline-specific, pseudo-ΔR shift from this specific regional open-ocean LSG OGCM estimate. However, extending this approach to Marine20 would be complex and reliant upon users identifying both the correct LSG OGCM regional baseline and the matching pseudo-ΔR shift.

Furthermore, for the non-polar sites illustrated in Figure 4, the relative shapes of the regional open-ocean MRAs estimated by the LSG OGCM are highly similar to both the LSG global-average estimate and also the mean BICYCLE global-average estimate. Changing to two alternative scenarios (CS or PD, Butzin et al. Reference Butzin, Köhler, Heaton and Lohmann2020 in this issue) within the LSG OGCM also has little effect on the relative shapes of the resultant regional open-ocean MRA estimates—the primary effect of such scenario changes are simply constant shifts, up and down respectively, in each regional open-ocean MRA estimate albeit with the size of these shifts potentially varying between regions.

Consequently, the difference in using constant pseudo-ΔR-type shifts from regional LSG OGCM curves, as opposed to constant ΔR shifts from the global-average curves, would be small for these sites. As a result, the proposed approach to radiocarbon age calibration described here, of using a global-average Marine20 radiocarbon age calibration curve based on BICYCLE with application of constant ΔR shifts, is consistent with an approach using regional LSG OGCM curves while also being much simpler and more efficient.

As discussed in more detail in Section 6, the output from the LSG OGCM matches that from BICYCLE at the global-average level suggesting that BICYCLE is a satisfactory model at this broader scale. However, by running much more quickly than the full LSG OCGM, BICYCLE provides greatly more flexibility which provides several key benefits for our calibration curve needs. In particular, LSG (or any other current OGCM) is too computationally expensive to carry out a rigorous exploration of the uncertainty in its output. The use of BICYCLE enables us to explore a much wider range of potential carbon cycle scenarios. Further, the LSG OGCM simulations were kept in a single scenario for the entire period from 0–55 cal kBP. The effect of a transient switch from a glacial scenario to an interglacial, more suitable for the Holocene, within a simulation was therefore unknown. With BICYCLE we can more easily investigate the effect of such scenario changes.

Our BICYCLE approach attempts to minimise potential inaccuracies in marine radiocarbon variability that stem from strictly unknown changes in the marine carbon cycle by exploring a wide range of possible carbon cycle states, providing a trade-off between expected accuracy and precision in resulting calibrated ages. As we describe later in Section 7, key future work lies in developing 3D OGCM models which can provide region-specific output further and understanding their inherent uncertainties. Work towards this goal would be invaluable in improving both marine calibration and our insight into the global carbon cycle.

Users calibrating open-ocean samples who wish to use individual regional MRA estimates provided by the LSG OGCM can find them on PANGAEA (https://doi.org/10.1594/PANGAEA.914500). This may also be relevant for users wishing to calibrate data outside Marine20’s intended range, i.e. where sea ice may have been present (north of 40/50°N or south of 40°S), especially if the LSG OGCM estimates exhibit significantly different relative changes. One approach to incorporate LSG OGCM MRA estimates, for example into OxCal, can be found via Alves et al. (Reference Alves, Macario, Urrutia, Cardoso and Bronk Ramsey2019).

Propagating Model Parameter Uncertainty through BICYCLE

While BICYCLE is a box model, it still contains a large number of parameters. To keep our approach practical, we focus our analysis of propagated model uncertainties on key processes and do not consider the uncertainty in all parameters. From sensitivity tests we identified that two processes are most important for simulated surface ocean Δ¹⁴C: piston velocity; and Atlantic meridional overturning circulation (AMOC), which is directly connected to the strength of North Atlantic Deep Water (NADW) formation.

Piston velocity is the rate of air-sea gas exchange and depends on near-surface wind speed. Thus, this parameter is important for the oceanic uptake of ¹⁴C from the atmosphere. In all previous applications, it was parametrized to 2.5 × 10⁻⁵ and 7.5 × 10⁻⁵ m s⁻¹, for the equatorial and high latitude surface boxes respectively (Heimann and Monfray Reference Heimann and Monfray1989). These values agree reasonably well with more recent estimates (Wanninkhof Reference Wanninkhof2014). To create Marine20, we maintain these values as the mean piston velocities but introduce uncertainty by placing a prior on each informed by a meta-analysis as described below.

The strength of the AMOC in the model is important for the transport of tracers from the surface to the deep ocean via the physical carbon pump. In our Marine20 setup, the AMOC, represented by strength of the NADW formation, is considered to have two states with some uncertainty on the value in each state. In interglacials, it is centred around 16 Sv (1 Sverdrup = 10⁶ m³ s⁻¹, e.g. Talley et al. Reference Talley, Pickard, Emery and Swift2011) while during glacials it is centred at 10 Sv. See below for more details.

In addition to uncertainties in piston velocity and AMOC, we aim to propagate uncertainty on the past levels of atmospheric Δ¹⁴C and CO₂. Below, we present a summary of how these various inputs were selected together with a description of their uncertainties. Those model parameters for which we consider uncertainty, and specify distributions on their values, are referred to as having uncertainty incorporated; while those for which we do not consider uncertainty are denoted as having uncertainty not incorporated. An overview of how these enter BICYCLE is given in Figure 2.

Piston Velocity:

Piston velocity is a key factor regulating the uptake of ¹⁴C by the ocean and hence largely determines simulated MRAs. Piston velocity depends on near-surface wind speed for which our knowledge is limited, in particular regarding values of the past. Therefore, piston velocities are perhaps the least well-constrained parameters in our model.

As stated above, BICYCLE operates with two piston velocities–dependent upon whether the surface ocean box falls in low or high-latitudes, which typically have either low or high wind speeds, respectively. Since BICYCLE has been tuned previously, and there is no common agreement on their absolute values, we do not change the mean values of these parameters as to do so would lead to a different carbon cycle baseline, potentially requiring further adjustments/analysis. We do however incorporate uncertainties around these means.

In specifying piston velocities and their uncertainties, we are not concerned with short-term, annual-scale, fluctuations. Short-term changes do not have a significant effect on the marine ¹⁴C reservoir since the slow air/ocean mixing process occurs over much longer time scales and so smooths such fine variations out. Instead, our interest is in quantifying the uncertainty in the long-term mean values. We model the piston velocities as constant over time and specify a ±15% (1σ) uncertainty on both the high-latitude and equatorial values:

$${\rm{Equatorial \ Boxes}}\,\,{k_{Eq}} \sim N\left( {2.5 \times {{10}^{ - 5}},{{\left( {0.375 \times {{10}^{ - 5}}} \right)}^2}} \right){\rm{m}} \,{{\rm{s}}^{ - 1}},$$

$${\rm{High - Latitude \ Boxes}}\,\,{k_{Hl}} \sim N\left( {7.5 \times {{10}^{ - 5}},{{\left( {1.125 \times {{10}^{ - 5}}} \right)}^2}} \right){\rm{m}} \,{{\rm{s}}^{ - 1}}.$$

Justification for Selection of Piston Velocity Uncertainties

To quantify historic piston velocity uncertainty, we first consider present day mean piston velocity with present day winds/climate. Naegler (Reference Naegler2009) provides five current day piston velocity estimates (with individual uncertainties) on comparable scales. We combine these values using a statistically rigorous meta-analysis (Figure 6) to give an estimate for the average present day piston velocity of $\overline k \sim N\left( {17,{{1.4}^2}} \right)$ cm hr^–1 (or equivalently N(4.7, 0.39²) 10^–5 m s^–1), i.e. the 1σ uncertainty on the level of piston velocity with present day climate is approximately ±10%.

Figure 6 Meta-analysis of present-day estimates of piston velocity from Naegler (Reference Naegler2009). Shown are, post-Naegler’s adjustments, the piston velocities (cm hr^–1) from 5 different studies together with their 95% confidence intervals (95%-CI). The grey diamonds represent the combined estimate, again with 95% intervals. Heterogeneity, as quantified by the measures ${I^2},{\tau ^2}$ and p (Cochran’s Q-statistic), assesses the extent to which the reported velocities differ between the different studies, see e.g. Borenstein et al. (Reference Borenstein, Hedges, Higgins and Rothstein2011) for more details. Here, the five reported velocities are seen to be consistent with a single shared underlying velocity.

Further uncertainty arises as a result of changing past climate, in particular different wind speeds. From the LSG OGCM (Butzin et al. Reference Butzin, Köhler and Lohmann2017, Reference Butzin, Köhler, Heaton and Lohmann2020 in this issue) we obtain (Table 1) gas exchange rates based on 10-year averaged monthly model results for three different climate scenarios (PD, GS and CS, Butzin et al. Reference Butzin, Prange and Lohmann2005).

Table 1 Piston velocities for ¹⁴CO₂ employed by the LSG OGCM according to various climate forcing scenarios, results are annual-mean values averaged over ocean areas considering the same latitude range as used by BICYCLE. The scenarios considered are PD (a present-day climate scenario which may also be considered as a surrogate for interstadials), and two glacial scenarios—CS based on the CLIMAP reconstruction (McIntyre et al. Reference McIntyre1981) and GS based on the GLAMAP climate reconstruction (Sarnthein et al. Reference Sarnthein, Gersonde, Niebler, Pflaumann, Spielhagen, Thiede, Wefer and Weinelt2003). The glacial climate scenarios involve some adjustment of atmospheric forcing fields in the tropics (CS) and in the Southern Ocean (both CS and GS; see Butzin et al. Reference Butzin, Prange and Lohmann2005 for further details).

We use these values to indicate the potential additional variability introduced by different past climate scenarios. These estimates suggest that, due to changes in climate, the mean piston velocity might change, within any particular basin, from its present day value by approximately another ±10% (at 1σ based upon the difference in the LSG OGCM estimates between the interglacial PD scenario and the glacial scenario CS). Finally, we combine the uncertainty on the present day value (±10%) with the LSG-based intra-basin climate variation (±10%) to get an approximate combined uncertainty of $\sqrt {{{0.1}^2} + {{0.1}^2}} \simeq 0.15,$ i.e. σ = 15%. This uncertainty is comparable with the ±20% proposed by Wanninkhof (Reference Wanninkhof2014) for the mean value of the present ocean.

Ocean Circulation

The strength of the AMOC, expressed by NADW formation, presumably varied between a weak glacial and a strong interglacial mode. In our model, the switch between both is a function of North Atlantic sea surface temperature and has been chosen such that the interglacial mode is reached at the onset of the Holocene. In each mode, we assume a constant, but unknown NADW value drawn from a suitable prior representing the two different strengths. Millennial-scale variations in the Atlantic meridional overturning circulation connected with the bipolar seesaw (e.g. Barker et al. Reference Barker, Diz, Vantravers, Pike, Knorr, Hall and Broecker2009) have been neglected for this update although are a topic identified for future work (see Section 7):

$$NADW\left( \theta \right) = \left\{ {\matrix{ {{\alpha _{Glacial}}\ if\ t \in Glacial \,mode,} \cr {{\alpha _{Inter}}\ if\ t \in Interglacial \,mode,} \cr } } \right.$$

where:

$$\matrix{ {{\alpha _{Glacial}} \sim N\left( {10,{2^2}} \right),} \cr {{\alpha _{Inter}} \sim N\left( {16,{1^2}} \right).} \cr } $$

Comparison with the LSG OGCM

The LSG model is a three-dimensional ocean general circulation model incorporating more processes (ocean physics) than a box model. Furthermore, it has the ability to provide location-specific MRAs. Consequently, one might expect the LSG OGCM to provide more detailed and accurate modeling of the global-average marine ¹⁴C reservoir than BICYCLE, albeit without the current ability to propagate model parameter uncertainty via Monte-Carlo.

As discussed in Section 3, the LSG OGCM was applied in three different climate scenarios: PD, CS and GS. Details are found in Table 1 and elsewhere (Butzin et al. Reference Butzin, Prange and Lohmann2005, Reference Butzin, Prange and Lohmann2012, Reference Butzin, Köhler and Lohmann2017, Reference Butzin, Köhler, Heaton and Lohmann2020 in this issue). Similarly to BICYCLE, the model was forced with the temporal changes in the atmospheric carbon records of CO₂ (Köhler et al. Reference Köhler, Nehrbass-Ahles, Schmitt, Stocker and Fischer2017) and IntCal20’s reconstruction of Δ¹⁴C (Reimer et al. Reference Reimer, Austin, Bard, Bayliss, Blackwell, Bronk Ramsey, Butzin, Cheng, Edwards, Friedrich, Grootes, Guilderson, Hajdas, Heaton, Hogg, Hughen, Kromer, Manning, Muscheler, Palmer, Pearson, van der Plicht, Reimer, Richards, Scott, Southon, Turney, Wacker, Adophi, Büntgen, Capano, Fahrni, Fogtmann-Schulz, Friedrich, Kudsk, Miyake, Olsen, Reinig, Sakamoto, Sookdeo and Talamo2020 in this issue), although using the pointwise means rather than individual realizations.

For each scenario, the LSG OGCM provided estimates of surface MRAs averaged over 50ºN–50ºS for the ocean surface above 50 m depth. These are shown in Figure 7B alongside the BICYCLE estimates with their 95% probability intervals. As can be seen, despite BICYCLE being a much simpler model, the BICYCLE estimates retain almost all of the temporal details present within the spatially averaged results from the LSG OGCM. BICYCLE’s probability interval also covers the range of LSG scenarios. Users should however be aware that the LSG scenario which has the highest agreement with the mean output from BICYCLE is that with present day climate (PD). This emphasises, that although comparable, climate change over our time window of interest leads to different implications in the two models when considered in detail. Nevertheless, the overall model comparison suggests that, despite its simplicity, BICYCLE is able to provide a robust and accurate global estimate.

Comparison with Present-Day Marine Observations

We also compare Marine20 against pre-bomb (0–200 cal BP) marine observations (Reimer and Reimer Reference Reimer and Reimer2001). In this comparison, it should be noted that these observations include regions of coastal upwelling and further, due to sampling biases, much of this data arise from the Western US coast, so may not be representative of the open ocean. Both data and model have not been corrected for fossil fuel contamination (Suess effect), thus should be comparable. Note that BICYCLE has previously been shown to track the Suess effect sufficiently accurately, even when applied without forced atmospheric Δ¹⁴C and CO₂ (Köhler et al. Reference Köhler, Knorr and Bard2014; Köhler Reference Köhler2016). We restrict our comparison to marine data from 50ºN–50ºS, to remove potential effects due to sea ice cover, and show the comparison in both Δ¹⁴C and radiocarbon age space (Figure 9, panels A and B). The previous Marine13 estimates are also shown. Figure 9C shows a zoomed-in plot of the estimated global-average MRA over this period.

Figure 9 Zoom-in on the pre-bomb time window (0–200 cal BP). IntCal20, Marine20, Marine13 and ¹⁴C determinations of observed marine samples (restricted to 50ºN–50ºS) in (A) Δ¹⁴C space, (B) Radiocarbon age space, and (C) illustrating the estimated global-average MRA. In (A, B) ¹⁴C determinations from observed marine samples are taken from the data base (http://calib.org/marine/), in (A) a subset of additional Δ¹⁴C data is also shown (Toggweiler et al. Reference Toggweiler, Druffel, Key and Galbraith2019). In (C) the global-average MRA from the LSG OGCM output (restricted to 50ºN–50ºS) for the scenario PD is also given for comparison. IntCal20 is shown with its 95% predictive intervals, while Marine13 and Marine20 their 95% probability intervals.

Figure 10 Plot of observed marine ¹⁴C ages (with no MRA or ΔR correction) against the Marine20 and IntCal20 curves 0–55 cal kBP. The datasets shown consist of: foraminifera (forams) from the Iberian and Pakistan Margins (Bard et al. Reference Bard, Ménot, Rostek, Licari, Böning, Edwards, Cheng, Wang and Heaton2013); corals from Tahiti and Barbados (Bard et al. Reference Bard, Hamelin, Fairbanks and Zindler1990, Reference Bard, Arnold, Hamelin, Tisnerat-Laborde and Cabioch1998, Reference Bard, Ménot-Combes and Rostek2004); corals from Vanuatu and Papua New Guinea (Cutler et al. Reference Cutler, Gray, Burr, Edwards, Taylor, Cabioch, Beck, Cheng and Moore2004); corals from Tahiti (Durand et al. Reference Durand, Deschamps, Bard, Hamelin, Camoin, Thomas, Henderson, Yokoyama and Matsuzaki2013); corals from Vanuatu, Kiritimati and Barbados (Fairbanks et al. Reference Fairbanks, Mortlock, Chiu, Cao, Kaplan, Guilderson, Fairbanks, Bloom, Grootes and Nadeau2005); corals from Papua New Guinea (Edwards et al. Reference Edwards, Beck, Burr, Donahue, Chappell, Bloom, Druffel and Taylor1993); and corals from Papua New Guinea and Vanuatu (Burr et al. Reference Burr, Beck, Taylor, Recy, Edwards, Cabioch, Correge, Donahue and O’Malley1998, Reference Burr, Galang, Taylor, Gallup, Edwards, Cutler and Quirk2004). Plotted are the 95% probability intervals for both the calendar age and radiocarbon age of each observation; together with the 95% probability/predictive intervals for Marine20 and IntCal20, respectively.

As we see, the modeled global-average MRA for Marine20 is significantly higher than in Marine13 (Figure 9C). This leads us to a Marine20 curve which consistently lies above Marine13 in terms of the radiocarbon age, or equivalently gives lower Δ¹⁴C values, for a particular calendar age θ.

In the Δ¹⁴C space, the weighted mean of the sampled marine observations covering 0–200 cal BP is –61 ± 16‰, agreeing with a further subset of samples published recently (Toggweiler et al. Reference Toggweiler, Druffel, Key and Galbraith2019) for 1940–1954 CE (i.e. 10 cal BP to –4 cal BP).

Although an over-time decreasing trend in Δ¹⁴C due to the Suess effect seems visible in the observed marine data from ca. 1900 CE onward (i.e. 50–0 cal BP) it is statistically non-significant (p-value 0.35). Results are nearly identical when restricted to 40ºN–40ºS in the data selection, although this leads to a weighted Δ¹⁴C mean of only –58‰. In the pre-bomb time window, the global-average MRA in Marine20 drops from slightly above 550 ¹⁴C yr at 200 cal BP to ~410 ¹⁴C at 0 cal BP, mostly due to the Suess effect seen in the atmospheric IntCal20 Δ¹⁴C. This most-recent trend is also very closely reproduced in the simulations with the LSG OGCM (Figure 9C).

The increase in the Δ¹⁴C deficit between the atmospheric and global-average surface ocean that we see with BICYCLE and Marine20, compared to the simpler box model of Marine13, results in a calibration curve which provides better agreement with our observed present-day observations. As seen in Figure 9A and 9B, the majority of these marine observations have an older radiocarbon age (i.e. their Δ¹⁴C is lower) than Marine13 and lie closer to Marine20. A further analysis suggested that the optimal offset to Marine20, in terms of maximising the statistical likelihood of the observed present-day data shown in Figure 9, was only 30 ¹⁴C yr. This offset was of sufficiently small size to lie within sampling uncertainty and so no adjustments of the BICYCLE simulations have been applied which might have led to a complete vanishing of the model/data offset.

Comparison with Corals and Foraminifera in IntCal20 Database

Figure 10 plots the complete Marine20 curve, in radiocarbon age space and from 0–55 cal kBP, against all the raw coral and foraminifera data available in the IntCal20 database, http://intcal.org. As described in the IntCal20 paper (Reimer et al. Reference Reimer, Austin, Bard, Bayliss, Blackwell, Bronk Ramsey, Butzin, Cheng, Edwards, Friedrich, Grootes, Guilderson, Hajdas, Heaton, Hogg, Hughen, Kromer, Manning, Muscheler, Palmer, Pearson, van der Plicht, Reimer, Richards, Scott, Southon, Turney, Wacker, Adophi, Büntgen, Capano, Fahrni, Fogtmann-Schulz, Friedrich, Kudsk, Miyake, Olsen, Reinig, Sakamoto, Sookdeo and Talamo2020 in this issue), corals which are older than 25 cal kBP have been excluded in view of potential diagenesis. Further, we do not include observations from the Cariaco Basin for comparison here since this basin appears to be a unique marine environment (Hughen and Heaton Reference Hughen and Heaton2020 in this issue) which no model is currently able to resolve, see Section 7 describing the limitations of our model.

We plot here the raw ¹⁴C ages of the observations, with no MRA or ΔR correction, along with their 95% probability intervals in both radiocarbon and calendar age. Since we have made no marine radiocarbon age correction for the observations in the plot, we would not expect the Marine20 curve to necessarily fit the raw data themselves. Rather our interest is in whether the offset from the curve (i.e. ΔR) is consistent for each dataset or location. Since marine data from older than 13.91 cal kBP does influence the construction of IntCal20, and hence also Marine20, the Marine20 curve is not entirely independent of the plotted data for this time period—although the large volume of other data informing IntCal20 still mean it is significantly so. From 0–13.91 cal kBP, the Marine20 curve is almost entirely independent of the plotted data since no marine data informed IntCal20 in this more recent period.

As with the present-day observations discussed in the previous section, most of this marine data is coastal and hence may not be representative of open-ocean. Most of the corals are seen to lie below the Marine20 curve suggesting that the locations of these data (Tahiti, Barbados, Vanuatu, Papua and Kiritimati) have a rather smaller region-specific MRA than the Marine20 global-average estimates, i.e. have a negative ΔR. For the period between 7 and 12 cal kBP (i.e. the early Holocene) we also see that the MRA for Tahiti is consistently smaller than for Papua New Guinea and Vanuatu, which is compatible with previous work (e.g. Table 1 in Reimer et al. Reference Reimer, Baillie, Bard, Bayliss, Beck, Blackwell, Bronk Ramsey, Buck, Burr, Edwards, Friedrich, Grootes, Guilderson, Hajdas, Heaton, Hogg, Hughen, Kaiser, Kromer, McCormac, Manning, Reimer, Richards, Southon, Talamo, Turney, van der Plicht and Weyhenmeyer2009). Beyond 12 cal kBP systematic differences are more difficult to detect because the analytical precision is lower.

Interestingly, between 15 and 17 cal kBP, the Iberian Margin data seem several centuries too old in ¹⁴C yr compared with the Marine20 curve. This age range corresponds to the Heinrich 1 stadial and is potentially compatible with the near collapse of AMOC during Heinrich events which are not incorporated into BICYCLE and therefore missing in Marine20. It is harder to detect if this effect is also present in older Heinrich events since the precision of ¹⁴C ages decreases through time, and some of these older events may have also had reduced intensity.