Using ice core measurements from Taylor Glacier, Antarctica, to calibrate in situ cosmogenic ¹⁴C production rates by muons

Cosmic rays entering the Earth's atmosphere produce showers of secondary particles such as protons, neutrons, and muons. The interaction of these particles with oxygen-16 (¹⁶O) in minerals such as ice and quartz can produce carbon-14 (¹⁴C). In glacial ice, ¹⁴C is also incorporated through trapping of ¹⁴C-containing atmospheric gases (¹⁴CO₂, ¹⁴CO, and ¹⁴CH₄). Understanding the production rates of in situ cosmogenic ¹⁴C is important to deconvolve the in situ cosmogenic and atmospheric ¹⁴C signals in ice, both of which contain valuable paleoenvironmental information. Unfortunately, the in situ ¹⁴C production rates by muons (which are the dominant production mechanism at depths of >6 m solid ice equivalent) are uncertain. In this study, we use measurements of in situ ¹⁴C in ancient ice (>50 ka) from the Taylor Glacier, an ablation site in Antarctica, in combination with a 2D ice flow model to better constrain the compound-specific rates of ¹⁴C production by muons and the partitioning of in situ ¹⁴C between CO₂, CO, and CH₄. Our measurements show that 33.7 % (±11.4 %; 95 % confidence interval) of the produced cosmogenic ¹⁴C forms ¹⁴CO and 66.1 % (±11.5 %; 95 % confidence interval) of the produced cosmogenic ¹⁴C forms ¹⁴CO₂. ¹⁴CH₄ represents a very small fraction (<0.3 %) of the total. Assuming that the majority of in situ muogenic ¹⁴C in ice forms ¹⁴CO₂, ¹⁴CO, and ¹⁴CH₄, we also calculated muogenic ¹⁴C production rates that are lower by factors of 5.7 (3.6–13.9; 95 % confidence interval) and 3.7 (2.0–11.9; 95 % confidence interval) for negative muon capture and fast muon interactions, respectively, when compared to values determined in quartz from laboratory studies (Heisinger et al., 2002a, b) and in a natural setting (Lupker et al., 2015). This apparent discrepancy in muogenic ¹⁴C production rates in ice and quartz currently lacks a good explanation and requires further investigation.

1 Introduction

1.1 Potential applications of ¹⁴C measurements in ice and in situ cosmogenic ¹⁴C production from ¹⁶O in Earth's surface minerals

As snow accumulates on ice sheets, it gradually densifies into firn and ice (Herron and Langway, 1980). During the firn-to-ice transition, the air in the interstitial space between the ice grains becomes trapped into bubbles within the ice matrix (Buizert, 2013). Included in the paleoatmospheric air trapped in the bubbles are ¹⁴C-containing atmospheric gases (¹⁴CO₂, ¹⁴CO, and ¹⁴CH₄; Fireman and Norris, 1982). ¹⁴C in ice is also produced through interactions of secondary cosmic rays with ¹⁶O directly in the lattice of the ice grains (i.e., in situ; Lal et al., 1990). Following the cosmogenic nuclear reactions, the hot ¹⁴C atom interacts with atoms in the surrounding ice lattice to produce ¹⁴CO₂, ¹⁴CO, and ¹⁴CH₄ (Lal et al., 1990; Petrenko et al., 2013).

Both the trapped atmospheric and in situ cosmogenic ¹⁴C signals in ice have unique applications. For example, the paleoatmospheric component of ¹⁴CH₄ in ice cores has been used to constrain past CH₄ emissions from old carbon reservoirs such as methane hydrates, permafrost, and geologic seeps (Dyonisius et al., 2020; Hmiel et al., 2020; Petrenko et al., 2009, 2017). Paleoatmospheric ¹⁴CO₂ can be potentially used for the absolute dating of ice core gases (Andree et al., 1984; Van De Wal et al., 1994) and to improve the radiocarbon calibration curve (Reimer et al., 2020; Hogg et al., 2020) in periods where tree ring data are not available. Measurements of ¹⁴CO in the modern atmosphere have been used to constrain the oxidative capacity of the atmosphere (Brenninkmeijer et al., 1992; Petrenko et al., 2021), and thus, paleoatmospheric ¹⁴CO in ice cores can be used for a similar application. The in situ cosmogenic component of ¹⁴CO at ice core sites can be potentially be used to reconstruct the past cosmic ray flux (BenZvi et al., 2019). Finally, measurements of the in situ cosmogenic component of ¹⁴CO₂ and ¹⁴CO can be used to constrain the accumulation/ablation rate of the ice core site (e.g., Lal et al., 1990; Lal and Jull, 1990). Unfortunately, the paleoatmospheric and in situ cosmogenic components of ¹⁴C in ice exist in a combined form and cannot be separated analytically (Petrenko et al., 2016). To separate these signals, it is important to have accurate estimates of the cosmogenic ¹⁴C production rates and the partitioning among the in situ produced ¹⁴C species (¹⁴CO₂, ¹⁴CO, and ¹⁴CH₄) in ice.

In situ cosmogenic ¹⁴C production in ice is analogous to the production in quartz because both minerals share the same target atom (¹⁶O). Measurements of in situ cosmogenic nuclides (³He, ¹⁰Be, ¹⁴C, ²¹Ne, ²⁶Al, and ³⁶Cl) in near-surface rocks are commonly used as tools to constrain various Earth surface processes such as the timing of glacial retreat and erosion rates (Gosse and Phillips, 2001; Balco, 2020). Due to its short half-life of 5700 ±30 years (Kutschera, 2019), ¹⁴C in quartz is uniquely suited to characterize surface processes on millennial timescales (e.g., Spector et al., 2019; Pendleton et al., 2019). In situ cosmogenic ¹⁴C measurements are also often paired with measurements of longer-lived nuclides such as ¹⁰Be and ²⁶Al (e.g., Hippe, 2017; Skov et al., 2019) to study complex surface processes such as subglacial erosion and millennial-scale glacier retreats and/or re-advances.

In situ cosmogenic ¹⁴C in Earth's surface minerals is produced from ¹⁶O by three nuclear reactions, namely (1) nucleon (neutron and proton)-induced spallation (Lal and Peters, 1967), (2) negative muon capture (Heisinger et al., 2002b), and (3) interactions with fast muons (Heisinger et al., 2002a). The depth dependence of the ¹⁴C production rate for each mechanism in ice is shown in Fig. 1. Nucleon-induced spallation dominates the ¹⁴C production at the surface but is quickly attenuated with depth. On the other hand, the relative contributions from the two muon mechanisms are lower near the surface but dominate at larger depths, as muons can penetrate deeper than nucleons (Fig. 1). Characterizing the in situ cosmogenic ¹⁴C production rates from muons is especially important for applications of cosmogenic surface exposure dating where the samples might be exposed to the subsurface cosmic ray flux for an extended period. One example of this would be a bedrock that is covered by a relatively thin (e.g., tens of meters) glacier.

Figure 1 (a) In situ cosmogenic ¹⁴C production rates scaled for the Taylor Glacier study site (77^∘44^′ S, 162^∘10^′ E; 526 m elevation) from the three nuclear mechanisms, i.e., nucleon-induced spallation (P_n), negative muon capture (P_neg), and fast muon interactions (P_fast). (b) Depth profiles of the three production mechanisms normalized to their respective surface production rates (the respective surface production rates are shown in the legend). For ¹⁴C production from neutron spallation, we used the surface production rate estimate from Young et al. (2014) with scaling from the Lifton et al. (2014) LSDn model. For the two muon mechanisms (negative and fast muons), we used the production rate model from Balco et al. (2008), which follows parameterizations by Heisinger et al. (2002a, b).

Understanding the muogenic ¹⁴C component is especially important for ¹⁴C studies in ice. Prior studies have shown that, at snow accumulation sites, most of the in situ ¹⁴C produced in the firn (including the majority of neutron-produced ¹⁴C) is lost to the atmosphere via gas movement in the firn open porosity (Petrenko et al., 2013; Van der Kemp et al., 2000; Wilson and Donahue, 1990). In situ cosmogenic ¹⁴C mainly starts to accumulate in deeper ice where gas exchange with the atmosphere no longer happens, and at these depths, the ¹⁴C production is entirely from the muon mechanisms. Thus, the in situ cosmogenic ¹⁴C signal in traditional deep ice cores is dominated by production from muons, and constraining the muogenic ¹⁴C production rates is critical to disentangle the in situ cosmogenic and atmospheric ¹⁴C signals in ice cores. Unfortunately, the in situ ¹⁴C production rates by muons in both ice and quartz are still highly uncertain (Hippe, 2017).

