Production of cosmogenic isotopes ⁷Be, ¹⁰Be, ¹⁴C, ²²Na, and ³⁶Cl in the atmosphere: Altitudinal profiles of yield functions

2.1 Formalism

A standard approach to model various effects of cosmic rays in the atmosphere, including production of cosmogenic isotopes, is based on the yield function formalism.

The yield function, Y(E,h), is defined as the production (the number of atoms per gram of air) of the isotope, at given atmospheric depth h, by primary particles of type i with the unit intensity (one primary particle with kinetic energy per nucleon E in the interplanetary space per steradian and cm²). The units of Y are (atoms g⁻¹ cm² sr). The production rate Q of cosmogenic isotope at time t is then defined as an integral of the product of the yield function and the energy spectrum of cosmic rays J_i(E,t) ((sr sec cm²)⁻¹), above the energy E_c corresponding to the local geomagnetic rigidity cutoff P_c:

(1)

where the summation is over different types of primary cosmic ray particles (protons, α particles, etc.). The relation between E_c,i and P_c (defined independently of the yield function computations) is

(2)

where Z_i and A_i are the charge and mass numbers of particles, respectively, E_r=0.938 GeV is the rest mass of a proton. For computations of the yield function we considered, as primary particles, only protons and α particles. Species heavier than helium can be effectively considered as scaled (by the nucleonic number) α particles [see Webber and Higbie, 2003]. An advantage of this approach is that the production rate can be calculated for any type of the energy spectrum beyond the standard modulated spectrum of galactic cosmic rays (GCR), for example, for solar energetic particle events or hypothetical nearby supernova explosions. Sometimes, production of cosmogenic isotopes is calculated directly, without the yield function [O'Brien, 1979; Masarik and Beer, 1999, 2009], but this is related to a prescribed cosmic ray spectrum and cannot be applied for a different or/and revised spectrum.

We note that the computational results are often given not as the strictly defined yield function Y (see above) but the so-called “production function” S, which gives production of the cosmogenic isotope per one primary particle impinging on the top of the atmosphere [e.g., Webber et al., 2007]. Here we show and tabulate the production function S which is a direct result of the simulations and which is related, for the isotropic flux of primary cosmic rays, to the true yield function as

(3)

The factor π appears as conversion between the flux on the top of the atmosphere and the CR intensity in the interplanetary space [cf. e.g., Grieder, 2001, Chapter 1.6.2].

2.2 Numerical Calculation of the Yield Function

Simulations of the nucleonic cascade in the atmosphere were made using direct Monte Carlo simulations by the general-purpose toolkit GEANT4 10.0 developed in CERN for modeling the particle transport and interactions [Agostinelli et al., 2003; Allison et al., 2006]. For our task we applied the embedded physics list QGSP_BIC_HP (Quark-Gluon String model for high-energy interactions + Geant4 Binary Cascade + High-Precision neutron package) [GEANT4 collaboration, 2013]. The Earth's atmosphere was modeled in a realistic manner as a set of spherical layers with homogeneous properties according to the empirical model of the atmosphere NRLMSISE-00 [Picone et al., 2002]. The top of the model was set at the altitude of 100 km, and the layers had the the thickness 1 g/cm² (for the top 20 g/cm²), 5 g/cm² (for the atmospheric depths 20–100 g/cm²), and 10 g/cm² (for the ones below 100 g/cm²) The total depth of the atmosphere was set to 1050 g/cm², and the soil was not included into the model.

The primary cosmic rays were modeled, in each simulation, as the monoenergetic (viz., with a fixed energy) isotropic flux impinging on the top of the atmosphere. We did a series of simulations for two types of cosmic rays, primary protons and α particles with fixed energies in the range from 20 MeV/nucleon to 100 GeV/nucleon with a quasi-logarithmically distributed values. We note that the computations are presented not for energy bins but for the fixed energies as denoted in the supporting information Table S1. From these simulations, the depth-energy distributions of the fluxes of cascade particles (protons, neutrons, and α particles, both primaries and secondaries) were stored as histograms with energy range 1 keV–100 GeV with logarithmic bins in energy (20 bins per decade).

We obtained sums of simulated secondary particles with their energy binned into energy bins of width centered at the energy , which have crossed a given horizontal level (atmospheric depth, h), and applied a weight of (where θ is the zenith angle of the secondaries) to account for the geometrical factor. The minimum value of was limited to 0.001 to avoid too high weights. These sums were divided by the energy bin width to correspond to the quantity which is defined as

(4)

where N_k and v_k are the concentration (in (MeV cm³]⁻¹)) and velocity of secondary particles of type k with energy at the atmospheric depth level h. Then the production function S (in units of atoms/g) at the given atmospheric level h is defined as

