The potential of in situ cosmogenic ¹⁴CO in ice cores as a proxy for galactic cosmic ray flux variations

Galactic cosmic rays (GCRs) interact with matter in the atmosphere and at the surface of the Earth to produce a range of cosmogenic nuclides. Measurements of cosmogenic nuclides produced in surface rocks have been used to study past land ice extent as well as to estimate erosion rates. Because the GCR flux reaching the Earth is modulated by magnetic fields (solar and Earth's), records of cosmogenic nuclides produced in the atmosphere have also been used for studies of past solar activity. Studies utilizing cosmogenic nuclides assume that the GCR flux is constant in time, but this assumption may be uncertain by 30 % or more. Here we propose that measurements of ¹⁴C of carbon monoxide (¹⁴CO) in ice cores at low-accumulation sites can be used as a proxy for variations in GCR flux on timescales of several thousand years. At low-accumulation ice core sites, ¹⁴CO in ice below the firn zone originates almost entirely from in situ cosmogenic production by deep-penetrating secondary cosmic ray muons. The flux of such muons is almost insensitive to solar and geomagnetic variations and depends only on the primary GCR flux intensity. We use an empirically constrained model of in situ cosmogenic ¹⁴CO production in ice in combination with a statistical analysis to explore the sensitivity of ice core ¹⁴CO measurements at Dome C, Antarctica, to variations in the GCR flux over the past ≈ 7000 years. We find that Dome C ¹⁴CO measurements would be able to detect a linear change of 6 % over 7 ka, a step increase of 6 % at 3.5 ka or a transient 100-year spike of 190 % at 3.5 ka at the 3σ significance level. The ice core ¹⁴CO proxy therefore appears promising for the purpose of providing a high-precision test of the assumption of GCR flux constancy over the Holocene.

1 Introduction

The galactic cosmic ray (GCR) flux at Earth is modulated by both the geomagnetic and the heliospheric (solar) magnetic fields. The heliospheric magnetic field strength is linked to solar activity and solar irradiance (e.g., Wu et al., 2018b; Steinhilber et al., 2009), with irradiance being a key climate driver. This has enabled the use of records of past cosmogenic nuclide production rates for studies of past solar variability (e.g., Adolphi et al., 2014; Bard et al., 2000; Steinhilber et al., 2009; Usoskin et al., 2016; Usoskin, 2023). The two main nuclides that have been used for these studies are ¹⁴C (mainly from tree rings, which record atmospheric ${}^{\mathrm{14}}\mathrm{C}{/}^{\mathrm{12}}\mathrm{C}$ ratio) and ¹⁰Be (from ice cores, which record the flux of ¹⁰Be at the snow deposition site). Measurements of cosmogenic nuclides have also been used extensively to study processes at the Earth's surface, such as the timing of glacial advance and retreat (e.g., Balco, 2020). In these applications, nuclides such as ¹⁰Be, ¹⁴C, ²⁶Al and ³⁶Cl that are produced in situ in surface rocks are of interest.

Cosmogenic nuclide-based reconstructions of past solar activity and ice extent have relied on the assumption that the GCR flux is constant in time (e.g., Balco, 2011; Gosse and Phillips, 2001; Muscheler, 2013). Measurements of cosmogenic radionuclides in meteorites provide arguably the strongest support for this assumption (e.g., Smith et al., 2019; Wieler et al., 2013). However, a number of important uncertainties are involved when interpreting these measurements, including meteoroid orbits, solar modulation of the GCR flux and break-up of meteoroids/fresh surface exposure on entry into the atmosphere. A review by Wieler et al. (2013) concluded that while overall the meteorite evidence indicates that the GCR flux is constant, this assumption is uncertain by 30 % or more. Measurements of cosmogenic nuclides in lunar rocks also indirectly indicate that the GCR flux could have been constant on million-year timescales, although there is still a confounding influence of solar modulation (e.g., Poluianov et al., 2018). Records of ${}^{\mathrm{10}}\mathrm{Be}{/}^{\mathrm{9}}\mathrm{Be}$ ratios in oceanic sediments and iron–manganese crusts (Willenbring and von Blanckenburg, 2010) have also been used to argue that the GCR flux is approximately constant on million-year timescales (Wieler et al., 2013). However, this approach also involves multiple confounding factors, such as solar and geomagnetic modulation of the GCR flux and ¹⁰Be transport, deposition, and oceanic cycling. Results from studies that have used cosmogenic ¹⁴C and ¹⁰Be to examine past solar activity also assume that there were no large changes in the GCR flux in the past few millennia (e.g., Knudsen et al., 2009; Steinhilber et al., 2012; Wu et al., 2018a). However, again, inferences about the GCR flux from such records are complicated by solar and geomagnetic modulation (e.g., Knudsen et al., 2008), carbon cycle (for ¹⁴C, e.g., Muscheler et al., 2007), and transport and deposition effects (for ¹⁰Be, e.g., Field et al., 2006).

Theoretical considerations also generally support the assumption that the GCR flux is constant, though small anisotropies are expected due to the effect of the nearest sources of GCRs and the diffusive propagation of cosmic rays in the galaxy (Erlykin and Wolfendale, 2006; Blasi and Amato, 2012; Ahlers and Mertsch, 2015; Mertsch and Funk, 2015). At energies above 100 GeV, the GCR flux at Earth today is isotropic to within 1 part in 1000, with the residual anisotropy characterized by a dipole plus statistically significant components on angular scales as small as 5° (e.g., Abeysekara et al., 2019, and references therein). The observations indicate that cosmic ray transport is dominated by diffusion in galactic magnetic fields, which should dampen the contributions of spatial and temporal point sources of cosmic rays. Nevertheless, significant GCR flux variations are in principle possible even on sub-millennial timescales. For example, Melott et al. (2017) and Thomas et al. (2016) consider the terrestrial effects of a supernova 50 parsecs from Earth and estimate that the production rate of atmospheric muons could increase by up to several orders of magnitude depending on how accelerated GCRs propagate through nearby galactic magnetic fields. While the predictions of such models should be understood to represent the extreme upper limit of possible effects, a number of observations suggest that supernova explosions in our galactic neighborhood do produce measurable effects on the local properties of GCRs. Such observations include the part-per-mille dipole anisotropy in the cosmic ray flux above 1 TeV (e.g., Abeysekara et al., 2019; Ahlers and Mertsch, 2015; Blasi and Amato, 2012; Erlykin and Wolfendale, 2006), the fluxes of positrons and antiprotons above 20 GeV and heavy nuclei above 1 TeV (e.g., Kachelriess et al., 2015), and the measurements of ⁶⁰Fe in ocean sediments (Wallner et al., 2016) and Antarctic snow (Koll et al., 2019). Thus, high-precision tests of GCR flux variations that are free of the confounding factors discussed above for meteorites and for cosmogenic ¹⁰Be and ¹⁴C produced in the atmosphere would be valuable.

2 Systematics of in situ cosmogenic ¹⁴CO in glacial ice

2.1 Overview of ¹⁴CO in glacial ice

We first provide an overview of the current understanding of the processes that control the abundance of in situ cosmogenic ¹⁴CO in glacial ice, which is needed to understand how the ice core ¹⁴CO proxy for GCR flux variations works. ¹⁴C in glacial ice originates from trapping of ¹⁴C-containing atmospheric gases such as carbon dioxide (CO₂), methane (CH₄) and carbon monoxide (CO) as well as from in situ cosmogenic production. In situ ¹⁴C is produced in glacial ice and firn via interactions of secondary cosmic ray neutrons and muons with ¹⁶O in the ice grains (Fig. 1a) (e.g., Lal et al., 1997; Petrenko et al., 2016; van der Kemp et al., 2002). Once produced, this ¹⁴C reacts rapidly to form predominantly ¹⁴CO₂ and ¹⁴CO, with a small amount of ¹⁴CH₄ and possibly other organics also being formed (e.g., Fang et al., 2021; Lal et al., 2000; Petrenko et al., 2013; van de Wal et al., 2007). ¹⁴C production rates are highest near the surface, where neutron-induced spallation of ¹⁶O is the main production mechanism. The neutron flux is attenuated rapidly with depth, however, and only affects the uppermost ≈ 20 m of the firn (or uppermost ≈ 10 m of solid ice) (e.g., Lal et al., 1987). Below these depths, production of ¹⁴C proceeds at lower rates and is dominated by negative muon capture as well as interactions with fast muons (Fig. 1b) (Petrenko et al., 2016; van der Kemp et al., 2002).

Figure 1 Introduction to production and loss of ¹⁴C in firn and ice. (a) Simplified schematic of the firn column, illustrating in situ cosmogenic production ¹⁴C by neutrons (n) and muons (µ) and loss in the upper, diffusive part of the firn as well as ¹⁴C production by muons below the lock-in depth where all of the ¹⁴C is retained. (b) ¹⁴CO production rates calculated as described in Sect. 2.2. Here $f_{\mathit{\mu }-}=\mathrm{0.0667}$ and f_μf=0.0722, which are mid-range values from the range constrained by Hmiel et al. (2024).

The concentration of in situ ¹⁴C in glacial ice at accumulation sites is controlled by the ¹⁴C production rates (site and depth-dependent), the snow accumulation rate and the retention of ¹⁴C in the firn. Sites at higher altitudes have less atmospheric shielding from cosmic rays, resulting in higher ¹⁴C production rates at the surface (e.g., Lifton et al., 2014). At sites with low accumulation rates, ice layers spend more time at relatively shallower depths, allowing for more in situ ¹⁴C to be produced. Further, prior work has shown that most of the in situ cosmogenic ¹⁴C that is produced in the firn is rapidly lost to the atmosphere (de Jong et al., 2004; Hmiel et al., 2024; Petrenko et al., 2013). Because of this, the majority of the in situ cosmogenic ¹⁴C in glacial ice is from production below the firn zone. Of the ¹⁴C-containing gases in glacial ice, ¹⁴CO has the highest ratio of in situ cosmogenic to trapped atmospheric ¹⁴C. This is due to (1) atmospheric ¹⁴CO concentrations being lower than those for ¹⁴CO₂ or ¹⁴CH₄ (mainly because global mean mole fractions for CO (≈ 80 nmol mol⁻¹) are much lower than those for CO₂ (≈ 420 µmol mol⁻¹) and CH₄ (≈ 1920 nmol mol⁻¹) (NOAA Global Monitoring Laboratory data viewer)) and (2) the relatively large fraction of in situ ¹⁴C that forms CO in ice (≈ 0.31; Dyonisius et al., 2023; Hmiel et al., 2024). This makes ¹⁴CO the best species for investigating the in situ cosmogenic component of ¹⁴C in ice.

In situ ¹⁴CO in glacial ice is present at very low concentrations (a few molecules per gram of ice is typical; see Figs. 2 and 3), making measurements very analytically challenging. Prior studies have either worked with relatively small (a few kg) ice samples available from a single shared ice core (e.g., van der Kemp et al., 2002), resulting in relatively large uncertainties, or required dedicated ice coring campaigns to obtain large ice amounts (100 kg or more) from multiple parallel ice cores for high-precision measurements (e.g., Dyonisius et al., 2023). Dry extraction of air from ice has been used for smaller ice samples (van der Kemp et al., 2002), and melt-extraction has been used for large samples (Dyonisius et al., 2023). CO in the extracted air is separated, it is combusted to CO₂, this CO₂ is subsequently converted to graphite, and then the ${}^{\mathrm{14}}\mathrm{C}{/}^{\mathrm{13}}\mathrm{C}$ or ${}^{\mathrm{14}}\mathrm{C}{/}^{\mathrm{12}}\mathrm{C}$ ratio is measured via accelerator mass spectrometry. A detailed description of the ice core ¹⁴CO measurement methodology can be found in Dyonisius et al. (2023).

Figure 2 Overview of ¹⁴CO results from Taylor Glacier. (a) Ice parcel back-trajectories for the deepest (72 m) Taylor Glacier ¹⁴CO sample. The solid black line shows the best-estimate flow trajectory, and the shaded envelope represents the 68 % confidence interval (CI). (b) Comparison of Taylor Glacier ¹⁴CO measurements with model predictions for accepted scenarios. (c) Accepted ranges of f_μ− and f_μf. Figures modified from Dyonisius et al. (2023).

Figure 3 Overview of ¹⁴CO results from Greenland Summit. (a) Measured ¹⁴CO content in ice grains and closed porosity along with a model fit. Horizontal error bars represent the depth range of firn and ice included in each sample. (b) Contour plot of the accepted ranges of the R₁ (initial retention) and L₁ (slow leakage) parameters in firn, together with the best-fit solution as well as solutions that result in maximum and minimum ¹⁴CO content in ice below the firn zone. (c) Contour plot of accepted ranges of f_μ− and f_μf from Taylor Glacier (dashed lines) and after further constraints from Greenland Summit measurements (solid lines).

2.2 Production of ¹⁴CO in glacial ice

Prior studies (Dyonisius et al., 2023; Hmiel et al., 2024) have presented detailed parameterizations of in situ cosmogenic ¹⁴CO production rates in glacial ice. This work uses the same parameterizations, which are described again here for the reader's convenience. The ¹⁴C production rate in ice via the neutron mechanism declines exponentially with depth, with the ¹⁴CO production rate calculated following Hmiel et al. (2024) as

$$\begin{array}{}\text{(1)}& P_n^{\mathrm{CO}}\left(h\right)={\mathrm{\Omega }}^{\mathrm{CO}}\cdot F_n\cdot S_n\cdot Q_c\cdot P_{n,\mathrm{SLHL}}^{\mathrm{Qtz}}\left(\mathrm{0}\right)\cdot e^{-h/{\mathrm{\Lambda }}_n}.\end{array}$$

In this equation, h is the mass-depth (in g cm⁻²), Ω^CO is the fraction of total in situ ¹⁴C that forms CO (we use 0.31, following Hmiel et al., 2024) and F_n is an adjustable dimensionless parameter that allows for tuning the neutron mechanism production rate within uncertainties (0.9–1.1 range). S_n is the site-specific dimensionless scaling factor which describes the ratio of ¹⁴C production rate at the site of interest to ¹⁴C production rate at a sea-level high-latitude reference site; S_n is determined using the model of Lifton et al. (2014). Q_c is a factor that translates ¹⁴C production rate from quartz to ice using the difference in oxygen atom density (atoms g⁻¹) between ice and quartz (Q_c=1.667). $P_{n,\mathrm{SLHL}}^{\mathrm{Qtz}}\left(\mathrm{0}\right)$ is the reference ¹⁴C production rate at the surface via the neutron mechanism in quartz at a sea-level high-latitude site; we use a value of 12.76 atoms g Qtz⁻¹ yr⁻¹ from the CRONUS-Earth project, which is defined for the 2001–2010 mean solar modulation and geomagnetic field conditions (Borchers et al., 2016). Λ_n is the absorption mean free path of neutrons in ice; we use a value of 150 g cm⁻² (Lal et al., 1987; van de Wal et al., 2007).

For ¹⁴CO production by the muon mechanisms, we use a model developed by Balco et al. (2008) (“Balco model”), which incorporates parameterizations of Heisinger et al. (2002a, b). The ¹⁴C production rate via negative muon capture in these parameterizations is calculated using

$$\begin{array}{}\text{(2)}& P_{\mathit{\mu }-}\left(h\right)=R_{\mathit{\mu }-}\left(h\right)\cdot f_{\mathrm{C}}\cdot f_{\mathrm{D}}\cdot f^{\ast },\end{array}$$

where R_μ−(h) is the stopping rate of negative muons (muons g⁻¹ yr⁻¹) at mass-depth h, f_C is the chemical compound factor representing the probability that the stopped muon is captured by one of the target atoms, f_D is the probability that the negative muon does not decay in the K shell before nuclear capture and f^∗ is the effective probability for production of the cosmogenic nuclide of interest after μ⁻ capture by the target nucleus. For production of ¹⁴C from ¹⁶O in ice, f_C=1, f_D=0.1828 and $f^{\ast }=\mathrm{0.137}$ (Heisinger et al., 2002a).

The ¹⁴C production rate via the fast muon mechanism is given by Heisinger et al. (2002b) as

$$\begin{array}{}\text{(3)}& P_{\mathit{\mu }\mathrm{f}}\left(h\right)={\mathit{\sigma }}_{\mathrm{0}}\cdot \mathit{\beta }\left(h\right)\cdot \mathit{\varphi }\left(h\right)\cdot \stackrel{\mathrm{‾}}{E}(h{)}^{\mathit{\alpha }}\cdot N,\end{array}$$

where σ₀ is the reference nuclear reaction cross section at a muon energy of 1 GeV (cm²), ϕ(h) is the total muon flux at mass-depth h (muons cm⁻² yr⁻¹ sr⁻¹), $\overline{E}\left(h\right)$ is the mean muon energy at mass-depth h (GeV), α is a power factor that describes the energy dependence of the cross section (we use α=0.75, consistent with Dyonisius et al., 2023, and Heisinger et al., 2002b) and N is the number of target nuclei per gram target mineral (Hmiel et al., 2024). We use σ₀=8.8 µb $=\mathrm{8.8}\times {\mathrm{10}}^{-\mathrm{30}}$ cm² (Heisinger et al., 2002b) and N= (6.022×10²³ atoms mol⁻¹)$/$(18.02 g mol⁻¹) for ice. β(h) is a unitless depth dependence factor (≈ 0.9 at our depths of interest, with only a slight dependence on depth), given by Heisinger et al. (2002b) as

$$\begin{array}{}\text{(4)}& \mathit{\beta }\left(h\right)={\displaystyle \frac{\stackrel{\mathrm{‾}}{E(h{)}^{\mathit{\alpha }}}}{{\stackrel{\mathrm{‾}}{E\left(h\right)}}^{\mathit{\alpha }}}}.\end{array}$$

Heisinger et al. (2002b) also provide an approximate function for β(h), which is used in the Balco model and hence in our model:

$$\begin{array}{}\text{(5)}& \mathit{\beta }\left(h\right)=\mathrm{0.846}-\mathrm{0.015}\ln \left(h+\mathrm{1}\right)+\mathrm{0.003139}(\ln \left(h+\mathrm{1}\right){)}^{\mathrm{2}}.\end{array}$$

The Balco model incorporates Eqs. (2) and (3) and also provides the muon fluxes and energies as a function of mass-depth h for a given site, using site atmospheric pressure as input.

2.3 Constraints on in situ ¹⁴CO production rates from measurements at Taylor Glacier, Antarctica

Recent studies at Taylor Glacier, Antarctica (an ice ablation site that exposes ancient ice at the surface), have provided measurements of ¹⁴C in ice older than 50 ka (Dyonisius et al., 2023; Petrenko et al., 2016). In such ice, any ¹⁴C from the snow accumulation site (from ¹⁴C-containing atmospheric gases trapped in air bubbles or from in situ cosmogenic production) has decayed away (¹⁴C half-life is 5700 years), and the only measurable ¹⁴C originates from relatively slow in situ cosmogenic production by muons as the glacier transports the ice at large depths and somewhat faster production as the ice gradually rises toward the surface via ablation. Due to the relatively fast ice ablation rate of ≈ 20 cm yr⁻¹, the ¹⁴C contribution from the neutron production mechanism is negligible for ice deeper than 6 m. This presented an opportunity to use ¹⁴CO measurements in Taylor Glacier to constrain the muogenic ¹⁴CO production rates in ice in a natural setting.

Dyonisius et al. (2023) presented measurements of ¹⁴CO in Taylor Glacier ice between the surface and 72 m depth. An ice flowline model for Taylor Glacier (Buizert et al., 2012b) was used to reconstruct the possible range of trajectories for the sampled ice parcels (Fig. 2a). The Balco model was used to calculate ¹⁴CO production via the muon mechanisms as ice parcels followed the trajectories. As prior work suggested that muogenic ¹⁴C production rates from Heisinger et al. (2002a, b) may be too high when applied to ice (Petrenko et al., 2016), Dyonisius et al. (2023) introduced production rate adjustment factors f_μ− and f_μf into production rate equations as follows:

$$\begin{array}{}\text{(6)}& {\displaystyle }{P}_{\mathit{\mu }-}^{\mathrm{CO}}\left(h\right)=f_{{\mathit{\mu }}^{-}}\cdot P_{{\mathit{\mu }}^{-}}^{\mathrm{Balco}}\left(h,P\right),\text{(7)}& {\displaystyle }{P}_{\mathit{\mu }\mathrm{f}}^{\mathrm{CO}}\left(h\right)=f_{\mathit{\mu }\mathrm{f}}\cdot P_{\mathit{\mu }\mathrm{f}}^{\mathrm{Balco}}\left(h,P\right).\end{array}$$

Here $P_{\mathit{\mu }}^{\mathrm{Balco}}\left(h,P\right)$ is the total ¹⁴C production rate (in atoms g⁻¹ yr⁻¹) in the Balco model for the respective muon mechanism at mass-depth h and surface pressure P. f_μ− and f_μf account for (1) the fraction of total ¹⁴C that forms ¹⁴CO (Ω^CO) and (2) the adjustment factor for total ¹⁴C production rate.

To define the best-estimate ¹⁴CO production rate adjustment factors f_μ− and f_μf, Dyonisius et al. (2023) used a grid search approach, as follows. Using the best-estimate ice parcel back-trajectory (Fig. 2a), an expected ¹⁴CO depth profile was calculated for each combination of f_μ− and f_μf between 0 and 0.2 at 0.001 resolution. The model results were then compared to ¹⁴CO measurements (Fig. 2b) with mean depths of 6.85 m or deeper (to avoid significant effects from the neutron mechanism), and a χ² metric was used to determine the goodness of fit. To define the possible range of f_μ− and f_μf, Dyonisius et al. (2023) used a Monte Carlo approach, as follows. First, 10 000 possible ice back-trajectories were generated by perturbing ablation rates along the glacier according to their uncertainties (Fig. 2a). Next, a wide prior distribution for f_μ− and f_μf was defined by starting with the best-estimate values and assuming a large and normally distributed 200 % uncertainty in these values. Thus, 100 000 Monte Carlo iterations of the model were then run, with each iteration randomly selecting a back-trajectory scenario and a pair of f_μ− and f_μf from the prior distribution described above. All pairs of f_μ− and f_μf that yielded ¹⁴CO depth profiles (Fig. 2b) that were within average measurement uncertainty (1σ or 2σ) from the best-fit solution were accepted (Fig. 2c).

2.4 Constraints on in situ ¹⁴CO retention and leakage in firn and production in ice at Greenland Summit

In situ cosmogenic ¹⁴C that is produced in the firn column above the lock-in depth can be lost to the atmosphere if it is able to leak out of the ice grains, resulting in low ¹⁴C retention into ice below the firn zone (e.g., de Jong et al., 2004; Petrenko et al., 2013, and references therein). Hmiel et al. (2024) used Greenland Summit to conduct the most comprehensive study to date of in situ cosmogenic ¹⁴C in the firn, with a focus on ¹⁴CO. This study measured ¹⁴CO in the ice grains in the firn matrix, in firn air, as well as in bubbly ice below the firn zone. Very large firn and ice samples (200–300 kg) were used for ¹⁴CO analysis, to provide sufficiently large numbers of ¹⁴C atoms for precise ¹⁴C measurements. Figure 3a shows the ¹⁴CO results for samples from the firn, firn-ice transition and bubbly ice below the firn zone. In the shallowest firn, ¹⁴CO increases rapidly with depth owing mainly to production by the neutron mechanism, reaching a peak in the 10–20 m depth range. Beyond 20 m, ¹⁴CO in the firn matrix declines gradually with depth in the diffusive part of the firn, reflecting leakage of in situ ¹⁴C from the ice grains. ¹⁴CO increases rapidly in the lock-in zone (≈ 70–80 m), reflecting addition of ¹⁴CO from trapped air. Below the lock-in zone, ¹⁴CO in the ice continues to increase gradually due to deeper production by the muon mechanisms.

To interpret the Greenland Summit ¹⁴CO results, Hmiel et al. (2024) employed a firn gas transport model that can also characterize trapped air in ice below the firn zone (Buizert et al., 2012a). Production of in situ ¹⁴C following the systematics described in Sect. 2.2 and tracking of ¹⁴C in ice grains and porosity was implemented in this model. With regard to ¹⁴C loss from ice grains in the firn, it was found that the model–data agreement was best if two separate loss processes were parameterized in the model: a fast process, with a timescale < 1 year and an additional slow process. This was described in the model using parameters R₁ and L₁. R₁ represents the fraction of in situ ¹⁴C in the ice grains that is initially retained. The fraction of in situ ¹⁴C in ice grains that is lost rapidly, given by 1−R₁, leaks out from the ice grains into the porosity at every model time step (0.5 year). L₁ represents the fraction of the initially retained ¹⁴C that is lost more slowly from the ice grains over the course of 1 year. Hmiel et al. (2024) used a grid search approach to constrain the possible ranges of R₁ and L₁ at Greenland Summit (Fig. 3b), showing that > 99 % of in situ ¹⁴C is lost rapidly from the ice grains, while the remaining ≈ 0.5 % (R₁) of in situ ¹⁴C continues to leak out slowly at a rate of ≈ 0.6 % yr⁻¹ (L₁). Hmiel et al. (2024) argued that the rapid loss is best explained by the process of gas diffusion through ice and suggested that the ≈ 0.5 % of ¹⁴C that is initially retained may be trapped in microbubbles or by impurities at dislocations or grain boundaries and is released via the process of recrystallization.

Greenland Summit ¹⁴CO measurements in ice below the firn zone also provided an opportunity to test muon mechanism ¹⁴CO production rate estimates from Taylor Glacier. For Greenland Summit ice samples, the contribution from trapped atmospheric ¹⁴CO is important (≈ 25 %–40 % of total), and uncertainties in the atmospheric ¹⁴CO history interfere with precise constraints on f_μ− and f_μf. Nevertheless, by trialing the Taylor Glacier sets of accepted f_μ−–f_μf pairs in combination with several possible atmospheric ¹⁴CO histories, Hmiel et al. (2024) were able to further narrow the possible ranges of f_μ− and f_μf (Fig. 3c).

3 In situ cosmogenic ¹⁴CO in ice cores as a possible proxy for GCR flux variability

3.1 Basic concept for using ¹⁴CO in ice cores as a GCR flux proxy

As the Greenland Summit ¹⁴CO results summarized above illustrate, the retention of in situ cosmogenic ¹⁴CO through the upper firn column is very low. This means that the majority of in situ ¹⁴CO found in ice below the firn zone originates from production by muons below the lock-in depth, where this ¹⁴CO can no longer escape to the atmosphere. If the firn layer is sufficiently thick (≈ 90–100 m actual depth or ≈ 65 m ice equivalent depth), the muons penetrating below the firn must have an energy of ≈ 15 GeV or greater at the surface (e.g., Rogers and Tristam, 1984). Such muons originate from primary GCR particles with energies of ≈ 100 GeV or greater (Gaisser et al., 2016). The part of the GCR flux possessing such energies is not affected appreciably by either the geomagnetic or the heliospheric magnetic fields. In situ cosmogenic ¹⁴CO content in ice cores drilled at such sites thus can serve as a proxy of variations in the primary GCR flux. This proxy is in principle free of the confounding effects discussed in the Introduction for other past GCR flux indicators.

Several considerations are important for site selection in order to increase the likelihood of success with this proxy. First, the in situ ¹⁴CO signal must be maximized to help with measurement sensitivity as well as to reduce interference from the trapped atmospheric ¹⁴CO component. Second, the site must have a thick firn column. This is needed to ensure that ¹⁴CO below the firn zone is produced only by muons originating from primary GCRs that are sufficiently energetic to be unaffected by solar magnetic field variations. Third, there should be as little in situ ¹⁴CO retained from the shallow firn as possible. ¹⁴CO produced in the shallow firn originates from neutrons or lower-energy muons that are affected by solar magnetic field variations and may complicate interpretation. Fourth, ideally the site must be glaciologically stable over time in terms of accumulation rate and lock-in depth. Large temporal variations in these parameters may introduce additional uncertainties in the interpretation, as they affect the predicted in situ ¹⁴CO content.

Considering the above, ice dome sites in the East Antarctica interior are most promising for attempting to examine past GCR flux variability using ¹⁴CO in ice cores. Low accumulation rates at such sites maximize cosmogenic exposure times and thereby the in situ ¹⁴CO signal. These sites also tend to have sufficiently thick firn columns (e.g., Buizert, 2013). The combination of low accumulation rate and thick firn column results in very long ice layer transit times through the firn, maximizing the chance that in situ ¹⁴CO produced by neutrons and low-energy muons in the shallow firn would be lost. Finally, dome sites are free of complications of upstream ice advection, and ice core water stable isotope records suggest that the interior East Antarctica climate has been stable over the last few thousand years (recent decades excepted) (e.g., Jouzel et al., 2001).

3.2 Using model predictions to explore Dome C, Antarctica, as a test case for the ¹⁴CO GCR flux proxy

Dome C, Antarctica, is a site that meets the criteria needed for the ¹⁴CO GCR flux proxy to be viable. It has been glaciologically very well characterized as a result of previous ice coring campaigns (e.g., EPICA Community Members, 2004) and has well-established logistical access owing to the permanent Concordia Station. Further, a scientific ice drilling expedition is planned for the near future to Dome C for the purpose of ¹⁴CO reconstruction at this site. We therefore use Dome C as an example site for more detailed model-based exploration of the ¹⁴CO past GCR flux proxy. We first applied the full firn-ice model mentioned above (Buizert et al., 2012a; Hmiel et al., 2024) to explore the (unwanted) contribution of ¹⁴CO originating from production in the shallow firn as well as trapped atmospheric ¹⁴CO to the overall ¹⁴CO signal in ice below the firn zone. In the model, we used an accumulation rate of 3.2 cm ice equivalent yr⁻¹ and the firn density profile from the FIRETRACC project (EU FIRETRACC Campaign participants, 2006), and we tuned the firn gas diffusivity profile based on a combination of available CO₂, CH₄, CFC-11, CFC-12, CFC-113, CH₃CCl₃, SF₆ and δ¹⁵N of N₂ measurements (EU FIRETRACC Campaign participants, 2006). For parameters relevant to in situ ¹⁴CO, we used F_n=1.03, R₁=0.44 % and L₁=0.45 % yr⁻¹ (see Sect. 2), which was the combination of values at Greenland Summit that maximized the amount of in situ ¹⁴CO produced in the shallow firn that is retained into ice below the firn zone (Fig. 3b). For muogenic ¹⁴CO production, we used f_μ− = 0.065 and f_μf=0.07, which are mid-range choices from the possible range of values that were consistent with both Taylor Glacier and Greenland Summit measurements (Fig. 3c). We used a constant concentration of 12 molecules cm⁻³ STP (standard temperature and pressure) for the atmospheric ¹⁴CO history, which is the average of the longest available Antarctic atmospheric ¹⁴CO record (Manning et al., 2005).

Figure 4a shows model-calculated ¹⁴CO content that represents the sum of ¹⁴CO in ice grains and closed porosity (this is what measurements done with a melt-extraction approach would provide). The solid black line shows results with both in situ and atmospheric ¹⁴CO included in the model. There is a sharp ¹⁴CO peak at ≈ 9 m depth that represents ¹⁴CO in ice grains and is driven by intense ¹⁴CO production by the neutron mechanism in near-surface firn. ¹⁴CO then declines to near zero by ≈ 70 m due to slow leakage out of ice grains (controlled by the L₁ parameter in the model). At depths > 80 m, the amount of closed porosity starts to increase, and this increases ¹⁴CO by trapping of ¹⁴CO from open porosity and by allowing more in situ ¹⁴CO to be retained. This process further accelerates at ≈ 95 m, which is the lock-in depth at Dome C. Below the close-off depth at Dome C (≈ 100 m), ¹⁴CO content continues to increase due to production by muons, rising to 8.5 ¹⁴CO molecules g⁻¹ ice at the deepest modeled level (110 m).

The dashed blue line shows the expected contribution to total ¹⁴CO from in situ ¹⁴CO originating only from the shallower part of the firn. This was assessed by setting the atmospheric ¹⁴CO history to zero and setting in situ production rates to zero for depths > 54 m ice equivalent. This contribution is < 0.5 ¹⁴CO molecules g⁻¹ ice and is due to ¹⁴CO that leaks out from ice grains in the shallow firn, diffuses into the deep firn and is subsequently trapped in air bubbles. The contribution from trapped atmospheric ¹⁴CO (dotted pink line; assessed by turning off in situ production in the model) is < 1.2 ¹⁴CO molecules g⁻¹ ice. ¹⁴CO originating from the sum of shallow in situ cosmogenic production and air trapping (solid red line) is < 1.6 ¹⁴CO molecules g⁻¹ ice at all depths below the firn zone.

We next examined the in situ cosmogenic ¹⁴CO component at Dome C arising from production by deep-penetrating muons, as well as its suitability for detecting changes in the past GCR flux. As this approach involved generating thousands of simulated data sets (see Sect. 3.3 below), we created a simple and computationally efficient ice-only model of in situ cosmogenic ¹⁴CO for this test of the proxy concept. This ice-only model has its starting (shallowest) depth in the lock-in zone and assumes an initial ¹⁴CO content of zero. ¹⁴CO production in the model is implemented following parameterizations described above in Sect. 2, with production rates within the range constrained by Taylor Glacier and Greenland Summit results. The model assumes that all of the in situ ¹⁴CO is retained and also includes ¹⁴C radioactive decay. The model defines annual ice layers and shifts these layers downward on an annual basis following the ice layer age scale for Dome C from Buizert et al. (2018). For the purposes of this test of the proxy concept, we set the deepest model depth at 300 m, as this is the practical limit for light ice coring projects that do not use drilling fluid and the deepest depth in the planned fieldwork. The exact starting depth of the model was chosen by comparing predictions of this ice-only model with predictions from the full firn-ice model in the 100–110 m depth range (below close-off depth) when using the same muogenic ¹⁴CO production rates and setting atmospheric ¹⁴CO history to zero; using 96.5 m for the starting depth yields the best match.

Figure 4 Predicted ¹⁴CO content at Dome C. (a) Predicted ¹⁴CO content in ice grains + closed porosity (bubbles) from the full firn-ice model considering all ¹⁴CO contributions as well as contributions from individual ¹⁴CO components. (b) Predicted in situ ¹⁴CO content in ice below the firn zone from a simple ice-only model. The legend indicates the combination of [f_μ−, f_μf] values (see Sect. 2) used in each model run, as well as whether the production rate was assumed to be constant (solid lines) or variable (dashed and dotted lines) in time. Markers illustrate what sample measurements might look like assuming 20 m depth averaging and a random 1σ measurement error of 3 %. (c) Time-variable production rate scenarios used to generate the corresponding depth–¹⁴CO profiles in panel (b).

Figure 4b illustrates predictions of the simple ice model for a few scenarios involving different combinations of f_μ− and f_μf as well as different production rate histories (representing past GCR flux variations), and Fig. 4c illustrates the time-variable production rate scenarios reflected in Fig. 4b. Because in situ cosmogenic ¹⁴CO production takes place at the full range of modeled depths (with production rate declining with depth as illustrated in Fig. 1b), the ¹⁴CO values at each depth represent a time integral of production rate minus the ¹⁴C decay rate. As expected, ¹⁴CO content increases most rapidly at the shallowest depths, followed by a broad peak in the 200–250 m depth range. For deeper ice, the rate of ¹⁴CO removal via radioactive decay exceeds the rate of production by muons, and ¹⁴CO values gradually decline. The modeled ice layers at Dome C span an age range of 7283 years, meaning that an ice core ¹⁴CO record reaching 300 m depth could offer information about past GCR flux variations for most of the Holocene.

Predicted ¹⁴CO content originating from deep-penetrating muons is between 20 and 30 molecules g ice⁻¹ for most of the modeled depth range. This means that the ¹⁴CO contribution arising from trapped atmospheric ¹⁴CO and ¹⁴CO production in shallow firn (< 1.6 molecules g ice⁻¹; Fig. 4a) would contribute < 8 % of total ¹⁴CO and is unlikely to interfere with ¹⁴CO signal interpretation. We note that there were three prior ¹⁴CO measurements that were made on the Dome C ice core in the depth range we are considering (de Jong et al., 2004). Those measurements were made on much smaller (1–2 kg) samples than Taylor Glacier and Summit measurements and thus were associated with very large uncertainties (50 %–100 %, considering only uncertainty reported for ¹⁴C activities). That said, de Jong et al. (2004) reported ¹⁴CO values in the 15–30 molecules g ice⁻¹ range for these three samples and concluded that there was no detectable in situ ¹⁴C retention from the firn, consistent with our model predictions.

Figure 4b illustrates that the absolute ¹⁴CO content in the ice as well as the depth of the ¹⁴CO peak depend on the balance of production rates from the negative muon capture and fast muon mechanisms (solid lines); this is controlled by the f_μ− and f_μf parameters in the model. In the modeled depth range, the fast muon mechanism is relatively more important (Fig. 1b), so maximizing f_μf at the expense of f_μ− (within the range in Fig. 3c) increases total ¹⁴CO and shifts the peak slightly deeper. Despite these differences, the shape of the depth–¹⁴CO curves remains largely similar. Figure 4b also illustrates a few scenarios where the production rates (controlled by GCR flux) vary in time (dashed and dotted lines). For time-variable production rates for the purpose of this illustration (Fig. 4c), we trialed (1) a scenario where production by each muon mechanism increases at a linear rate over the entire duration of the model run, reaching 15 % higher rates by the end of the run (black dotted line), (2) a scenario where there is a 15 % step increase in production rates halfway through the model run (pink dash-dot line) and (3) a scenario where there is a 2 times transient step increase in production rates that takes places between 3000 and 3100 years in the model run (dashed green line). As can be seen, all of these types of variations produce depth–¹⁴CO profiles that are distinct in their shape from the constant production rate scenarios.

3.3 Analysis of sensitivity of ice core ¹⁴CO measurements at Dome C to variations in past GCR flux

We compare several time-varying scenarios to the baseline model of a constant GCR flux with muonic production rates (f_μ−, f_μf) that are consistent with ice core ¹⁴CO measurements in both Dyonisius et al. (2023) and Hmiel et al. (2024). As shown in Fig. 4b, deviations from the baseline model can be produced by temporal variations in the GCR flux. However, in the presence of a steady-state flux, uncertainties in the muonic production rates also create deviations from the baseline model. While the normalization of the depth–¹⁴CO profile is affected by both the temporal variations in the GCR flux and the production rates, the shape of the profile is more sensitive to temporal variations in the flux. Therefore, we develop an analysis that is sensitive to the shape of the ¹⁴CO profile as a function of depth.

To discriminate the steady-state GCR scenario H₀ from the time-varying scenario H₁, we construct a test statistic using a Bayes factor (Jeffreys, 1998; Kass and Raftery, 1995):

$$\begin{array}{}\text{(8)}& B_{\mathrm{01}}={\displaystyle \frac{P\left(\mathit{c}\mathrm{|}{H}_{\mathrm{0}}\right)}{P\left(\mathit{c}\mathrm{|}{H}_{\mathrm{1}}\right)}}={\displaystyle \frac{\int \mathrm{d}{\mathit{\theta }}_{\mathrm{0}}P(\mathit{c}\mathrm{|}{\mathit{\theta }}_{\mathrm{0}},H_{\mathrm{0}})P\left({\mathit{\theta }}_{\mathrm{0}}\mathrm{|}{H}_{\mathrm{0}}\right)}{\int \mathrm{d}{\mathit{\theta }}_{\mathrm{1}}P(\mathit{c}\mathrm{|}{\mathit{\theta }}_{\mathrm{1}},H_{\mathrm{1}})P\left({\mathit{\theta }}_{\mathrm{1}}\mathrm{|}{H}_{\mathrm{1}}\right)}}.\end{array}$$

In this expression, c={c(h_j)} is a ¹⁴CO profile measured as a function of discrete depths $h_j=\left\{h_{\mathrm{1}},\mathrm{\dots },h_N\right\}$. The Bayes factor computes the ratio of the marginal probabilities of measuring ¹⁴CO profile c given the steady-state and time-varying scenarios H₀ and H₁. If the data provide greater evidence for the steady-state hypothesis H₀, B₀₁>1, and if the time-varying hypothesis is supported, B₀₁<1.

During the calculation of B₀₁, each marginal probability, P(c|H_i), can be factorized into two terms: a conditional probability P(c|θ_i,H_i), where θ_i lists the free parameters of GCR flux model H_i, and a “prior” probability distribution P(θ_i|H_i). To complete the calculation, we integrate over the possible values of the parameters θ_i in each model H_i. The prior probabilities specify the allowed ranges of parameters θ_i in model H_i and allow us to weight the calculation toward more probable values of the model parameters. Note that we are free to choose the functional form of the prior distributions using theoretical considerations, past measurements or even our subjective degree of belief in the most likely values of the parameters of a given model. In this work, we use non-informative (or flat) prior distributions that do not favor any particular values of the model parameters, beyond restricting their ranges to physically motivated regions.

In the sensitivity calculation, the muonic production rates (f_μ−, f_μf) are nuisance parameters folded into both θ₀ and θ₁. Using the confidence intervals on (f_μ−, f_μf) from Dyonisius et al. (2023) and Hmiel et al. (2024), we can factorize the prior probability P(θ_i|H_i) for model i into a joint prior $P(f_{\mathit{\mu }-},f_{\mathit{\mu }\mathrm{f}}\mathrm{|}{H}_i$) and a set of independent priors dependent on the parameters of the model. For example, the joint prior $P(f_{\mathit{\mu }-},f_{\mathit{\mu }\mathrm{f}}\mathrm{|}{H}_i$) is given by the “2σ acceptance from Summit ice” contour in Fig. 3c. If we wish to test a cosmic ray model H₁ with a flux that varies linearly in time, the model includes an additional free parameter a representing the rate of change of the flux as a function of time. In the calculation of the Bayes factor, the prior distribution of a is a uniform probability density function:

$$\begin{array}{}\text{(9)}& P\left(a\mathrm{|}{H}_{\mathrm{1}}\right)=\left\{\begin{array}{l}\frac{\mathrm{1}}{a_{\max }-a_{\min }},a\in [a_{\min },a_{\max }];\\ \mathrm{0},\mathrm{otherwise}.\end{array}\right.\end{array}$$

Here a_min and a_max represent the allowed range of values we consider for the rate of change of the flux. We use a uniform distribution for P(a|H₁) because it is unbiased, giving equal weight to all values of a between a_min and a_max.

Table 1 Simulated sensitivity to temporal changes in the GCR flux. We report the magnitude of GCR flux changes in time-varying models required to produce a 3σ or 5σ detection in at least 50 % of simulated data sets, assuming 2 % (3 %) relative uncertainties in the ¹⁴CO measurements. For example, to produce a 3σ detection of a linearly increasing or decreasing GCR flux, the rate of change of the flux must be at least 6 % (7 %) over 7 ka.

Our calculation assumes the ice core ¹⁴CO measurements are depth-averaged over 20 m, and each measurement has independent Gaussian uncertainties with relative sizes of 2 % at 1σ. The 20 m depth averaging is assumed because this would provide the needed amount of ice for high-precision ¹⁴CO measurements (≈ 140 kg with a 10 cm diameter ice core). Recent improvements in analytical techniques for ice core and atmospheric ¹⁴CO measurements (Petrenko et al., 2023, 2021) make 2 % uncertainties achievable, although for completeness we also repeat the calculations assuming 3 % 1σ uncertainties. The conditional probability of observing a ¹⁴CO profile c given GCR model H_i with parameters θ_i is

$$\begin{array}{}\text{(10)}& \begin{array}{rl}& P\left(\mathit{c}\mathrm{|}{\mathit{\theta }}_i,H_i\right)=\\ & {\prod }_{j=\mathrm{1}}^N{\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\pi }}{\mathit{\sigma }}_j}}\exp \left\{-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left({\displaystyle \frac{c_j-\widehat{c}(h_j\mathrm{|}{\mathit{\theta }}_i,H_i)}{{\mathit{\sigma }}_j}}\right)}^{\mathrm{2}}\right\}.\end{array}\end{array}$$

