A New Full 3-D Model of Cosmogenic Tritium 3 H Production in the Atmosphere (CRAC:3H)
aa r X i v : . [ a s t r o - ph . E P ] S e p manuscript submitted to JGR: Atmospheres
A new full 3D model of cosmogenic tritium Hproduction in the atmosphere (CRAC:3H)
S. V. Poluianov , , G. A. Kovaltsov , I. G. Usoskin , Sodankyl¨a Geophysical Observatory, University of Oulu, Finland Space Physics and Astronomy Research Unit, University of Oulu, Finland Ioffe Physical-Technical Institute, St. Petersburg, Russia
Key Points: • A new CRAC:3H model of cosmogenic tritium ( H) production in the atmosphereis presented. • For the first time, it provides 3D production, also explicitly treating particles heav-ier than protons. • This model provides a useful tool for the use of H as a tracer of atmospheric andhydrological circulation.
Corresponding author: Stepan Poluianov, [email protected] –1–anuscript submitted to
JGR: Atmospheres
Abstract
A new model of cosmogenic tritium ( H) production in the atmosphere is presented. Themodel belongs to the CRAC (Cosmic-Ray Atmospheric Cascade) family and is namedas CRAC:3H. It is based on a full Monte-Carlo simulation of the cosmic-ray induced at-mospheric cascade using the Geant4 toolkit. The CRAC:3H model is able, for the firsttime, to compute tritium production at any location and time, for any given energy spec-trum of the primary incident cosmic ray particles, explicitly treating, also for the firsttime, particles heavier than protons. This model provides a useful tool for the use of Has a tracer of atmospheric and hydrological circulation. A numerical recipe for practi-cal use of the model is appended.
Tritium ( H) is a radioactive isotope of hydrogen with the half-life time of approx-imately 12.3 years. As an isotope of hydrogen, it is involved in the global water cycle andforms a very useful tracer of atmospheric moisture (e.g., Sykora & Froehlich, 2010; Juh-lke et al., 2020) or hydrological cycles (Michel, 2005). In the natural environment, tri-tium is mostly produced by galactic cosmic rays (GCR) in the atmosphere, as a sub-productof the induced nucleonic cascade, and is thus a cosmogenic radionuclide. On the otherhand, tritium is also produced artificially in thermonuclear bomb tests. Before the nuclear-test ban became in force, a huge amount of tritium had been produced artificially andrealised into the atmosphere, leading to an increase of the global reservoir inventory oftritium by two orders of magnitude above the natural level (e.g., Sykora & Froehlich,2010; Cauquoin et al., 2016). Thus, the cosmogenic production of tritium was typicallyneglected as being too small against anthropogenic one. However, as nearly 60 years havepassed since the nuclear tests, its global content has reduced to the natural pre-bomblevel (Palcsu et al., 2018) and presently is mostly defined by the cosmogenic production.Accordingly, natural variability of the isotope production can be again used in atmospherictracing, water vapour transport, dynamics of the stratosphere-troposphere exchanges overAntarctica (Cauquoin et al., 2015; Fourr´e et al., 2018; Palcsu et al., 2018; Juhlke et al.,2020; L´aszl´o et al., 2020). Moreover, a combination of the H data with other tracerslike atmospheric Be, which is also produced by cosmic-ray spallation reactions, but whosetransport is different, can be a very powerful research tool. For this purpose, a reliableproduction model is needed, which is able to provide a full 3D and time variable pro-duction of tritium in the atmosphere.Some models of tritium production by cosmic rays (CR) in the atmosphere weredeveloped earlier. First models (Fireman, 1953; Craig & Lal, 1961; Nir et al., 1966; Lal& Peters, 1967; O’Brien, 1979) were based on simplified numerical or semi-empirical meth-ods of modelling the cosmic-ray induced atmospheric cascade. Later, a full Monte-Carlosimulation of the cosmogenic isotope production in the atmospheric cascade had beendeveloped (Masarik & Beer, 1999) leading to higher accuracy of the results. However,that model had some significant limitations: (1) were considered only GCR protons (heav-ier GCR species were treated as scaled protons); (2) the energy spectrum of GCR wasprescribed; (3) only global and latitudinal zonal mean productions were presented, im-plying no spatial resolution. That model was slightly revisited by Masarik & Beer (2009),but the methodological approach remained the same. A more recent tritium productionmodel developed by Webber et al. (2007) is also based on a full Monte-Carlo simulationof the atmospheric cascade and was built upon the yield-function approach which allowsdealing with any kind of the cosmic-ray spectrum. However, only columnar (for the en-tire atmospheric column) production was provided by those authors, making it impos-sible to model the height distribution of isotope production. Moreover, that model wasdealing with CR protons only, while the contribution of heavier species to cosmogenicisotope production can be as large as 40% (see section 3). –2–anuscript submitted to
JGR: Atmospheres
Here we present a new model of cosmogenic tritium production in the atmosphere,that is based on a full simulation of the cosmic-ray induced atmospheric cascade. Thismodel belongs to the CRAC (Cosmic-Ray Atmospheric Cascade) family and is namedas
CRAC:3H . The CRAC:3H model is able, for the first time, to compute tritium pro-duction at any location and time, for any given energy spectrum of the primary incidentCR particles, explicitly treating, also for the first time, particles heavier than protons.This model provides a useful tool for the use of H as a tracer of atmospheric and hy-drological circulation.
The local production rate q of a cosmogenic isotope, in atoms per second per gramof air, at a given location with the geomagnetic rigidity cutoff P c and the atmosphericdepth h can be written as q ( h, P c ) = X i Z ∞ E c ,i J i ( E ) · Y i ( E, h ) · dE, (1)where J i ( E ) is the intensity of incident cosmic-ray particles of the i -th type (character-ized by the charge Z i and atomic mass A i numbers) in units of particles per (s sr cm GeV), Y i ( E, h ) is the isotope yield function in units of (atoms sr cm g − , — see sec-tion 2.1 for details), E is the kinetic energy of the incident particle in GeV, h is the at-mospheric depth in units of (g/cm ), E c,i = q ( Z i · P c /A i ) + E − E is the energycorresponding to the local geomagnetic cutoff rigidity for a particle of type i , and thesummation is over the particle types. E = 0 .
938 GeV is the proton’s rest mass. Thegeomagnetic rigidity cutoff P c quantifies the shielding effect of the geomagnetic field andcan be roughly interpreted as a rigidity/energy threshold of primary incident chargedparticles required to imping on the atmosphere (see formalism in Elsasser, 1956; Smartet al., 2000). Here we computed the tritium production function in a way similar to our previ-ous works in the framework of the CRAC-family models (e.g., Usoskin & Kovaltsov, 2008;Kovaltsov & Usoskin, 2010; Kovaltsov et al., 2012; Poluianov et al., 2016), viz. by ap-plying a full Monte-Carlo simulation of the cosmic-ray induced atmospheric cascade, asbriefly described below. Full description of the nomenclature and numerical approachis available in Poluianov et al. (2016).The yield function Y i ( E, h ) (see equation 1) of a nuclide of interest provides thenumber of atoms produced in the unit (1 g/cm ) atmospheric layer by incident parti-cles of type i (e.g., cosmic ray protons, α -particles, heavier species) with the fixed en-ergy E and the unit intensity (1 particle per cm per steradian). The yield function shouldnot be confused with the so-called production function S i ( E, h ), which is defined as thenumber of nuclide atoms produced in the unit atmospheric layer per one incident par-ticle with the energy E . In a case of the isotropic particle distribution, these quantitiesare related as Y = πS, (2)where π is the conversion factor between the particle intensity in space and the parti-cle flux at the top of the atmosphere (see, e.g., chapter 1.6.2 in Grieder, 2001).The production function in units (atoms cm /g) can be calculated, for the isotropicflux of primary CR particles of type i , as S i ( E, h ) = X l Z E η l ( E ′ ) · N i,l ( E, E ′ , h ) · v l ( E ′ ) · dE ′ , (3) –3–anuscript submitted to JGR: Atmospheres S pe c i f i c σ ( m b ) Energy (MeV) a) N(n,x) H N(p,x) HO(n,x) H O(p,x) H 10 -4 -3 A gg r ega t e η ( c m / g ) Energy (MeV) b) protonsneutrons
Figure 1.
Specific σ (panel a) adopted from Nir et al. (1966); Coste et al. (2012) and aggre-gate η ( E ) (panel b) cross-sections for production of tritium as a function of the particle’s energy. where summation is over types l of secondary particles of the cascade (can be protons,neutrons, α -particles), η l is the ‘ aggregate ’ cross-section (see below) in units (cm /g), N i,l ( E, E ′ , h )and v l ( E ′ ) are concentration and velocity of the secondary particles of type l with en-ergy E ′ at depth h . The aggregate cross-section η l ( E ′ ) is defined as η l ( E ′ ) = X j κ j · σ j,l ( E ′ ) , (4)where j indicates the type of a target nucleus in the air (nitrogen and oxygen for tritium), κ j is the number of the target nuclei of type j in one gram of air, σ j,l ( E ′ ) is the totalcross-section of nuclear reactions between the l -th atmospheric cascade particle and the j -th target nucleus yielding the nuclide of interest. Atmospheric tritium is produced byspallation of target nuclei of nitrogen and oxygen, which have the values of κ N = 3 . · g − and κ O = 8 . · g − , respectively. The reactions yielding tritium are causedmainly by the cascade neutrons and protons and include: N(n,x) H; N(p,x) H; O(n,x) H;O(p,x) H. The cross-sections used here were adopted from Nir et al. (1966) and Costeet al. (2012), as shown in Figure 1a. We assumed that cross-sections of the neutron-inducedreactions are similar to those for protons above the energy of 2 GeV. For reactions causedby α -particles, N( α ,x) H and O( α ,x) H, the cross-sections were assessed from proton onesaccording to Tatischeff et al. (2006). These reactions are induced mostly by α -particlesfrom the primary CRs and are, hence, important only in the upper atmospheric layers.The tritium aggregate cross-sections η (equation 4) are shown in Figure 1b. Althoughproduction efficiencies of protons and neutrons are similar at high energies, they differsignificantly in the <
500 MeV range. Because of the lower energy threshold and highercross-sections for neutrons in this energy range, comparing to protons, tritium produc-tion is dominated by neutrons in a region where the cascade is fully-developed, viz., inthe lower part of the atmosphere.The quantity N i,l ( E, E ′ , h ) · v l ( E ′ ) describing the cascade particles (equation 3)was computed using a full Monte Carlo simulation of the cascade induced in the atmo-sphere by energetic cosmic-ray particles. The general computation scheme was similarto that applied by Poluianov et al. (2016). The simulation code was based on the Geant4toolkit v.10.0 (Agostinelli et al., 2003; Allison et al., 2006). In particular, we used thephysics list QGSP BIC HP (Quark-Gluon String model for high-energy interactions; Geant4Binary Cascade model; High-Precision neutron package) (Geant4 collaboration, 2013),which was shown to describe the cosmic ray cascade with sufficient accuracy (e.g., Mesick –4–anuscript submitted to JGR: Atmospheres et al., 2018). We simulated a real-scale spherical atmosphere with the inner radius of 6371km, height of 100 km and thickness of 1050 g/cm . The atmosphere was divided intohomogeneous spherical layers with the thickness ranging from 1 g/cm (at the top) to10 g/cm near the ground. The atmospheric composition and density profiles were takenaccording to the atmospheric model NRLMSISE-00 (Picone et al., 2002). Cosmic rayswere simulated as isotropic fluxes of mono-energetic protons and α -particles, while heav-ier species were considered as scaled α − particles (see section 2.2). The simulations wereperformed with a logarithmic grid of energies between 20 MeV/nuc and 100 GeV/nuc.The number of simulated incident particles was set so that the statistical accuracy of theisotope production should be better than 1% in any location. This number varied from1000 incident particles for α − particles with the energy of 100 GeV/nucleon to 2 · simulations for 20-MeV protons. The results were extrapolated to higher energies, upto 1000 GeV/nuc, by applying a power law. The yield of the secondary particles (pro-tons, neutrons and α -particles) at the top of each atmospheric layer was recorded as his-tograms with the spectral (logarithmic) resolution of 20 bins per one order of magnitudein the range of the secondary particle’s energy between 1 keV and 100 GeV. The primaryCR particles were also recorded in the same way.The production functions S i ( E, h ) were subsequently calculated from the simula-tion results, using equation (3), for a prescribed grid of energies and atmospheric depthsand are tabulated in the Supporting Information. Some examples of the tritium produc-tion function are shown in Figure 2 for primary CR protons. One can see in Figure 2athat the efficiency of tritium atom production grows with the energy of the incident par-ticles because of larger atmospheric cascades induced. Contributions of different com-ponents to the total production are shown in Figure 2b for low (0.1 GeV) and medium(1 GeV) energies of the primary proton. The red curve for the 0.1 GeV incident protonsdepicts a double-bump structure: the bump in the upper atmospheric layers ( h <
10 g/cm )is caused by spallation reactions caused mostly by the primary protons (as indicated bythe red dotted curve), while the smooth curve at higher depths is due to secondary neu-trons (red dashed curve). Overall, production of tritium at depths greater than 10 g/cm is very small for the low-energy primary protons. On the other hand, higher-energy (1GeV, blue curves in Figure 2b) protons effectively form a cascade reaching the ground,where the contribution of secondary neutrons dominates below ≈
50 g/cm depths.This type of the depth/altitude profiles or the tritium production function was notstudied in earlier works, where only columnar functions, viz. integrated over the full at-mospheric column, were presented (Webber et al., 2007). Therefore, in order to compareour results with the earlier published ones, we also calculated the columnar productionfunction S C ( E ) = Z h sl S ( E, h ) · dh, (5)where h sl = 1033 g/cm is the atmospheric depth at the mean sea level or the thick-ness of the entire atmospheric column. The columnar production function is tabulatedin the Supporting Information and depicted in Figure 3 along with the earlier results pub-lished by Webber et al. (2007) for incident protons. No results for incident α -particleshave been published earlier, and the production function of cosmogenic tritium by cosmic-ray α -particles is presented here for the first time. One can see that, while the produc-tion functions for incident protons generally agree between our work and the results byWebber et al. (2007), there are some small but systematic differences. In particular, ourresult is lower than that of Webber et al. (2007) in the low-energy range below 100 MeV.It should be noted that the contribution of this energy region to the total productionof tritium is negligible because of the geomagnetic shielding in such a way that low-energyincident particles can impinge on the atmosphere only in spatially small polar regions.In the energy range above 200 MeV, the tritium production function computed here ishigher than that from Webber et al. (2007). The difference is not large, ≈ –5–anuscript submitted to JGR: Atmospheres -8 -7 -6 -5 -4 -3 -2 -1
10 100 100010 -8 -7 -6 -5 -4 -3 -2 -1 S ( c m / g ) Depth h (g/cm ) a) 0.1 GeV 0.2 GeV0.5 GeV 1 GeV 2 GeV 5 GeV 10 GeV b) 0.1 GeV: sum n p 1 GeV: sum n p Depth h (g/cm ) Figure 2.
Production function S = Y /π of tritium by primary protons. Panel a: productionfunction S by primary protons with energies between 0.1–10 GeV, as denoted in the legend.Panel b: contribution of protons (p) and secondary neutrons (n) to the production function(sum) for 0.1 GeV (red) and 1 GeV (blue) primary protons. cade simulation (FLUKA vs. Geant4). Overall, our model predicts slightly higher pro-duction of tritium than the one by Webber et al. (2007), for the same cosmic-ray flux. The first term J i ( E ) in equation (1) refers to the spectrum of differential intensityof the incident cosmic-ray particles. A standard way to model the GCR spectrum forpractical applications is based on the so-called force-field approximation (Gleeson & Ax-ford, 1967; Caballero-Lopez & Moraal, 2004; Usoskin et al., 2005), which parameterizesthe spectrum with reasonable accuracy even during disturbed periods, as validated bydirect in-space measurements (Usoskin et al., 2015). In this approximation, the differ-ential energy spectrum of the i -th component of GCR near Earth (outside of the Earth’smagnetosphere and atmosphere) is parameterized in the following form: J i ( E, t ) = J LIS ,i ( E + Φ i ( t )) E ( E + 2 E )( E + Φ i ( t ))( E + Φ i ( t ) + 2 E ) , (6)where J LIS ,i is the differential intensity of GCR particles in the local interstellar medium,often called the local interstellar spectrum (LIS), E is the particle’s kinetic energy pernucleon, E is the rest energy of a proton (0.938 GeV), and Φ i ( t ) ≡ φ ( t ) · Z i /A i is themodulation parameter defined by the modulation potential φ ( t ) as well as the charge ( Z i )and atomic ( A i ) numbers of the particle of type i , respectively. The spectrum at any mo-ment of time t is fully determined by a single time-variable parameter φ ( t ), which hasthe dimension of potential (typically given in MV or GV) and is called the modulationpotential. The absolute value of φ makes no physical sense and depends on the exact shapeof LIS (see discussion in Usoskin et al., 2005; Herbst et al., 2010, 2017; Asvestari et al.,2017).In this work, we made use of a recent parameterization of the proton LIS (Vos &Potgieter, 2015), which is partly based on direct in situ measurements of GCR: J LIS ( E ) = 0 . E . β (cid:18) E + 0 . . (cid:19) − . , (7) –6–anuscript submitted to JGR: Atmospheres -6 -4 -2 S C Energy (GeV/nuc) protons, this work α -particles, this work protons, Webber et al., 2007 Figure 3.
Columnar production function S C = Y C /π (number of atoms per primary incidentnucleon) of tritium by incident protons (blue line) and α -particles (red line). Tabulated valuesare available in the Supporting Information. Circles indicate the production function for protonsfrom Webber et al. (2007). where J LIS ( E ) is the differential intensity of GCR protons in the local interstellar mediumin units of particles per (s sr cm GeV), E and β = v/c are the particle’s kinetic en-ergy (in GeV) and the velocity-to-speed-of-light ratio, respectively. Following a recentwork (Koldobskiy et al., 2019) based on a joint analysis of data from the space-borneexperiment AMS-02 (Alpha Magnetic Spectrometer) and from the ground-based neutron-monitor network, we assumed that LIS (in the number of nucleons) of all heavier ( Z ≥ P c , and at the time moment t it is defined as Q C ( P c , t ) = Z h sl q ( h, P c , t ) · dh. (8)The global production rate Q global is the spatial average of Q C ( P c ) over the globe, whilethe integral of Q over the globe yields the total production of tritium.Production of tritium by GCR, which always bombard the Earth’s atmosphere, isdescribed above. Production by solar energetic particles (SEP) can be computed in asimilar way, with the SEP energy spectrum entering directly in equation (1). Using the production function computed here (section 2.1) and applying equations (1)and (8), we calculated the mean production rate Q of tritium in the atmosphere for dif-ferent levels of solar modulation (low, moderate and high), for the entire atmosphere andonly for the troposphere. The results are shown in Table 1. The modeled local produc-tion rates q ( h, P c ) (equation 1) used for the computation can be found in a tabular formin the supporting information. –7–anuscript submitted to JGR: Atmospheres
Table 1.
Tritium production rates (in atoms/(s cm )) averaged globally (see also Figure 5)and over the polar regions (geographical latitude 60 ◦ –90 ◦ ), separately in the entire atmosphereand only the troposphere for different levels of solar activity: low, medium and high ( φ =400, 650and 1100 MV, respectively). The values of the modulation potential correspond to the formal-ism described in section 2.2. The geomagnetic field is taken according to IGRF (InternationalGeomagnetic Reference Field, Th´ebault et al., 2015) for the epoch 2015. The tropopause heightprofile is adopted from Wilcox et al. (2012). Solar activity Entire atm. TroposphereGlobal Polar Global PolarLow 0.41 0.92 0.12 0.16Moderate 0.345 0.72 0.11 0.14High 0.27 0.51 0.09 0.10The global production rate of tritium for a moderate solar activity ( φ = 650 MV),which is the mean level for the modern epoch (Usoskin et al., 2017), is 0.345 atoms/(scm ). This value can be compared with earlier estimates of the global production rateof tritium. We performed a literature survey and found that the estimates performed be-fore 1999 were based on different approximated approaches and vary by a factor of 2.5,between 0.14–0.36 atoms/(s cm ) (Craig & Lal, 1961; Nir et al., 1966; O’Brien, 1979; Masarik& Reedy, 1995). Modern estimates, based on full Monte-Carlo simulations, are more con-strained. The early value of the global production rate of 0.28 atoms/(s cm ) publishedby Masarik & Beer (1999) was revised by the authors to 0.32 atoms/(s cm ) in Masarik& Beer (2009). Our value is very close to that, despite the different computational schemesand assumptions made. The computed global production rate also agrees with the es-timates obtained from reservoir inventories (e.g. Craig & Lal, 1961), that are, however,loosely constrained within a factor of about four, between 0.2–0.8 atoms/(s cm ). Wenote that heavier-than-proton primary incident particles contribute about 40% to theglobal production of tritium, in the case of GCR, and thus, it is very important to con-sider these particles explicitly.Geographical distribution of the columnar production rate Q C ( P c ) of tritium is shownin Figure 4. It is defined primarily by the geomagnetic cutoff rigidity (e.g., Smart & Shea,2009; Nevalainen et al., 2013) and varies by an order of magnitude between the high-cutoffspot in the equatorial west-Pacific region and polar regions.Dependence of the global production rate of tritium on solar activity quantified viathe modulation potential φ is shown in Figure 5, both for the entire atmosphere and forthe troposphere. The tropospheric contribution to the global production is about 31%on average, ranging from 30% (solar minimum) to 34% (solar maximum).Even though the production rate is significantly higher in the polar region, its con-tribution to the global production is not dominant, because of the small area of the po-lar regions. Figure 6 (upper panel) presents the production rate of tritium in latitudi-nal zones (integrated over longitude in one degree of geographical latitude) as a func-tion of geographical latitude and atmospheric depth. It has a broad maximum at mid-latitudes (40–70 ◦ ) in the stratosphere (10–100 g/cm of depth) and ceases both towardsthe poles and ground. The bottom panel of the Figure depicts the zonal mean contri-bution (red curve) of the entire atmospheric column into the total global production. Itillustrates that the distribution with a maximum at mid-latitudes shape is defined bytwo concurrent processes: the enhanced production (green curve) and reduced zonal area –8–anuscript submitted to JGR: Atmospheres La t i t ude ( deg ) Longitude (deg) -90-60-300306090 -135 -90 -45 0 45 90 135 1800.1 . . . . . . . . . . . . . . . . . . . . . . . . . Figure 4.
Geographical distribution of the columnar production rate Q C (atoms/(s cm )) oftritium by GCR corresponding to a moderate level of solar activity ( φ =650 MV). The geomag-netic cutoff rigidities were calculated using the eccentric tilted dipole approximation (Nevalainenet al., 2013) for the IGRF model (epoch 2015). Other model parameters are as described above.The background map is from Gringer/Wikimedia Commons/public domain. (blue curve) from the equator to the pole. The zonal contribution is proportional to theproduct of these two processes.The altitude profile of the tritium production rate by GCR for the moderate levelof solar activity is shown in Figure 7. The maximum of the globally averaged produc-tion is located at about 40 g/cm or 20 km of altitude in the stratosphere, correspond-ing to the Regener-Pfotzer maximum where the atmospheric cascade is most developed.The maximum of production is somewhat higher in the polar region because of the re-duced geomagnetic shielding there, so that lower-energy CR particles can reach the lo-cation.Figure 8 depicts temporal variability of the global tritium production for the pe-riod 1951–2018, computed using the model presented here. To indicate the solar cycleshape, the sunspot numbers are also shown in the bottom. The contribution from GCRis shown by the blue curve and computed using the modulation potential reconstructedfrom the neutron-monitor network (Usoskin et al., 2017). Red dots consider also addi-tional production of tritium by strong SEP events, identified as ground-level enhance-ment (GLE) events ( http://gle.oulu.fi ). This is negligible on the long run but maycontribute essentially on the short-time scale. Overall, the production of tritium is mostlydriven by the heliospheric modulation of GCR as implied by obvious anti-correlation withthe sunspot numbers. A new full model CRAC:3H of tritium cosmogenic production in the atmosphereis presented. It is able to compute the tritium production rate at any location in 3D andfor any type of the incident particle energy spectrum/intensity — slowly variable galac- –9–anuscript submitted to
JGR: Atmospheres Q g l oba l ( a t o m s / ( s c m )) Modulation potential φ (MV)total atm.troposphere Figure 5.
Global columnar production Q global of tritium, in the entire atmosphere and only inthe troposphere, as a function of solar activity quantified via the heliospheric modulation poten-tial. The shaded area denotes the range of a solar cycle modulation for the modern epoch. Thegeomagnetic field corresponds to the IGRF for the epoch 2015. The tropopause height profileis adopted from Wilcox et al. (2012). The values of the modulation potential correspond to theformalism described in section 2.2. tic cosmic rays or intense sporadic events of solar energetic particles. The core of the modelis the yield/production function, rigorously computed by applying a full Monte-Carlosimulation of the cosmic-ray induced atmospheric cascade with high statistics and is tab-ulated in the Supporting Information. Using this tabulated function, one can straight-forwardly and easily calculate the production of tritium for any conditions in the Earth’satmosphere (see Appendix A), including solar modulation of GCR, sporadic SEP events,changes of the geomagnetic field, etc. The columnar and global production of tritium,computed by the new model, is comparable with most recent estimates by other groups,but is significantly higher than the results of earlier models, published before 2000. Italso agrees well with empirical estimates of the tritium reservoir inventories, consider-ing large uncertainties of the latter. Thus, for the first time, a reliable model is devel-oped that provides a full 3D production of tritium in the atmosphere. These results canbe used as an input for atmospheric transport models or for direct comparison with tri-tium observations that are important for the study of solar activity in link with the hy-drological cycle or for evaluation of the atmospheric dynamics in models. Appendix A Calculation of tritium production: Numerical algorithm
Using the production function S ( E, h ) presented here in the Supporting Informa-tion, one can easily compute tritium production at any given location (quantified by thelocal geomagnetic rigidity cutoff P c and atmospheric depth h ), and time t , following thenumerical algorithm below.1. For a given moment of time t , the intensity of incident primary particles can beevaluated, in case of GCR, using equations (6) and (7) for the independently knownmodulation potential φ (e.g., as provided at http://cosmicrays.oulu.fi/phi/phi.html ). These formulas can be directly applied for protons, while the contri- –10–anuscript submitted to JGR: Atmospheres D ep t h h ( g / c m ) Zonal mean contribution C zonal Columnar production rate Q C Cosine of latitude C z ona l ( deg - ) Latitude (deg)
Figure 6.
Upper panel: Tritium zonal production rate by GCR ( φ =650 MV, geomagneticfield IGRF epoch 2015) as a function of the atmospheric depth and northern geographical lati-tude. The color scale (on the right) is given in units of atoms per second per degree of latitudeper gram/cm . Bottom panel: zonal mean contribution C zonal (red curve, per degree of latitude)to the tritium global production rate (a columnar integral of the distribution shown in the upperpanel), normalized so that its total integral over all latitudes is equal to unity. Green dot-dashedand blue dashed lines represent the columnar production rate and cosine of latitude, respectively(both in arbitrary units), and C zonal is directly proportional to their product.–11–anuscript submitted to JGR: Atmospheres D ep t h h ( g / c m ) Production rate q (10 -3 atoms/(g s)) global polar 20 km Figure 7.
Altitude profile of the tritium differential production q (equation 1) by GCR for themoderate solar activity level ( φ =650 MV). The red solid and blue dash lines represent the globaland polar (60 ◦ –90 ◦ ) production rates, respectively. The horizontal marks on the right indicatethe approximate altitude, which depends on the exact atmospheric conditions. Q g l oba l ( a t o m s / ( s c m )) Years GCRGCR+SEP SSN SS N Figure 8.
Monthly means of the global production rates Q global of tritium computed herefor the period 1951–2018. The blue curve is for the GCR production (modulation potential andgeomagnetic field were adopted from Usoskin et al. (2017) and IGRF, respectively). The reddots indicate periods of GLE events ( http://gle.oulu.fi ) with additional production of tritiumby SEPs as computed using the spectral parameters adopted from Raukunen et al. (2018). Thegrey-shaded curve in the bottom represents the sunspot number (right-hand side axis) adoptedfrom SILSO ( , Clette & Lef`evre, 2016).–12–anuscript submitted to JGR: Atmospheres bution of heavier species ( Z ≥
2) can be considered, using the same formulas, butapplying the scaling factor of 0.353 for LIS, which is given in number of nucleons,and considering kinetic energy per nucleon. Thus, the input intensities of the in-cident protons J p ( E, t ) and heavier species J α ( E, t ), the latter effectively includ-ing all heavier species, can be obtained. Energy should be in units of GeV, and J ( E ) in units of nucleons per (sr cm s GeV). The energy grid is recommendedto be logarithmic (at least 10 points per order of magnitude).2. The production function S i ( E, h ) for the given atmospheric depth h can be ob-tained, for both protons S p and heavier species S α , from the Supporting Informa-tion in units of (cm /g). The yield function is defined as Y = π · S , in units of(sr cm /g), also separately for protons and heavier species. The product of theyield function and the intensity of incident particles is called the response func-tion F i ( E, h ) = Y i ( E, h ) · J i ( E ), separately for protons and heavier species.3. As the next step, the local geomagnetic rigidity cutoff P c , which is related to thelower integration bound in equation (1), needs to be calculated for a given loca-tion and time. A good balance between simplicity and realism is provided by theeccentric tilted dipole approximation of the geomagnetic field (Nevalainen et al.,2013). The value of P c in this approximation can be computed using a detailednumerical recipe (Appendix A in Usoskin et al., 2010). This approach works wellwith GCR, but is too rough for an analysis of SEP events, where a detailed com-putations of the geomagnetic shielding is needed (e.g., Mishev et al., 2014).4. Next, the response function F i should be integrated above the energy bound de-fined by the geomagnetic rigidity cutoff P c , as specified in equation (1) separatelyfor the protons and α − particles (the latter effectively includes also heavier Z > F ( E ) whose val-ues are defined at grid points E and E as F and F , respectively, be approx-imated by a power law between these grid points. Then its integral on the inter-val between these grid points is Z E E F ( E ) · dE = ( F · E − F · E ) · ln ( E /E )ln ( F /F ) + ln ( E /E ) . (A1)The final production rate at the given location, atmospheric depth and time is thesum of the two components (protons and α − particles).5. In a case when not only the very local production rate of tritium is required, butspatially integrated or averaged, the columnar production function (equation 8)can be used. The spatially averaged/integrated production can be then obtainedby averaging/integration over the appropriate area considering the changes in thegeomagnetic cutoff rigidity P c . Acknowledgments
The yield/production functions of tritium, obtained in this work, are available in the Sup-porting Information to this paper. The used cross-section data can be found in Nir etal. (1966) and Coste et al. (2012). The toolkit Geant4 (Agostinelli et al., 2003; Allisonet al., 2006) is freely distributed under Geant4 Software License at . This work used publicly available data for SEP events from the GLE database ( http://gle.oulu.fi ), sunspot number series from SILSO ( ,Clette & Lef`evre, 2016), heliospheric modulation potential series provided by the Oulucosmic ray station ( http://cosmicrays.oulu.fi/phi/phi.html ). S.P. acknowledgesthe International Joint Research Program of ISEE, Nagoya University, and thanks NaoyukiKurita from Nagoya University for valuable discussion. This work was partly supported –13–anuscript submitted to
JGR: Atmospheres by the Academy of Finland (Projects ESPERA no. 321882 and ReSoLVE Centre of Ex-cellence, no. 307411).
References
Agostinelli, S., Allison, J., Amako, K., Apostolakis, J., Araujo, H., Arce, P., . . .Zschiesche, D. (2003). Geant4 - a simulation toolkit.
Nucl. Instr. Meth.Phys. A , (3), 250–303. Retrieved from doi: http://dx.doi.org/10.1016/S0168-9002(03)01368-8Allison, J., Amako, K., Apostolakis, J., Araujo, H., Dubois, P., Asai, M., . . .Yoshida, H. (2006). Geant4 developments and applications. Nuclear Science,IEEE Transactions on , (1), 270-278. doi: 10.1109/TNS.2006.869826Asvestari, E., Gil, A., Kovaltsov, G., & Usoskin, I. (2017). Neutron Monitorsand Cosmogenic Isotopes as Cosmic Ray Energy-Integration Detectors: Effec-tive Yield Functions, Effective Energy, and Its Dependence on the Local In-terstellar Spectrum. J. Geophys. Res. Space Phys. , (10), 9790-9802. doi:10.1002/2017JA024469Caballero-Lopez, R. A., & Moraal, H. (2004). Limitations of the force field equa-tion to describe cosmic ray modulation. J. Geophys. Res.: Space Phys. , (A1),A01101. Retrieved from http://dx.doi.org/10.1029/2003JA010098 doi: 10.1029/2003JA010098Cauquoin, A., Jean-Baptiste, P., Risi, C., Fourr´e, ´E., & Landais, A. (2016). Mod-eling the global bomb tritium transient signal with the AGCM LMDZ-iso: Amethod to evaluate aspects of the hydrological cycle. J. Geophys. Res. (Atmos.) , (21), 12,612-12,629. doi: 10.1002/2016JD025484Cauquoin, A., Jean-Baptiste, P., Risi, C., Fourr´e, ´E., Stenni, B., & Landais, A.(2015). The global distribution of natural tritium in precipitation simulated withan Atmospheric General Circulation Model and comparison with observations. Earth Planet. Sci. Lett. , , 160-170. doi: 10.1016/j.epsl.2015.06.043Clette, F., & Lef`evre, L. (2016). The New Sunspot Number: Assembling All Correc-tions. Solar Phys. , , 2629-2651. doi: 10.1007/s11207-016-1014-yCoste, B., Derome, L., Maurin, D., & Putze, A. (2012). Constraining galac-tic cosmic-ray parameters with z ≤ Astron. Astrophys. , , A88.Retrieved from https://doi.org/10.1051/0004-6361/201117927 doi:10.1051/0004-6361/201117927Craig, H., & Lal, D. (1961). The Production Rate of Natural Tritium. Tellus Ser. A , (1), 85-105. doi: 10.1111/j.2153-3490.1961.tb00068.xElsasser, W. (1956). Cosmic-Ray Intensity and Geomagnetism. Nature , , 1226-1227. doi: 10.1038/1781226a0Fireman, E. L. (1953). Measurement of the (n, H ) Cross Section in Nitrogen andIts Relationship to the Tritium Production in the Atmosphere. Phys. Rev. , (4),922-926. doi: 10.1103/PhysRev.91.922Fourr´e, E., Landais, A., Cauquoin, A., Jean-Baptiste, P., Lipenkov, V., & Petit,J. R. (2018). Tritium Records to Trace Stratospheric Moisture Inputs in Antarc-tica. J. Geophys. Res. (Atmos.) , (6), 3009-3018. doi: 10.1002/2018JD028304Geant4 collaboration. (2013). Physics reference manual (ver-sion geant4 9.10.0) [Computer software manual]. (available from http://geant4.cern.ch/support/index.shtml )Gleeson, J. J., & Axford, W. I. (1967). Cosmic rays in the interplanetary medium. Astrophys. J. , , L115. Retrieved from http://adsabs.harvard.edu/abs/1967ApJ...149L.115G doi: 10.1086/180070Grieder, P. K. F. (2001). Cosmic Rays at Earth . Amsterdam: Elsevier Science.Herbst, K., Kopp, A., Heber, B., Steinhilber, F., Fichtner, H., Scherer, K., & –14–anuscript submitted to
JGR: Atmospheres
Matthi, D. (2010). On the importance of the local interstellar spectrum for the so-lar modulation parameter.
J. Geophys. Res.: Atmos. , (D1), D00I20. Retrievedfrom http://dx.doi.org/10.1029/2009JD012557 doi: 10.1029/2009JD012557Herbst, K., Muscheler, R., & Heber, B. (2017). The new local interstellar spectraand their influence on the production rates of the cosmogenic radionuclides Beand C. Journal of Geophysical Research: Space Physics , (1), 23-34. doi:10.1002/2016JA023207Juhlke, T., S¨ultenfuß, J., Trachte, K., Huneau, F., Garel, E., Santoni, S., . . . vanGeldern, R. (2020). Tritium as a hydrological tracer in Mediterranean pre-cipitation events. Atmos. Chem. Phys. , (6), 3555-3568. doi: 10.5194/acp-20-3555-2020Koldobskiy, S. A., Bindi, V., Corti, C., Kovaltsov, G. A., & Usoskin, I. G. (2019).Validation of the Neutron Monitor Yield Function Using Data From AMS-02 Ex-periment, 2011-2017. J. Geophys. Res.: Space Phys. , (4), 2367-2379. doi:10.1029/2018JA026340Kovaltsov, G. A., Mishev, A., & Usoskin, I. G. (2012). A new model of cosmogenicproduction of radiocarbon C in the atmosphere.
Earth Planet. Sci. Lett. , ,114-120. doi: 10.1016/j.epsl.2012.05.036Kovaltsov, G. A., & Usoskin, I. G. (2010). A new 3D numerical model of cosmogenicnuclide Be production in the atmosphere.
Earth Planet. Sci. Lett. , , 182-188.doi: 10.1016/j.epsl.2010.01.011Lal, D., & Peters, B. (1967). Cosmic ray produced radioactivity on the earth. InK. Sittle (Ed.), Handbuch der physik (Vol. 46, pp. 551–612). Berlin: Springer.L´aszl´o, E., Palcsu, M., & Leel˝ossy, ´A. (2020). Estimation of the solar-induced nat-ural variability of the tritium concentration of precipitation in the Northern andSouthern Hemisphere.
Atmos. Envir. . doi: 10.1016/j.atmosenv.2020.117605Masarik, J., & Beer, J. (1999). Simulation of particle fluxes and cosmogenic nuclideproduction in the Earth’s atmosphere.
J. Geophys. Res. , , 12099-12112. doi:10.1029/1998JD200091Masarik, J., & Beer, J. (2009). An updated simulation of particle fluxes and cos-mogenic nuclide production in the Earth’s atmosphere. J. Geophys. Res. , ,D11103. doi: 10.1029/2008JD010557Masarik, J., & Reedy, R. C. (1995). Terrestrial cosmogenic-nuclide production sys-tematics calculated from numerical simulations. Earth Planet. Sci. Lett. , , 381-395. doi: 10.1016/0012-821X(95)00169-DMesick, K. E., Feldman, W. C., Coupland, D. D. S., & Stonehill, L. C. (2018).Benchmarking Geant4 for Simulating Galactic Cosmic Ray Interactions WithinPlanetary Bodies. Earth Space Sci. , (7), 324-338. doi: 10.1029/2018EA000400Michel, R. (2005). Tritium in the hydrological cycle. In P. Aggarwal, J. Gat, &K. Froehlich (Eds.), Isotopes in the Water Cycle. Past, Present and Future of aDeveloping Science (p. 53-66). Dordtrecht: Springer.Mishev, A. L., Kocharov, L. G., & Usoskin, I. G. (2014). Analysis of the groundlevel enhancement on 17 may 2012 using data from the global neutron moni-tor network.
J. Geophys. Res.: Space Phys. , (2), 670–679. Retrieved from http://dx.doi.org/10.1002/2013JA019253 doi: 10.1002/2013JA019253Nevalainen, J., Usoskin, I., & Mishev, A. (2013). Eccentric dipole approximation ofthe geomagnetic field: Application to cosmic ray computations. Adv. Space Res. , (1), 22-29. doi: http://dx.doi.org/10.1016/j.asr.2013.02.020Nir, A., Kruger, S. T., Lingenfelter, R. E., & Flamm, E. J. (1966). Natural Tritium. Rev. Geophys. Space Phys. , , 441-456. doi: 10.1029/RG004i004p00441O’Brien, K. (1979). Secular variations in the production of cosmogenic iso-topes in the earth’s atmosphere. J. Geophys. Res. , , 423-431. doi: 10.1029/JA084iA02p00423Palcsu, L., Morgenstern, U., S¨ultenfuss, J., Koltai, G., L´aszl´o, E., Temovski, M., . . . –15–anuscript submitted to JGR: Atmospheres
Jull, A. J. (2018). Modulation of Cosmogenic Tritium in Meteoric Precipitationby the 11-year Cycle of Solar Magnetic Field Activity.
Sci. Rep. , , 12813. doi:10.1038/s41598-018-31208-9Picone, J. M., Hedin, A. E., Drob, D. P., & Aikin, A. C. (2002). Nrlmsise-00 empir-ical model of the atmosphere: Statistical comparisons and scientific issues. J. Geo-phys. Res.: Space Phys. , (A12), SIA 15-1–SIA 15-16. Retrieved from http://dx.doi.org/10.1029/2002JA009430 (1468) doi: 10.1029/2002JA009430Poluianov, S. V., Kovaltsov, G. A., Mishev, A. L., & Usoskin, I. G. (2016). Pro-duction of cosmogenic isotopes Be, Be, C, Na, and Cl in the atmosphere:Altitudinal profiles of yield functions.
Journal of Geophysical Research (Atmo-spheres) , , 8125-8136. doi: 10.1002/2016JD025034Raukunen, O., Vainio, R., Tylka, A. J., Dietrich, W. F., Jiggens, P., Heynderickx,D., . . . Siipola, R. (2018). Two solar proton fluence models based on ground levelenhancement observations. Journal of Space Weather and Space Climate , (27),A04. doi: 10.1051/swsc/2017031Smart, D., & Shea, M. (2009). Fifty years of progress in geomagnetic cutoff rigiditydeterminations. Adv. Space Res. , (10), 1107–1123. Retrieved from doi: http://dx.doi.org/10.1016/j.asr.2009.07.005Smart, D. F., Shea, M. A., & Fl¨uckiger, E. O. (2000). Magnetospheric Modelsand Trajectory Computations. Space Sci. Rev. , , 305-333. doi: 10.1023/A:1026556831199Sykora, I., & Froehlich, K. (2010). Radionuclides as Tracers of Atmospheric Pro-cesses. In K. Froehlich (Ed.), Environmental Radionuclides: Tracers and Timersof Terrestrial Processes (Vol. 16, p. 51-88). Amsterdam: Elevier.Tatischeff, V., Kozlovsky, B., Kiener, J., & Murphy, R. J. (2006). Delayed X- andGamma-Ray Line Emission from Solar Flare Radioactivity.
Astrophys. J. Suppl. , , 606-617. doi: 10.1086/505112Th´ebault, E., Finlay, C. C., Beggan, C. D., Alken, P., Aubert, J., Barrois, O., . . .Zvereva, T. (2015). International Geomagnetic Reference Field: the 12th genera-tion. Earth Planet. Space , , 79. doi: 10.1186/s40623-015-0228-9Usoskin, I. G., Alanko-Huotari, K., Kovaltsov, G. A., & Mursula, K. (2005). Helio-spheric modulation of cosmic rays: Monthly reconstruction for 1951–2004. J. Geo-phys. Res. , , A12108. doi: 10.1029/2005JA011250Usoskin, I. G., Gil, A., Kovaltsov, G. A., Mishev, A. L., & Mikhailov, V. V. (2017).Heliospheric modulation of cosmic rays during the neutron monitor era: Calibra-tion using PAMELA data for 2006-2010. J. Geophys. Res. , , 3875-3887. doi:10.1002/2016JA023819Usoskin, I. G., & Kovaltsov, G. A. (2008). Production of cosmogenic Be isotope inthe atmosphere: Full 3-D modeling.
J. Geophys. Res. , (D12), D12107. doi: 10.1029/2007JD009725Usoskin, I. G., Kovaltsov, G. A., Adriani, O., Barbarino, G. C., Bazilevskaya, G. A.,Bellotti, R., . . . Zverev, V. G. (2015). Force-field parameterization of the galac-tic cosmic ray spectrum: Validation for Forbush decreases. Adv. Space Res. , ,2940-2945. doi: 10.1016/j.asr.2015.03.009Usoskin, I. G., Mironova, I. A., Korte, M., & Kovaltsov, G. A. (2010). Regional mil-lennial trend in the cosmic ray induced ionization of the troposphere. J. Atmos.Solar-Terrestr. Phys. , , 19-25. doi: 10.1016/j.jastp.2009.10.003Vos, E. E., & Potgieter, M. S. (2015). New Modeling of Galactic Proton Modulationduring the Minimum of Solar Cycle 23/24. Astrophys. J. , , 119. doi: 10.1088/0004-637X/815/2/119Webber, W. R., Higbie, P. R., & McCracken, K. G. (2007). Production of the cos-mogenic isotopes H, Be, Be, and Cl in the Earth’s atmosphere by solar andgalactic cosmic rays.
Journal of Geophysical Research (Space Physics) , (A11), –16–anuscript submitted to JGR: Atmospheres
A10106. doi: 10.1029/2007JA012499Wilcox, L. J., Hoskins, B. J., & Shine, K. P. (2012). A global blended tropopausebased on ERA data. Part I: Climatology.
Q. J. R. Meteorol. Soc. , (664), 561-575. doi: 10.1002/qj.951(664), 561-575. doi: 10.1002/qj.951