Reviewing 39 Ar and 37 Ar underground production in shallow depths with implications for groundwater dating

total depth-dependent production rate is estimated with Monte Carlo simulations assuming a uniform distribution of the parameter uncertainties. This work aims to provide a comprehensive framework for interpreting 39 Aractivitiesintermsofgroundwaterresidencetimesandforexposureagedatingofrocks.Theproductionof 37 Ar is also addressed since this isotope is relevant as a proxy for 39 Ar production, for the timing of river-groundwater exchanges, and in the context of on-site inspections (OSI) within the veri ﬁ cation framework of the Comprehensive Nuclear-Test-Ban Treaty (CTBT). In this perspective, we provide an interactive web-based application for the calculation of 37 Ar and 39 Ar production rates in rocks.

the water underground is calculated based on the radioactive decay from a known and relatively constant atmospheric activity of 16.4 ± 0.7 mBq m −3 air [atomic ratio: 39 Ar/Ar tot = 8.0 (± 0.6) × 10 −16 Benetti et al. (2007); Gu et al. (2021)].Dissolved 39 Ar activities larger than the atmospheric value have been observed in groundwater studies for many decades (Andrews et al., 1986;Loosli et al., 1989).Until recently, these underground-produced activities were attributed solely to neutron activation reactions of lithogenic potassium [ 39 K(n,p) 39 Ar] with neutrons originating from natural radioactivity in the rock [U and Th decay chains; Fabryka-Martin, 1988;Šrámek et al., 2017].Production resulting from cosmogenic neutron interactions was largely discarded because of the short penetration depth of these particles in regard to the typical depth of 39 Ar groundwater dating experiments and the short residence time of groundwater at shallow depths.Biases in the ages concluded from dissolved activities occur when the subsurface production processes are not considered properly.
Argon-37 ( 37 Ar; t 1/2 = 34.9days) activities in the atmosphere are low [0.5-10 mBq m −3 air , Purtschert et al., 2017] and therefore, infiltrating surface water contains virtually no 37 Ar.In the shallow underground (< 5 m), 37 Ar production is dominated by reactions with cosmogenic neutrons [ 40 Ca (n,α) 37 Ar], leading to depth-dependent concentrations in soil air and/or shallow groundwater (Riedmann and Purtschert, 2011).Thereby, the activity observed at a given depth encompasses information about the residence time and the percolation depth of the freshly recharged fraction of the groundwater.This can for instance be used for studying surface water -groundwater (SW-GW) interactions (Schilling et al., 2017;Delottier et al., 2022).In this context, 37 Ar closes the dating gap between 222 Rn (days -weeks) and 3 H/ 3 He (months -years).Assessing the natural 37 Ar background in soils and rocks is also crucial for distinguishing between artificial and natural sources for nuclear explosion monitoring in the frame of the Comprehensive Nuclear-Test-Ban Treaty [CTBT; Carrigan et al., 1996;Carrigan and Sun, 2014;Guillon et al., 2016].
Muons are cosmogenic secondary particles penetrating deeper than neutrons and capable of inducing nuclear reactions up to depths of several hundreds of meters (Aglietta et al., 1999;Braucher et al., 2003;Mei and Hime, 2006;Malgin and Ryazhskaya, 2008;Alfimov and Ivy-Ochs, 2009;Malgin, 2017;Kneißl et al., 2019).Muons were recently identified as significant contributors to 39 Ar underground production by the muon capture reactions on potassium [ 39 K(μ − ,γ) 39 Ar] in the context of Dark Matter research in deep underground laboratories (Mei and Hime, 2006).The 37 Ar production processes related to muon capture [ 39 K(μ − ,2n) 37 Ar] were also reported in early literature (Spannagel and Fireman, 1972).However, the importance of these processes for shallower groundwater dating applications was never recognized until now.
Thereby, 39 Ar and 37 Ar production reactions by neutrons take place on different target atoms (K and Ca respectively) whereas production by muon capture partly occurs on the same targets (Table 1).Thus, neutron and muon induced production leads to distinct 37 Ar/ 39 Ar production ratios.In the case of water transport from one production regime to another, additional evidence about the origin and mechanism of underground production can be derived due to their different half-lives.LLC allows for the measurement of those two radioargon isotopes in the same sample.
In the present paper, we review the most relevant production channels for 39 Ar and 37 Ar in typical rocks [as proposed by Fabryka-Martin, 1988] at depths between 0 and 200 meters below the surface (m.b.s).For each rock formation, the elemental weight fractions are summarized in Table A.1 in the supplementary material.The depth range considered is relevant for groundwater dating applications, as well as for the assessment of 37 Ar background profiles in the context of an on-site inspection (OSI) in the frame of the CTBT.The inclusion of muon-induced processes aims to provide a framework for the understanding of observed over-modern 39 Ar activities in groundwater, in particular in low U/Th rocks, as well as for the use of radioargon isotopes for studying erosion rates and exposure ages.
We distinguish the production processes associated with slow-muon capture reactions (Section 3) or neutron interactions (Section 4).For the processes involving neutrons, the calculation of neutron production rates and conversion to neutron fluxes is required.The energy dependency of the interaction process could be considered owing to recent numerical simulations of the neutron energy spectra with the code FLUKA (Grieger et al., 2020).At a given depth, the total radioargon production rate was calculated as the sum of the contributions from all the channels (Section 5).An interactive web-based application for the calculation of the production rate profiles for a rock with a custom elemental rock composition is available.Finally, the uncertainties of these calculations were addressed with Monte Carlo simulations (Section 6).