Here $\widehat{c}(h_j\mathrm{|}{\mathit{\theta }}_i,H_i$) is the expected ¹⁴CO profile at depth h_j, which we compute using the model, while c_j=c(h_j) is the observed ¹⁴CO concentration at depth h_j. Since the measurement uncertainties are independent and Gaussian, the probability is the product of N independent Gaussian probability density functions over the N measurements in the depth profile. Multiplying this probability by the prior distributions of the nuisance parameters (f_μ−, f_μf) and the allowed prior ranges of the model parameters (such as the slope of the linear change in the GCR flux) allows us to account for both systematic and statistical uncertainties in the measurement.

We calculate our sensitivity to a given GCR scenario as follows:

1. We produce 5×10⁶ random realizations of the ¹⁴CO profile at Dome C, assuming a constant production rate but accounting for the systematic uncertainties in (f_μ−, f_μf). The profiles are generated with depth averaging of 20 m, and relative measurement uncertainties of ${\mathit{\sigma }}_j/c_j=\mathrm{2}$ % and 3 % are both investigated.

2. For each time-varying model under consideration, we compute a distribution of Bayes factors B₀₁ using the random constant-flux data sets. This provides us with a distribution of the Bayes factor when the null constant-flux hypothesis H₀ is true.

3. We next produce a large set of independent ¹⁴CO profiles assuming the alternative time-varying hypothesis H₁ is true and compute the Bayes factor $B_{\mathrm{01}}^{\ast }$ for each simulated data set. We expect that $B_{\mathrm{01}}^{\ast }$ will be much smaller than B₀₁, on average, since P(c|H₁)>P(c|H₀) when the alternative hypothesis H₁ is true.

