A Short Note on Marine Reservoir Age Simulations Used in IntCal20

ABSTRACT Beyond ~13.9 cal kBP, the IntCal20 radiocarbon (14C) calibration curve is based upon combining data across a range of different archives including corals and planktic foraminifera. In order to reliably incorporate such marine data into an atmospheric curve, we need to resolve these records into their constituent atmospheric signal and marine reservoir age. We present results of marine reservoir age simulations enabling this resolution, applying the LSG ocean general circulation model forced with various climatic background conditions and with atmospheric radiocarbon changes according to the Hulu Cave speleothem record. Simulating the spatiotemporal evolution of reservoir ages between 54,000 and 10,700 cal BP, we find reservoir ages between 500 and 1400 yr in the low- and mid-latitudes, but also more than 3000 yr in the polar seas. Our results are broadly in agreement with available marine radiocarbon reconstructions, with the caveat that continental margins, marginal seas, or tropical lagoons are not properly resolved in our coarse-resolution model.


INTRODUCTION
The spatiotemporal variability of marine radiocarbon ( 14 C) records is superimposed by a systematic isotopic depletion with respect to the atmosphere. This effect is frequently expressed as marine reservoir age (MRA) and has to be taken into account when 14 C calibration considers marine archives (Alves et al. 2018). Prebomb MRAs have ranged from~400 yr in subtropical oceans to more than 1000 yr in polar seas (Key et al. 2004). Prior to the continuous tree-ring record, MRAs are poorly constrained and have to be inferred from ad-hoc assumptions or through modeling. IntCal20 aims to reconstruct the atmospheric 14 C concentration at mid-latitudes in the Northern Hemisphere. From 0 to~13.9 cal kBP, IntCal20 is constructed from dendrochronologically dated tree rings providing direct atmospheric observation. However, further back in time, the lack of available trees means that a variety of alternative archives are required including macro-fossils, speleothems, corals and planktic foraminifera. Each of these archives have their own specific characteristics. In particular, the use of data from marine environments (such as corals and foraminifera) is complicated by the MRA meaning such observations are offset, probably floating in space and time, from the atmospheric 14 C we wish to reconstruct. To reliably reconstruct atmospheric 14 C, we therefore need to be able to resolve such records into their constituent atmospheric and MRA components.
To provide such a resolution for IntCal20, we have taken the following approach. Firstly, we construct a preliminary estimate for atmospheric 14 C based solely upon the Hulu cave record (Southon et al. 2012;Cheng et al. 2018). This Hulu-Cave-based atmospheric 14 C reconstruction is then fed into an ocean general circulation model (Butzin et al. 2017) to provide temporal MRA estimates at all of our marine sites contributing to IntCal. These MRA estimates are then fed back into the construction of the main IntCal20 atmospheric curve allowing us to compile the MRA-adjusted marine records alongside the other archives, see Heaton et al. (2020 in this issue) for further details. Such an approach can be seen somewhat analogously to the first step in an iterative backfitting procedure (or a single update of the MRA within a wider Markov Chain Monte Carlo method). While ideally one may wish to iterate this procedure several times to improve MRA estimation, we believe a single run will provide a reliable first-order approximation. Here, we give an overview of the applied simulation scenarios to create these MRA estimates given the Hulu-Cave-based curve and briefly discuss the results.