Overview
In the subsurface, radioargon isotopes ( 37 Ar and 37 Ar) are mainly produced by negative muon-capture and neutron activation reactions on potassium and calcium atoms (Table 1).The total production rate in the rock P Ar [atoms g −1 yr −1 ] as a function of lithospheric depth z [m] is thus the sum of the contributions from direct muon capture (P Ar,μ À ) and reactions with neutrons (P Ar,n ): Neutrons and muons are cosmogenic particles and their atmospheric flux depends on the altitude, latitude, and solar magnetic activity (Gordon et al., 2004).In the following, we use the reference scaling: a flat surface at high latitude (> 60°) and sea level altitude (SLHL) with a midsolar modulation.Deviations from these conditions can be addressed by scaling factors (Dunai, 2000).Particle fluxes in the underground are conventionally expressed in equivalent water depth [m.w.e] or [g cm −2 ] because particles are attenuated according to the total mass of matter traversed (1 m.w.e = 100 g cm −2 ).We use the rock densities ρ r and water-filled porosities φ provided by Fabryka-Martin (1988) to calculate the rock bulk densities ρ (in Table A.1), which controls the attenuation of the neutron flux in the underground.The conversion from the depth z [m] to h [m.w.e], is given by h ¼ z Á ρ=ρ water .The most relevant depth range for groundwater dating applications is 0 < z < 200 m of lithospheric depth since the production rates decrease by several orders of magnitude over this interval.
Muons have been extensively studied (i) in the context of underground deep laboratories (> 1 km.w.e; e.g. for dark matter research) because of the background they induce (Aglietta et al., 1995;Wang et al., 2001;Mei and Hime, 2006;Agafonova et al., 2011) and (ii) for geomorphic applications at shallow depths [0-100 m.w.e, (Stone et al., 1998;Schaller et al., 2002Schaller et al., , 2004;;Alfimov and Ivy-Ochs, 2009;Nishiyama et al., 2019;Lechmann et al., 2021)].On average, 100 μ m −2 s −1 reach the earth's surface at SLHL, with a mean energy of 4 GeV (Nakamura and Particle Data Group, 2010).As they travel through matter, they slow down and are ultimately stopped.The negative fraction of the muon flux can then be captured by the nucleus Coulomb field.Captured muons quickly cascade to the innermost electron orbit (1 s), where they either decay or interact with a proton of the nucleus (p + μ − → n + ν μ ).The resulting nucleus de-excites then by the emission of X neutrons, together with γ-rays.The total process is summarized by: where ν μ is a neutrino and X = 0, 1, 2, etc. is the number of neutrons released, which determines the product nucleus.This process is illustrated for radioargon products in Fig. 1.These muon-induced neutrons further contribute to radioargon underground production (Section 4.1.2).
Neutrons underground originate from lithogenic or cosmogenic sources.Lithogenic sources comprise neutrons from spontaneous fission and (α,n) reactions on light elements as part of U and Th decay chains (Šrámek et al., 2017).Cosmogenic sources include spallation and evaporation neutrons produced in the atmosphere by cosmic ray nucleons (the so-called "hadronic component of cosmic rays") and underground muon-induced processes.Here, we use the term "primary cosmogenic neutrons" to refer to the neutrons from the hadronic component of cosmic rays.Although the term "primary" is not strictly correct (because those are secondary particles in the atmosphere), it is opposed to muon-induced neutrons, which are also cosmogenic but result from secondary processes in the subsurface.Neutron interactions with the elements of the rock matrix are either high-energy spallation or nuclear capture of lower-energy neutrons (Fabryka-Martin, 1988;Smith, 1989).