The production rates of cosmogenic nuclides are usually determined from calibration sites where independent controls on exposure history are available, such as ¹⁴C dating from organic materials (e.g., Lifton et al., 2015) or argon (${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$) dating from lava flows (e.g., Balbas and Farley, 2020; Fenton et al., 2019). However, the commonly used estimates of muogenic ¹⁴C production rates (for both negative muon capture and fast muon reactions) were derived through the laboratory irradiation of artificial target compounds (Heisinger et al., 2002a, b). To our knowledge, there is only one prior study (Lupker et al., 2015) that provided estimates of total muogenic in situ ¹⁴C production rates based on measurements in a natural setting. Using ¹⁴C measurements from a 15.5 m deep quartzite core from Leymon High, Spain, Lupker et al. (2015) estimated a sea level high latitude (SLHL) surface production rate of 3.34 ($+\mathrm{0.43}/-\mathrm{1.07}$) ¹⁴C atoms per gram of quartz per year for negative muon capture and 0 ($+\mathrm{0.42}/-\mathrm{0.00}$) ¹⁴C atoms per gram of quartz per year for fast muon interactions (1σ uncertainties). The large uncertainties in the ¹⁴C production rates (especially the production rate from fast muons) estimated by Lupker et al. (2015) were due to relatively large measurement uncertainty for their deepest samples and a small contribution to the ¹⁴C signal from fast muons. Balco (2017) also refitted the Leymon High data and obtained similar results regarding the magnitude of the SLHL ¹⁴C production rate from negative muon capture.

Petrenko et al. (2016) used ¹⁴C measurements (¹⁴CO, ¹⁴CO₂, and ¹⁴CH₄) in >50 ka ice for the 2–20 m depth range from Taylor Glacier, Antarctica, to constrain the ¹⁴C production rates in ice. The old age of the ice ensured that all in situ cosmogenic and paleoatmospheric ¹⁴C inherited from the ice accumulation site had decayed away. Unfortunately, Petrenko et al. (2016) were unable to accurately constrain the total ¹⁴C production rates because of the high uncertainty resulting from the melt–extraction technique used to obtain their ¹⁴CO₂ measurements (see Sect. 1.3).

1.2 Overview of ¹⁴C production from muons

Following Heisinger et al. (2002b), the production rate of ¹⁴C (atoms per gram per year) by negative muon capture (P_neg) as a function of lithospheric depth (h; typically in g cm⁻²) is given by the following:

$$\begin{array}{}\text{(1)}& {\displaystyle }{P}_{\mathrm{neg}}\left(h\right)=R_{\mathit{\mu }-}\left(h\right)\cdot f_{\mathrm{tot}}\text{(2)}& {\displaystyle }{f}_{\mathrm{tot}}=f_{\mathrm{C}}\cdot f_{\mathrm{D}}\cdot f^{*},\end{array}$$

where R_μ−(z) is the stopping rate of negative muons (muons per gram per year) at lithospheric depth h, and f_tot is the overall probability of ¹⁴C production in ice from a stopped negative muon (unitless). The stopping rate of negative muons at the given depth R_μ−(h) has been empirically determined from measurements at deep underground laboratories (Heisinger et al., 2002b). The lithospheric depth (h) is a product of actual depth (z) and density (ρ) of the target mineral (ρ_ice=0.92 g cm⁻³).

The total probability (f_tot) of ¹⁴C production from negative muon capture is expressed by the product of the chemical compound factor (f_C), representing the probability that the stopped muon is captured by one of the target atoms (¹⁶O in case of ¹⁴C production), the probability that the negative muon does not decay in the K shell before nuclear capture (f_D), and the effective probability for the production of cosmogenic nuclides after μ- capture by the target atom (f^*; Eq. 2; Heisinger et al., 2002b; Lupker et al., 2015). All probability (f) terms in Eq. (2) are unitless. From experiments involving laboratory irradiation of artificial targets, the overall probability (f_tot) for ¹⁴C production in ice from negative muon was estimated to be 0.025±0.002 (Heisinger et al., 2002b).

An expression for the production rate of nuclides by fast muon interactions (P_fast) as a function of lithospheric depth (h) is given by Heisinger et al. (2002a):

$$\begin{array}{}\text{(3)}& {\displaystyle }{P}_{\mathrm{fast}}\left(h\right)={\mathit{\sigma }}_{\mathrm{0}}\cdot \mathit{\beta }\left(h\right)\mathit{\varphi }\left(h\right)\cdot \stackrel{\mathrm{‾}}{E}\left(h\right)\mathit{\alpha }\cdot N\text{(4)}& {\displaystyle }\mathit{\beta }\left(h\right)=\mathrm{0.846}-\mathrm{0.015}\ln (h+\mathrm{1})+\mathrm{0.003139}\left(\ln \right(h+\mathrm{1}){)}^{\mathrm{2}},\end{array}$$

where ϕ(h) is the total muon flux at depth z (muons per centimeter squared per year per steradian), σ₀ is the reference nuclear reaction cross section at muon energy of 1 GeV (millibarn, mb), β(h) is the unitless parameterized depth-dependence factor (Eq. 4), $\stackrel{\mathrm{‾}}{E}\left(h\right)$ is the mean muon energy at depth h (GeV), α is a power factor that describes the energy dependence of the cross section (unitless; α=0.75), and N is the number of target nuclei per gram of the target mineral. The overall production rate of ¹⁴C from fast muons provided by Heisinger et al. (2002a) has a high (±50 %) uncertainty because of the uncertainty in the reference nuclear reaction cross section σ₀ (${\mathit{\sigma }}_{\mathrm{0}}=\mathrm{0.0088}\pm \mathrm{0.0049}$ mb). Following Lupker et al. (2015), in this study we used f_tot and σ₀ as tuning parameters for the two muogenic production mechanisms in a cosmogenic nuclide production model (Sect. 3.2) to fit our ¹⁴C measurements.

1.3 Gas extraction methods for ice core ¹⁴C analysis

Common methods to liberate gas trapped in ice core bubbles include melting (wet extraction; e.g., Sowers et al., 1992; Mitchell et al., 2011) and mechanical destruction of the ice lattice (dry extraction; e.g., Bereiter et al., 2013; Ahn et al., 2009; Zumbrunn et al., 1982). Dry extraction is generally preferable for CO₂ analysis because the presence of liquid water in a wet extraction introduces extraneous CO₂ from the carbonate–acid reaction between the meltwater and impurities in the ice (e.g., Delmas et al., 1980; Raynaud et al., 1982). Multiple studies of ¹⁴CO₂ in ice have used dry extraction methods (e.g., Van De Wal et al., 1994, 2007; Smith et al., 2000; Van der Kemp et al., 2000). However, dry extraction systems (e.g., Lüthi et al., 2008) can potentially introduce biases in the CO₂ mole fraction [CO₂] due to incomplete gas extraction (Bereiter et al., 2015). Considering that the in situ cosmogenic production of ¹⁴C occurs directly in the ice lattice (Lal et al., 1990), it has been argued that dry extraction may also not liberate all of the ¹⁴C from the ice (e.g., van Roijen et al., 1994).

Other studies of ¹⁴C in ice (e.g., Lal et al., 1990, 1997, 2001; Jull et al., 1994) have used wet extraction methods. These wet extraction studies involved an addition of acid to drive off all dissolved CO₂ from the meltwater (Lal et al., 1990, 1997, 2001; Jull et al., 1994). The acidification process may have resulted in an additional CO₂ release from impurities in the ice (e.g., carbonate dust). In dust-rich Greenland ice, the presence of liquid water in a wet extraction produced excess CH₄ artifacts (Lee et al., 2020). It is thus possible that a wet extraction approach for ¹⁴C analysis may also result in additional C release from organics in the ice which are not ¹⁴C free.

A third method to liberate gases trapped in ice cores is sublimation under vacuum (e.g., Wilson and Donahue, 1989, 1990; Wilson and Long, 1997; Siegenthaler et al., 2005; Schmitt et al., 2011). Sublimation can occur when the pressure and temperature on the surface of the ice are below the triple point of the water-phase-change diagram. In addition to being free of problems associated with wet extraction methods, sublimation guarantees 100 % gas extraction efficiency (Schmitt et al., 2011; Bereiter et al., 2013, 2015), which includes any ¹⁴C trapped in the ice lattice. Therefore, sublimation is likely an optimal method for ¹⁴CO₂ measurements in ice.

This study presents new ¹⁴C measurements in three gas species (¹⁴CO, ¹⁴CO₂, and ¹⁴CH₄) in ancient (>50 ka) ice from the ablation zone of the Taylor Glacier, Antarctica, to constrain the compound-specific ¹⁴C production rates in ice by muons. Ice at this location does not contain a significant amount of ¹⁴C inherited from the accumulation site (Petrenko et al., 2016), and the ¹⁴C content is entirely due to production by muons during transport within the glacier. We improved on the earlier work by Petrenko et al. (2016) by (1) using a newly developed ice sublimation extraction device for ¹⁴CO₂ measurements (see Sect. 2.3.2), (2) collecting deeper samples to ∼72 m to better characterize the ¹⁴C production rate from the fast muon mechanism, and (3) using a more realistic 2D ice flow model from Buizert et al. (2012) to account for the flow trajectory and exposure history of the samples (see Sect. 3.1).

2 Field sampling and analytical methods

2.1 Site description

The blue ice area of the Taylor Glacier (Fig. 2) provides access to near-unlimited amounts of well-dated ancient ice (Baggenstos et al., 2017; Bauska et al., 2016; Menking et al., 2019; Schilt et al., 2014; Shackleton et al., 2020). This allows ice from Taylor Glacier to be measured for ultra-trace-gas species that require a very large amount of ice (Dyonisius et al., 2020; Petrenko et al., 2016, 2017; Buizert et al., 2014). In this study, we used the same site as Petrenko et al. (2016; 77^∘43.699^′ S, 161^∘43.179^′), where ice >50 ka in age at the surface has been previously identified.

2.2 Field sampling

Approximately 1000 kg of ice is needed to obtain both the necessary CH₄-derived and CO-derived C mass for ¹⁴C analyses. Because of this large sample requirement, and to avoid post-coring in situ ¹⁴C production at the surface, the melt extraction for ¹⁴CH₄ and ¹⁴CO samples was performed on site using the large-volume melter apparatus and technique described in Petrenko et al. (2016). The liberated air was transferred to 34.9 L electropolished stainless steel canisters and shipped to our laboratories for processing and analyses. Similar to other studies using this large-volume ice melter (e.g., Dyonisius et al., 2020; Petrenko et al., 2016, 2017), four procedural blanks (two with modern ¹⁴CH₄ standard gas and two with ¹⁴C dead ¹⁴CH₄ standard gas) were collected in the field. These field procedural blanks allow us to characterize the addition of extraneous ¹⁴C to the samples. The standard gases used in the field procedural blanks were passed through a Sofnocat 423 reagent, which removes CO (and thus ¹⁴CO) but leaves CH₄ (and ¹⁴CH₄) intact.

Figure 2 Map of the Taylor Glacier study site. The sampling location is marked by a red cross on the map. The orange star in the inset map shows the location of the Taylor Glacier relative to the Antarctic continent. The map was made using ArcGIS Pro, with imagery layers from Esri and Earthstar Geographics.

The sampling scheme for this study is shown in Fig. S1 in the Supplement. We used the 9.5 in. (241 mm) diameter Blue Ice Drill (BID; Kuhl et al., 2014) to collect seven large-volume samples during the 2015/2016 austral summer field season for ¹⁴CO and ¹⁴CH₄ analyses. The surface sample was collected from 21×1.5 m deep shallow cores, each with an average mid-depth of ∼0.75 m. Six additional deep samples with mid-depths of 19.5, 30, 40.5, 51, 61.5, and 72 m were also collected by combining ice from three ∼78 m deep boreholes. Each of the deep large-volume samples spanned approximately 10.5 m depth. Continuous sticks of ice subsamples (3×3 cm; spanning the whole length of the core) were taken from one of the three ice core boreholes (TGDeep3) for age control (see Sect. S3 in the Supplement). The continuous sample sticks were measured for the CH₄ mole fraction [CH₄] using the continuous flow analysis (CFA) system described in Rhodes et al. (2013) at Oregon State University (OSU).

In addition to the large-volume samples, we collected 26 smaller subsamples (∼1.5–2 kg) from 13 depth levels and two boreholes for ¹⁴CO₂ measurements. Each depth level contained a pair of replicates; however, only 9 out of the 13 replicate pairs were true replicates (i.e., collected from the same borehole and cut from the same depth interval; Fig. S1). Collecting same depth-adjacent samples below 50 m depth from a single borehole was challenging because of the reduced core quality (i.e., more fractures in the ice), and thus, the replicates had to be collected from a different borehole. Immediately after removal from the borehole, ice samples become exposed to a more intense cosmic ray bombardment (post-coring in situ cosmogenic ¹⁴C production). Five artificial bubble-free ice (BFI) samples were manufactured in the field, following methods from Mitchell et al. (2011), but upscaled to produce 1.5–2 kg samples. The field-produced BFI samples were shipped together with the collected glacial ice samples to characterize the effects of the post-coring in situ cosmogenic ¹⁴CO₂ production in the samples.

2.3 Laboratory analytical methods

2.3.1 Large-volume samples for ¹⁴CO and ¹⁴CH₄ measurements

The detailed approach for sample processing, measurements, and associated procedural corrections for the large-volume samples have been previously described in detail (Petrenko et al., 2016). In this section, we only provide a brief overview and highlight the differences between our methods and those of Dyonisius et al. (2020). First, the δ¹³CH₄ measurements were conducted at the Institute of Arctic and Alpine Research (INSTAAR), following the methods described by Miller et al. (2002; Table S1 in the Supplement). The δ¹³CH₄ measurements were not corrected for gravitational (Sowers et al., 1992) and diffusive isotopic fractionation (Buizert et al., 2013) because these corrections are only necessary to reconstruct the paleoatmospheric δ¹³CH₄ signal. In this study, the δ¹³CH₄ values are only used to normalize and calculate the absolute ¹⁴CH₄ abundance (in molec. per gram of ice).

The large-volume samples and field procedural blanks were measured for [CH₄] using a gas chromatograph–multidetector (GC-MD) system (Prinn et al., 2008; Table S2). Pressure in the sample canisters was measured using a Paroscientific, Inc., Digiquartz^® series 740 absolute pressure transducer at the Scripps Institution of Oceanography (SIO) for total air content (TAC) determination (Table S3). Two of the field procedural blanks were also measured for the $\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$, $\mathrm{Xe}/{\mathrm{N}}_{\mathrm{2}}$, and $\mathrm{Xe}/\mathrm{Kr}$ ratio (Table S4) at the Scripps Institution of Oceanography (SIO), following the procedures described in Bereiter et al. (2018). The noble gas ratios were used to constrain the degree of gas solubility during the melt extraction. The large-volume samples were measured for CO mole fraction [CO] using a Picarro G2401 analyzer (Table S5) and again for pressure at the University of Rochester (UR; Table S4).

The CH₄ in the large-volume samples and blanks was combusted to CO₂, cryogenically separated, and flame-sealed in glass ampules using the air processing line at the University of Rochester (Dyonisius et al., 2020). We also processed 3×100 µg of CH₄-derived C samples each from the modern ¹⁴CH₄ standard gas and ¹⁴C dead standard gas used for the field procedural blanks. The sample air that remained after CH₄ processing (∼10 L standard temperature pressure, STP) was diluted with a gas containing 10.02±0.26 µmol mol⁻¹ (95 % confidence interval, CI) of ¹⁴C-depleted CO (¹⁴CO $=\mathrm{0.19}\pm \mathrm{0.08}$ pMC; 95 % CI) to increase the CO-derived C mass for the accelerator mass spectrometry (AMS) measurements. The dilutant gas was measured for δ¹³CO using methods described in Vimont (2017; δ¹³CO $=-\mathrm{23.36}\pm \mathrm{0.2}$ ‰; 95 %CI).

The CO- and CH₄-derived CO₂ was graphitized using the Australian Nuclear Science and Technology Organization (ANSTO) micro-furnaces, following Yang and Smith (2017). We used the ¹⁴C activity measured on the 100 µg samples as the true ¹⁴C activity of the standard gases (Table S6). Because of the larger sample size, the effect of extraneous C introduced by graphitization on these 100 µg samples is assumed to be negligible. Using a mass balance approach described in Petrenko et al. (2017), the total extraneous C mass for the ¹⁴CH₄ samples was determined to be 0.63±0.28 µgC, and the corresponding ¹⁴C activity for the extraneous C was 16.7±10.2 pMC (95 % CI).

Table 1 Measured ¹⁴CO₂, ¹⁴CO, and ¹⁴CH₄ after all associated corrections and calculated total ¹⁴C, ${}^{\mathrm{14}}{\mathrm{CH}}_{\mathrm{4}}{/}^{\mathrm{14}}\mathrm{CO}$ ratios, ¹⁴CO₂, and ¹⁴CO fractions. All errors presented indicate the 95 % CI. NA – not available.

^a Data from Petrenko et al. (2016). ^b The total ¹⁴C value for the 10 m sample was scaled by a factor of 1.003±0.003 (95 % CI) to account for the lack of ¹⁴CH₄ measurements (Sect. 4.2).

Figure 3 (a) Measured ¹⁴CO molec. per gram of ice after all corrections. (b) Measured ¹⁴CH₄ molec. per gram of ice after all corrections. (c) Measured ¹⁴CO₂ molec. per gram of ice after all corrections. (d) Total measured ¹⁴C atoms per gram of ice. This represents the sum of ¹⁴CO, ¹⁴CH₄, and sublimation-based ¹⁴CO₂ measurements. All error bars shown in this figure are 95 % CI.

In prior studies (e.g., Dyonisius et al., 2020; Petrenko et al., 2017), ¹⁴CO measurements from the field procedural blanks were used to characterize the effects of extraneous ¹⁴C addition from sample extraction, handling, storage, transport, and processing (including the graphitization step). For this study, the field procedural blanks were still used to characterize the effects from in situ production of ¹⁴CO in the sample air canisters by cosmic rays during storage and transport. However, to better characterize the effects from the addition of extraneous C during the graphitization process, we used a linear empirical correction from 10 commensurately sized ¹⁴C standards and blanks at ANSTO (see Fig. S2a; Table S7), following Petrenko et al. (2021). This approach has the benefit of bracketing the effects of extraneous C from graphitization at ANSTO with low and high ¹⁴C standards, similar to the approach for the ¹⁴CH₄ samples. The ¹⁴CO blank for this sample set is 22.45±3.24 molec. ¹⁴CO per cc of air STP (95 % CI), which is higher than the ¹⁴CO blanks reported in Dyonisius et al. (2020). This is mainly because there was an extra year between the retrieval and processing of the samples (thus, there was more in situ ¹⁴CO production in sample canisters during storage). The amount of ¹⁴C  per gram of ice for ¹⁴CO, ¹⁴CH₄, and ¹⁴CO₂ (Table 1) is calculated using the same method as in Petrenko et al. (2016) and is consistent with the Hippe and Lifton (2014) formulations for in situ ¹⁴C concentrations. ¹⁴CH₄ and ¹⁴CO measurements in our samples after all associated corrections, in addition to earlier Taylor Glacier results from Petrenko et al. (2016), are shown in Table 1 and Fig. 3.

2.3.2 Sublimation and processing of samples for ¹⁴CO₂ measurements

CO₂ was liberated from ice samples using a newly developed ice sublimation device at the University of Rochester (Hmiel, 2020), roughly following the design of Schmitt et al. (2011). To briefly summarize the procedure, 1.5–2 kg ice samples were loaded into a vacuum glass vessel. The vessel was then evacuated, and the ice was sublimated at vacuum with six infrared emitters (Emitted Energy, USA) for 8–10 h. We did not sublimate 100 % of the samples because, as the ice sublimates away, impurities such as dust and organics start to accumulate on the surface. The aggregation of impurities on the sublimation front might enhance unwanted chemical reactions that produce extraneous carbon (Schmitt et al., 2011). Furthermore, towards the end of the extraction, the sublimation became less efficient as less surface area was available to absorb radiation. Approximately 1 kg of ice was sublimated in 8–10 h. However, the incomplete sublimation does not compromise the 100 % extraction efficiency, as all the gases trapped in the ice that is sublimated away are still released (Schmitt et al., 2011).

The liberated CO₂ was cryogenically trapped with liquid nitrogen, and the air was also cryogenically trapped with 5Å molecular sieve (Sigma-Aldrich, USA) under liquid nitrogen. After the sublimation was completed, the trapped CO₂ and air were expanded into separate volume-calibrated manometers in which pressure measurements were taken to calculate the [CO₂]. Finally, the isolated CO₂ was cryogenically transferred to and flame-sealed into a Pyrex glass ampule. The CO₂ was graphitized at ANSTO using the micro-furnaces (Yang and Smith, 2017), and the graphitized samples were measured for ¹⁴C activity at the ANTARES AMS facility (Smith et al., 2010). One ¹⁴CO₂ sample (replicate for 30 m depth sample) was unfortunately lost during sublimation because the ice fractured under vacuum during the evacuation step.

A ∼50–75 g ice subsample was taken from every ¹⁴CO₂ sample and shipped to OSU. The aliquots were measured for [CO₂], following Ahn et al. (2009), and [CH₄] and TAC, following Mitchell et al. (2013; Table S8). Five field-produced bubble-free ice (BFI) samples and nine laboratory-produced BFI samples were also sublimated along with the glacial ice samples. During the sublimation of the BFI samples, a standard gas with known ¹⁴CO₂ activity and [CO₂] was introduced into the bottom of the glass sublimation vessel at 0.15 scc min⁻¹ flow rate for 8–10 h. The flow rate was set to mimic the rate of air liberation from glacial ice samples and the processing time also mimicked the amount of time needed to sublimate glacial ice samples. We used a standard gas with dead ¹⁴CO₂ activity for four laboratory-produced BFI samples and a standard gas with modern ¹⁴CO₂ activity for the other five laboratory-produced BFI samples. The CO₂ was cryogenically trapped downstream, processed, and measured for ¹⁴C activity, following the same methods as the ice samples. In combination with the OSU [CO₂] and TAC measurements, the BFI samples were used to constrain the amount of extraneous carbon and ¹⁴C introduced by sample transport, storage, and processing (see Sect. S1; Table S9). Finally, 11 commensurately sized ¹⁴C standards and blanks (14–16 µgC) with known ¹⁴C activities (in 0–135 pMC range) were prepared, graphitized, and measured at ANSTO, concurrently with all the samples (Table S7), to characterize the effects from the addition of extraneous C during the graphitization process.

The detailed corrections for the ¹⁴CO₂ samples are discussed in the Supplement. We correct for the effects of extraneous C from graphitization and other ANSTO processing using a linear empirical correction from the commensurately sized ¹⁴C standards (Sect. S1.1; Fig. S2b). The effects of extraneous carbon from ice sublimation/CO₂ extraction are calculated from the difference in measured ¹⁴C activity of the laboratory-produced BFI samples relative to the measured ¹⁴C activity of the standard gases with a mass balance approach (Sect. S1.2; Tables S8 and S9). Finally, the samples were corrected for the effects of post-coring in situ ¹⁴CO₂ production in ice using results from the field-produced BFI samples (Sect. S1.3; Table S10). The ¹⁴CO₂ measurements in our samples after all associated corrections with their error-propagated uncertainties are shown in Table 1 and Fig. 3.

Figure 4 (a) ${}^{\mathrm{14}}{\mathrm{CH}}_{\mathrm{4}}{/}^{\mathrm{14}}\mathrm{CO}$. (b) ¹⁴CO/total ¹⁴C and ¹⁴CO₂/total ¹⁴C ratio. The solid lines represent the mean, and the dashed lines represent 2 standard deviations of the ${}^{\mathrm{14}}{\mathrm{CH}}_{\mathrm{4}}{/}^{\mathrm{14}}\mathrm{CO}$ ratio, ¹⁴CO₂, and ¹⁴CO fractions for samples deeper than 6.85 m where production by muons dominates. The ratios for rejected samples (Sects. 4.1, S2) are not shown. All error bars shown in this figure are 95 % CI.

An in-depth discussion about the analytical uncertainty in the ¹⁴CO₂ measurements obtained from the sublimation method (which is important for the interpretation of the data because it is the largest source of uncertainty for total ¹⁴C) is provided in Sect. S1.4. In brief, we used the pooled standard deviation of replicate pairs (±7.8 ¹⁴CO₂ molec. per gram of ice; 2σ) as the uncertainty for all ¹⁴CO₂ measurements, except the 2.25 m sample pair (where we used the error-propagated uncertainties instead; Table 1). For the rest of the paper, we refer to the sum of measured ¹⁴CO, ¹⁴CO₂, and ¹⁴CH₄ as the total ¹⁴C. The ${}^{\mathrm{14}}{\mathrm{CH}}_{\mathrm{4}}{/}^{\mathrm{14}}\mathrm{CO}$, ¹⁴CO/total and ¹⁴CO₂/total ¹⁴C ratio of the samples are shown in Fig. 4.

2.4 Sample integrity

Several samples were excluded from the data analysis; the detailed reasoning for rejecting these samples is discussed in Sects. S2 and S3. The surface samples (0.75 m depth) for all three ¹⁴C species (¹⁴CO, ¹⁴CH₄, and ¹⁴CO₂) are rejected because of ambient air contamination from abundant fractures in the ice sample (due to thermal stresses in near-surface ice) and the likely chemical and/or biological alteration of CH₄, CO, and CO₂. The 19.5 and 30 m ¹⁴CO and ¹⁴CH₄ samples from the 2015/2016 season were rejected because of anomalous alterations in TAC, [CO], and [CH₄]. The 30 m ¹⁴CO₂ sample was also rejected due to alteration in [CO₂].

3 Estimating the muogenic ¹⁴C production rates

3.1 Ice flow model to constrain sample exposure history

We used a 2D ice flow model from Buizert et al. (2012) to generate flow trajectories for the ice parcels corresponding to each sample depth (Fig. 5). The model first computes the 2D steady-state ice flow velocity field based on the observed surface velocities, ablation rates, and glacier thickness and then generates an ice parcel back-trajectory using a 2D linear interpolation of the ice flow velocity field (Kavanaugh et al., 2009a, b; Kavanaugh and Cuffey, 2009; Bliss et al., 2011). The largest source of uncertainty for the trajectories are the ablation rates (Buizert et al., 2012), which are based on measurements of 163 poles initially planted in 2002/2003 (Kavanaugh et al., 2009b; Bliss et al., 2011). All survey poles were measured a year later, providing 1-year average ablation rate estimates (Kavanaugh et al., 2009b) and again in the 2006/2007 season. Additionally, 17 poles were remeasured in the 2009/2010 and 2010/2011 seasons (Buizert et al., 2012). The 4-year average ablation rates were 4.7 cm yr⁻¹ higher than the 1-year average (Kavanaugh et al., 2009b); following Buizert et al. (2012), the ablation rate uncertainty for each pole was calculated by dividing 4.7 cm yr⁻¹ with $\sqrt{N}$, where N is the length of the observation period in years (N=1, 4, 7, or 8). Figure S3 shows the ablation rates along the glacier and their uncertainties inferred from the survey pole data.

The ice flow model used a bedrock profile from Kavanaugh et al. (2009a); however, the bedrock profile only extends to 72 km away from the glacier terminus, a point which we refer to as the glacier head (Fig. 5). This bedrock profile corresponds to 5–6 kyr of ice flow history and approximately one ¹⁴C half-life. Beyond the constraints from the bedrock profile, we had to make assumptions about the depth of the long-term transport (z_deep). Morse et al. (1998) provided a radar-based bedrock profile that includes the Taylor Glacier snow accumulation area (Baggenstos et al., 2018) north of the Taylor Dome – approximately 60 km upstream from the glacier head, where the Kavanaugh et al. (2009a) bedrock profile ends. Based on the bedrock profile from Morse et al. (1998), at the Taylor Glacier accumulation area, the depth of ∼80 kyr ice (which corresponds to our 72 m sample) is ∼575 m. We thus assumed that the depth of long-term transport (z_deep) for the 72 m sample under the best-estimate ablation rate scenario (which we define as the reference sample) is 575 m. For other ice parcel trajectories (i), we scaled the depth of long-term transport (z_deep) as follows:

$$\begin{array}{}\text{(5)}& z_{\mathrm{deep}}\left(i\right)=\mathrm{575}-(z_{\mathrm{ref}}-z_{\mathrm{head}}(i\left)\right),\end{array}$$

where z_ref represents the depth of the 72 m reference sample in the model at the glacier head under the best-estimate flowline (z_ref is 699 m), and z_head represents the depth of the ice parcel of interest at the glacier head. We assumed that the difference in depth between the reference sample and the sample of interest (i) at the glacier head and during long-term transport within the glacier is the same.

Figure 5 Example of the ice parcel back-trajectory and associated uncertainties. For the Monte Carlo estimate of uncertainties (Sect. 5.4.2), for each given sample depth (72 m in this figure), 10 000 back-trajectories are generated. Each back-trajectory corresponds to a different ablation rate scenario (the ablation rates are perturbed within their experimental measurement uncertainties to generate the scenarios). The shaded region represents the 68 % CI uncertainty envelope of the flow trajectory.

Figure 6 (a) Comparison between total ¹⁴C measurements with modeled best-estimate σ₀ and f_tot parameters from this study and Heisinger et al. (2002a, b). (b) Comparison between total ¹⁴C measurements with modeled best-estimate σ₀ and f_tot parameters from this study and modeled total ¹⁴C with best-fit f_tot when σ₀ is forced to be zero. (c) Comparison between ¹⁴CO measurements with modeled best-estimate f_neg and f_fast parameters from this study and modeled ¹⁴CO with best-fit f_neg when f_fast is forced to be zero. (d) Comparison between total ¹⁴C measurements with modeled best-estimate σ₀ and f_tot parameters from this study and modeled total ¹⁴C from the sensitivity analyses when we assume 25 % contribution from organics and high ablation rate scenario (Fig. S9). The thin colored lines represent the 95 % CI envelope of the model results (corresponding to the contour plot of Fig. 7a for Fig. 6a and b and the contour plot of Fig. 7b for Fig. 6c). The error bars shown in the data are 95 % CI. In Fig. 6d, the solid blue line represents the ¹⁴C profile from the modeled best-estimate f_tot and σ₀ under the best-estimate ablation rate. The solid red line represents the ¹⁴C profile from the sensitivity analysis when f_tot and σ₀ are tuned to fit the total ¹⁴C data that are scaled by 25 % to account for contribution from organics (red dots) under high ablation rate scenario (Fig. S9). The dashed red line represents the ¹⁴C profile from the sensitivity analysis when f_tot is kept constant at 0.021 (which is the minimum value provided by Heisinger et al., 2002b) and σ₀=0 (which provides the best-fit against the data).

3.2 ¹⁴C production in sample ice parcel

We used the model for in situ cosmogenic nuclide production by muons from Balco et al. (2008) and Balco (2017; model 1A in Balco, 2017), with all relevant parameters adjusted for ice (Fig. 1). This model in turn uses the Heisinger et al. (2002a, b) parameterizations described above and the additional altitude scaling of the muon fluxes described in Balco et al. (2008). We then used a forward model that numerically integrates the total ¹⁴C in the ice sample along its flow path. For the initial condition, we assumed that, at the depth of long-term transport (z_deep), the ¹⁴C concentration in the ice parcel is at the steady state, as follows:

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\mathrm{d}C}{\mathrm{d}t}}\left(\mathrm{at}{z}_{\mathrm{deep}}\right)=\mathrm{0}=P_{\mathrm{neg}}\left(z_{\mathrm{deep}}\right)+P_{\mathrm{fast}}\left(z_{\mathrm{deep}}\right)-C_{\mathrm{0}}\mathit{\lambda }.\end{array}$$