4. For each $B_{\mathrm{01}}^{\ast }$, we compute the tail probability, or p value, that gives the probability that a constant flux model will produce a Bayes factor smaller than the time-varying model purely by a chance statistical fluctuation, i.e.,

   $$\begin{array}{}\text{(11)}& p=P\left(B_{\mathrm{01}}<B_{\mathrm{01}}^{\ast }\mathrm{|}{H}_{\mathrm{0}}\right).\end{array}$$

The reported sensitivity of a given model is the value of the model parameter(s) in which at least 50 % of simulated data sets yield $p<{\mathrm{10}}^{-\mathrm{3}}$ (3σ evidence against the steady-state model). We also report the value of the model parameter(s) yielding $p<\mathrm{3}\times {\mathrm{10}}^{-\mathrm{7}}$, corresponding to a 5σ discovery of a time-varying flux. This “calibration” of the Bayes factor accounts for the chance probability that a steady-state flux could produce a false positive report of a time-varying flux.

We investigated scenarios involving (1) a linear GCR flux increase over the entire duration of the record, (2) a step-like increase at approximately the mid-point of the record and (3) a brief (100-year) burst in the GCR flux. The scenarios start at a common point in model year 0 (≈7 ka) and diverge during the model run (toward present day). The results are reported in Table 1. For a scenario H₁ where the GCR flux increases linearly with time and assuming 2 % (3 %) relative uncertainties in the measured ¹⁴CO profile, a flux increase a=6 % (7 %) over 7 ka is required to produce a 3σ evidence of a non-steady flux in at least 50 % of simulated data sets. For a 5σ detection, the rate of change of the flux must be at least a=8 % (10 %). We also investigated and found similar sensitivities for a scenario involving a step-like increase in the GCR flux at 3.5 ka. Much larger GCR flux changes are required for detection in the impulsive burst scenario: 190 % (240 %) for 3σ evidence. This is likely due to the large amount of temporal averaging (≈ 700 years) that is imposed by the 20 m depth averaging for the measurements and the fact that the ¹⁴CO content at each depth level represents a time integral of production rates. We further note that improving the relative uncertainty in the ¹⁴CO measurement from 3 % to 2 % has a minor effect on the sensitivity to linear and step-like increases in the GCR flux, but the change in sensitivity to burst-like increases in the flux is substantial.