Capture reactions of slow-negative muons
The production of radioargon isotopes by muon capture P Ar,μ-is the product of the local muon stopping rate R μ -(z) [μ − g −1 yr −1 ] and the yield of atoms produced per stopped negative muon For practical applications R μ -(z) can be approximated by the exponential depth-dependence (Eq.(B.6) in supplementary material).The rate of stopping muons at SLHL is R μ -(0) = 190 g −1 yr −1 (Heisinger et al., 2002a).The absorption-free path for stopped negative muons Λ μ [m] depends on the muon energy, which varies with depth.This effect is accounted for by considering k exponents with attenuation coefficients Λ k,μ [m] and relative intensities p k,μ [−] (Table B.1) for the change in the shape of R μ -(z), as proposed by Schaller et al. (2004) in Eq. ( 5) (Alfimov and Ivy-Ochs, 2009).
Once a muon stops, it can be captured by the atoms composing the matter in its vicinity.The yield of atoms of argon produced per stopped negative muon Y Ar μ À depends on four factors [Eq.( 6)], which are depicted in Fig. 1 and discussed in the following.The abundance of the target isotope of element i is referred to as f a in Eq. ( 6).
The chemical compound factor f c,i is the fraction of stopped muons reaching the 1 s muonic level of the target element i, which is one constituent of a compound.The factor f c,i [Eq.( 7)] is determined as the ratio of the probability to be captured on the target element i relative to the (arbitrary) reference element oxygen [W i , Eq. ( 8)], in relation to the probability to be captured by any other element in the material (von Egidy and Hartmann, 1982;Alfimov and Ivy-Ochs, 2009).In Eq. ( 7), N i , resp.N k [atoms g −1 ] are the elemental concentration of the target element i, or any other element k, in the sample.
Fig. 1.Slow negative muon capture reaction processes and factors for radioargon production.Slowed-down muons ultimately stop in the subsurface [with a stopping rate R μ -(z)] and are captured by the atoms composing the media (f c,i ).After the capture, they cascade into the innermost orbital where they can be captured by the nucleus (f d ).This reaction yields nuclides via the formation of an intermediate excited compound, whose de-excitation is susceptible to producing neutrons (f r ).
The nuclear capture probability f d is the probability for the muon in the 1 s level to be captured in the nucleus before decaying [Eq.( 9)].The probability for a muon to be removed from the 1 s level (λ S ) is the sum of the probability that it decays λ d and the probability that it is captured λ c .The decay probability of a muon in an atom (λ d = 1/ τ d , where τ d is the mean lifetime in the atom) is slightly reduced compared to the decay probability of a free muon (λ 0 = 1/τ 0 = 1/2200 ns): λ d = Q•λ 0 with the Huff Factor Q being in the range 0.8-1 (Suzuki et al., 1987).This marginal difference is neglected here (i.e.λ d ≈ λ 0 ).The factors λ S and λ c are calculated from the values compiled by Suzuki et al. (1987) and differ marginally (< 5 %) from those computed by Alfimov and Ivy-Ochs (2009).
The nuclear capture process produces a high-energy neutron, which either leaves the nucleus, ejects a particle from the nucleus (through direct interactions), or transfers its energy to other nucleons, thereby putting the nucleus in an excited state (Fabryka-Martin, 1988).Since the release of charged particles is impeded by the Coulomb barrier, most of the time the nucleus de-excites by the emission of one or more neutrons together with γ-rays (Gosse and Phillips, 2001).
The probability of reaction or neutron multiplicity distribution f r is then the probability for the excited (Z − 1, A)* nucleus to yield the desired (Z − 1, A − X) nucleus.In other words, this is the effective probability for the production of the nuclide i after the negative muon capture in the target nucleus (Heisinger et al., 2002a).This neutron multiplicity depends on the target nucleus through the momentum of the proton capturing the muon [Eq.( 2) and Eq. ( 3)].For our calculations, we refer to the compilation of reaction probabilities summarized by Fabryka-Martin (1988).
The yield of radioargon atoms produced per stopped muons Y Ar μ − [Eq. ( 6)] was calculated for each reaction channel for 39 Ar and 37 Ar (Table 1).The total production yield for each isotope, Y 39 μ − or Y 37 μ − , is the sum of the contributions from all channels.The factors f a , f r , and f d are specific to the target element and independent of the rock composition.The factor f c,i however depends on the elemental rock composition [Eq.( 6) and Eq. ( 7)] .Table C.1 in supplementary material summarizes these numbers for sandstone as an example.Table 2 shows the total production yields Y 39 μ − and Y 39 μ − , as well as the fraction of contribution of each reaction channel to the total Y Ar μ À for typical rocks.Depending on the composition of the rock, this contribution can vary substantially.Generally, for 39 Ar the production following muon capture (P Ar,μ ) is dominated by reactions on 39 K, except for ultramafic and carbonate rocks, where the reaction 40 Ca(μ − , p) 39 Ar accounts for 74 %, resp.64 % of the production by muon capture.For 37 Ar, both 39 K and 40 Ca reactions need to be accounted for in most of the rocks.

Reactions with neutrons from different origins
The radioargon production rate by neutrons P Ar,n [at g −1 yr −1 ] is: where N tg [at g −1 ] is the target element concentration, σ(E n ) [cm 2 ] is the energy-dependent reaction cross section, and φ n (E n ) [n cm −2 yr −1 ] is the energy-dependent neutron flux.The production of 39 Aris dominated by reactions with high-energy neutrons on 39 K (Fig. 2A).The production on 42 Ca is negligible due to the small reaction cross section at lower energies and low 42 Ca concentrations compared to 39 K.The isotope 37 Ar is chiefly produced on 40 Ca with fast neutrons ≥1 MeV (Fig. 2B).Unlike 39 Ar, 37 Ar might also be produced by thermal neutrons (note that the thermal cross section of the 40 Ca reaction is poorly constrained).In watersaturated media, 37 Ar can also be produced on 36 Ar dissolved in airequilibrated water (AEW).However, the reaction on 40 Ca is dominating in most cases due to the low 36 Ar abundance [ 36 Ar/Ar = 0.33 %] (Musy et al., 2022).
The total radioargon underground production rate induced by neutrons P Ar,n [Eq.( 11)] is the sum of the contributions from interactions with neutrons from (i) the hadronic component of cosmic rays P Ar(n,f) , (ii) slow negative muon-capture reactions P Ar(n,μ − ) , (iii) fast-muon induced reactions P Ar(n,μf) , and (iv) from spontaneous fission and (α,n) reactions on light elements, with α's from U\ \Th decay chains P Ar(n,U/Th) :

Neutron production rates
The overall production rate of neutrons at a given lithospheric depth z, P n (z) [n g −1 yr −1 ] is the sum of neutrons from all origins:

