A New Full 3-D Model of Cosmogenic Tritium $^3$H Production in the Atmosphere (CRAC:3H)

A new model of cosmogenic tritium ($^3$H) production in the atmosphere is presented. The model belongs to the CRAC (Cosmic-Ray Atmospheric Cascade) family and is named as CRAC:3H. It is based on a full Monte-Carlo simulation of the cosmic-ray induced atmospheric cascade using the Geant4 toolkit. The CRAC:3H model is able, for the first time, to compute tritium production at any location and time, for any given energy spectrum of the primary incident cosmic ray particles, explicitly treating, also for the first time, particles heavier than protons. This model provides a useful tool for the use of $^3$H as a tracer of atmospheric and hydrological circulation. A numerical recipe for practical use of the model is appended.


Introduction
Tritium ( 3 H) is a radioactive isotope of hydrogen with the half-life time of approximately 12.3 years. As an isotope of hydrogen, it is involved in the global water cycle and forms a very useful tracer of atmospheric moisture (e.g., Sykora & Froehlich, 2010;Juhlke et al., 2020) or hydrological cycles (Michel, 2005). In the natural environment, tritium is mostly produced by galactic cosmic rays (GCR) in the atmosphere, as a sub-product of the induced nucleonic cascade, and is thus a cosmogenic radionuclide. On the other hand, tritium is also produced artificially in thermonuclear bomb tests. Before the nucleartest ban became in force, a huge amount of tritium had been produced artificially and realised into the atmosphere, leading to an increase of the global reservoir inventory of tritium 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 typically neglected as being too small against anthropogenic one. However, as nearly 60 years have passed since the nuclear tests, its global content has reduced to the natural pre-bomb level (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 atmospheric tracing, water vapour transport, dynamics of the stratosphere-troposphere exchanges over Antarctica (Cauquoin et al., 2015;Fourré et al., 2018;Palcsu et al., 2018;Juhlke et al., 2020;László et al., 2020). Moreover, a combination of the 3 H data with other tracers like atmospheric 10 Be, which is also produced by cosmic-ray spallation reactions, but whose transport is different, can be a very powerful research tool. For this purpose, a reliable production model is needed, which is able to provide a full 3D and time variable production of tritium in the atmosphere.
Some models of tritium production by cosmic rays (CR) in the atmosphere were developed 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 methods of modelling the cosmic-ray induced atmospheric cascade. Later, a full Monte-Carlo simulation of the cosmogenic isotope production in the atmospheric cascade had been developed (Masarik & Beer, 1999) leading to higher accuracy of the results. However, that model had some significant limitations: (1) were considered only GCR protons (heavier GCR species were treated as scaled protons); (2) the energy spectrum of GCR was prescribed; (3) only global and latitudinal zonal mean productions were presented, implying no spatial resolution. That model was slightly revisited by Masarik & Beer (2009), but the methodological approach remained the same. A more recent tritium production model developed by Webber et al. (2007) is also based on a full Monte-Carlo simulation of the atmospheric cascade and was built upon the yield-function approach which allows dealing with any kind of the cosmic-ray spectrum. However, only columnar (for the entire atmospheric column) production was provided by those authors, making it impossible to model the height distribution of isotope production. Moreover, that model was dealing with CR protons only, while the contribution of heavier species to cosmogenic isotope production can be as large as 40% (see section 3).
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. This model belongs to the CRAC (Cosmic-Ray Atmospheric Cascade) family and is named as CRAC:3H. The CRAC:3H model is able, for the first time, to compute tritium production at any location and time, for any given energy spectrum of the primary incident CR particles, explicitly treating, also for the first time, particles heavier than protons. This model provides a useful tool for the use of 3 H as a tracer of atmospheric and hydrological circulation.

Production model
The local production rate q of a cosmogenic isotope, in atoms per second per gram of air, at a given location with the geomagnetic rigidity cutoff P c and the atmospheric depth h can be written as where J i (E) is the intensity of incident cosmic-ray particles of the i-th type (characterized by the charge Z i and atomic mass A i numbers) in units of particles per (s sr cm 2 GeV), Y i (E, h) is the isotope yield function in units of (atoms sr cm 2 g −1 , -see section 2.1 for details), E is the kinetic energy of the incident particle in GeV, h is the atmospheric depth in units of (g/cm 2 ), E c,i = (Z i · P c /A i ) 2 + E 2 0 − E 0 is the energy corresponding to the local geomagnetic cutoff rigidity for a particle of type i, and the summation is over the particle types. E 0 = 0.938 GeV is the proton's rest mass. The geomagnetic rigidity cutoff P c quantifies the shielding effect of the geomagnetic field and can be roughly interpreted as a rigidity/energy threshold of primary incident charged particles required to imping on the atmosphere (see formalism in Elsasser, 1956;Smart et al., 2000).

Production function
Here we computed the tritium production function in a way similar to our previous works in the framework of the CRAC-family models (e.g., Usoskin & Kovaltsov, 2008;Kovaltsov et al., 2012;Poluianov et al., 2016), viz. by applying a full Monte-Carlo simulation of the cosmic-ray induced atmospheric cascade, as briefly described below. Full description of the nomenclature and numerical approach is available in Poluianov et al. (2016).
The yield function Y i (E, h) (see equation 1) of a nuclide of interest provides the number of atoms produced in the unit (1 g/cm 2 ) atmospheric layer by incident particles of type i (e.g., cosmic ray protons, α-particles, heavier species) with the fixed energy E and the unit intensity (1 particle per cm 2 per steradian). The yield function should not be confused with the so-called production function S i (E, h), which is defined as the number of nuclide atoms produced in the unit atmospheric layer per one incident particle with the energy E. In a case of the isotropic particle distribution, these quantities are related as where π is the conversion factor between the particle intensity in space and the particle flux at the top of the atmosphere (see, e.g., chapter 1.6.2 in Grieder, 2001).
The production function in units (atoms cm 2 /g) can be calculated, for the isotropic flux of primary CR particles of type i, as -3-manuscript submitted to JGR: Atmospheres 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 2 /g), N i,l (E, E ′ , h) and v l (E ′ ) are concentration and velocity of the secondary particles of type l with energy E ′ at depth h. The aggregate cross-section η l (E ′ ) is defined as 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 total cross-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 by spallation of target nuclei of nitrogen and oxygen, which have the values of κ N = 3.22· 10 22 g −1 and κ O = 8.67·10 21 g −1 , respectively. The reactions yielding tritium are caused mainly by the cascade neutrons and protons and include: N(n,x) 3 H; N(p,x) 3 H; O(n,x) 3 H; O(p,x) 3 H. The cross-sections used here were adopted from Nir et al. (1966) and Coste et al. (2012), as shown in Figure 1a. We assumed that cross-sections of the neutron-induced reactions are similar to those for protons above the energy of 2 GeV. For reactions caused by α-particles, N(α,x) 3 H and O(α,x) 3 H, the cross-sections were assessed from proton ones according to Tatischeff et al. (2006). These reactions are induced mostly by α-particles from 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. Although production efficiencies of protons and neutrons are similar at high energies, they differ significantly in the <500 MeV range. Because of the lower energy threshold and higher cross-sections for neutrons in this energy range, comparing to protons, tritium production is dominated by neutrons in a region where the cascade is fully-developed, viz., in the 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 atmosphere by energetic cosmic-ray particles. The general computation scheme was similar to that applied by Poluianov et al. (2016). The simulation code was based on the Geant4 toolkit v.10.0 (Agostinelli et al., 2003;Allison et al., 2006). In particular, we used the physics list QGSP BIC HP (Quark-Gluon String model for high-energy interactions; Geant4 Binary Cascade model; High-Precision neutron package) (Geant4 collaboration, 2013), which was shown to describe the cosmic ray cascade with sufficient accuracy (e.g., Mesick et al., 2018). We simulated a real-scale spherical atmosphere with the inner radius of 6371 km, height of 100 km and thickness of 1050 g/cm 2 . The atmosphere was divided into homogeneous spherical layers with the thickness ranging from 1 g/cm 2 (at the top) to 10 g/cm 2 near the ground. The atmospheric composition and density profiles were taken according to the atmospheric model NRLMSISE-00 (Picone et al., 2002). Cosmic rays were simulated as isotropic fluxes of mono-energetic protons and α-particles, while heavier species were considered as scaled α−particles (see section 2.2). The simulations were performed 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 the isotope production should be better than 1% in any location. This number varied from 1000 incident particles for α−particles with the energy of 100 GeV/nucleon to 2 · 10 7 simulations for 20-MeV protons. The results were extrapolated to higher energies, up to 1000 GeV/nuc, by applying a power law. The yield of the secondary particles (protons, neutrons and α-particles) at the top of each atmospheric layer was recorded as histograms with the spectral (logarithmic) resolution of 20 bins per one order of magnitude in the range of the secondary particle's energy between 1 keV and 100 GeV. The primary CR particles were also recorded in the same way.
The production functions S i (E, h) were subsequently calculated from the simulation results, using equation (3), for a prescribed grid of energies and atmospheric depths and are tabulated in the Supporting Information. Some examples of the tritium production function are shown in Figure 2 for primary CR protons. One can see in Figure 2a that the efficiency of tritium atom production grows with the energy of the incident particles because of larger atmospheric cascades induced. Contributions of different components 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 protons depicts a double-bump structure: the bump in the upper atmospheric layers (h <10 g/cm 2 ) is caused by spallation reactions caused mostly by the primary protons (as indicated by the red dotted curve), while the smooth curve at higher depths is due to secondary neutrons (red dashed curve). Overall, production of tritium at depths greater than 10 g/cm 2 is very small for the low-energy primary protons. On the other hand, higher-energy (1 GeV, blue curves in Figure 2b) protons effectively form a cascade reaching the ground, where the contribution of secondary neutrons dominates below ≈ 50 g/cm 2 depths. This type of the depth/altitude profiles or the tritium production function was not studied in earlier works, where only columnar functions, viz. integrated over the full atmospheric column, were presented (Webber et al., 2007). Therefore, in order to compare our results with the earlier published ones, we also calculated the columnar production function where h sl = 1033 g/cm 2 is the atmospheric depth at the mean sea level or the thickness of the entire atmospheric column. The columnar production function is tabulated in the Supporting Information and depicted in Figure 3 along with the earlier results published by Webber et al. (2007) for incident protons. No results for incident α-particles have been published earlier, and the production function of cosmogenic tritium by cosmicray α-particles is presented here for the first time. One can see that, while the production functions for incident protons generally agree between our work and the results by Webber et al. (2007), there are some small but systematic differences. In particular, our result 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 production of tritium is negligible because of the geomagnetic shielding in such a way that low-energy incident 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 is higher than that from Webber et al. (2007). The difference is not large, ≈ 30%, but systematic and can be related to the uncertainties in the cross-sections or details of the cas- cade simulation (FLUKA vs. Geant4). Overall, our model predicts slightly higher production of tritium than the one by Webber et al. (2007), for the same cosmic-ray flux.