The steady-state assumption means that, at z_deep, the rate of radioactive decay (C₀λ) is balanced by production from negative muon capture (P_neg) and fast muon reaction (P_fast). We use a ¹⁴C decay constant $\mathit{\lambda }=\mathrm{1.216}\times {\mathrm{10}}^{-\mathrm{4}}$ yr⁻¹ for all our calculations (including the conversion of ¹⁴C units in the previous section), which corresponds to the recommended ¹⁴C half-life of 5700 years (Kutschera, 2019). For each ice parcel, we calculated the steady state, the initial ¹⁴C concentration (C₀) from Eq. (6), and then used the following differential equation,

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\mathrm{d}C}{\mathrm{d}t}}=P_{\mathrm{neg}}\left(z\left(t\right)\right)+P_{\mathrm{fast}}\left(z\left(t\right)\right)-C\mathit{\lambda },\end{array}$$

to numerically integrate the ¹⁴C concentration of the ice parcel along the flow trajectory. To avoid interference from spallogenic ¹⁴C, we only considered samples with mean depths greater than 6.85 m depth.

We sampled the parameter space in a grid search approach to obtain the best-estimate values for muogenic ¹⁴C production parameters σ₀ and f_tot. Using the best-estimate flow trajectory, we calculated the expected ¹⁴C in the samples corresponding to all combinations of σ₀ and f_tot, with each of the parameters ranging between 0 % and 100 % of the values from Heisinger et al. (2002a, b). To save computational time, we first conducted the grid search at a coarse resolution of 10 % increments (Fig. S4a). The goodness of the fit (χ2) for each simulation was calculated as follows:

$$\begin{array}{}\text{(8)}& {\mathit{\chi }}^{\mathrm{2}}=\sum {\displaystyle \frac{(C_{\mathrm{obs}\left(z\right)}-C_{\exp }(z){)}^{\mathrm{2}}}{C_{\exp }\left(z\right)}},\end{array}$$

where C_obs(z) is the measured total ¹⁴C, and C_exp(z) is the total ¹⁴C (¹⁴CO₂ + ¹⁴CO + ¹⁴CH₄; Fig. 3d) calculated by the forward model at sample depth z. To find more precise best-estimate σ₀ and f_tot, we conducted the grid search again at a higher resolution in 0.2 % increments from the Heisinger et al. (2002a, b) values near the χ2 minimum, between 0 and 0.0352 mb for σ₀ and 0 to 0.01 for f_tot (Fig. S4b).

To estimate the uncertainties in σ₀ and f_tot, we used a Monte Carlo sampling of model parameters. We assumed that the ablation rate uncertainties (Fig. S3) represent 2σ normally distributed uncertainties. We then perturbed the ablation rates within their uncertainties and generated a pool of 10 000 possible flow trajectories for each sample depth. However, in 69 out of 10 000 flow scenarios, the ice parcel back-trajectories hit the bedrock and became unphysical afterwards. These unphysical trajectories were removed from the pool of possible ice flow trajectories. Next, we started with the best-estimate σ₀ and f_tot and assumed a normally distributed and large 200 % (1σ) error for each parameter (Fig. S5a) as the prior distribution for the Monte Carlo method. We removed σ₀ and f_tot values that are below zero from the prior distribution because they are unphysical and conducted 100 000 Monte Carlo simulations using the forward ¹⁴C production model. For each Monte Carlo simulation, we randomly picked one of the previously generated possible ice flow trajectories and a random pair of σ₀ and f_tot from the generated prior distributions (Fig. S5a). We then calculated the expected ¹⁴C concentrations for each sample depth using the forward model and compared the model–data fit. We accept all pairs of σ₀ and f_tot values that produce the model-calculated total ¹⁴C within the 95 % CI (7.8 ¹⁴C atoms per gram of ice) and 67 % CI (3.9 ¹⁴C atoms per gram of ice) analytical uncertainty of the best-fit, model-calculated total ¹⁴C (black line in Fig. 6). The ranges of accepted σ₀ and f_tot pairs are shown in Fig. 7a as contours. The discussion about the selection of acceptance criteria for estimating σ₀ and f_tot uncertainties is provided in Sect. S1.4.