Hadronic component of cosmic rays
In the atmosphere and the shallow underground, neutrons are produced by interactions between galactic cosmic rays (GCRs, i.e. protons and secondary neutrons) and oxygen and nitrogen nuclei (Fabryka-Martin, 1988;Zanini et al., 2003).These neutrons are emitted with an average energy of 1 MeV but can range up to 10 MeV or more (Fabryka-Martin, 1988).The fast neutron flux from the hadronic component of cosmic rays drops sharply with depth withan attenuation length of Λ n,f = 160 g cm −2 ≈ 0.61 m in standard rock (Gosse and Phillips, 2001).The neutron production rate induced by the moderated atmospheric neutron flux, P n,f [n g −1 yr −1 ], as a function of lithospheric depth z [m], is then: Table 2 Yields of 39 Ar (Y 39 μ − Þ and 37 Ar ðY 37 μ − ) isotopes produced per stopped μ − [atoms (μ − ) −1 ] [Eq.( 6)] for typical rock formations and compositions (Table A.1) compiled by Fabryka-Martin (1988).The contribution of each production channel to the total yield is also indicated; is the average number of neutrons produced per muon captured in the rock [Eq.( 15)]; The neutron production rate from natural rock radioactivity P n,U/Th [n g −1 yr −1 ] was calculated with the code proposed by Šrámek et al. (2017).P n,0 = 2525 n g −1 yr −1 is the neutron production rate at the surface, normalized for SLHL (Heisinger et al., 2002a).It becomes insignificant compared to neutrons induced by muons or from natural radioactivity at depths >2-3 m of standard rock [∼5 m.w.e, Fig. 3A, Esch, 2001;Bogdanova et al., 2006;Garrison, 2014].

Neutrons from slow negative muon-capture reactions
The neutron production rate induced by muon-capture reactions in a given rock, P n,μ − (z) [n g −1 yr −1 ] [Eq.( 14)] is the product of the muon stopping rate R μ − [Eq.( 5)] and the average number of neutrons produced per negative muon stopped in the rock Y n μ À .After being captured by the atom of the target element i (f c,i ), and interacting with the nucleus (f d,i ), each muon yields an average ).The yield y n i is calculated from the probability distribution f r (X) for the emission of X neutrons following the muon capture reaction [Eq.(2); Charalambus, 1970].Most of the time, the capture reaction energy is dissipated by the emission of X ≤ 2 neutrons (Fabryka-Martin, 1988;Malgin and Ryazhskaya, 2008).The average neutron yield Y n μ À is obtained by summing up over all elements i contributing to the specific rock type (Eq.( 15), Table 2).The average energy of the emitted neutrons is 4-8 MeV (Fabryka-Martin, 1988;Malgin and Ryazhskaya, 2008;Grieger et al., 2020).

Fast muon-induced neutrons
Neutrons are also produced underground by fast muons via the following processes: a) muon interaction with nuclei via a virtual photon producing nuclear disintegration, b) muon elastic scattering with neutrons bound in nuclei, c) photo-nuclear reactions associated with electromagnetic showers generated by muons, d) (n,n) reactions with secondary neutron production following any of the above processes, (Wang et al., 2001;Niedermann, 2002;Malgin and Ryazhskaya, 2008;Palermo, 2016).
The fast muon-induced neutron production rate P n,μf (z) [n g −1 yr −1 ] [Eq.( 16)] can be calculated from the total muon flux ϕ μ (z) [cm −2 s −1 ] (Eq.(B.1) in supplementary material) and the average muon energy E μ (z) [GeV] [Eq.(D.1)] in standard rock. 1 We use the approach proposed by Heisinger et al. (2002b) in Eq. ( 16), which is based on a fit to the experimental data from Wang et al. (2001) to the number of produced neutrons per muon.In Eq. ( 16), α ≈ 0.75 and β(z) is an empirical factor that depends weakly on depth and can be approximated by Eq. (D.2).
1 Mass Number: A = 22; rock density: ρ r = 2.65 g cm −3 (Malgin, 2015) Even though most fast muon-induced neutrons are emitted with energies in the MeV range, their energy spectrum extends up to hundreds of GeV (Wang et al., 2001;Garrison, 2014;Palermo, 2016;Malgin, 2017).Generally, the energy distribution of these neutrons is uncertain since measuring energy spectra over a wide range is a complex task, and simulating the physical processes involved in neutron production from muons requires sophisticated models.Grieger et al. (2020) propose for the first time an assessment of the muon-induced neutron spectra in a shallow underground laboratory (140 m.w.e) based on a measured muon flux (Ludwig et al., 2019).We used the outputs from these simulations for the calculation of radioargon production in the following sections.

Neutrons from U and Th decay chains
In the uranium ( 238,235 U) and thorium ( 232 Th) decay series, α-particles are produced and interact with light elements of the host rock.These (α, n) reactions lead to the emission of neutrons, which contribute to  39 Ar and (C) 37 Ar production rate profiles; (D) 37 Ar/ 39 Ar production rate ratio (solid line) and 37 Ar/ 39 Ar normalized by the target concentration ratio Ca/K (dashed line) as a function of lithospheric depth z for all neutron and muon induced channels in typical clay-shales, sandstone, and carbonates rock formations (Table A.1 in supplementary material).radioargon production.Neutrons are also released during spontaneous fission of 238 U.The resulting total neutron production rate P n,U/Th [n g −1 yr −1 ] is calculated with the code developed by Šrámek et al. (2017) using the composition of typical rock types (Table A.1).The results are given in Table 2.
Fig. 3A summarizes the neutron production profiles P n (z) [Eq.( 12)] in the clay-shales, sandstone, and carbonate rock formations as examples.Except for the in-situ produced neutrons P n,U/Th , the production rates depend strongly on the depth of the rock overburden.Neutron production from the hadronic component of cosmic rays P n,f and fast-muon induced neutrons (P n,μ f ) are essentially independent of the rock elemental composition.Cosmogenic neutrons (P n; f , P n;μ − , P n,μ f ) dominate up to depths of about 25-50 m in the sandstone and carbonate formations.Since in-situ neutron production is almost an order of magnitude higher in the clay-shales (Table 2), this source prevails in depths deeper than ~5 m for this formation.