4 Conclusions

¹⁴CO in ice cores at low-accumulation sites such as Dome C, Antarctica, has a good potential to provide a test of the assumption of GCR flux constancy over the Holocene and to serve as a proxy for past variations in the GCR flux on timescales of a few thousand years. ¹⁴CO measurements in the proposed approach would be most sensitive to gradual linear or step-like changes in the GCR flux, in principle allowing us to test the assumption of GCR flux constancy to 7 % or better. This would represent a large improvement over the ≈ 30 % uncertainty associated with constraints from meteorite measurements. Because our approach involves a large amount of temporal averaging, sensitivity to short-lived GCR bursts is much worse. However, such bursts (if present) would have been captured by high-resolution records of other cosmogenic nuclides such as ice core ¹⁰Be and tree-ring ¹⁴C.

We note that cosmogenic nuclides produced in the atmosphere such as ¹⁰Be are primarily sensitive to the GCR flux below 10 GeV, while the ¹⁴CO proxy discussed here is sensitive to the flux above 100 GeV. The extent to which temporal variations in the GCR flux above 100 GeV would produce proportional changes below 10 GeV, while beyond the scope of this paper, is an interesting question to consider, as the answer depends on the origin of the temporal variations. Since the diffusion length of cosmic rays increases with energy, it is reasonable to expect that a constant GCR flux at high energy is likely to imply a constant GCR flux below 10 GeV, while a time-varying flux above 100 GeV could still be consistent with a constant flux at or below 10 GeV.

For the most precise results, the ¹⁴CO proxy approach requires an ice dome site that is glaciologically stable (accumulation rate, lock-in depth) over the duration of the GCR flux reconstruction. Although our work indicates that the ¹⁴CO GCR flux proxy is likely to provide useful results for most of the Holocene, we expect that GCR flux reconstructions beyond the Holocene with this approach would be more challenging, owing to (1) the need for drilling fluid to obtain ice below ≈ 300 m, which would greatly increase logistical requirements and introduce added challenges of CO contamination from the drilling fluid; (2) glaciological changes beyond the Holocene; and (3) reduced ¹⁴CO signal at greater depths due to ¹⁴C radioactive decay.