METHOD
We employ the Hamburg Large Scale Geostrophic (LSG) ocean general circulation model (Maier-Reimer et al. 1993) with a horizontal resolution of 3.5°and a vertical resolution of 22 unevenly spaced levels. The LSG model has been constantly improved by including a bottom boundary layer scheme (Lohmann 1998), a sophisticated numerical advection scheme (Schäfer-Neth and Paul 2001;Prange et al. 2003) and has demonstrated its skills in 14 C simulations (Butzin et al. 2005(Butzin et al. , 2012(Butzin et al. , 2017. In the current setup, the model can simulate the entire radiocarbon timescale (~60 kyr) within a few days of CPU time on state-of-theart computer systems. The model is forced with monthly fields of recent and glacial wind stress, surface air temperature, and freshwater flux derived in previous climate simulations (Lohmann and Lorenz 2000;Prange et al. 2004). In the same way as in our precursor study (Butzin et al. 2017), we consider three climate forcing scenarios to assess the impact of past ocean-climate variability on marine 14 C records. To summarize, scenario PD employs present-day climate background conditions approximating the Holocene and interstadials. Glacial scenario GS aims at representing the Last Glacial Maximum (LGM), featuring a shallower Atlantic meridional overturning circulation (AMOC) weakened by about 30% compared to PD. A second glacial climate scenario (CS) mimics cold stadials with further AMOC weakening by about 60%. A thorough discussion of these scenarios can be found in Butzin et al. (2005).
Radiocarbon is simulated on-line as F 14 R oce-atm (the 14 C enrichment of the ocean relative to the contemporaneous atmosphere, Soulet et al. 2016) following Toggweiler et al. (1989). That is, we do not consider 14 C and 12 C separately but directly simulate their fractionation-corrected ratio. This approach neglects biological effects which are one order of magnitude smaller than the effects of ocean circulation and radioactive decay on F 14 R oce-atm (Fiadeiro 1982). Oceanic uptake of 14 C is calculated according to Sweeney et al. (2007), using atmospheric CO 2 concentrations as compiled in the spline of Köhler et al. (2017) and prescribed timeinvariant concentrations of dissolved inorganic carbon in surface water as simulated by Hesse et al. (2011). As described in the introduction, atmospheric 14 C forcing was based solely upon the Hulu Cave speleothem record with a Dead Carbon Fraction (DCF) of about (450 ± 70) 14 C yr (Southon et al. 2012;Cheng et al. 2018). To create this Hulu-Cave-based atmospheric curve, the same Bayesian spline errors-in-variables statistical methodology was used as in construction of the final IntCal20 curve (Heaton et al. 2020 in this issue). Within each individual speleothem, the DCF was considered constant but taking potentially different levels for each speleothem in the cave. Priors for the DCF of each speleothem were taken from the relevant papers of Southon et al. (2012) and Cheng et al. (2018). The spline-based Markov Chain Monte Carlo (MCMC) was run for 250,000 iterations (with the first half being discarded as burn-in) to obtain a posterior mean estimate for the Hulubased 14 C atmospheric reconstruction along with upper and lower 95% credible (2-sigma) pointwise intervals (see Figure 1).
We combine all three ocean climates (PD, GS, and CS) and all three 14 C forcing scenarios (posterior mean, upper and lower 95% credible intervals) to create nine numerical experiments. Each experiment is spun up over 20,000 yr with fixed atmospheric values of 14 C and CO 2 at 54 cal kBP (Δ 14 C atm = 73-504‰ and pCO 2 = 219 μatm), before the model is run transiently forward in time until 10.7 cal kBP forced by time-variable 14 C and pCO 2 values. Marine reservoir ages are calculated afterwards in 14 C-yr as MRA = -8033 × ln (F 14 R oce-atm ), based on simulation results for the upper 50 m of the ocean corresponding to the typical habitat depth of corals and foraminifera analyzed in IntCal20. The impact of different ocean surface depth ranges on MRA has been analyzed in Butzin et al. (2017) and is not investigated any further here. Figure 1 shows the simulated history of global-mean MRAs. The envelope spans all MRA simulations and may be interpreted as total uncertainty range. A closer look at the MRA uncertainty of any specific ocean climate scenario, e.g. scenario GS, indicates that the contribution of the atmospheric 14 C forcing uncertainty is small between 10.7 and~44 cal kBP. That is, the uncertainty of simulated MRAs is largely due to the difference in the three climate scenarios, with the upper or lower bound corresponding to scenario CS or PD, respectively.

RESULTS AND DISCUSSION
In the following discussion we focus on scenario GS, as (most) MRA estimates in IntCal20 are based on this scenario, which also corresponds to the ensemble median of simulated MRAs when forced with the mean atmospheric 14 C record. Disregarding high-latitude regions where sea ice inhibits oceanic 14 CO 2 uptake, we find global mean MRA variations between 500 and 1400 14 C yr. The lowest reservoir ages are simulated during the last deglaciation with two minima around 11 and 14 cal kBP caused by the combination of small atmospheric 14 C values and rising concentrations of atmospheric CO 2 , the latter allowing a higher oceanic 14 CO 2 uptake. At the LGM (19-23 cal kBP) the global mean MRA was about 800 14 C yr. Prior to the LGM our simulations show variations in global mean MRA of several hundreds of years with positive anomalies associated with the Mono Lake (~34 cal kBP, Lund et al. 2017) and Laschamps (~42 cal kBP, Lascu et al. 2016) geomagnetic excursions. In our precursor study (considering scenario GS and atmospheric 14 C according to IntCal13, black line in Figure 1) the amplitude of these variations was slightly larger between 27 cal kBP and 34 cal kBP, but more than 200 14 C yr smaller around the Laschamps excursion. This is due to the revised atmospheric 14 C forcing during this period. Note that the global mean MRA peak at 42 cal kBP precedes the atmospheric 14 C maximum by about 3000 cal yr. Instead, the global mean MRA peak coincides with an inflection point of the atmospheric 14 C forcing curve after which the slowly adjusting ocean is able to catch up with the 14 C excursion in the atmosphere.
Since the climate scenarios applied here are identical to those used in our previous work, we also find-as expected-similar results: the youngest surface waters are always found in the subtropical oceans and the oldest surface waters in the polar seas, where ice cover can boost MRAs to more than 3000 14 C yr. During the LGM the general spatial pattern of the MRA is in line with the few data-based MRA reconstructions (Skinner et al. 2017) available for this time slice (Figure 2).
Here, the lowest MRAs are simulated for the subtropical Northeast Pacific where 14 C paleorecords are particularly scarce. However, in no scenario we arrive at such low MRAs close to preindustrial values as contained in the reconstruction for the subtropical Northwest Atlantic (Skinner et al. 2017). To assess potential model biases not only for a single date but for the entire simulation period, we compare reconstructed marine Δ 14 C records contributing to IntCal20 with marine Δ 14 C simulated at nearby grid cells, keeping in mind that none of these sites is representative for open-ocean conditions for which the LSG model has been designed. This is shown as scatterplot in Figure 3, in which unbiased model results and reconstructions should align along a straight line with slope = 1 and intercept = 0. Our simulations frequently underestimate marine Δ 14 C prior to~30 cal kBP where the uncertainties in the reconstructions are large. Simulated Δ 14 C is also low in the Cariaco Basin prior to~19 cal kBP.
The underestimate of simulated Δ 14 C suggests that the MRAs derived for these data points may be too high or that much of the coral data older than 30 ka BP has undergone diagenesis. The Cariaco Basin Δ 14 C data could be reconciled with our simulations by assuming unrealistically large variations in the dead carbon fraction of the Hulu Cave 14 C record. Instead, IntCal20 has used another approach (Bayesian adaptive learning, see Heaton et al. 2020 in this issue) to estimate past MRAs in this specific marine environment. The LSG model resolution of~380 km in the region of the Cariaco Basin is too coarse to capture its small depression (160 km long, 60 km wide) on the continental shelf off Venezuela, which was decoupled from the open Caribbean during the last glacial sea level lowstand (Peterson et al. 1991). First MRA simulations utilizing a global ocean general circulation model with enhanced resolution down to~20 km in marginal seas (FESOM2 developed by Danilov et al. 2017) show promising improvements for the Cariaco region but long-term simulations, which would be necessary for a potential application of FESOM2 within this calibration effort here, are not yet available. Figure 3 Comparison of simulated and reconstructed marine Δ 14 C for locations which contain records contributing to IntCal20 (Bard et al. 1990(Bard et al. , 1998(Bard et al. , 2004a(Bard et al. , 2004bBurr et al. 1998Burr et al. , 2004Hughen et al. 2000Hughen et al. , 2004Hughen et al. , 2006Cutler et al. 2004;Fairbanks et al. 2005;Durand et al. 2013;Heaton et al. 2013). Simulation results are for ocean climate scenario GS forced with the mean atmospheric Δ 14 C record where error bars indicate the combined uncertainty of atmospheric Δ 14 C and scenarios CS and PD. Four data points are not shown (two data points with reconstructed Δ 14 C < -150‰ and simulated Δ 14 C~20‰, plus two data points with reconstructed Δ 14 C > 1500‰ and simulated Δ 14 C~200‰).