Neutron fluxes and energy spectra
The reaction cross sections for 37 Ar and 39 Ar production vary considerably with the neutron energy (Fig. 2A and Fig. 2B).The energy distribution of the subsurface neutron flux must therefore be considered when calculating radioargon production rates as a function of depth [Eq.( 10)].Depending on their origin, neutrons are produced with primary energies up to several tens of MeV.When they travel through matter (rock and water), they lose their energy mainly due to elastic scattering.They are finally thermalized at an energy of approximately 0.025 eV and ultimately absorbed by the surrounding material (Gosse and Phillips, 2001).The characterization of the resulting moderated secondary neutron energy spectrum is particularly difficult in the case of fast-muon-induced neutrons because of the complexity of their production channels.Secondary subsurface neutron energy spectra and fluxes can be estimated using Monte Carlo simulations (Fassò et al., 2000;Empl et al., 2014;Zhang and Mei, 2014;Grieger et al., 2020).
In the following, we use the neutron spectra in the Felsenkeller Tunnel IV in Germany (140 m.w.e) simulated by Grieger et al. (2020) with the code FLUKA (Fassò et al., 2000).These results include the muoninduced neutrons and the neutrons from U/Th decay chains [i.e.(μ,n) and (α,n) neutrons].A similar neutron energy distribution is assumed for the primary cosmogenic and the (μ,n) neutrons.Simulations were performed for three locations (Fig. 4A): two measurement chambers (MK1 and MK2) and a workshop room (WS).An averaged spectrum from MK1 and WS locations is used for our calculations (Fig. 2C; Section 6.2.2).In the outputs of the simulations, the (μ,n) neutrons include neutrons from fast muon and slow muon capture reactions.Both (α,n) and (μ,n) neutrons show a similar spectral shape from thermal energies up to 1 MeV, where the (μ,n) flux peaks while the (α,n) flux decreases.Above this energy, the (μ,n) spectrum extends to hundredths of MeV while the (α,n) flux drops.
The neutron flux φ n [n cm −2 yr −1 ] relates to the neutron production rates P n [n g −1 yr −1 ] following Eq.( 17) (Gosse and Phillips, 2001;Malgin, 2015).The parameter Λ n is the attenuation length for isotropic neutron flux in rocks.This factor depends weakly on the chemical composition of the media but rather on the initial neutron energy, the rock density, and its pore water content.For primary cosmogenic neutrons, Λ n,f = 160 g cm −2 is used.For muon-induced neutrons, we use the value Λ n,μ = 35 g cm −2 proposed by Malgin (2015).
The neutron production rate P n and thus, the total energy-integrated neutron flux depends on the chemical composition of rocks.Conversely, the neutron energy distribution (i.e. the shape of the neutron spectrum) is assumed to be relatively independent of the elemental composition (Lowrey, 2013).Implications of this assumption are discussed in Section 6.2.2.The neutron spectra for (μ,n) and (α,n) reactions simulated in the Felsenkeller Tunnel (Fig. 2C) are normalized by the total energyintegrated neutron flux at their respective location (i.e.MK1 or WS), resulting in dimensionless neutron spectra (i.e.[cm −2 s −1 ][cm −2 s −1 ] −1 ).In a subsequent step, these dimensionless spectra are multiplied with our calculated total neutron fluxes [Eq.( 17)] from all origins in the subsurface.The (μ,n) dimensionless spectrum is scaled with the neutron fluxes φ n, calculated with P n = P n,μf , P n,μ −, or P n,f .The (α,n) spectrum shape is used for the neutron flux calculated from P n = P n,U/Th [Eq.( 17)].The resulting φ n (z,E) [n cm −2 yr −1 ] are then folded with the production cross sections (Fig. 2A and Fig. 2B), resulting in the 39 Ar and 37 Ar production rates [Eq.( 10)].Among the reaction channels with neutrons (Table 1), only the reactions 39 K(n,p) 39 Ar and 40 Ca(n,α) 37 Ar are considered since the others contribute to < 1 % of the total production.

