Toward Reconciling Radiocarbon Production Rates With Carbon Cycle Changes of the Last 55,000 Years

Since it is currently not understood how changes in 14 C production rate ( Q ), and in the carbon cycle, can be combined to explain the reconstructed atmospheric Δ 14 C record, we discuss possible reasons for this knowledge gap. Reviewing the literature, we exclude that changes in the content of atoms in the atmosphere, which produce cosmogenic 14 C after being hit by galactic cosmic rays, might be responsible for parts of the observed differences. When combining Q with carbon cycle changes, one needs to understand the changes in the atmospheric 14 C inventory, which are partially counterintuitive. For example, during the Last Glacial Maximum, Δ 14 C was ∼400‰ higher compared with preindustrial times, but the 14 C inventory was 10% smaller. Some pronounced changes in atmospheric Δ 14 C do not correspond to any significant changes in the atmospheric 14 C inventory, since CO 2 was changing simultaneously. Using two conceptually different models (BICYCLE-SE and LSG-OGCM), we derive hypothetical Q s by forcing the models with identical atmospheric CO 2 and Δ 14 C data. Results are compared with the most recent data-based estimates of Q derived from cosmogenic isotopes. Millennial-scale climate change connected to the bipolar seesaw is missing in the applied models, which might explain some, but probably not all, of the apparent model-data disagreement in Q . Furthermore, Q based on either data from marine sediments or ice cores contains offsets, suggesting an interpretation deficit in the current data-based approaches.


Introduction
Radiocarbon is a valuable tracer in paleoclimatology and paleoceanography.It is used for chronological purposes, as for example, in the most recent international calibration efforts of IntCal20, Marine20, and SHCal20 (e.g., Heaton et al., 2020;Hogg et al., 2020;Reimer et al., 2020), and provides the main method to determine the age of carbon-bearing samples over the last ∼55 kyr.In addition, 14 C is one of the three carbon isotopes, along with 12 C and 13 C, used to disentangle past changes in the carbon cycle.Carbon cycle models, initially restricted to box models, have a long history of helping in the interpretation of 14 C data (e.g., Bard et al., 1994;Broecker et al., 1990).Carbon cycle models have played an integral part in international calibration efforts since the IntCal project began in 1998, originally only for the Holocene (Stuiver, Reimer, Bard, et al., 1998, Stuiver, Reimer, & Braziunas, 1998) but now, with the 2020 curve updates, for the whole time window, which is reachable with 14 C (Butzin et al., 2020;Heaton et al., 2020).Due to the coastal locations of most marine samples, as these carbon cycle models become more highly resolved spatially, our model-based estimates of the marine reservoir ages (MRA) should improve further (e.g., Lohmann et al., 2020).However, there is still one major unresolved challenge in the interpretation of the 14 C cycle.Independent evidence on 14 C production rates (Q) and model-based carbon cycle changes do not combine to successfully reproduce the atmospheric Δ 14 C record.A transient modeling approach (Köhler et al., 2006) indicated that 14 C production rates inferred from 10 Be records available at the time (Muscheler et al., 2004(Muscheler et al., , 2005) ) did not agree with those inferred from the geomagnetic field records (Laj et al., 2004).Further, neither could explain the reduction in atmospheric Δ 14 C from its high values of more than 400‰ around the Last Glacial Maximum (LGM) to lower than 100‰ at the onset of the Holocene (Figure 1b).Simulated carbon cycle changes, which were able to track reconstructed deglacial changes in atmospheric CO 2 and δ 13 C (Köhler et al., 2005) reasonably well, were also only able to explain less than a quarter (or 110‰) of that deglacial decline in atmospheric Δ 14 C. (Galbraith et al., 2015), as similarly shown already in Bard (1998).This suggests that simulations are prevented from meeting the Δ 14 C target successfully if CO 2 is not simulated reasonably well.Importantly, the different estimates of the 14 C production rates inferred from the various 10 Be and paleomagnetic records seem to be too diverse.Neither a single approach alone, nor any mean value calculated so far, was able to successfully reproduce Δ 14 C. Dinauer et al. (2020) tested five different time series for Q (one based on 10 Be in ice cores, the others on various paleomagnetic reconstructions) in the Bern3D model and all led to systematically lower glacial atmospheric Δ 14 C than the IntCal estimates.(Köhler et al., 2017) and δ 13 C of CO 2 (Eggleston et al., 2016).(b) IntCal20 Δ 14 C (Reimer et al., 2020).The calculated Δ 14 C (magenta line) is based on an atmospheric 14 C inventory, which is held constant at 60,000 mol (1σ = 3,000 mol).(c) Back calculation of the atmospheric 14 C inventory.The impact of CO 2 on 14 C inventory is tested with two additional calculations in which CO 2 is either kept constant at preindustrial (277 ppm) or LGM (182 ppm) values.(d) Residual of temporal changes in atmospheric 14 C inventory and the estimate of the 14 C decay D in the atmosphere both as 400-year running mean (rm400y).Blue vertical bands mark Heinrich stadials 1-5 (H1-H5) and the Younger Dryas (YD), similarly to Reimer et al. (2020).During these periods, dramatic changes in Atlantic Ocean circulation may have occurred.
10.1029/2021PA004314 3 of 21 We therefore believe it is time to take a different view on the 14 C cycle.We briefly review what is known on the problem before investigating potential causes for the disparity.First, we evaluate if potential changes in the relative abundance of the precursor atoms, from which cosmogenic radionuclides are produced, could have had an influence on 14 C production rates (Section 2.1).Second, we discuss the changes in the total atmospheric 14 C inventory.This provides an alternative perspective, since the total 14 C inventory is influenced by overall atmospheric CO 2 levels as well as the 14 C production rate (Section 2.2).Third, in Section 2.3, we compile and discuss the most recent data-based estimates of Q.Finally, we use two conceptually different carbon cycle models (Sections 3 and 4) to calculate the 14 C production rate by inverting reconstructed Δ 14 C under different carbon cycle scenarios.Here, we use the carbon cycle box model BICYCLE-SE (Köhler & Munhoven, 2020) and the ocean general circulation model LSG (Butzin et al., 2020).While both models can be similarly forced by the changes in the atmospheric carbon reservoirs, BICYCLE-SE can additionally simulate a carbon cycle internally, which closely matches both the atmospheric CO 2 reconstructed from ice cores and deep ocean carbonate ion concentrations from sediment cores.From comparison of our model-based 14 C production rates with the most recent independent 10 Be-and paleomagnetic-based estimates, we derive where shortcomings in the data and/or models might hinder progress.The results from BICYCLE-SE also provide insight into how mass conservation in the carbon cycle, land carbon storage, ocean circulation, or solid Earth processes might be of relevance for the 14 C cycle.