(5)

where κ_j is the content of the target nuclei in one gram of air (atoms/g), σ_j,k is the cross-section of the corresponding nuclear reactions, and summation is over the type of secondary particle k and the type j of the target nucleus.

Radiocarbon ¹⁴C is produced mostly via capture of secondary neutrons by atmospheric nitrogen which composes 78% of the atmosphere by volume. The ⁷Be and ¹⁰Be isotopes are produced by spallation of atmospheric nitrogen and oxygen (forming together about 99% of the atmosphere by volume). The ²²Na and ³⁶Cl isotopes are produced by spallation of atmospheric argon (about 1% of the atmosphere by volume) which is much less abundant than nitrogen and oxygen. In computations we adopted the cross-sections from Reyss et al. [1981], Lange et al. [1995], Jull et al. [1998], Webber and Higbie [2003], Tatischeff et al. [2006], Beer et al. [2012], and also from the Experimental Nuclear Reaction Database (EXFOR/CSISRS) http://www.nndc.bnl.gov/exfor/exfor00.htm. We note that cross-sections we used to compute the production of ⁷Be, ¹⁰Be, ¹⁴C, and ²²Na are the same as in our previous works [Usoskin, 2008; Kovaltsov and Usoskin, 2010; Kovaltsov et al., 2012; Leppänen et al., 2012]. Transport of neutrons with energy below 1 keV, for production of ¹⁴C, was calculated in a way similar to the work of Kovaltsov et al. [2012].

The number of simulated cascades was set to assure the statistical accuracy of the computed columnar isotope production to be better than 1%. It varied with the type of primaries and their initial energy, ranging from 1000 simulated cascades (for α particles with the energy of 100 GeV/nucleon) to 2·10⁷ cascades (for 20 MeV protons).

Production functions computed in this way are tabulated in the supporting information, for different isotopes and atmospheric depths. We emphasize that the computational results for α particles are given per nucleon but not per the entire α particle.

An example of the altitude (depth) profile of the production of ³⁶Cl is shown in Figure 1 for primary cosmic ray protons at several selected energies ranging from 0.1 GeV to 3 GeV. The total production is shown as solid lines with dots, while dashed curves depict contribution from secondary neutrons (the difference between the two is due to protons). One can see that for lower energy range Figure 1a, the production of ³⁶Cl is dominated by direct spallation of atmospheric argon nuclei by primary protons in the upper atmospheric layer of a few tens of g/cm², while contribution of secondary particles is much smaller but becomes dominant at greater depths, where the cascade is fully developed. This is because low-energy primary particles have insufficient energy to initiate a developed nucleonic cascade. For energies of the primary particles greater than 1 GeV Figure 1b, the cascade is well developed and the production curve is smooth with nearly exponential attenuation with depth.

3.1 Columnar Production

Columnar (viz., integrated within the entire atmospheric column) production of cosmogenic isotopes is often discussed [Webber et al., 2007; Kovaltsov et al., 2012]:

(6)

where H = 1033 g/cm² is the atmospheric depth at the mean sea level. Accordingly, we present here columnar productions for the purpose of comparison with earlier results (see Figures 1).