Total radioargon production rates
Fig. 3 summarizes the total production rates [P Ar , Eq. ( 1)] for 39 Ar (Fig. 3B), resp. 37Ar (Fig. 3C), resulting from muon capture reactions Fig. 4. Neutron energy spectra simulated with FLUKA in the Felsenkeller Tunnel Laboratory (Grieger, 2021).(A) (μ,n) and (α,n) neutrons energy spectra at the WS, MK1, and MK2 locations; (B) Bands between the energy spectra of the WS and MK1 locations for neutrons from (μ,n) [in purple] and (a,n) reactions [in pink].The spectra are normalized by the total energy-integrated neutron fluxes at each location.The solid and dashed bolder lines represent the mean fluxes.
[P Ar,μ-, Eq. ( 4)] and interactions with neutrons from all origins discussed previously [P Ar,n , Eq. ( 10)].The general pattern is the following: -The nucleonic component of cosmic rays is the dominant production channel in the first few meters of the lithosphere, regardless of the rock chemical composition (Fabryka-Martin, 1988).-At intermediate depths, muon-induced reactions become the major production path.This includes slow muon capture and reactions with muon-induced neutrons (Heisinger et al., 2002a;Heisinger et al., 2002b).-Ultimately, production with neutrons originating from the U and Th decay chains in rocks is the prevailing process (Šrámek et al., 2017).
For both 39 Ar and 37 Ar, the production with cosmogenic particles is dominating up to 70-200 m in the sandstone and carbonate formations.In contrast, in the clay-shales, the production with in-situ neutrons is the most important channel from 25 to 30 m because of the high natural radioactivity of this formation (U = 3.6 ppm, Th = 22 ppm).
In the range where muon-induced reactions dominate (from 5 m up to 25-200 m), the relative importance of production by direct muon capture (P Ar,μ À ; pink line in Fig. 3) or with the neutrons emitted following this process (P Ar n,μ À ð Þ ; orange line in Fig. 3) depends on the rock chemical composition.In most cases, these two production rates are comparable.The situation is yet different in the carbonate formation, where the Ca/K = 112.0 is 3-10 times larger than in the other formations (Table A.1).In this case, P Ar,μ À dominates up to 50 m because 39 Ar is produced by muon capture on both K and Ca targets.
Among the radioargon production with cosmogenic particles, the reaction with fast-muon-induced neutrons (P Ar n,μ f ð Þ ) is the most penetrating in the subsurface.For example in the sandstone formation, the number of neutrons produced by fast muons dominates up to depths of ∼ 50 m depth, where U/Th neutrons start to prevail (Fig. 3A).In greater depths, the fewer fast muon-induced neutrons have higher energies than in-situ neutrons.Because of the high reaction cross section at high neutron energies, the reactions with fast-muon-induced neutrons (P Ar n,μ f ð Þ ) prevail up to depths of 130 m for 39 Ar (Fig. 3B), resp.200 m for 37 Ar (Fig. 3C).
Overall, these findings are in agreement with Smith (1989), who concluded that radioargon production is dominated by cosmogenic secondaries up to ∼ 100 m depth.This contrasts with Mei et al. (2010) who stated that (i) 39 Ar is produced predominantly by muon-capture reactions over the reactions with muon-induced neutrons at all depths (in rock with Ca/K = 1.4), and (ii) cosmogenic particles are dominating 39 Ar production up to 2 km.w.e [i.e.∼ 750 m of standard rock, Fig. 4 in Mei et al., 2010].This discrepancy is likely linked to the different empirical parametrization for the depth-dependency of the total muon flux, as well as the other parameters that enter the calculation (e.g.muon energy, stopping rate ratio, etc.).The uncertainties on these parameters are discussed in Section 6 for our calculations.
The production rates of 39 Ar and 37 Ar following the reactions with neutrons are largely determined by the target concentration ( 39 K for 39 Ar, resp. 40Ca for 37 Ar).This is visible in the 37 Ar / 39 Ar production rate ratio (Fig. 3D) in the Ca-rich carbonate formation (Ca = 30 %wt, K = 0.03% wt), where the 37 Ar production rate is one to two orders of magnitude higher than 39 Ar (i.e. 37Ar / 39 Ar ≥ 10).In contrast, in the clay-shales (Ca = 2.5 % wt, K = 2.3 %wt) and sandstone (Ca = 4%wt, K = 1.1 % wt) formations, the production rates for 37 Ar and 39 Ar are comparable for the whole depth range (i.e. 37Ar / 39 Ar ≈ 1).When the production rate ratios are normalized by the neutron target element concentrations (i.e.Ca for 37 Ar and K for 39 Ar), the profiles (dashed lines in Fig. 3D) converge toward a value of 37 Ar / 39 Ar ∼ 0.3 at the surface and below 50 m for all rock types.At intermediate depths (i.e. between 4 and 50 m), the important production of 39 Ar by muon capture on Ca target reduces the normalized ratio in the Ca-rich formations.In the clay-shales formation the 37 Ar / 39 Ar norm is on the contrary increased over this intermediate depth range.This is linked with the preferential 37 Ar by the direct muon-capture reaction on K rather than Ca in this formation (Table 2).
An interactive web-based application allowing for online calculations of the production rates of neutrons, 39 Ar, and 37 Ar in rocks is provided in conjunction with this research paper.Details about the application, as well as related instructions, can be accessed online: https://www.climate.unibe.ch/radioargon