Figure 7 (a) The 68 % and 95 % CI contours of accepted σ₀ and f_tot values for total ¹⁴C. (b) The 68 % and 95 % CI contours of accepted f_neg and f_fast values for ¹⁴CO (see Sect. 5.4.3). For comparison, the σ₀ and f_tot values from Heisinger et al. (2002a, b) are shown as a blue star. The best-fit values for σ₀, f_tot, f_neg, and f_fast are shown as a red star in both figures. The best-fit values for σ₀ and f_tot are when α (energy scaling in Eq. 3).

3.3 ¹⁴CO production model in sample ice parcel

The in situ cosmogenic ¹⁴CO production rates in ice are of specific interest, as discussed in Sect. 1.1. To characterize the ¹⁴CO production rates, we introduced additional scaling factors f_neg and f_fast for negative muon and fast muon mechanisms, respectively, as tuned model parameters. The differential equation of Eq. (7) is modified into the following:

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{\mathrm{d}{(}^{\mathrm{14}}\mathrm{CO})}{\mathrm{d}t}}=f_{\mathrm{neg}}{P}_{\mathrm{neg}}\left(z\left(t\right)\right)+f_{\mathrm{fast}}{P}_{\mathrm{fast}}\left(z\left(t\right)\right)-{(}^{\mathrm{14}}\mathrm{CO})\mathit{\lambda }.\end{array}$$