The Process of 14 C Production and Potential Changes in the Precursor Material for Radiocarbon
Cosmogenic radionuclides are produced by cosmic radiation from GCRs, solar energetic particles (SEPs), or by the outcome of more exotic sources like a nearby supernova explosion hitting atoms in Earth's stratosphere.The vast majority of newly produced 14 C results from GCRs, while SEPs are on average responsible for only about 0.25% (Kovaltsov et al., 2012).Within the energy range most important for the production of the cosmogenic radionuclides (4 ⋅ 10 2 − 8 ⋅ 10 3 MeV, Herbst et al. (2017)), the unmodulated intensity of GCRs entering the solar system is uncertain but generally assumed to have been constant over time.
In addition to the rate at which GCRs enter the solar system, the global 14 C production rate by GCRs depends on two further parameters, the solar magnetic activity, and the Earth's geomagnetic field.These are quantified via the Sun's modulation potential ϕ and Earth's dipole moment M, respectively.The reconstruction of the modulation potential after 1700 CE found 10-year averages between 100 and 800 MV (Alanko-Huotari et al., 2007).At present-day values of Earth's dipole moment M of 7.8 ⋅ 10 22 Am 2 , this would correspond to a range in Q of nearly a factor of two-between 2.7 and 1.5 atoms/(cm 2 s) (Kovaltsov et al., 2012).Similarly, the simulated change in Earth's dipole moment across the last 100 kyr was according to Panovska et al. (2019) between 2.5 ⋅ 10 22 and 9.9 ⋅ 10 22 Am 2 .For an assumed constant preindustrial ϕ of ∼400 MV, this corresponds to a range in Q between 3.4 and 1.7 atoms/(cm 2 s) (Kovaltsov et al., 2012).The solar modulation potential ϕ largely follows the 11-year cycle of sun spots, but is also connected to longer term variations (McCracken & Beer, 2007).Before the year 1953 CE, observations are missing; however, ϕ can be reconstructed from radionuclides (Muscheler et al., 2007;Steinhilber et al., 2012;Brehm et al., 2021).One of the most recent preindustrial estimates for Q, averaged over the years 1750-1900 CE, is 1.88 atoms/(cm 2 s) or 504 mol/yr; while for modern time, Q is 1.64 atoms/(cm 2 s) or 440 mol/yr (Kovaltsov et al., 2012), which was slightly revised in Poluianov et al. (2016) and Herbst et al. (2017).
The bulk (99.999%) of the 14 C production rate Q from GCR occurs from thermal neutrons hitting 14 N atoms, which then subsequently lose one proton and transform into 14 C nuclei (Kovaltsov et al., 2012).Some 4.2% of the thermal neutrons are lost due to other reactions not leading to 14 C.These other potential chain reactions also include nitrogen as well as oxygen and argon ( 14 N into 15 N; 16 O into 17 O; 18 O into 19 O; 40 Ar into 41 Ar).In the following, we discuss whether a variable amount of nitrogen in the past atmosphere might have influenced the 14 C production rates arising from GCRs.We also roughly approximate whether the atmospheric amounts of oxygen and argon have significantly varied over glacial time scales and therefore, if they might, via a change in the neutron sink, indirectly affect Q.The tiny (0.001%) contribution to Q from 17 O hit by a neutron and the negligible contribution from further spallation reactions (Kovaltsov et al., 2012) are ignored here.An overview on potential changes in atmospheric gases of interest for Q is found in Table 1.

10.1029/2021PA004314
4 of 21 Nitrogen dominates our present-day atmosphere with 78.08% of the volumetric share of dry air (Schlesinger, 1997).With a total of 1.77 × 10 20 mol for the whole dry air atmosphere (Headly & Severinghaus, 2007), the present-day atmosphere contains ∼138,200 Pmol of N 2 .In the ocean, the present-day nitrogen inventory is mainly (∼94%) found in form of N 2 (357 Pmol N 2 ≈ 10 4 PgN) (Gruber, 2008).This estimate has an uncertainty of about ±10%.The inventory of other forms of nitrogen in the ocean, which are relevant for the biological carbon pumps, is at least one order of magnitude smaller than the amount of N 2 in the ocean.Thus, for our purposes, changes in the marine biology, including nitrification and denitrification, which would consume or provide N 2 , can be ignored.The amount of nitrogen on land is 140 PgN (Batjes, 1996), two orders of magnitude smaller than the oceanic inventory, and is therefore negligible.In total, the present-day nitrogen inventory is distributed between atmosphere (∼99.7%),ocean (∼0.3%), and land (<0.001%), respectively.
Due to temperature-dependent solubility, a colder glacial ocean should store more nitrogen (Ritz et al., 2011;Weiss, 1970).The solubility of gases is related to mean ocean temperature (Headly & Severinghaus, 2007).A data-based approach estimates that during the LGM, the mean ocean temperature would have been around 3 K lower than at preindustrial times (Shakun et al., 2015).This would increase the solubility of N 2 by less than 10%.We can therefore reasonably assume that the oceanic N 2 inventory might increase by the same order of magnitude (∼36 Pmol N 2 ).In addition, the effects of a glacial increase in the ocean salinity by ∼3% and of a lower glacial sea level of about 120 m need to be considered.The first effect reduces nitrogen solubility (Ritz et al., 2011;Weiss, 1970), but by an order of magnitude less than the related temperature-based increase.The second effect leads to a larger partial pressure of N 2 at a sea level of ∼1.6% (Headly & Severinghaus, 2007) that would according to Henry's Law, linearly increase the oceanic N 2 by a similar amount (6 Pmol N 2 ).These two effects (sea level and salinity) nearly cancel each other.Altogether, atmospheric N 2 at the LGM might therefore have been smaller by 0.3‰ than at modern times.
In the present day, oxygen's and argon's volumetric share of dry atmosphere are 20.95% and 0.93%, respectively (Schlesinger, 1997).Besides the atmosphere and ocean, there are no further significant pools for argon (Ar).Following similar arguments as for N 2 , during the LGM, a colder, more saline glacial ocean with a lower sea level might have contained about 10%-15% more Ar than at present.This would have reduced atmospheric Ar-which contains, at ∼1,650 Pmol, two orders of magnitude more Ar than the 21 Pmol Ar in the ocean (Sarmiento & Gruber, 2006)-by the relative fraction of 2‰.Calculations for oxygen are more difficult.The present atmosphere contains two orders of magnitude more O 2 than the ocean (∼37,100 vs. 170 Pmol O 2 ).Again, during the LGM, the temperature-, salinity-, and pressure-related effects would increase glacial oceanic O 2 by ∼10%.However, oxygen is a fundamental part of the biological carbon cycle.Using the BICYCLE-SE carbon cycle model (Köhler & Munhoven, 2020), we estimate that the glacial ocean stored 70 Pmol less O 2 than the preindustrial ocean.This value lies roughly in the middle of two other estimates: Cliff et al. (2021) (Plattner et al., 2002).Both marine and terrestrial carbon cycle processes have large uncertainties, but more or less compensate each other.At the LGM, any relative overall change in atmospheric O 2 is less than 2‰ compared to present-day levels.Finally, ice core records suggest that the δ 18 O of atmospheric O 2 (the target of one possible sink for the GCR neutrons) at the LGM was 1‰ smaller than at present day (Severinghaus et al., 2009).This suggests that subsequent effects on the 14 C production rate Q would be small.Their precise quantification, including nonlinearities ignored here, would have to consider the whole nuclear cascade caused by GCRs with a glacial atmosphere and reduced sea level.
In summary, we find that atmospheric N 2 changed by less than 0.3‰ during the last 55 kyr.When considering the total amount of nitrogen in the atmosphere, these relative changes in N 2 on glacial/interglacial timescales have a larger effect than those in N 2 O, another N-containing gas of the atmosphere.The atmospheric concentration of N 2 O was 190 ppb at LGM and rose to 270 ppb at preindustrial times (Schilt et al., 2010).This corresponds to a relative change in total atmospheric nitrogen of less about 0.1‰.Other relative changes in the amount of the GCRs' target material over glacial timescales are smaller than 1‰.Since these might affect less than 5% of the alternative thermal neutron sinks, their overall impact on Q would be less than 0.05‰.The interactions of GCRs with the various nuclei involved in the radionuclide production rates are considered to occur mostly (about two thirds) in the stratosphere (Masarik & Beer, 1999).However, since they do not change the amount of atoms available as GCR targets, no further assumptions on specific aspects of gases in the stratosphere, such as the chemistry related to the ozone cycle, are necessary.We therefore exclude any significant impact of changes in the amount of the precursor (GCR target) material on the 14 C production rate.Furthermore, we can safely state that nearly all discussed processes would point toward a slightly smaller amount of target material at the LGM and are therefore unable to explain any of the increased atmospheric Δ 14 C found in the reconstructions.