Uncertainties
Estimating the uncertainties in the calculated production rates as a function of depth is a complex task because of the numerous, not wellconstrained factors involved.This includes the rock chemical composition, the energy-dependent reaction cross sections, and the parametrization of complex nuclear reactions (Šrámek et al., 2017).The uncertainties of these parameters propagate through the production calculations.In the last decades, models have been proposed and discussed in the literature for the description of the flux and energy spectra of muons and neutrons in the subsurface for instance (Stone et al., 1998;Braucher et al., 2003;Mei and Hime, 2006;Malgin and Ryazhskaya, 2008;Malgin, 2017;Ludwig et al., 2019).
In the following, we estimate uncertainty bands for the calculated production profiles based on the assumptions that (a) the chemical composition measured in rock samples is representative of the whole formation (which is likely an optimistic perspective); and (b) that the production relevant parameters are independent and uncorrelated with uniformly distributed uncertainties.The uniform distribution was preferred over the normal one because the uncertainties are mainly based on the spread of the different models proposed in the literature (e.g. for the muon stopping rate).
The total uncertainty at a given depth is calculated with Monte Carlo iterations.The values of the parameter are randomly selected within the uniform distribution bounded by the uncertainty range (i.e. the value ± its uncertainty Δ).For each profile, 10′000 simulations were performed and used to compute the mean and standard deviation as a function of depth.Table 3 summarizes the uncertainties (Δ) for all parameters that enter into the calculation of the 39 Ar and 37 Ar production rate profiles.Justifications and explanations for the assumed ranges can be found in Sections 6.1 and 6.2.The total uncertainties along a depth profile P Ar(z) depend on the weight of the different production channels in the depth z and are discussed in Section 6.3.

Uncertainties of the muon-capture reaction
The two main factors involved in the muon-capture reactions are the muon stopping rate R μ − (z) and the yield Y Ar μ À of radioargon atoms produced per stopped muon for each reaction listed in Table 1.The uncertainty of Y Ar μ − was discussed before in the context of 36 Cl underground production by Alfimov and Ivy-Ochs (2009), which concluded an overall uncertainty of ±10 %.The stopping rate R μ -(z) was estimated by different models in the literature.Fig. E3 in supplementary material shows the comparison between the parametrizations of R μ -(z) from Schaller et al. (2004), which was used in the present study, and previously published approaches.Most models agree within 30 % up to 50 m depths, which is the relevant depth range in this study.Therefore, an uncertainty of 30 % was used for the muon-stopping rate.Further comparison of these parametrizations can be found in Supplementary Material in Section E.