We note that P_neg and P_fast in Eq. (9) are the total ¹⁴C production rates calculated from the Balco et al. (2008) model. The scaling factors f_neg and f_fast each encompass two terms, i.e., one that adjusts the total ¹⁴C production rates and another that accounts for the ¹⁴CO fraction of total ¹⁴C. The determination of best-estimate f_neg and f_fast and their uncertainties was similar to the approach for σ₀ and f_tot described above. The χ2 grid search was conducted with all combinations of f_neg and f_fast values ranging from 0 to 0.2 at 0.001 resolution (Fig. S4c). Similar to the total ¹⁴C data, we used the average analytical uncertainty in the ¹⁴CO sample set as the acceptance criteria for the Monte Carlo simulations to estimate the uncertainties of f_neg and f_fast. We accepted all sets of f_neg and f_fast from the 100 000 Monte Carlo simulations that yielded model-predicted ¹⁴CO within 1.2 ¹⁴CO molec. per gram of ice (95 % CI uncertainty) and 0.6 ¹⁴CO molec. per gram of ice (68 % CI uncertainty) from the best-fit model (Fig. 8). Figure 7b shows (as contours) the accepted sets of f_neg and f_fast values.

Figure 8 (a) Comparison between measured total ¹⁴C from Scharffenbergbotnen, expected total ¹⁴C using the production rates inferred in this study, and expected total ¹⁴C using the Heisinger et al. (2002a, b) production rates. (b) Comparison between measured ¹⁴CO from Scharffenbergbotnen and expected ¹⁴CO, using the production rates inferred from Taylor Glacier. The colored lines in both figures represent the 95 % CI envelope of the model results. At the depths plotted in this figure (deeper than 5 m), the production from neutron-induced spallation is negligible.

3.4 Comparison with Scharffenbergbotnen ablation site

Van der Kemp et al. (2002) measured ¹⁴CO₂ and ¹⁴CO in ice from the Scharffenbergbotnen ice ablation site, Antarctica. Using a 1D ablation model, we examined how the estimates of muogenic ¹⁴C production rates from Taylor Glacier compare to the Scharffenbergbotnen data. We assumed that the measured ¹⁴CO₂ plus ¹⁴CO from Scharffenbergbotnen are comparable to our measurements of total ¹⁴C in Taylor Glacier (since our data show that less than 0.3 % of the total ¹⁴C from muon production forms ¹⁴CH₄; Sect. 4.1). We then used the ¹⁴C concentration from the deepest Scharffenbergbotnen sample (45 m) as the initial condition. In the 1D ablation model, the Scharffenbergbotnen ice parcel moves upward at a rate ($\mathrm{d}z/\mathrm{d}t$) equal to the ablation rate from stake measurements (Eq. 10; $a=\mathrm{16}\pm \mathrm{4}$ cm yr⁻¹).

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{\mathrm{d}z}{\mathrm{d}t}}=-a\text{(11)}& {\displaystyle \frac{\mathrm{d}C}{\mathrm{d}t}}=P_n^{\prime }\left(z\left(t\right)\right)+{P^{\prime }}_{\mathrm{neg}}\left(z\left(t\right)\right)+{P^{\prime }}_{\mathrm{fast}}\left(z\left(t\right)\right)-C\mathit{\lambda }.\end{array}$$

The expected ¹⁴C concentration in the ice is given by the differential equation (Eq. 11), where $P_n^{\prime }$ is the ¹⁴C spallogenic production rate from Young et al. (2014; 12.0±0.9 atoms per gram of quartz per year at the surface), first scaled to the SLHL production rate in ice (20.0±1.5 atoms per gram of ice per year at the surface) accounting for the number of ¹⁶O atoms per gram of ice vs. quartz (variable N; Eq. 3; Petrenko et al., 2016) and then to the production rate at the Scharffenbergbotnen site (1173 m above sea level or m a.s.l.) using the Lifton et al. (2014) LSDn nuclide-specific model ($P_n^{\prime }=\mathrm{71.2}\pm \mathrm{3.6}$ atoms per gram of ice per year at the surface). $P_{\mathrm{neg}}^{\prime }$ and $P_{\mathrm{fast}}^{\prime }$ are the muogenic production rates inferred from the Taylor Glacier data scaled to the elevation of Scharffenbergbotnen (1173 m a.s.l.), using the altitude scaling factors from Balco (2017). We also repeated this calculation for ¹⁴CO only to compare the muogenic ¹⁴CO production rates with the ¹⁴CO data from Scharffenbergbotnen.

4 Results and discussion

4.1 Measured ¹⁴C values and partitioning of ¹⁴CO₂, ¹⁴CO, and ¹⁴CH₄

Table 1 and Fig. 3a–c show the depth profiles of ¹⁴CO, ¹⁴CH₄, and ¹⁴CO₂ after all corrections. For the ¹⁴CO₂ measurements, comparison with prior results that used a wet extraction approach (Fig. S6) confirms the caveats discussed by Petrenko et al. (2016) in that their ¹⁴CO₂ measurements were uncertain and represent the upper bound. The ${}^{\mathrm{14}}{\mathrm{CH}}_{\mathrm{4}}{/}^{\mathrm{14}}\mathrm{CO}$ ratios from the new samples (0.0074±0.0004; 95 % CI; n=4; from all samples below 19.5 m) appear to be constant within uncertainties (Fig. 4a), which is in agreement with earlier results (0.0076±0.0004; 95 % CI; n=4) from Petrenko et al. (2016). This confirms that the two muon reactions produce ¹⁴C in a constant ${}^{\mathrm{14}}{\mathrm{CH}}_{\mathrm{4}}{/}^{\mathrm{14}}\mathrm{CO}$ ratio. The ¹⁴CO and ¹⁴CO₂ fractions of total ¹⁴C are also relatively constant at depth (Fig. 4b), suggesting that the two muon reactions produce all three ¹⁴C species in constant ratios.

For samples deeper than 6.85 m, on average 33.7 % (±11.4 %; 95 % CI) of the produced cosmogenic ¹⁴C becomes ¹⁴CO and 66.1 % (±11.5 %; 95 % CI) of the produced cosmogenic ¹⁴C becomes ¹⁴CO₂ (Fig. b). The uncertainties of ¹⁴CO and ¹⁴CO₂ fractions on the deepest samples (72 m depth) are relatively large because of the small ¹⁴CO₂ signal (11.8 to 13.6 ¹⁴CO₂ molec. per gram of ice) relative to the uncertainty in our measurements (±7.8 ¹⁴CO₂ molec. per gram of ice; 95 % CI). The ¹⁴CO₂ fraction in samples that are deeper than 6.85 m (0.66±0.12; 95 % CI) is also in agreement with prior reported ¹⁴CO₂ fraction of 0.69 from the Scharffenbergbotnen ablation site (Van der Kemp et al., 2002). Finally, the shallow samples (<6 m ice equivalent) show higher ${}^{\mathrm{14}}{\mathrm{CH}}_{\mathrm{4}}{/}^{\mathrm{14}}\mathrm{CO}$ ratios (Fig. 4a) and ¹⁴CO₂ $/$ total ¹⁴C ratios. This may indicate that spallation produces higher amounts of ¹⁴CH₄ and ¹⁴CO₂ relative to ¹⁴CO (Petrenko et al., 2016) or that CO (and ¹⁴CO) is not well-preserved in near-surface ice of the Taylor Glacier.

At depths where production from muons dominates (>6 m), less than 0.3 % of the produced cosmogenic ¹⁴C in ice forms ¹⁴CH₄ (Table 1; Fig. 4). Although the ¹⁴CH₄ measurement from the 10 m depth sample is not available (Petrenko et al., 2016), we still include the 10 m data point in the total ¹⁴C dataset used to infer σ₀ and f_tot values and their uncertainties. The contribution from ¹⁴CH₄ (which would have been on the order of ∼0.2 ¹⁴CH₄ molec. per gram of ice; Fig. 3b) is insignificant compared to the uncertainty in total ¹⁴C. We account for the lack of ¹⁴CH₄ measurement at this depth by scaling the total ¹⁴C of the 10 m sample by a factor of 1.003±0.003 (95 % CI; Table 1).

Table 2 Probability of the ¹⁴C production from stopped negative muons (f_tot), reference nuclear reaction cross section for production via fast muon interactions (σ₀), and total ¹⁴C production rates in ice at the surface from the two muon reactions rescaled to SLHL (sea level high latitude), using the Lifton et al. (2014) LSDn scaling. All errors shown represent the 95 % confidence interval.