Cosmic-ray spectrum
The first term J i (E) in equation (1) refers to the spectrum of differential intensity of the incident cosmic-ray particles. A standard way to model the GCR spectrum for practical applications is based on the so-called force-field approximation (Gleeson & Axford, 1967;Caballero-Lopez & Moraal, 2004;Usoskin et al., 2005), which parameterizes the spectrum with reasonable accuracy even during disturbed periods, as validated by direct in-space measurements (Usoskin et al., 2015). In this approximation, the differential energy spectrum of the i-th component of GCR near Earth (outside of the Earth's magnetosphere and atmosphere) is parameterized in the following form: 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 per nucleon, E 0 is the rest energy of a proton (0.938 GeV), and Φ i (t) ≡ φ(t)·Z i /A i is the modulation 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 moment of time t is fully determined by a single time-variable parameter φ(t), which has the dimension of potential (typically given in MV or GV) and is called the modulation potential. The absolute value of φ makes no physical sense and depends on the exact shape of LIS (see discussion in Usoskin et al., 2005;Herbst et al., 2010Herbst et al., , 2017Asvestari 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:  where J LIS (E) is the differential intensity of GCR protons in the local interstellar medium in units of particles per (s sr cm 2 GeV), E and β = v/c are the particle's kinetic energy (in GeV) and the velocity-to-speed-of-light ratio, respectively. Following a recent work (Koldobskiy et al., 2019) based on a joint analysis of data from the space-borne experiment AMS-02 (Alpha Magnetic Spectrometer) and from the ground-based neutronmonitor network, we assumed that LIS (in the number of nucleons) of all heavier (Z ≥2) GCR species can be represented by the LIS for protons scaled with a factor of 0.353 for the same energy per nucleon.
The integral production rate in the entire atmospheric column is called the columnar production rate. For a given location, characterized by the geomagnetic cutoff rigidity P c , and at the time moment t it is defined as The global production rate Q global is the spatial average of Q C (P c ) over the globe, while the integral of Q over the globe yields the total production of tritium.
Production of tritium by GCR, which always bombard the Earth's atmosphere, is described above. Production by solar energetic particles (SEP) can be computed in a similar way, with the SEP energy spectrum entering directly in equation (1).

Results
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 different levels of solar modulation (low, moderate and high), for the entire atmosphere and only for the troposphere. The results are shown in Table 1. The modeled local production rates q(h, P c ) (equation 1) used for the computation can be found in a tabular form in the supporting information. Table 1. Tritium production rates (in atoms/(s cm 2 )) averaged globally (see also Figure 5) and over the polar regions (geographical latitude 60 • -90 • ), separately in the entire atmosphere and only the troposphere for different levels of solar activity: low, medium and high (φ=400, 650 and 1100 MV, respectively). The values of the modulation potential correspond to the formalism described in section 2.2. The geomagnetic field is taken according to IGRF (International Geomagnetic Reference Field, Thébault et al., 2015) for the epoch 2015. The tropopause height profile is adopted from Wilcox et al. (2012) The global production rate of tritium for a moderate solar activity (φ = 650 MV), which is the mean level for the modern epoch , is 0.345 atoms/(s cm 2 ). This value can be compared with earlier estimates of the global production rate of tritium. We performed a literature survey and found that the estimates performed before 1999 were based on different approximated approaches and vary by a factor of 2.5, between 0.14-0.36 atoms/(s cm 2 ) (Craig & Lal, 1961;Nir et al., 1966;O'Brien, 1979;Masarik & Reedy, 1995). Modern estimates, based on full Monte-Carlo simulations, are more constrained. The early value of the global production rate of 0.28 atoms/(s cm 2 ) published by Masarik & Beer (1999) was revised by the authors to 0.32 atoms/(s cm 2 ) in Masarik & Beer (2009). Our value is very close to that, despite the different computational schemes and assumptions made. The computed global production rate also agrees with the estimates 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 2 ). We note that heavier-than-proton primary incident particles contribute about 40% to the global production of tritium, in the case of GCR, and thus, it is very important to consider these particles explicitly.
Geographical distribution of the columnar production rate Q C (P c ) of tritium is shown in 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-cutoff spot in the equatorial west-Pacific region and polar regions.
Dependence of the global production rate of tritium on solar activity quantified via the modulation potential φ is shown in Figure 5, both for the entire atmosphere and for the 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 contribution to the global production is not dominant, because of the small area of the polar regions. Figure 6 (upper panel) presents the production rate of tritium in latitudinal zones (integrated over longitude in one degree of geographical latitude) as a function of geographical latitude and atmospheric depth. It has a broad maximum at midlatitudes (40-70 • ) in the stratosphere (10-100 g/cm 2 of depth) and ceases both towards the poles and ground. The bottom panel of the Figure depicts the zonal mean contribution (red curve) of the entire atmospheric column into the total global production. It illustrates that the distribution with a maximum at mid-latitudes shape is defined by two concurrent processes: the enhanced production (green curve) and reduced zonal area . Geographical distribution of the columnar production rate QC (atoms/(s cm 2 )) of tritium by GCR corresponding to a moderate level of solar activity (φ=650 MV). The geomagnetic cutoff rigidities were calculated using the eccentric tilted dipole approximation (Nevalainen et 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 the product of these two processes.
The altitude profile of the tritium production rate by GCR for the moderate level of solar activity is shown in Figure 7. The maximum of the globally averaged production is located at about 40 g/cm 2 or 20 km of altitude in the stratosphere, corresponding 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 reduced geomagnetic shielding there, so that lower-energy CR particles can reach the location. Figure 8 depicts temporal variability of the global tritium production for the period 1951-2018, computed using the model presented here. To indicate the solar cycle shape, the sunspot numbers are also shown in the bottom. The contribution from GCR is shown by the blue curve and computed using the modulation potential reconstructed from the neutron-monitor network . Red dots consider also additional production of tritium by strong SEP events, identified as ground-level enhancement (GLE) events (http://gle.oulu.fi). This is negligible on the long run but may contribute essentially on the short-time scale. Overall, the production of tritium is mostly driven by the heliospheric modulation of GCR as implied by obvious anti-correlation with the sunspot numbers.

Conclusion
A new full model CRAC:3H of tritium cosmogenic production in the atmosphere is presented. It is able to compute the tritium production rate at any location in 3D and for any type of the incident particle energy spectrum/intensity -slowly variable galac- tic cosmic rays or intense sporadic events of solar energetic particles. The core of the model is the yield/production function, rigorously computed by applying a full Monte-Carlo simulation of the cosmic-ray induced atmospheric cascade with high statistics and is tabulated in the Supporting Information. Using this tabulated function, one can straightforwardly and easily calculate the production of tritium for any conditions in the Earth's atmosphere (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. It also agrees well with empirical estimates of the tritium reservoir inventories, considering large uncertainties of the latter. Thus, for the first time, a reliable model is developed that provides a full 3D production of tritium in the atmosphere. These results can be used as an input for atmospheric transport models or for direct comparison with tritium observations that are important for the study of solar activity in link with the hydrological 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 Information, one can easily compute tritium production at any given location (quantified by the local geomagnetic rigidity cutoff P c and atmospheric depth h), and time t, following the numerical algorithm below.
1. For a given moment of time t, the intensity of incident primary particles can be evaluated, in case of GCR, using equations (6) and (7) for the independently known modulation potential φ (e.g., as provided at http://cosmicrays.oulu.fi/phi/ phi.html). These formulas can be directly applied for protons, while the contri- bution of heavier species (Z ≥2) can be considered, using the same formulas, but applying 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 incident protons J p (E, t) and heavier species J α (E, t), the latter effectively including all heavier species, can be obtained. Energy should be in units of GeV, and J(E) in units of nucleons per (sr cm 2 s GeV). The energy grid is recommended to 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 obtained, for both protons S p and heavier species S α , from the Supporting Information in units of (cm 2 /g). The yield function is defined as Y = π · S, in units of (sr cm 2 /g), also separately for protons and heavier species. The product of the yield function and the intensity of incident particles is called the response function 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 the lower integration bound in equation (1), needs to be calculated for a given location and time. A good balance between simplicity and realism is provided by the eccentric tilted dipole approximation of the geomagnetic field (Nevalainen et al., 2013). The value of P c in this approximation can be computed using a detailed numerical recipe (Appendix A in Usoskin et al., 2010). This approach works well with GCR, but is too rough for an analysis of SEP events, where a detailed computations 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 defined by the geomagnetic rigidity cutoff P c , as specified in equation (1) separately for the protons and α−particles (the latter effectively includes also heavier Z >2 species). Since the response function is very sharp, the use of standard methods of numerical integration, such as trapezoids, Gauss, etc., may lead to large uncertainties. For numerical integration of equation (1), the piecewise power-law approximation is recommended, as described below. Let function F (E) whose values are defined at grid points E 1 and E 2 as F 1 and F 2 , respectively, be approximated by a power law between these grid points. Then its integral on the interval between these grid points is E2 E1 F (E) · dE = (F 2 · E 2 − F 1 · E 1 ) · ln (E 2 /E 1 ) ln (F 2 /F 1 ) + ln (E 2 /E 1 ) .
The final production rate at the given location, atmospheric depth and time is the sum of the two components (protons and α−particles). 5. In a case when not only the very local production rate of tritium is required, but spatially integrated or averaged, the columnar production function (equation 8) can be used. The spatially averaged/integrated production can be then obtained by averaging/integration over the appropriate area considering the changes in the geomagnetic cutoff rigidity P c .
by the Academy of Finland (Projects ESPERA no. 321882 and ReSoLVE Centre of Excellence, no. 307411).