Uncertainties of the production with neutrons
Radioargon production with neutrons depends on the magnitude of the different neutron sources (P n ; Eq. ( 12)), the primary neutron energy, the shape of the moderated energy spectra (for (μ,n) or (α,n) neutrons), and the energy-dependent reaction cross sections.The uncertainty of the latter is difficult to estimate since the simulated values provided in the nuclear data libraries do not include uncertainties.Šrámek et al. (2017) concluded an uncertainty of 30 % by comparing the experimental data available on the EXFOR database (https://www-nds.iaea.org/exfor)with values of integrated cross sections.Accordingly, this estimate of 30 % was adopted for the whole energy range in our Monte Carlo iterations.The target element concentrations (N tg ) measured by ICP-MS are assumed to be precise values.

Neutron production rates
The neutron flux and production rate associated with the hadronic component of cosmic rays (P n,f ) is well understood (Fabryka-Martin, 1988;Riedmann and Purtschert, 2011;Guillon et al., 2016) and becomes negligible after the first few meters of the lithosphere, which corresponds to the unsaturated zone in the vast majority of cases.Inaccuracies in the calculation of P n,f would have a minor impact on the activities measured in groundwater and were therefore not considered.The average neutron yield per captured stopped muon Y n μ − (Eq.( 15)) which controls P n,μ-has an uncertainty of about 10 % similar to Y Ar μ − : The neutron production rate induced by fast-muons P n,μf is a controversial topic in the literature (Palermo, 2016).Alternative parametrizations to the solution proposed by Heisinger et al. (2002b) for P n,μf (z) [Eq.( 16)] are discussed in Section E in Supplementary Material.An agreement within 50 % is observed between the different models and was used in our Monte Carlo iterations.
The uncertainty on neutron production rates from U/Th decay chains and spontaneous fission P n,U/Th (Table 2) concluded by Šrámek et al. ( 2017) is ±30 %.

Neutron energy spectra
The neutron energy distributions for the (μ,n) and (α,n) neutrons depend on the localization of the simulation in the Felsenkeller Tunnel (140 m.w.e ≈ 45 m depth in standard rock), as well as the rock water content.Grieger et al. (2020) performed simulations for three different sites (Fig. 4A), which were differently shielded from the surrounding hornblende rock: (i) Measurement chamber 1 (MK1) is shielded by serpentinite rock and pre-bomb steel with an estimated 160 g cm −2 total areal density; (ii) Measurement chamber 2 (MK2) is surrounded by steel, iron-pellets and lead with an estimated total areal density of 210 g cm −2 and (iii) the workshop (WS) is shielded by a 24 cm thick brick wall.In the simulations, the energy-integrated (α,n) neutron flux dominates over (μ,n) neutrons at WS and MK1 (Table 5.7 in Grieger, 2021).This is in agreement with the results obtained in this study at 45 m depth (Fig. 3A).The high-purity detector present in MK2 forms a significant target for (μ,n) neutrons.In this case, (μ,n) neutrons dominate over (α,n) neutrons.
Since the neutron spectra are normalized by the total energy-integrated neutron fluxes, the absolute shielding differences between the measurement sites are not so important for our application.However, the materials and detectors present in the tunnel cavities and the wall compositions affect the neutron energy distributions.The high-density shielding at MK2 favors the production of high-energy neutrons from (μ,n) reactions and leads to a distribution peaking at higher energy in comparison to the other locations (Fig. 4A), which are more representative of spectra shapes in rocks (Hebel, 2010;Šrámek et al., 2017;Purtschert et al., 2021).Hence, the spectrum from MK2 was not used in this study and instead, we used an average of the normalized fluxes at the WS and MK1 locations.The area between the curves (shaded area) defines the uncertainty range in the Monte Carlo iterations (Fig. 4B).
These spectra were simulated assuming a water content of 3 wt% while this value can go up to 12 wt% in aquifers [assuming fully saturated conditions, with φ = 30 %; Grieger et al., 2020].Higher water content has a minor effect on the total energy-integrated neutron flux.However, it would lead to a more effective neutron moderation, and thus, a shift of the spectra distribution toward lower energies.The production in the higher energy range, where the cross sections for both 39 Ar and 37 Ar are peaking (Fig. 2C) would hence be reduced.The fraction of production from neutron interactions (P Ar,n ) to the total radioargon production (P Ar ) would also decrease.More detailed simulations are required to parametrize this effect.

Total uncertainty of the production rates
The total uncertainties on radioargon production rates determined by the Monte Carlo simulations are shown in Fig. 5A and B. The full lines represent the mean (μ) of the 10'000 iterations and the ranges depict the 95 % (± 2σ) confidence interval.The uncertainty range changes with depth because the contributions of different production channels also change.At intermediate depth (25 < z < 150 m) the uncertainty range is larger because of the dominance of the less well-constrained muon-induced processes.This is particularly true in the sandstone formation, where the production with fast-muon induced neutrons is the prevailing production channel up to depths of 130 m for 39 Ar (Fig. 3B), resp.200 m for 37 Ar (Fig. 3C).The maximal relative uncertainty (2σ/μ) is 67 % for 39 Ar, resp.74 % for 37 Ar (right panel in Fig. 5).Our calculations should therefore be considered as a robust estimate of the order of magnitude of radioargon underground production.

Conclusions
Nuclear reactions producing radio noble gases and other isotopes in the underground have been investigated and used for a long time in earth science.For example, at the earth's surface, cosmogenic exposure dating is used to determine timescales of glacial advances and retreats.At greater depths, where "old" groundwater is expected to circulate, in-situ production of isotopes in the water and rocks was recognized as a disturbing interference or a formation-specific diagnostic tool for dating and flow path reconstruction purposes.In this study, we focus on the transition zone between these two situations, where both cosmogenic and in-situ production processes are relevant.We reevaluated the reaction channels producing 39 Ar and 37 Ar at intermediate depths where groundwater resides often only during the recharge process.This depth range has received little attention so far because the cosmogenic production was assumed to take place only at shallower depths, where it wouldn't substantially impact the initial 39 Ar activity.However, this might not always be the case with the consideration of slow infiltration and deep-reaching production reactions involving muons.In this context, our recompilation of the significant 39 Ar and 37 Ar underground production processes aims to provide production rate profiles in typical aquifer types including both cosmogenic and insitu channels.
In most rocks and up to depths of ~25-50 m, the neutron flux in the subsurface originates mainly from cosmogenic sources (i.e.hadronic component of cosmic rays, and muon-induced neutrons).Below this depth, lithogenic neutron sources dominate (i.e.natural U/Th decay chains in the rocks).The cross sections of the main 37 Ar and 39 Ar producing neutron reactions (i.e. 40Ca(n,α) 37 Ar and 39 K(n,p) 39 Ar) are energy-dependent, and therefore, the moderated neutron spectra of all sources must be determined.This was considered here by using the spectra shapes simulated with an integrated Monte Carlo code (FLUKA, Grieger et al., 2020) for (μ, n) and (α,n) reactions.
Our results show that muon-induced neutrons contribute significantly to 39 Ar and 37 Ar underground production in rocks.Radioargon production profiles calculated for clay-shale, carbonate, and sandstone formations showed a significant contribution of muon-induced reactions at depths between 3 m and 150-200 m.The occurrence of muon-capture reactions and the rock formations where it prevails can be constrained by the measurement of 37 Ar, which reflects local conditions due to its short halflife.The 37 Ar/ 39 Ar production ratios are different for muon-capture and neutron reactions when normalized by the main target elemental concentrations.
The computation of complex reactions and processes in environmental systems is inherently fraught with uncertainties.These were addressed by a careful literature review and Monte Carlo simulations, which resulted in relative uncertainties on the order of up to 70 %.This conservative estimate may not be sufficient for precise "groundwater exposure age dating" yet but provides an important basis for understanding the significance of the processes involved in specific field cases.In the future, more simulations in different rock formations and with various water contents will further improve our knowledge about subsurface production processes affecting radioargon groundwater dating.

Fig. 3 .
Fig. 3. (A) Neutron (B)39 Ar and (C) 37 Ar production rate profiles; (D)37 Ar/ 39 Ar production rate ratio (solid line) and 37 Ar/ 39 Ar normalized by the target concentration ratio Ca/K (dashed line) as a function of lithospheric depth z for all neutron and muon induced channels in typical clay-shales, sandstone, and carbonates rock formations (TableA.1 in supplementary material).

Table 1
Production channels for 39 Ar and 37 Ar in the underground by neutron interactions and slow negative muon capture.The * reactions are of minor importance and not considered in the underground production calculations in the following.

Table 3
Uncertainties Δ of the parameters used to compute the 39 Ar and 37 Ar production rate profiles.