^* Adjusted to ice, assuming that the chemical compound factor (f_c) of ice is 1.0, and f_c for quartz is 0.704 (Heisinger et al., 2002b).

4.2 Inferred muogenic ¹⁴C production rates in ice and comparison with production rates in quartz

Assuming that the majority of in situ cosmogenic ¹⁴C in ice forms ¹⁴CO₂, ¹⁴CO, and ¹⁴CH₄, the muogenic ¹⁴C production parameters from Heisinger et al. (2002a, b; f_tot for negative muon capture and σ₀ for fast muon reaction) are well outside the confidence intervals of our measurements (Table 2; Figs. 6a and 7a). Using the larger uncertainty for ¹⁴CO₂ measurements obtained from step-by-step error propagation (Sect. S1.4; Fig. S7) does not change this conclusion. We calculated factors of 5.7 (3.6–13.9; 95 % CI) and 3.7 (2.0–11.9; 95 % CI) lower probability of the negative muon capture reaction (f_tot) and the reference cross section for the fast muon mechanism (σ₀) compared to the values given by Heisinger et al. (2002a, b).

One possible explanation for the disagreement between our results and those of Heisinger et al. (2002a, b) is that our ¹⁴C measurements (mostly either ¹⁴CO or ¹⁴CO₂, as ¹⁴CH₄ only constitutes <0.3 % of total ¹⁴C) might be incorrect. However, in the following we thoroughly explore this possibility and argue that it is very unlikely. Our ¹⁴CO measurements used a well-established analytical technique (e.g., Dyonisius et al., 2020; Hmiel et al., 2020; Petrenko et al., 2013, 2017, 2021). With regards to ¹⁴CO measurements in air, this analytical technique (Petrenko et al., 2021) yields comparable results to independent, atmospheric ¹⁴CO measurements from other research groups (e.g., Manning et al., 2005; Mak and Southon, 1998). We also have no reason to believe that there is a systematic loss of ¹⁴CO during the ice melting process. The ice melting (wet extraction) ensures that all ¹⁴C and CO are liberated. CO is not very soluble in water (the dissolved CO fraction at equilibrium in our system is on the order of 1 %), and we used the measured $\mathit{\delta }\mathrm{Xe}/\mathrm{Kr}$ (Table S4) to correct for the solubility effects for both [CO] and [CH₄]. The on-site field extraction within hours of the sample retrieval ensures that there is minimal post-coring gas loss. Finally, ice core and firn air ¹⁴CO measurements at Summit, Greenland, are consistent within uncertainties with muogenic ¹⁴CO production rates inferred from Taylor Glacier (Hmiel et al., 2020; Hmiel, 2020).

The sublimation technique used for our ¹⁴CO₂ measurements also ensures the 100 % extraction efficiency of gases (Schmitt et al., 2011). We tested the system performance by sublimating BFI (bubble-free ice) samples while adding standard gases with known ¹⁴CO₂ activities. As mentioned in Sect. 2.3.2, we used two standard gases with known ¹⁴CO₂ activities, i.e., one with modern ¹⁴CO₂ activity and the other with dead ¹⁴CO₂ activity, when sublimating the BFI samples. There is neither a significant alteration in the CO₂ mole fraction and ¹⁴CO₂ activity of both standard gases (Table S8) nor in the CO₂ mole fraction of the ice samples (Table S9), which suggests that the processes of sublimating ice and flowing gas through the system components do not result in loss of ¹⁴CO₂. Finally, we can also rule out the possibility of post-coring ¹⁴CO₂ loss. In a separate measurement campaign (Hmiel, 2020), we brought the sublimation system to Summit, Greenland, and sublimated the ice samples on site (within days from the time when the ice core was drilled). We compared the ¹⁴CO₂ from the on-site field sublimation with the ¹⁴CO₂ from depth-adjacent replicates sublimated at the University of Rochester laboratory and found that they are indistinguishable within uncertainty.

We argue that another strong indication that our measurements are robust is the good agreement with independent results from Van der Kemp et al. (2002). Van der Kemp et al. (2002) measured ¹⁴CO₂ and ¹⁴CO in ice from the Scharffenbergbotnen ice ablation site using a dry extraction technique. The total measured ¹⁴C values were significantly lower than the expected values, based on the stake-measured ablation rates and muogenic production rates based on laboratory irradiations of quartz targets (Heisinger et al., 2000a, b). Van der Kemp et al. (2002) initially hypothesized that the low extraction efficiency of dry mechanical extraction (which can result in an incomplete release of the in situ produced ¹⁴C from the ice grains) might be responsible for this discrepancy. However, we used a sublimation method for our ¹⁴CO₂ measurements and a melt extraction method for our ¹⁴CO measurements; both methods guarantee that all in situ cosmogenic ¹⁴C in the ice lattice is released. Figure 8 shows that the Scharffenbergbotnen data are consistent with the expected total ¹⁴C and ¹⁴CO from the Taylor-Glacier-derived production rates.

The good agreement in the ratio of ¹⁴C compounds (¹⁴CO₂ fraction is 0.66±0.12 in this study and 0.69 in Van der Kemp et al., 2002) suggests that our extraction methods and analytical techniques were not systematically losing ¹⁴CO or ¹⁴CO₂ (which would then bias the ¹⁴CO₂ and ¹⁴CO fraction). It is theoretically possible that both our measurements and Van der Kemp et al. (2002) are wrong. However, to produce the same ratio of ¹⁴C compounds, it would require all three analytical systems from these studies to be systematically wrong in the same direction and by the same magnitude, which is highly unlikely. The good agreement between the Taylor Glacier and Scharffenbergbotnen data suggests that dry mechanical extraction used by Van der Kemp et al. (2002) is a valid technique for extracting ¹⁴CO₂ and ¹⁴CO from bubbly, non-clathrated ice cores. One possible explanation is that, after production, in situ ¹⁴CO₂ and ¹⁴CO quickly migrates from the ice matrix to the air bubbles. This result is consistent with previous observations that the retention of in situ cosmogenic ¹⁴C in firn grains is very low (Petrenko et al., 2013; Van der Kemp et al., 2000; Wilson and Donahue, 1990).

One way to reconcile our measurements and those of Van der Kemp et al. (2002) with the Heisinger et al. (2002a, b) production rates would be to have much higher (factor of 3 or greater) long-term (hundreds to thousands of years) ablation rates at Taylor Glacier and Scharffenbergbotnen compared to recent ablation stakes measurements. The ablation rate over the last 100 years is especially important with regards to in situ ¹⁴C production rate from negative muon capture (Fig. S9). The ablation rate at blue ice areas is controlled by climate via a combination of temperature, insolation, and wind (mainly katabatic; e.g., Bintanja, 1999). To obtain much higher long-term ablation rates at both Taylor Glacier and Scharffenbergbotnen (which are on opposite sides of Antarctica), we would need either the temperatures to have dropped sharply in the last couple of decades, the winds to have slowed dramatically, or for insolation to have sharply decreased. However, the two deep ice cores nearest to Taylor Glacier, namely Taylor Dome (Steig et al., 2000) and RICE (Roosevelt Island Climate Evaluation; Bertler et al., 2018) have not shown large climate changes in the region over the last 1000 years. The European Project for Ice Coring in Antarctica (EPICA) Dronning Maud Land (EDML) ice core drilled nearby Scharffenbergbotnen blue ice area also has not shown large climate variability in the region over the Holocene period (EPICA Community Members, 2010). The glaciological survey of Taylor Glacier also indicated that the glacier is approximately at steady state, given the stake-measured ablation rates (Kavanaugh et al., 2009a, b). The 30-year record of weather observations over the McMurdo Dry Valleys area show that the recent climate in this region has been stable (Obryk et al., 2020). Finally, a 14-year-long observation study (Sinisalo et al., 2003) over the Scharffenbergbotnen blue ice area has also shown no significant recent change in the ablation rates. We thus argue that a large decrease in ablation rates in recent years, as compared to the long-term average (over the last hundreds to thousands of years), does not seem to be a realistic explanation.

Table 3 ¹⁴CO-specific surface production rates in ice from the two muon mechanisms normalized to the SLHL (sea level high latitude) site using the Lifton et al. (2014) LSDn scaling. All errors shown represent the 95 % confidence interval.

For direct comparison with other studies, we used the scaling factors from the Lifton et al. (2014) LSDn model to calculate the corresponding sea level high latitude (SLHL) total ¹⁴C and ¹⁴CO-specific production rates in ice (Tables 2 and 3). Our estimates of the ¹⁴CO-specific production rates agree with those of Petrenko et al. (2016), within errors (Table 3). Compared with the results from Petrenko et al. (2016), we also calculated a slightly smaller uncertainty in the ¹⁴CO-specific production rate by negative muon capture (Table 3). We converted the Lupker et al. (2015) estimates of f_tot in quartz into f_tot for ice (Table 2), using the chemical compound factors (f_C) for quartz and ice from Heisinger et al. (2002b). With regards to the negative muon capture, the Lupker et al. (2015) estimate of f_tot is in close agreement with Heisinger et al. (2002b; Table 2). This result is supported by Balco (2017), who refitted the Leymon High data and obtained similar f_tot estimates. However, the high f_tot in Lupker et al. (2015), as compared to our result, was offset by their best σ₀ estimate of zero (lower than our result). For a direct comparison with results from Lupker et al. (2015), we fit our data, while forcing σ₀ (and hence ¹⁴C production from fast muons) to be zero (Fig. 6b) and cannot find a scenario with reasonable model–data agreement.