One can see that the columnar production curves computed here are in good agreement (within 5–10% for the energy above 100 MeV/nucleon) with earlier results, except for one case. While the present results for ¹⁰Be are fully consistent with those published earlier by Kovaltsov and Usoskin [2010] for energies above 100 MeV/nucleon, the discrepancy with the results by Webber et al. [2007] is more systematic, by a factor of 1.5–1.7 (Figure 3a). We note that we used, for this isotope, the same cross-sections as Webber and Higbie [2003]. We have no clear explanation for this discrepancy, especially taking into account that the result for another beryllium isotope ⁷Be (Figure 5a), which is very similar to ¹⁰Be in the sense of production, is in good agreement between the two models.

We note that the columnar production is shown only for illustration, while advanced studies should include also complex atmospheric transport which can be modeled only using the altitude profiles of the production functions S.

3.2 Global Production Rate

Also, for the purpose of illustration we depict the mean global production rate, due to GCR, of the five isotopes discussed here. The global production rate of an isotope is the averaged over the Globe columnar production rate, defined as

(7)

where Q is given by equation 1 and integration is over the entire atmospheric column h (as in the columnar production) and over the Earth's surface (longitude and latitude) Ω. The time dependence is included into variability of the modulation potential ϕ and in the slow changes of the local geomagnetic cutoff rigidity P_c at each location.

The modulation potential ϕ parameterizes the differential intensity of GCR in the vicinity of Earth, J(E,ϕ(t)) (see equation 1). It varies in accord with solar activity being low (higher GCR flux) and high (low GCR flux) around solar minima and maxima, respectively. The formalism of the modulation potential approach is described in detail elsewhere [e.g., Webber and Higbie, 2003; Vainio et al., 2009]. Here we use the modulation potential defined as in Usoskin et al. [2011] with the local interstellar spectrum, according to Burger et al. [2000]. We note that the exact value of ϕ makes sense only for a fixed reference local interstellar spectrum of GCR [Usoskin et al., 2005; Herbst et al., 2010].

The large-scale geomagnetic field provides additional shielding of the Earth from charged cosmic ray particles. The shielding is often parameterized in terms of the local effective geomagnetic rigidity cutoff P_c [Cooke et al., 1991]. In the first approximation, the value of P_c at each location is determined by the geomagnetic latitude and the geomagnetic dipole moment M, which slowly varies on the centennial-millennial time scale in the range (6–12) ·10²² A m² [e.g., Genevey et al., 2008; Nilsson et al., 2014].

The global production rates are shown in Figure 7 as a function of the modulation potential ϕ for several values of the geomagnetic dipole moment M. It qualitatively resembles other similar plots shown earlier [e.g., Masarik and Beer, 1999; Webber et al., 2007]. One can see that both parameters play an important role. The global production rates are also shown in Table 1 for the modern value of the geomagnetic dipole moment M = 7.8·10²² A m², for the cosmic ray modulation levels corresponding to the mean, minimum, and maximum of solar cycles [Usoskin et al., 2011].

We note different levels of the production rates for different isotopes: from <10⁻⁴ for ²²Na to a few atoms per cm² per second for ¹⁴C. This is defined by cross-sections of the processes (neutron capture for ¹⁴C is much more likely than spallation reaction with a high energy threshold for other isotopes) and by abundance of the target atoms (argon, which is the target for ²²Na, is a factor ≈100 less abundant in the atmosphere than nitrogen and oxygen, which are targets for other elements).

Global production rates calculated here for ⁷Be, ¹⁰Be, ¹⁴C and ²²Na are very close, within 5%, to those published by us earlier [Usoskin, 2008; Kovaltsov and Usoskin, 2010; Kovaltsov et al., 2012; Leppänen et al., 2012] using an older version of the CRAC model. Accordingly, comparisons with other simulation results and direct data that are discussed in great detail in those works, are valid also here. The only new for us isotope is ³⁶Cl, whose global rate agrees within 15% with the value calculated by Masarik and Beer [2009].