Radiocarbon Inventory in the Atmosphere
The inventory of 14 C (in atoms or mol) can be calculated as with CO 2 in ppm and f = 176.734⋅ 10 12 mol/ppm being a constant for unit conversation.The relationship between Δ 14 C and δ 14 C is given by (Mook, 1980), and is the modern 14 C/ 12 C standard, when the half-life of 14 C of 5730 ± 40 years used in IntCal20 is considered.Note that we do not use the most recent estimate of the half-life of 14 C of 5700 ± 30 years (Audi et al., 2003) as this was not implemented in IntCal20.For more details on the choice of this standard see, for example, Trumbore et al. (2016).The denominator in Equation 2 is often considered uniform and neglected (e.g., Stuiver & Polach, 1977).We calculate the inventory of 14 C in the atmosphere from estimates of atmospheric Δ 14 C, CO 2 , and δ 13 C of CO 2 .We neglect throughout this study the 14 C contribution from methane (CH 4 ) since the atmospheric concentration of CH 4 is three orders of magnitude smaller than that of CO 2 .
For most of the last 55,000 years, the obtained 14 C inventory in the atmosphere (Figure 1c) stays within ±10% of its preindustrial value of 60 kmol.The exceptions are an overshoot between 42 and 37 kcal BP, covering the Laschamps geomagnetic excursion, and slightly lower values between 51 and 43 kcal BP.The atmospheric 14 C inventory during the LGM is up to 10% smaller than at preindustrial times.This change is opposite to the higher Δ 14 C in the LGM compared to preindustrial times.The large and abrupt drop in Δ 14 C by 40‰ in 200 years at the end of Heinrich stadial 1 is clearly associated with the simultaneous rise in atmospheric CO 2 at the onset of the Bølling/Allerød around 14.7 kcal BP (Marcott et al., 2014)-we see hardly any changes in the atmospheric 14 C inventory.However, at the same time, the high-resolution 10 Be record (Adolphi et al., 2017) indicates a drop in Q (albeit at the end of a long-term rise, Figure 2c).Both paleo reconstructions and climate models indicate an overshoot in the Atlantic meridional overturning circulation (AMOC) at the end of Heinrich stadial 1 (Barker et al., 2010;Skinner et al., 2021), pointing to massive changes in the oceanic carbon cycle and suggesting a rapid transfer of old 14 C-depleted carbon from the deep ocean to the surface ocean and the atmosphere (Rae et al., 2018).However, some studies also point toward a contribution to this prominent CO 2 peak from terrestrial sources  et al., 2018, 2019) for three different values of ϕ from ice core ( 10 Be, 36 Cl) data (Adolphi et al., 2018), normalized to the last 2 kyr, and from sediment core ( 10 Be) data (Simon et al., 2020), normalized to a dipole moment with solar modulation potential ϕ = 400 MV for comparison.(d) Related dipole moments M for two of the approaches.Shaded areas, whenever available, give 1σ uncertainties as denoted in the original studies.Some highly fluctuating records are shown as 400-year running mean (rm400y).Blue vertical bands mark Heinrich stadials 1-5 (H1-H5) and the YD, similarly to Reimer et al. (2020).During these periods, dramatic changes in Atlantic Ocean circulation may have occurred.(Bauska et al., 2016;Winterfeld et al., 2018).The combination of all information indicates that at this point in time, all the processes influencing Δ 14 C were changing at the same time.The impact of δ 13 C on the derived 14 C inventory is negligibly small.
We now use this information of a nearly constant atmospheric 14 C inventory to provide an alternative way to illustrate the importance of CO 2 on atmospheric Δ 14 C. Here, we calculate a hypothetical Δ 14 C for an atmospheric 14 C inventory that stayed constant at its preindustrial value, but with a relative σ of 5%.This value of σ is consistent with a 95% confidence interval that corresponds to the ±10% range around the mean (Figure 1b).Both IntCal20 and the hypothesized Δ 14 C agree within their uncertainties after 37 kcal BP.However, some structure within Δ 14 C can now be more clearly associated with CO 2 .For example, while simulations by Muscheler et al. (2004) have shown that the fine structure in atmospheric Δ 14 C during the Holocene can be explained by the changes in Q, they did not consider the gradual rise in CO 2 after 8 kcal BP.This gradual rise might help to improve the match of simulated Δ 14 C, based on independent estimates of Q, with Δ 14 C estimated directly from observations, as also discussed in Roth and Joos (2013).
This perspective on the 14 C cycle clearly shows that models which poorly represent or even neglect changes in atmospheric CO 2 are of limited applicability.For models that consider transient changes in the carbon cycle, it is of course also important to simulate CO 2 changes for the right reasons, that is, with carbon fluxes having the correct isotopic signatures.If CO 2 is fixed at either LGM or preindustrial levels, the atmospheric 14 C inventory (based on IntCal20's estimate of Δ 14 C) no longer stays within a ±10% range around its preindustrial value, but decreases by ∼30% over the last 25 kyr (Figure 1c).This implies that carbon cycle simulations that fail to reconstruct CO 2 would-if forced with independent (data-based) Q-contain a very different atmospheric 14 C inventory from those scenarios that meet CO 2 data closely.On the other hand, from 42 kcal BP until the onset of the Heinrich stadial 1 at 17.5 kcal BP, the atmospheric 14 C inventory inferred from Δ 14 C when using a constant glacial CO 2 level of 182 ppm stays within a ±10% range of the preindustrial value.This suggests that when applied in glacial modes (including reduced atmospheric CO 2 ), models might be useful for understanding the glacial 14 C cycle even if they do not represent Holocene CO 2 changes well.When models are used to inversely calculate Q from atmospheric Δ 14 C, a rough fit of the model-internal CO 2 to data seems necessary, but the exact value is only of minor importance.These considerations are discussed further with model output below.
The time series of atmospheric carbon can be investigated one step further.Calculating the time derivative of Equation 1gives Neglecting temporal changes in 13 C, the changes in δ 14 C are similar to those in Δ 14 C.The residual change in the atmospheric 14 C inventory, identically to this time derivative of Equation 4, is therefore directly available from the data, as is the loss due to radioactive decay (Figure 1d).Both are intrinsic features of 14 C cycle and on the order of 10 mol/yr.Interestingly, the two unknowns of Equation 5, the 14 C production rate Q (see Section 2.3 and Figure 2) and the carbon cycle sinks are a factor of 50 larger than the residual and the loss due to decay.This tells us that we only have precise information about the difference of two large fluxes, both changing the atmospheric 14 C inventory.Small uncertainties in one or the other of these two fluxes, when independently estimated, might quickly lead to problems in successfully closing the 14 C mass balance.

Previous Estimates of Changes in the 14 C Production Rate
Different approaches exist for the estimate of changes in the 14 C production rate over time.We briefly discuss the most important ones in this subsection.Variations in other cosmogenic nuclides, most prominently 10 Be, but also 36 Cl, indicate how the GCR flux reaching the Earth might have varied over time.However, specific differences between 10 Be and 14 C need further attention.While 10 Be is produced by spallation reactions, thermal neutron capture is responsible for the production of 14 C (e.g., Herbst et al., 2017).Furthermore, climate impacts on 10 Be deposition, which would not have any imprint on 14 C, need to be considered (Heikkilä et al., 2013).
In 30-year-long deglacial simulations, 10 Be is dominated by the 11-year solar cycle (production-driven), while higher frequencies are climate-driven (Heikkilä et al., 2014).How glacial climate and multi-millennial climate variability influences 10 Be is not yet completely understood.By first-order principles, one can expect transport and washout effects analogous to other aerosols deposited on the polar ice sheets (e.g., Alley et al., 1995;Schüpbach et al., 2018), albeit likely of smaller magnitude since 10 Be originates above the cloud layer and is thus only affected by sub-and in-cloud scavenging for a fraction of its transport path.Still, any uncertainties in past precipitation rates over the ice sheet (Gkinis et al., 2014;Rasmussen et al., 2013) will directly impact our ability to reconstruct the past cosmogenic radionuclide production rates from ice core data.Further climate influences on 10 Be mixing and deposition might possibly originate from changes in the tropopause height (Rind et al., 2001) and the Brewer-Dobson circulation (Fu et al., 2020).
One of the most recent estimates of relative changes in Q was given in Adolphi et al. (2018).This was based on a stack of different 10 Be and 36 Cl records from ice cores, extended to modern time with data from Muscheler et al. (2016).After the usage of this ice core-based estimate of Q in Dinauer et al. ( 2020), we discovered that the inclusion of 36 Cl samples that had likely been contaminated caused a bias in the stack, affecting parts of the Holocene.Therefore, in the version of the stack that is used here, the Holocene 36 Cl data have been removed.
The recently published NEEM 10 Be data set (Zheng et al., 2021) is not included here.However, this omission does not affect our conclusions, since the NEEM data are largely consistent with the other records in the stack.This ice core-based approach gives only a relative change in Q.For the calculation of absolute numbers, we use a Q of 440 mol/yr at 0 cal BP.This value is also used in the BICYCLE-SE model and leads to stable atmospheric Δ 14 C in model validation experiments run with constant climatic boundary conditions.In order to avoid a dependency of the whole time series of ice core-based Q on the data point at 0 cal BP, our final time series is obtained by first dividing the data by its mean values calculated for the last 2 kyr.Second, data scatter is reduced by using a 400-year running average (Figure 2c).We note that this approach to estimating Q contains some empirical climate corrections, since 10 Be and 36 Cl fluxes have been compared to other climate proxies (e.g., δ 18 O or aerosols) with the common signal being removed by linear regression.We plot our ice core-based estimate for Q on the GICC05 age scale (Andersen et al., 2006;Svensson et al., 2006), even though the link to U/Th-dated speleothems found age offsets to GICC05 of up to 500 years (Adolphi et al., 2018).We here refrain from an age correction (which would affect the ice core-based Q through its influence on the ice accumulation rates used for flux calculation).Applying this correction would lead to discontinuities in the 10 Be fluxes, since the age correction is based on very few tie points in the glacial.Furthermore, the effect would be less than +2% for the period older than 22 kyr-where the largest discrepancies between ice core-based Q and Δ 14 C occur.
Alternatively, the dipole moment M of Earth's magnetic field itself has been used to reconstruct Q.In Figure 2c, we plot the most recent approach based on M from Panovska et al. (2018Panovska et al. ( , 2019) ) using dependencies of Q from ϕ, M, and the local interstellar spectrum (Herbst et al., 2017;Poluianov et al., 2016).We show the estimates for three different but constant values of ϕ, hence neglecting solar variations.Differences between the Q based on cosmogenic radionuclides from ice cores and the M-based estimates of Q are either caused by unconsidered temporal changes in the GCR flux, long-term solar activity changes, or are 10 Be-specific.They might furthermore point to biases in the M-based approaches, since reconstructions of M differ widely dependent upon the records on which they are based (e.g., Dinauer et al., 2020;Panovska et al., 2019).Recently, 10 Be from marine sediment records has been used to reproduce M across the Laschamps geomagnetic excursion (Simon et al., 2020).The dipole moments from Panovska et al. (2019) and Simon et al. (2020) generally agree between 60 and 55 kcal BP and 45-20 kcal BP, but differ considerably at the local maximum in M in-between (Figure 2d).This suggests that the estimates of Q based on M alone might be biased.Such approaches are therefore not followed any further here.
An alternative approach to estimating Q can be provided by carbon cycle models.These models can be used to inversely calculate the Q needed to generate simulated Δ 14 C time series that agree with the data-based reconstructions (Muscheler et al., 2005;Stuiver & Braziunas, 1993).This model-based approach has also been applied using the Bern3D model for the Holocene (Roth & Joos, 2013) and has recently been extended to the last 50 kyr (Dinauer et al., 2020).Dinauer et al. (2020) simulated CO 2 that was close to the data-based CO 2 reconstructions at glacial times, but contained offsets of up to 50 ppm during the deglaciation.This CO 2 offset affects the modelbased Q as we demonstrate in simulations further below.

Methods: The Models BICYCLE-SE and LSG-OGCM
We will discuss model-based reconstructions of Q using either the carbon cycle box model BICYCLE-SE (Köhler & Munhoven, 2020) or the ocean general circulation model LSG (Butzin et al., 2020).Both models can be forced similarly-with data-reconstructed changes in atmospheric CO 2 (Köhler et al., 2017) and in atmospheric Δ 14 C (Reimer et al., 2020).An overview on all simulation scenarios analyzed is given in Table 2.The table includes information on which of the three variables (atmospheric Δ 14 C, CO 2 , and Q) is prescribed or whether they are internally calculated.Note.Either atmospheric Δ 14 C or the 14 C production rate (Q) is prescribed.The other variable is then derived from the model simulations.Additionally, atmospheric CO 2 is either calculated internal within the model or prescribed from ice core data.

Table 2 Simulation Scenarios Investigated in This Paper
the model, the carbon content and the isotopic signatures (δ 13 C, δ 14 C) in all boxes are followed.In the ocean, the marine carbonate system is explicitly calculated, tracing dissolved inorganic carbon and alkalinity from which the other variables of the carbon cycle (dissolved CO 2 , bicarbonate ion, carbonate ion, and pH) are calculated as a function of temperature, salinity, and pressure.Marine export production is a function of the macronutrients  PO 3− 4 at the surface.For preindustrial times, export production is restricted to 10 PgC/yr in agreement with more complex ocean general circulation models that have been validated with field data (e.g., Sarmiento et al., 2002;Schlitzer, 2002).Especially in the surface Southern Ocean  PO 3− 4 is then partly unused, which is conceptually understood due to a lack of the micronutrient iron.All  PO 3− 4 in the Southern Ocean surface box might be used for export production if the prescribed dust fluxes are above its value at 18 kyr BP, indicating sufficient glacial iron supply.Marine biology is then assumed to be no more limited by iron and global export production might be as high as 13 PgC/yr.
All temporal changes of climate are prescribed externally.In previous versions of BICYCLE (e.g., Köhler et al., 2006), the carbonate compensation was mimicked by a response function to variability in deep ocean carbonate ion.Recently, solid Earth processes (a process-based sediment module for exchange fluxes of CaCO 3 between deep ocean and sediment, volcanic CO 2 outgassing, silicate and carbonate weathering as riverine input of bicarbonate, and coral reefs as another shallow water sink of CaCO 3 ) have been implemented (version SE used here).A more complete model description of this version, and of the forcing, is found in Köhler and Munhoven (2020).In all applications considered here, the model runs transiently over two glacial cycles (the last 210 kyr) in order that the initial conditions of the sediments have a minimal effect on the simulated results.Changing boundary conditions (e.g., temperature, sea level, sea ice, ocean circulation, and terrestrial and marine biology) and the process-based sediment-ocean exchange, volcanic CO 2 outgassing, weathering fluxes, and growth of coral reefs are considered as implemented in the most recent model extension BICYCLE-SE (Köhler & Munhoven, 2020).The model can alternatively be forced to use prescribed CO 2 , which then overwrites model internal CO 2 .However, since prescribing CO 2 violates mass conservation of carbon, we mainly discuss scenarios of BICYCLE-SE in which model-internal CO 2 is used.Here, mass conservation is guaranteed, all carbon fluxes are model consistent, and the ice core CO 2 is met reasonably well (Figure 3a).
When prescribing atmospheric Δ 14 C, Q is calculated for every time step in BICYCLE-SE (using Equation 2 for conversion) by The subscripts denote atmosphere a, ocean o, terrestrial biosphere b, volcanism v, and weathering w.C a describes the atmospheric carbon pool, f a carbon flux with the subscript denoting the related process, 14 C a the atmospheric 14 C inventory, and λ the 14 C decay rate.Δ(δ 14 C a ) is the difference between the prescribed data and the model-internal calculation.The fluxes between atmosphere and ocean, and between atmosphere and terrestrial biosphere, consist of several subfluxes depending on the number of relevant boxes involved (one atmospheric box, five surface ocean boxes, and seven terrestrial boxes).
While analyzing results from BICYCLE in Heaton et al. (2020), it was stated that the numerics which are used to solve the ordinary differential equations prevented us from calculating Q as in Equation 6.This statement was incorrect and erroneously given because during the analysis of the model results, we picked the wrong variable.However, no other conclusions from Heaton et al. (2020) need to be revised.

The LSG Ocean General Circulation Model
LSG, the Hamburg Large Scale Geostrophic ocean general circulation model (Maier-Reimer et al., 1993), is applied here with a horizontal resolution of 3.5° and vertically 22 unevenly spaced levels.Model improvements contain a bottom boundary layer scheme (Lohmann, 1998) and a sophisticated numerical advection scheme (Schäfer-Neth & Paul, 2001;Prange et al., 2003).LSG-OGCM has been used for 14 C simulations before (Butzin et al., 2005(Butzin et al., , 2012(Butzin et al., , 2017(Butzin et al., , 2020) ) and the results discussed here are based on simulations performed for Heaton et al. (2020).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 & Lorenz, 2000;Prange et al., 2004).Millennial-scale climate variability connected with the bipolar seesaw that involves transient changes in AMOC are not implemented in any of the discussed simulations.They have been addressed in sensitivity studies for these two models before.They led to an increased atmospheric Δ 14 C during times of reduced AMOC (Butzin et al., 2005;Köhler et al., 2006), similarly to the response of reduced ocean ventilation during the Younger Dryas in different C cycle models (e.g., Meissner, 2007;Muscheler et al., 2000Muscheler et al., , 2008;;Ritz et al., 2008;Singarayer et al., 2008).We note however that Matsumoto and Yokoyama ( 2013) obtained a decrease in atmospheric Δ 14 C to a reduction in AMOC, an opposite response to that seen in other models.It can therefore be summarized that the response of the carbon cycle to the bipolar seesaw is highly model dependent (Gottschalk et al., 2019).
The model-based Qs drop to small values around 53 kcal BP in both models in scenario PD (LSG-OGCM) even to nonphysical, negative values (Figure 2a).These negative values indicate that the climate state at that time was different from the interglacial-like PD scenario, but may also suggest that the reconstruction of Δ 14 C beyond 50 kcal BP needs further revisions.Furthermore, with the current knowledge of M, very small values of Q at 53 kcal BP are only possible if one invokes rather massive changes in the local interstellar spectrum.Since this would also be seen in 10 Be, and available data give no indication of such 10 Be changes, we manually adjusted the values of Q before 50 kcal BP in scenario V2 to avoid such small numbers.Thus, in V2, Q now stays constant at 660 mol/yr before 70 kcal BP and then decreases linearly to the values found in scenario V1 at 50 kcal BP.After 50 kcal BP, both scenarios V1 and V2 are identical in Q.Before 50 kcal BP, the simulated Δ 14 C in scenario V2 stays within the large 95% confidence interval of IntCal20, but does not directly meet its mean values.
Since the 14 C cycle has a long memory, production rates Q inferred by BICYCLE-SE are compared for the last 75 kcal BP.After 20,000 simulated years, results based on different initial conditions are reasonably similar (Köhler et al., 2006).Note that glacial/interglacial changes as contained in BICYCLE-SE not only meet atmospheric CO 2 , but also in general terms, agree with deep ocean reconstructions: the LGM-to-preindustrial change in a global mean deep ocean reservoir age of 644 ± 101 14 C years in the model is in agreement with 689 ± 53 14 C years derived from 256 deep ocean 14 C samples (Skinner et al., 2017); and basin-wide deep ocean  CO 2− 3 changes agree with reconstructions (Köhler & Munhoven, 2020;Yu et al., 2013).
The Q calculated by both BICYCLE-SE and LSG-OGCM share the same fine structure with results from LSG-OGCM being 50-100 mol/yr smaller than those from BICYCLE-SE (Figure 2b).In LSG-OGCM, the two different glacial-like climate states (GS and CS) lead to nearly identical results, while the interglacial scenario (PD) generally leads to smaller Q.These lower Qs in the PD scenario correspond to a higher surface ocean Δ 14 C.This is caused by higher global-average gas exchange rates in PD, partly related to reduced sea ice cover (Butzin et al., 2005;Heaton et al., 2020).For a better understanding of the difference between both models, we perform a further experiment in which the terrestrial carbon cycle in BICYCLE-SE is switched off and CO 2 is prescribed (V1CO2-TB).Such a scenario should be suited better for a comparison to results obtained with the "ocean-only" LSG-OGCM.The difference in Q between this no-carbon-cycle scenario to V1CO2 (a scenario with prescribed CO 2 and Δ 14 C) is around 50 mol/yr or up to half of the overall difference between BICYCLE-SE and LSG-OGCM (Figure 2b).Terrestrial carbon is quite young, for example, the difference of mean Δ 14 C in land carbon to the atmosphere in BICYCLE-SE is 10-50‰ most of the time.However, whether the additional load of 1,500-2,200 PgC in the land carbon cycle with its small 14 C depletion is considered or ignored has an effect on Q.The rest of the model differences is probably caused by higher, on average, gas exchange rates in the nonpolar regions of LSG-OGCM compared to BICYCLE-SE (Heaton et al., 2020).
When comparing a model simulation that internally calculate atmospheric CO 2 (scenario V1), against a simulation where CO 2 is prescribed (scenario V1CO2), one learns about the effect of violating mass conservation (Figure 3).Prescribing CO 2 implies in our setups that positive or negative artificial CO 2 fluxes with the same 14 C signature as the atmosphere are generated in order to keep the internally simulated CO 2 in agreement with the observation.These artificial fluxes therefore generate or destroy 14 C and have an impact on the 14 C inventory.Internally calculated and prescribed simulations differ in CO 2 especially during abrupt changes, for example, within Termination I, and consequentially also in the 14 C inventory.For example, around the beginning of Heinrich stadial 1, the rise in CO 2 is offset between simulations by about 1 kyr, as is the rise in atmospheric 14 C inventory (Figure 3a).The contribution to Q based solely on the forced (artificial) carbon fluxes needed to keep CO 2 in agreement with the data (as in scenario V1CO2) is surprisingly small, but not zero (<50 mol/yr, blue line in Figure 3b).This implies that although one has to be careful when calculating Q from Δ 14 C and prescribed CO 2 , its relevance for the overall pattern of Q is relatively small.Even when CO 2 is not altered (scenario PRECO2), while the deduced changes in Q are nonnegligible, they remain small (Figure 3b).However, the resulting atmospheric 14 C inventories differ from inventories based on reasonable changes in the carbon cycle (Figure 3c).This illustrates that while setups with a preindustrial CO 2 might indeed be used for a model-based calculation of Q, they lead to a completely different 14 C cycle whose interpretation, especially in the deep ocean, might be difficult.
The model-based Qs contain some millennial-scale variability, which is contained in IntCal20, but not seen in the sediment-based (scenario SED) version of Q, and only present to a limited extent in the ice core-based (scenario ICE) version (Figure 4b).This suggests that these changes are mainly related to the carbon cycle and not to Q.We therefore additionally use a smoothed version of the Hulu Cave Δ 14 C (Cheng et al., 2018;Southon et al., 2012).This Hulu Δ 14 C curve has previously been employed by LSG-OGCM when calculating the MRA for the marine records contributing to IntCal20 (Butzin et al., 2020;Heaton et al., 2020) (Figure 4b).The smoothed Hulu Cave Δ 14 C curve agrees with IntCal20 within their 95% confidence intervals for all times with the exception of a few wiggles between Heinrich stadials 2 and 1.The difference between the model-based Qs obtained from the IntCal20 and Hulu Cave Δ 14 C records (HULU-V2) indicates that the millennial-scale variability within IntCal20 might be related to the climate change connected with the bipolar seesaw (Figure 4c).In detail, differences seemed to be not related specifically to Heinrich stadials, but rather more generally to all Greenland stadials.Changes connected with the Heinrich events are probably long and large enough that they are also included in the smoothed Δ 14 C version of Hulu Cave.According to 231 Pa/ 230 Th data from the Bermuda rise (Böhm et al., 2015;Henry et al., 2016;Lippold et al., 2019;McManus et al., 2004), all Greenland stadials of the last 50 kyr, no matter whether they contain a Heinrich event or not, are connected to abrupt changes in the AMOC (Figures 4d and 4e).If this is correct, the lack of abrupt ocean circulation changes in the transient simulations performed with LSG-OGCM and BICYCLE-SE introduces several biases.First, during times of AMOC shutdown, the simulated LSG-OGCM MRAs for the marine records used in the construction of IntCal20 would be too large (Butzin et al., 2012).This has the potential to alter, in periods where marine data provide a substantial contribution to the curve, the millennial-scale variability contained in IntCal20.We note that MRA ages for Cariaco basin are unaffected since they are not based on LSG-OGCM (Heaton et al., 2020).We also note that the unsmoothed Hulu Cave Δ 14 C (Cheng et al., 2018) already contains the millennial-scale variability.To estimate the size of this AMOC shutdown effect, one would need to build a revised IntCal20 without these marine records.Second, depending on the effect of AMOC shutdown on atmospheric Δ 14 C-either an increase as in Butzin et al. (2005) and Köhler et al. (2006) or a decrease as in Matsumoto and Yokoyama (2013)-the model-based Q would differ from the values deduced here that have not incorporated such AMOC changes.This second bias would also be present during the times of AMOC shutdown, where the Hulu Cave and IntCal20 Δ 14 C agree with each other, for example, the Heinrich stadials.Third, the Marine20 curve (Heaton et al., 2020) may also miss the related variabilities.This would occur since Marine20, an estimate of the mean nonpolar surface ocean MRA, is created using IntCal20's atmospheric Δ 14 C as a forcing to an older version of the BICYCLE model without abrupt climate changes connected to the bipolar seesaw.
The most robust feature in Q occurs where the reconstructions based on ice core and on sediment core data agree.Together, they point to smaller values across the Laschamps geomagnetic excursion than the model-based estimates.This suggests that models which depict the right processes at the right time (including AMOC shutdowns during the Greenland stadials 10 and 11, which are dated to be roughly synchronous with the Laschamps geomagnetic excursion, Figure 4b) should simulate higher atmospheric Δ 14 C for an AMOC shutdown.Such models should subsequently calculate a lower model-based Q that lies in greater agreement with the data-based approaches here.Recently available atmospheric Δ 14 C from New Zealand kauri trees across the Greenland stadial 11 (Cooper et al., 2021), including the highest pre-bomb Δ 14 C values known at >800‰, challenges our understanding even further.There is no indication, from ice core and sediment data, that this Δ 14 C peak has been caused by changes in Q alone.
Due to the nature of the archive, the marine sediment cores contain a smoother version of 10 Be than the ice cores.However, there is nevertheless an unexplained offset between the Qs based on the two 10 Be data sets.The ice core-based Q has consistently lower values than the sediment core-based version (Figure 4b).The current ice core estimate suggests similar values of Q during the LGM as those in preindustrial times (this is not the case for earlier versions of ice core-based Q (Köhler et al., 2006)).Such a lack of difference between the preindustrial and LGM implies that all deglacial changes in Δ 14 C need to be explained by the carbon cycle.At the same time,  (Reimer et al., 2020) and an added extension from new kauri-based data around 42 kcal BP (Cooper et al., 2021).We also show a smoothed Δ 14 C estimate obtained from Hulu Cave alone (Cheng et al., 2018;Southon et al., 2012).This Hulu estimate was used previously by the LSG-OGCM to calculate prior estimates for the MRA of those marine records incorporated in the construction of IntCal20.it illustrates a fundamental uncertainty of ice core 10 Be-flux records: any uncertainty in ice-accumulation rates is directly transferred to 10 Be fluxes.In particular, on long timescales (glacial/interglacial), different methods to reconstruct accumulation rates lead to differences of up to 30% (2σ) (Rasmussen et al., 2013).We note that marine 10 Be/ 9 Be records are also not free of slowly evolving biases, of up to 15%, due to changes in climate (Savranskaia et al., 2021).
Times with abrupt changes in CO 2 during the Heinrich stadial 1-Bølling/Allerød-Younger Dryas-Early Holocene sequence of events in Termination I need special attention.Here, model-based Q and data-based Q not only disagree, but also show the opposite behavior (Figure 5b).We hypothesize that this disagreement is caused by omitted processes in the models related to dynamic changes in the AMOC as discussed before and potentially also by old carbon mobilized from permafrost thaw (Meyer et al., 2019;Winterfeld et al., 2018).Note that there are also unresolved differences between different Greenland ice core 10 Be records during the Younger Dryas (Adolphi et al., 2014).
Although this detailed comparison of the different Qs already highlights some shortcomings in the applied carbon cycle models, we briefly discuss potential implications for the simulated 14 C cycle.We input the data-based estimates for Q into BICYCLE-SE and calculate the resultant simulated Δ 14 C (Figure 5).To allow a comparison, the independent data-based estimates for Q obtained from sediment cores (scenario SED) and ice cores (scenario ICE) are linearly extended from 60 kcal BP to give a Q at 70 kcal BP that is 50% larger than at preindustrial times.This matches the approach taken to Q in scenario V2 (Figure 5b).In doing so, for the ICE scenario, we overwrite the very small values of Q before 60 kcal BP (Figure 2b), which would otherwise due to the long-term memory effect in 14 C, prevent the model output from coming even close to IntCal20's Δ 14 C prior to the Laschamps geomagnetic excursion (∼42 kcal BP).In simulations with Q based on sediment data (SED), the resultant atmospheric Δ 14 C is in agreement with IntCal20 until ∼43 kcal BP, but differs widely thereafter.Between Heinrich stadial 5 (∼47 kcal BP) and the early Holocene (9 kcal BP), simulated Δ 14 C based on Q from ice cores (ICE) is clearly smaller than the reconstructed IntCal20 Δ 14 C.If the climate and the carbon cycle are held constant at preindustrial levels, the simulated Δ 14 C in the glacial is 50-100‰ smaller than when CO 2 roughly follows the data (as shown by scenarios ICEPRE vs. ICE, Figure 5c).The difference in simulated atmospheric Δ 14 C between the ice core-and sediment core-based Q is as large as that between the sediment core-based Q and IntCal20's Δ 14 C.Even with the shortcomings in the applied model in mind, this points to some missing processes in the independent reconstructions of Q.However, it is not straightforward to envision a Q that leads to a perfect atmospheric Δ 14 C and which makes sense with the available 10 Be data.

Discussing the Wider 14 C Cycle Within BICYCLE-SE
The contributions of different physical and biological processes to changes in atmospheric Δ 14 C in an older version of the BICYCLE model have been analyzed in detail before (Köhler et al., 2006).Briefly, it was found that for a constant preindustrial Q, individual processes contributed less than 30‰ each to changes in atmospheric Δ 14 C. Reduced glacial Southern Ocean vertical mixing, reduced glacial gas exchange due to sea ice cover, and reduced glacial AMOC contributed most to the enhanced glacial values; while a reduced glacial land carbon storage led to smaller atmospheric Δ 14 C at LGM.Here, we have applied a revised model version, BICYCLE-SE, that also contains solid Earth processes (Köhler & Munhoven, 2020) whose contributions to Δ 14 C need to be identified.However, a complete switch off of any of the solid Earth processes (volcanic CO 2 outgassing, weathering, or sediment-ocean interaction) quickly leads to runaway trends in the carbon cycle simulated with BICY-CLE-SE.We therefore analyze the contributions of these processes by manipulating the 14 C signature of these fluxes and then calculate anomalies, from the control setup, obtained with these manipulated scenarios: 1. Volcanic CO 2 outgassing-Volcanic CO 2 outgassing in BICYCLE-SE is a function of sea level change and varies between 7 and 15 Tmol/yr.This would also include 14 C-free CO 2 outgassing from hydrothermal activities at mid ocean ridges and island volcanoes as suggested elsewhere (e.g., Hasenclever et al., 2017;Ronge et al., 2016;Stott et al., 2019).In scenario V14, we assume that the volcanic CO 2 is not 14 C-free, but has the same Δ 14 C signature as the atmosphere.2. Weathering-In the control run (V1), carbonate weathering is held constant at 12 Tmol/yr from which 50% of the carbon entering the surface ocean as bicarbonate ion is of 14 C-free rock origin.In scenario W14, none of the carbon contributing to carbonate weathering is considered to be 14 C-free.3. Sediment-ocean interaction-Depending on the deep ocean carbonate ion concentration, CaCO 3 in the sedimentary mixed layer might dissolve.This would lead to DIC and alkalinity fluxes into the deep ocean with carbon isotopic signatures as tracked in the sedimentary mixed layer that have a Δ 14 C as low as −630‰.In scenario M14, any dissolved CaCO 3 is assumed to contain the same Δ 14 C as the ocean box through which these fluxes enter the deep ocean.
The contributions of these solid Earth processes to atmospheric Δ 14 C are as follows: We find with less than −5‰ a marginal contribution of CaCO 3 dissolution (M14) to atmospheric Δ 14 C; about −20‰ from carbonate weathering (W14); and about −35‰ from 14 C-free volcanic emissions (V14).All these contributions show only small trends over time (Figure 6).
Additionally, it has been suggested that the release of about 700 PgC from permafrost over the last deglaciation might have contributed significantly to the deglacial decline in atmospheric Δ 14 C (Ciais et al., 2012;Winterfeld et al., 2018).We have also calculated, in an additional scenario P14, the maximum contribution of a 14 C-free release from such inert carbon.Here, this permafrost carbon release is assumed to be constant over 10 kyr (from 17 to 7 kcal BP) and accompanied by a similar large additional uptake of carbon by the terrestrial carbon cycle.Thus, we are able to isolate, in scenario P14, the pure effect of the permafrost carbon on atmospheric Δ 14 C without additional perturbation of the carbon cycle.We find a contribution of less than −15‰ to the deglacial decline in atmospheric Δ 14 C (Figure 6).
Altogether, these additional processes focusing on the solid Earth and on permafrost, which have been previously neglected, contribute little to changes in the 14 C cycle.In particular, they do not help in explaining the reconstructed changes in Δ 14 C seen in IntCal20.Although different in detail from our transient simulations, we find no evidence for an important contribution of sediments to the 14 C cycle as suggested by Dinauer et al. (2020).

Conclusions
We investigated various avenues toward answering why, when used in carbon cycle models, the independent estimates of the 14 C production rate are still unable to provide simulated atmospheric Δ 14 C time series, which are in agreement with data-based reconstructions.We find no evidence for significant changes in the inventories of precursor material, which when hit by GCRs in the atmosphere, produces 14 C. We calculate, from independent data-based estimates of Q, the atmospheric 14 C inventory.This informs when changes in CO 2 dominated changes in Δ 14 C and highlights the importance of model-internal CO 2 for any simulated 14 C cycle.Models that poorly represent, or even neglect, changes in atmospheric CO 2 have a restricted applicability.When based solely on the dipole moment of Earth's magnetic field, data-based reconstructions of Q seem to be biased.Furthermore, estimates of Q from cosmogenic isotopes in ice cores partly disagree with those from 10 Be in sediment cores.This leads to significantly different simulated atmospheric Δ 14 C and points to some lack of process understanding in the data-based Qs.
Two conceptually different carbon cycle models calculate Q based on atmospheric Δ 14 C from IntCal20.These model-based estimates show similar patterns.However, in the application of both models, millennial-scale climate variability connected to the bipolar seesaw and abrupt changes in the AMOC are ignored.The same models, in similar setups as used here, were used in the construction of IntCal20.Thus, these missing AMOC changes might have introduced a small bias in the marine reservoir age for marine records contributing to IntCal20.This lack of abrupt AMOC changes in the applied models helps somewhat to explain why the model-based estimates of Q are larger than the independent data-based reconstructions for Q.However, the new maximum in atmospheric Δ 14 C at 42 kcal BP, seen in recently discovered New Zealand kauri trees, poses a challenge for our understanding.This peak seems to be purely related to, as yet unexplained, carbon cycle changes.
The abrupt climate change (the bipolar seesaw) needs to be included in future modeling applications that cover the last 55 kyr.This is challenging due to the current model dependencies in the response of the carbon cycle to such changes.Here, improvements in our understanding of the end of the Heinrich 1 stadial seem to be particularly crucial, since all major players in the 14 C cycle (CO 2 , 10 Be, AMOC, and atmospheric Δ 14 C) changed at the same time.Our investigations were focused on atmospheric 14 C.This is prone to uncertainties since its change is mainly based on the difference of two large, but insufficiently constrained, fluxes-the source Q and the carbon cycle-related sinks.In future applications, the inclusion of marine 14 C data, and of more highly spatially -resolved models, should help to improve our understanding of the radiocarbon cycle further.

Figure 1 .
Figure 1.Atmospheric carbon time series.(a) CO 2(Köhler et al., 2017) and δ 13 C of CO 2(Eggleston et al., 2016).(b) IntCal20 Δ 14 C(Reimer et al., 2020).The calculated Δ 14 C (magenta line) is based on an atmospheric 14 C inventory, which is held constant at 60,000 mol (1σ = 3,000 mol).(c) Back calculation of the atmospheric 14 C inventory.The impact of CO 2 on 14 C inventory is tested with two additional calculations in which CO 2 is either kept constant at preindustrial (277 ppm) or LGM (182 ppm) values.(d) Residual of temporal changes in atmospheric 14 C inventory and the estimate of the 14 C decay D in the atmosphere both as 400-year running mean (rm400y).Blue vertical bands mark Heinrich stadials 1-5 (H1-H5) and the Younger Dryas (YD), similarly toReimer et al. (2020).During these periods, dramatic changes in Atlantic Ocean circulation may have occurred.

Figure 2 .
Figure 2. Comparing different 14 C production rates (Q) (a-c) and related dipole moments M (d).(a) Model-based Qs from BICYCLE-SE and LSG-OGCM.(b) Differences in Q between scenario V2 and other model setups in BICYCLE-SE.(c) Data-based Qs calculated from M(Panovska et al., 2018(Panovska et al., , 2019) ) for three different values of ϕ from ice core ( 10 Be, 36 Cl) data(Adolphi et al., 2018), normalized to the last 2 kyr, and from sediment core ( 10 Be) data(Simon et al., 2020), normalized to a dipole moment with solar modulation potential ϕ = 400 MV for comparison.(d) Related dipole moments M for two of the approaches.Shaded areas, whenever available, give 1σ uncertainties as denoted in the original studies.Some highly fluctuating records are shown as 400-year running mean (rm400y).Blue vertical bands mark Heinrich stadials 1-5 (H1-H5) and the YD, similarly toReimer et al. (2020).During these periods, dramatic changes in Atlantic Ocean circulation may have occurred.
Three climate forcing scenarios are distinguished.Present-day climate background conditions are contained in scenario PD.Glacial conditions (scenario GS) represent the LGM and feature a shallower AMOC weakened by about 30% compared to PD.An alternative glacial climate setup (scenario CS) contains a further AMOC weakening by about 60% potentially representing stadials.No temporal changes of climate within a scenario are considered.See Butzin et al. (2005) for further details.Radiocarbon is simulated by the 14 C enrichment of the ocean relative to the contemporaneous atmosphere by directly simulating their fractionation-corrected 14 C/ 12 C ratio neglecting biological effects.Oceanic uptake of 14 C follows Sweeney et al. (2007) and uses atmospheric CO 2 concentrations as compiled in the spline of Köhler et al. (2017), atmospheric 14 C/ 12 C following IntCal20

Figure 3 .
Figure 3.The effect of (a) CO 2 on model-based Q (b) and atmospheric 14 C inventory (c, right y-axis) in BICYCLE-SE dependent upon whether CO 2 is either internally calculated (V1, PRE) or prescribed from data (V1CO2, PRECO2).In (b), we also plot the contributions to Q based solely on the forced (artificial) carbon fluxes caused by prescribing CO 2 .In (c) the prescribed Δ 14 C from IntCal20 (left y-axis) is additionally plotted.Blue vertical bands mark Heinrich stadials 1-5 (H1-H5) and the YD, similarly to Reimer et al. (2020).During these periods dramatic changes in Atlantic Ocean circulation may have occurred.

Figure 4 .
Figure 4. Radiocarbon production rates Q versus climate change.(a) Atmospheric Δ 14 C from IntCal20(Reimer et al., 2020) and an added extension from new kauri-based data around 42 kcal BP(Cooper et al., 2021).We also show a smoothed Δ 14 C estimate obtained from Hulu Cave alone(Cheng et al., 2018;Southon et al., 2012).This Hulu estimate was used previously by the LSG-OGCM to calculate prior estimates for the MRA of those marine records incorporated in the construction of IntCal20.(b) Q from model-based (scenarios V2, KAURI, HULU in BICYCLE-SE) and data-based (ICE, SED) reconstructions.We plot on the left y axis Q, as atoms/(cm 2 s), and on the right axis (c) changes in Q, in mol/yr, as in all other figures.(d) NGRIP δ 18 O on the GICC05 timescale (NGRIP Members, 2004).(e) 231 Pa/ 230 Th from the Bermuda rise as compiled in Lippold et al. (2019), including data from McManus et al. (2004); Böhm et al. (2015); Henry et al. (2016).(f) WDC δ 18 O on the WD2014 timescale (WAIS Divide Project Members, 2015).(g) Atmospheric CO 2 estimate used in this work (spline based on multi-records as in Köhler et al. (2017)), compared with a recent compilation from WDC (Bauska et al., 2021).Vertical bands mark Heinrich stadials and the YD (light blue, following Waelbroeck et al. (2019), as also shown in Reimer et al. (2020)) and other non-Heinrich Greenland stadials (GS) (pink, Rasmussen et al. (2014)) with their numbers shown as labels on the top.
Figure 4. Radiocarbon production rates Q versus climate change.(a) Atmospheric Δ 14 C from IntCal20(Reimer et al., 2020) and an added extension from new kauri-based data around 42 kcal BP(Cooper et al., 2021).We also show a smoothed Δ 14 C estimate obtained from Hulu Cave alone(Cheng et al., 2018;Southon et al., 2012).This Hulu estimate was used previously by the LSG-OGCM to calculate prior estimates for the MRA of those marine records incorporated in the construction of IntCal20.(b) Q from model-based (scenarios V2, KAURI, HULU in BICYCLE-SE) and data-based (ICE, SED) reconstructions.We plot on the left y axis Q, as atoms/(cm 2 s), and on the right axis (c) changes in Q, in mol/yr, as in all other figures.(d) NGRIP δ 18 O on the GICC05 timescale (NGRIP Members, 2004).(e) 231 Pa/ 230 Th from the Bermuda rise as compiled in Lippold et al. (2019), including data from McManus et al. (2004); Böhm et al. (2015); Henry et al. (2016).(f) WDC δ 18 O on the WD2014 timescale (WAIS Divide Project Members, 2015).(g) Atmospheric CO 2 estimate used in this work (spline based on multi-records as in Köhler et al. (2017)), compared with a recent compilation from WDC (Bauska et al., 2021).Vertical bands mark Heinrich stadials and the YD (light blue, following Waelbroeck et al. (2019), as also shown in Reimer et al. (2020)) and other non-Heinrich Greenland stadials (GS) (pink, Rasmussen et al. (2014)) with their numbers shown as labels on the top.

Figure 5 .
Figure 5. BICYCLE-SE model output.(a) Data (Köhler et al., 2017), and model-derived estimates of atmosphere CO 2 concentrations.(b) The value of Q assumed in different model scenarios (V1, V2, ICE, and SED).Before 60 kcal BP, Q in scenarios ICE and SED are linearly extended.See text for details.(c) Simulated atmospheric Δ 14 C and the IntCal20 (Reimer et al., 2020) data-based reconstruction.Scenario ICEPRE is identical to ICE with respect to Q; however climate and hence CO 2 are fixed at preindustrial levels.

Figure 6 .
Figure6.The contributions of selected processes to changes in simulated atmospheric Δ 14 C in BICYCLE-SE.The carbon fluxes are identical to scenario V1, but the 14 C-signature of individual processes has been modified in these scenarios.Plotted anomalies are calculated with respect to V1. See Subsection 4.2 for details.

Table 1
also used a model and suggested that the LGM ocean stored 33 Pmol less O 2 , while Anderson et al. (2019) derived a reduction of 100 Pmol O 2 based on global extrapolation of deep Pacific proxy data.At the same time, a reduction in photosynthesis during the LGM would have led to a reduction of 41-115 Pmol in atmospheric O 2 .This estimate is based upon a reduction of between 450 and 1,250 PgC in terrestrial glacial carbon content (Jeltsch-Thömmes et al., 2019) and the typical relationship of 1.1 mol O 2 production per 1 mol CO 2 uptake during photosynthesis Amount of Atmospheric Gases of Interest for the 14 C Production Rate and Their Suggested Increase at the LGM Compared to Modern-Day Values