Because of the relatively large uncertainty in the ¹⁴CO₂ measurements, the total ¹⁴C data still allow σ₀ to be close to zero, given a sufficiently large f_tot (Fig. 7a). However, our ¹⁴CO data (which have much lower relative uncertainties and use a more established measurement technique) unambiguously show that σ₀ and ¹⁴C production from fast muons cannot be zero (Figs. 8, 7b). As discussed in Lupker et al. (2015) and Balco (2017), the ¹⁴C data from the 15.5 m Leymon High quartzite core might not cover the depth range in which production from fast muons dominates. In contrast, when integrated over the whole flow history, production from fast muons represents the dominant source of ¹⁴C in our deeper samples. The very high-end estimate of our reference nuclear reaction cross section σ₀ (for ¹⁴C production from fast muons) is still within the large uncertainty of σ₀ from both Heisinger et al. (2002a) and Lupker et al. (2015; Table 2). However, our estimated total probability of negative muon capture (f_tot; and hence the ¹⁴C production rate from negative muon) is well outside the confidence intervals of f_tot reported by both Heisinger et al. (2002b) and Lupker et al. (2015; Table 2).

One caveat to our estimated in situ muogenic ¹⁴C production rates in ice (and that of Van der Kemp et al., 2002) is that the total ¹⁴C from the gas species we measured (¹⁴CO, ¹⁴CO₂, and ¹⁴CH₄) might not account for all the muogenic in situ ¹⁴C. Although ¹⁴CO₂ and ¹⁴CO likely constitute the large majority (Lal et al., 1997, 2000), a small amount of in situ ¹⁴C can also form ¹⁴C-bearing organic materials. Measurements of ¹⁴C in organic carbon from alpine ice for the purpose of radiocarbon dating have shown elevated ¹⁴C values attributed to in situ cosmogenic production (Fang et al., 2021; Hoffmann, 2016). A laboratory irradiation experiment of glacier ice with an artificial neutron flux showed that 11 %–25 % of produced ¹⁴C form organic compounds (Hoffmann, 2016). Earlier work involving the irradiation of ice samples to produce ¹⁴C (e.g., Roessler et al., 1984) also found that organics accounted for a minor fraction of the total ¹⁴C. However, we are not aware of any existing studies that specifically investigated production of ¹⁴C-bearing organic materials in ice from muons. Measuring ¹⁴C in organic compounds is unfortunately beyond the scope of this study, as it requires an entirely different analytical setup.

Another possible explanation for the disagreement is that the 2D ice flow model (and thus the time-dependent exposure history of the ice parcels) might be inaccurate. In the following, we conducted sensitivity analyses to combine both uncertainties by using the +2.4σ maximum ablation rate scenario (which corresponds to the deepest physically possible ice trajectory; Fig. S8) and scaling our total ¹⁴C upward by 25 % to account for the in situ produced ¹⁴C in organics (red dots in Fig. 6d). First, we kept f_tot as a constant (f_tot=0.021, which is the minimum f_tot from the reported uncertainty in Heisinger et al., 2002b) and tuned σ₀ to fit the measurements (dashed red line in Fig. 6d) under the high ablation rate scenario. We find that the best-fit σ₀ (and in situ ¹⁴C production from fast muon) is zero. The modeled total ¹⁴C under this scenario underestimates the total ¹⁴C at lower depths where production from fast muon dominates (>20 m depth) and overestimates the total ¹⁴C at depths where production from negative muon capture dominates (<20 m depth). This means that, even under these extreme scenarios, the ¹⁴C production rate from negative muon capture has to be lower than the lower-bound estimate of Heisinger et al. (2002b), and some production from fast muons is needed to compensate for the lower production rate from negative muon capture to improve the fit. We then repeated our grid search approach (Sect. 3.2) to find the best-fit f_tot and σ₀ that correspond to the 25 % higher total ¹⁴C and maximum ablation rate scenario (solid red line in Fig. 6d). The best-fit f_tot is 0.0055, which is 22 % of the Heisinger et al., 2002b, value) and the best-fit σ₀ is 0.0040 mb (which is 46 % of the Heisinger et al., 2002a, value). We note that these values are within uncertainties in our original best-fit f_tot and σ₀ (Table 2) and still cannot be reconciled with values from Heisinger et al. (2002a, b). We conclude that additional uncertainties from the ice flow history and ¹⁴C contribution from organics likely cannot reconcile the difference between the negative muon capture ¹⁴C production rate inferred by our data and that of Heisinger et al. (2002b).

In their experimental determination of the ¹⁴C production rate by fast muons, Heisinger et al. (2002a) used a single muon energy of 190 GeV (σ(E)). The reference nuclear reaction cross section at 1 GeV (σ₀) was then scaled using the following equation:

$$\begin{array}{}\text{(12)}& \mathit{\sigma }\left(E\right)={\mathit{\sigma }}_{\mathrm{0}}{\stackrel{\mathrm{‾}}{E}}^{\mathit{\alpha }},\end{array}$$

where α is a power factor that describes the energy dependence of the cross section (unitless). However, the mean muon energy ($\stackrel{\mathrm{‾}}{E}$) of 190 GeV used by Heisinger et al. (2002a) and the muon flux intensity were much higher than those expected in the first few hundred meters of ice in natural settings (for the top 200 m of Taylor Glacier, $\stackrel{\mathrm{‾}}{E}=\mathrm{32}$ GeV; Fig. S10). It may be possible that the power factor α of 0.75 is incorrect. Balco (2017) has tried fitting the calibration data with α=1 (which simplifies Eq. 12 into a linear relationship between σ and E). Following Balco (2017), we conducted a sensitivity analysis, fitting our data with similar methods to those described above but with α=1. Increasing α from 0.75 to 1, while keeping σ₀ constant, reduces the overall ¹⁴C production rate from fast muons (Fig. S13). To compensate for the lower production rate from fast muons, the new best-fit σ₀ is now 0.0032 mb, which is 29.8 % higher than the best-fit σ₀ when α=0.75 and 36.6 % that of Heisinger et al. (2002b) value. The best-fit f_tot (overall probability of ¹⁴C production from negative muon capture) becomes 0.0051, which is 17.2 % higher than the best-fit f_tot when α=0.75 and 20.4 % that of Heisinger et al. (2002a) value (Fig. 7a). These values are within the uncertainties in the f_tot and σ₀ derived with α=0.75 that we presented (Table 2; Fig. 7a) and still cannot be reconciled with original Heisinger et al. (2002a, b) values.

5 Conclusions

This study presents ¹⁴CO₂ measurements in ablating ice obtained via a new ice sublimation technique, combined with ¹⁴CO and ¹⁴CH₄ measurements obtained from a well-established large-volume melt–extraction method to estimate the species-specific and total in situ muogenic ¹⁴C production rates in ice. Under the assumption that the majority of in situ ¹⁴C in ice exists as ¹⁴CO, ¹⁴CO₂, and ¹⁴CH₄, we estimated lower muogenic in situ ¹⁴C production rates in ice (by a factor of 5.7 (3.6–13.9) and 3.7 (2.0–11.9) with 95 % CI for negative muon capture and fast muon interactions, respectively) compared to values inferred from laboratory irradiation experiments (Heisinger et al., 2002a, b) and measurements in quartz (Lupker et al., 2015; Balco 2017). Prior ice core measurements from Scharffenbergbotnen (Van der Kemp et al., 2002) obtained with an independent technique also appear to be consistent with these lower muogenic ¹⁴C production rates in ice. This comparison with the Van der Kemp et al. (2002) results also suggests that a dry extraction technique appears to release essentially all in situ ¹⁴C in bubbly (non-clathrated) ice.

At present, there does not appear to be a way to reconcile our Taylor Glacier ice core results and the independent ice core measurements from Scharffenbergbotnen (Van der Kemp et al., 2002) with the muogenic ¹⁴C production rates determined in quartz (Heisinger et al., 2002a, b; Lupker et al., 2015). This is a problem that needs further investigation, and we recommend that future studies address this via laboratory muon irradiation experiments involving both ice and quartz targets and studies that include a quantification of the organic fraction of muogenic ¹⁴C in ice.

Finally, the constraints on muogenic ¹⁴C production rates in ice and the partitioning between the in situ produced ¹⁴C-bearing gas species provided by this study will allow for future measurements of ¹⁴C-containing gases in other ice cores to be used for several applications, including using ¹⁴CO₂ measurements for the absolute dating of the bubbles in ice cores (Andree et al., 1984; Van De Wal et al., 1994) and using ¹⁴CO measurements to either constrain the past oxidative capacity of the atmosphere (Brenninkmeijer et al., 1992; Petrenko et al., 2021) or reconstruct the past cosmic ray flux (BenZvi et al., 2019).