We note that the mean global production makes sense only for ¹⁴C, which is globally mixed in the atmosphere. Recent direct in situ measurements of stratospheric radiocarbon during 2002–2005 imply the global production of (2.2 ± 0.6)·10²⁶ atoms of ¹⁴C per year [Kanu et al., 2016]. Our model (Figure 7a) predicts the global ¹⁴C production of 2.36·10²⁶ atoms/yr for the period 2002–2005 (the mean ϕ = 802 MV [Usoskin et al., 2011]), which is in good agreement with the measurements. Since the regional atmospheric transport and deposition prevents global mixing for other isotopes [Beer et al., 2012], it makes little sense to consider the globally averaged production rate for them, especially for short-living ones such as ⁷Be.

For illustration, we show in Table 1 also (geomagnetically) polar and equatorial production rates of the isotopes, for the modern value of the geomagnetic dipole moment. One can see that the geomagnetic field shields cosmic rays effectively (equatorial production rates are about 15% of those in polar regions), while the solar cycle variability is strongest in polar regions (a factor of 2 versus 10–15% at the geomagnetic equator).

4.1 Solar Cycle in ¹⁰Be

As an example of an application of the approach presented here, we have computed the deposition flux of ¹⁰Be in the northern polar region and compared it with the measurements in the North Greenland Ice Core Project (NGRIP) ice core [Berggren et al., 2009]. The deposition flux was calculated in two steps. First, a 3-D time varying pattern of the isotope atmospheric production was calculated using the yield functions presented here and applying the reconstruction of the modulation potential ϕ based on data from the global neutron monitor network since 1951 [Usoskin et al., 2011]. For the atmospheric transport we used a parameterization by Heikkilä et al. [2009, Table 3], applying the mean latitudinal height profile of the tropopause. Finally, we calculated the deposition flux of ¹⁰Be in the Northern polar region. A 1 year delay due to the transport was applied. The calculated ¹⁰Be flux is in very good agreement with the real data, especially for the period 1951–1970 (see Figure 8). A minor (about 5%) discrepancy after 1970 is most likely related to the postdeposition effects in firn and/or to the regional climate variability on annual-decadal time scale [Pedro et al., 2006, 2012]. The good agreement with data validates the yield function for ¹⁰Be.

4.2 The Period Circa 775 A.D.

The strongest known solar energetic particle (SEP) event took place in 775 A.D. [Miyake et al., 2012; Usoskin et al., 2013; Mekhaldi et al., 2015]. Precise measurements of cosmogenic isotopes, in particular ¹⁰Be in ice cores, have been made for that particular period [Miyake et al., 2015; Sigl et al., 2015; Mekhaldi et al., 2015]. Not discussing the event itself, we compared the mean levels of the ¹⁰Be deposition flux for a few decades after the event with the prediction of our model. In Figure 9 we compare the modeled curves for ¹⁰Be deposition flux with that measured in four polar ice cores in Antarctica and Greenland (see Figure 9), averaged over the period 780–800 A.D. (the ice core dating corrections applied according to Sigl et al. [2015]), i.e., after the event and defined by GCR, not SEPs. These measured mean levels of the ¹⁰Be flux are shown by the horizontal hatched strip. The vertical hatched bar corresponds to the (cycle-averaged) value of the modulation potential for that period defined independently using the ¹⁴C record [Usoskin et al., 2016]. The two black curves show the modeled deposition flux (production rate according to the present model and the atmospheric transport and deposition according to the parameterization by Heikkilä et al. [2009], the geomagnetic dipole moment M = 10²³ A m² as reconstructed for that epoch by Licht et al. [2013]), for the northern and southern polar regions. One can see that the measured fluxes of ¹⁰Be are directly reproduced by the model within the possible uncertainties without any ad hoc adjustment or normalization which is typically applied to ¹⁰Be data.

Accordingly, we conclude that with the new yield function, we are able to quantitatively model the production of the cosmogenic radioisotopes in the atmosphere.