The Effect of Grain Size on Porewater Radiolysis

The radiolysis of porewaters by uranium, thorium, and potassium in mineral grains is a recognized source of molecular hydrogen in rock‐hosted and sediment‐hosted fluids. This radiolytic hydrogen is of geomicrobiological interest as a potential energy source (electron donor) for microbial metabolism, especially in energy‐limited settings such as the marine deep biosphere or the subsurface of Mars. Previous efforts to predict the production of radiolytic hydrogen from columns of rock and sediment have tended to rely upon analytic models that cannot account for the attenuation of mineral radiation by grains larger than ∼30 µm. To address this, we have developed a Monte Carlo method to simulate the physics of mineral radiation and evaluate the production of H2 as a function of mineral grain size and radioisotope composition. The results confirm that grain size is a major control on radiolytic H2 yield. For example, using the standard geological classification of grain sizes, we find that clay can produce up to an order of magnitude more H2 per unit time than sand. The magnitude of this effect is illustrated using compositional data from real geological units in order to demonstrate the dependence of radiolytic hydrogen flux on natural radionuclide concentration and bulk porosity.

• Porewater radiolysis in rocks and minerals yields molecular hydrogen (H 2 ), an energy source for microbial life on Earth and perhaps Mars • Our Monte Carlo simulations of mineral radiation reveal how grain size controls H 2 yield for a given radioisotope composition and porosity • For some realistic scenarios, claysized grains can produce up to an order of magnitude more H 2 per unit time than sand-sized grains Hitherto, studies of porewater radiolysis have not properly accounted for the attenuation of each radiogenic particle or ray within its mineral grain (clast or crystal) of origin. This attenuation can be safely ignored when the grain diameter is much smaller than the stopping distance of the relevant particle in the mineral medium, as for example, in clay. However, alpha particles in the kinetic energy range characteristic of mineral radiation (∼5.6 MeV on average) typically have stopping distances in SiO 2 (and similar materials) of ∼23 μm (National Institute of Standards and Technology, 2021a). Consequently, if such an alpha particle is produced near the edge of a grain of diameter ≳20 μm (i.e., coarser than a medium silt) and travels toward the interior of the grain, it will not have enough kinetic energy to escape from the other side, resulting in a null contribution to the radiolysis of the surrounding porewater. Similar considerations apply to beta particles-although their stopping distances are substantially larger-given a characteristic mean kinetic energy of ∼0.33 MeV, equating to a stopping distance of ∼452 μm in SiO 2 (National Institute of Standards and Technology, 2021b). Furthermore, the radiolysis contribution due to gamma (photon) radiation requires special treatment given their chargeless, massless character. Here, we quantitatively account for the variability of grain size in order to expand the range of environments in which radiolytic H 2 production and its potential contribution to microbial productivity can be estimated. The method detailed in this paper is generalizable to the extent that it can be adapted for a variety of geological contexts and sediment or rock types with any level of radioactivity, whose constituent grains can be of arbitrary size and shape. In this particular study, we seek to model radiolytic H 2 production from porous rocks and sediments with grain sizes ≳20 μm for the first time (to our knowledge). One motivation for doing so is the fact that high concentrations of radioactive minerals on Earth commonly occur in relatively coarse sedimentary rocks, for example, sandstone-type uranium deposits. In addition, relatively coarse granular materials can accommodate more microorganisms, which may tend to be excluded by the restrictive pore throat sizes of finer-grained materials, even at the same porosity (e.g., Krumholz et al., 1997).

Materials and Methods
We have written a Fortran program that numerically simulates the physics of natural radioactive decay within a matrix of water and mineral grains using a Monte Carlo (MC) method. (The algorithmic foundation of this program is given in Supporting Information S1.) MC models of alpha attenuation in mineral grains have previously been deployed in some geochronological studies (e.g., Farley et al., 1996;Hourigan et al., 2005). Although analytical models are sometimes easier to implement than numerical models, an MC approach offers maximum flexibility since it can ultimately be adapted for grains of any size and shape, with any distribution of 10.1029/2021EA002024 3 of 21 radionuclides, any arrangement of pore space, etc. (Ketcham et al., 2011). In the present study, we aim to address the effect of grain size independent of other factors. For simplicity, grains are therefore represented as perfect spheres composed of pure SiO 2 . Real mineral grains typically deviate from a spherical geometry (especially clay minerals) but modeling the geometry in this way leverages the rotational symmetry of the system, which eases the conceptual and computational aspects of the MC method.
The primary goal of our simulations is to score the total amount of kinetic energy transferred to the surrounding water layer due to the decay of radioactive isotopes within the SiO 2 grains. Figure 1 shows the geometry of a spherical SiO 2 grain and a composite of sample thorium decay events occurring within it.
The program leverages the MC method and the physics of radioactive decay to drive a given simulation toward a well-averaged endpoint. The steps the simulation uses are as follows: 1. Choose a random SiO 2 grain whose radius is within the desired interval of clay, silt, or sand (Wentworth, 1922), or a logarithmic distribution of all three; 2. Choose a random K, Th, or U constituent (weighted by the proportion of that element in the composition; see Section 2.2); 3. Choose a random point ( , , ) within the grain at which the decay of the constituent plus its daughters will proceed; Figure 1. A randomly selected SiO 2 grain of radius 0.059 cm (590 µm) showing a sample thorium decay against , , and coordinate axes (black arrows with gradations). Alpha, beta, and gamma decay events register as red, green, and blue, respectively. These depth vectors locate exit points that correspond to solutions of the sphere 2 = 2 + 2 + 2 .
4. Examine the entire decay chain in order, choosing random exit vectors, and applying the appropriate physics of alpha, beta, and gamma propagation in SiO 2 ; 5. Repeat Steps 1-4 many times (typically 10 5 trials) in order to ensure the simulation converges while keeping a running total of the amount of kinetic energy (in MeV) deposited in the surrounding water layer; 6. Output radiolytic H 2 production from the volume of rock or sediment; 7. Using Step 6 as a baseline, scale H 2 production by volume-corrected porosity (see Section 2.4). Athy's law can be used to scale radiolytic H 2 production in the case where porosity can be estimated for a column of rock or sediment as a function of depth.
Point 2 assumes that all radioactive nuclides are uniformly distributed throughout each grain, that all grains are the same in U, Th, and K content. This is a useful simplification for examining the effect of grain size as an independent variable, but in real rocks radionuclides may be heterogeneously distributed within grains and most concentrated in particular minerals (e.g., zircon) that have a different grain size from the rest of the sediment or rock. For example, radioactivity is commonly associated with secondary diagenetic minerals such as iron oxides and clays that adsorb U and Th (Murray et al., 2014). Point 6 assumes that all radiolytic H 2 is generated from water and not from the radiolysis of hydrated or hydroxylated minerals, a less well understood phenomenon (Smetannikov, 2011;Westbrook et al., 2015). These considerations are not taken further in the present study but suggest possible avenues for future study.
The "length" of the simulation is controlled by the number of K, Th, and U decays + daughters that are tracked (read: trials), with the understanding that a large number of trials will enable the simulation to converge in a manner dictated by the physics of radioactive decay. For the purposes of this study, our simulations output (a) H 2 production in units of moles per [Earth] year cm 3 as a function of sediment depth in meters, or (b) H 2 production as a function of grain size (logarithmically distributed) in cm.

Radioactive Decay
This section describes how the components of natural radioactivity interact with the medium through which they are propagating. This is also referred to as stopping power or linear energy transfer (LET). Our simulations treat alpha, beta, and gamma decay piecewise, wherein each decay daughter in a given chain is assigned a random exit vector. The entire decay chains for 238 U, 235 U, 232 Th, and 40 K were hard-coded into the program using data from the Interactive Chart of Nuclides provided by Brookhaven National Laboratory (2021). (The average energies required for each transition within a given decay chain were implemented rather than the "endpoint" or maximum energies, as previous studies have done.) These vectors are then used to propagate each decay product through a geometrically computed depth of SiO 2 material before depositing the remainder of its kinetic energy in the surrounding water layer. This strategy is advantageous since it avoids the assumption of alpha, beta, and gamma dominance through the use of lump energy sums (Blair et al., 2007). Overall, alpha, beta, and gamma emission resulting from natural radioactive decay will occur with average kinetic energies of 5.63 ± 1.68, 0.334 ± 0.267, and 0.227 ± 0.639 MeV, respectively.

Alpha Particles
The interaction of a heavy charged projectile with the electrons in a stopping medium (absorber) is described by the Bethe-Bloch formula (Bethe, 1930). The amount of kinetic energy d imparted to the target medium per unit path length d is given as where 2 ∕4 0 is the coupling strength of electromagnetism, is the atomic number of the charged particle or ion projectile, e 2 is the rest mass-energy of the electron, A∕ is the electron density of the absorber with effective atomic number and effective atomic mass (when referring to compounds and mixtures), and is the mean ionization potential of the target medium. The dimensionless number ≡ ∕ = √ 1 − [ e 2 ∕ ( + e 2 ) ] 2 from special relativity. The last term in the square brackets ∕2 is the density effect correction that accounts for the phenomenon that, in a dense medium, the field which perturbs electrons far from the projectile track is modified by the dielectric polarization of the atoms between the distant electrons and the projectile (Fermi, 1940). The quantity d ∕d in Equation 1 is sometimes preceded by a negative sign to indicate that kinetic energy is transferred to the stopping medium.
Alpha particles (helium nuclei, = 2) that are emitted during the course of natural radioactive decay possess an average kinetic energy that is relatively low, or ∼5.6 MeV. This low kinetic energy necessitates a correction to the stopping power due to the phenomenon of electron capture, which accounts for the fact that the bare nuclear charge of the projectile is reduced and can thus be replaced by an effective projectile charge (Weaver & Westphal, 2002). Once corrected, Equation 1 is used to numerically propagate alpha particles through a known amount of SiO 2 with characteristics A∕ and . From a radiolysis standpoint, the production of molecular hydrogen per unit path length is greatest at the end of an alpha particle's range, a region commonly referred to as the Bragg peak in radiation physics.
A charged particle or ion projectile with atomic number of rest mass significantly larger than that of the electron (charge ) is considered a "heavy" charged particle. This would include mesons, protons, alpha particles, and higher nuclei. Electrons and positrons are considered "light" charged particles; the physics of electron propagation through matter is discussed in the next section.

Beta Particles
While the Bethe-Bloch formula can be used in principle to calculate the stopping power of any charged particle passing through matter, it must be modified in the case of beta particles (electrons). There are several reasons for this: 1. one must always consider special relativity when handling electron-electron interactions; 2. an electron that acts as a projectile can lose nearly all of its kinetic energy in a single interaction with a target electron; 3. quantum mechanics states that we will not be able to distinguish between the projectile and the target electrons after the collision (they are identical particles); 4. for a high energy electron that traverses a high stopping medium, radiative processes such as bremsstrahlung (from the German: "breaking radiation") will be significant.
The total stopping power of an electron as a function of kinetic energy in a specified target material is the sum of the collision plus bremsstrahlung contributions (Bethe & Heitler, 1934): where ⟨d ∕d ⟩c is the stopping power due to electron collisions, where is the fine-structure constant. Equations 2-4 are used to numerically propagate beta particles through a known amount of SiO 2 with characteristics similar to those described by Equation 1.
The simulation of radioactive decay in a mineral medium requires that Equations 1-4 be applied to the propagation of alpha and beta particles inside spherical SiO 2 grains. As a means of developing and validating our numerical method, it would be well if the algorithm given in Supporting Information S1 were checked against a pure geometric argument. It has been found analytically that if particles with stopping distance are emitted at random vectors from points inside a spherical grain of diameter Θ , the proportion (or percent fraction) of these particles that escape from the grain rather than terminating (or stopping) within it is given by Equation 5 has been obtained previously (e.g., Farley et al., 1996); a derivation is provided in Text S2 in Supporting Information S1. Figure 2 shows a graphical comparison of these two methods in the case of alpha particles, along with an extension to beta particles using the same algorithm given in Supporting Information S1.
The comparison shown in Figure 2 corresponds to an alpha particle test case with stopping distance = 19.6 μm (kinetic energy: 4.93 MeV). These curves disagree on average by 0.2% over 100,000 trials per micron radius. The good agreement between the core algorithm of the simulation and our geometric analysis partially validates our MC approach in this study.

Gamma Rays
Unlike energetic charged particles which have a finite range when traversing matter, gamma rays (photons) of a given energy do not possess such a discrete range, but rather are absorbed exponentially as a function of distance in the absorbing medium. This phenomenon is referred to as attenuation. (Alpha and beta particles, strictly speaking, do not attenuate; we use this term contextually to service the reader.) For a beam of 0 monoenergetic photons, the rate of absorption of the photons as a function of depth in the absorber is where is the absorption cross section of each atomic scattering center (or absorbing molecule), and is the volume density of atomic scattering centers. (The negative sign above represents absorption.) We define ≡ , where is referred to as the linear attenuation coefficient. We can then integrate Equation 6 with respect to distance to obtain Figure 2. Escape efficiency of alpha particles with stopping distance = 19.6 μm as a function of spherical grain radius, where the methods of MC simulation and geometric analysis are compared. The escape efficiency of beta particles is also shown using the same algorithm.
Similarly, if we consider depth in terms of not distance, but of areal density (in g/cm 2 ), we can define the mass attenuation coefficient such that where is the mass density (in g/cm 3 ), is the atomic molar mass (in g/mol) of the absorbing medium, and A is Avogadro's number (in molecules/mol), respectively. Equation 7 now becomes where the values of m are tabulated.
Using Equation 9, a characteristic attenuation (or absorption) length can be formulated as an analog to the stopping distance of alpha and beta particles; for a single photon, this is taken as the distance into a material when the probability of absorption has dropped to ∕ 0 = 1∕ , or ∼37%. Explicitly, this gives a propagation distance that is equal to the inverse of the mass attenuation coefficient, or 1∕ m . A typical attenuation length in this study is on the order of centimeters.
There are a number of different mechanisms by which photons can interact with matter. However, only three of these mechanisms-the photoelectric effect, Compton scattering, and electron/positron pair production-dominate, while the remainder usually make only a negligible contribution. The total absorption cross section in Equation 8 is the sum of the cross sections for each type of photon interaction: where PE is the cross section for photoelectric interactions, C is the cross section for Compton interactions, and PP is the cross section for pair production. (A cross section in nuclear physics can informally be read as an energy-dependent probability for interaction.) Since ∝ from Equation 8, there is a corresponding mass attenuation coefficient for each cross section.

Photoelectric Effect
Free electrons cannot fully absorb photons without violating conservation of energy and conservation of momentum, and for that reason, free electrons are not subject to the photoelectric effect. Indeed, the photoelectric effect should be thought of as interaction between an incident photon and an entire atom, not just one of its electrons (Einstein, 1905).
Photoelectric interactions are more likely to eject electrons from the more tightly bound inner electron shells of the atom than the less tightly bound outer shells, since interactions involving inner shell electrons occur closer to the center of the atom where it is easier to transfer an incident photon's linear momentum to the atom. This is because nearly all the mass of the atom is concentrated in its nucleus.
The cross section for the photoelectric effect is: where is the fine-structure constant, is the atomic number of the absorbing medium, e is the classical electron radius, the quantity e 2 is the rest mass-energy of the electron, and finally is the energy of the incident photon.
Equation 11 should be applied if the energy of the incident photon is between 0.1 and 0.35 MeV, or where the photon energy is above the K-absorption edge.

Compton Scattering
The interaction of photons with "nearly" free electrons is described by a number of mechanisms: 1. Photon interactions with "truly" free electrons is referred to as Thomson scattering; 2. Coherent (or elastic) photon interactions with "nearly" free electrons is referred to as Rayleigh scattering; 3. Incoherent (or inelastic) photon interactions with "nearly" free electrons is referred to as Compton scattering.

Conservation of energy for the Compton scattering interaction gives
] 2 from special relativity, with and ′ the energies of the incident and scattered photon, respectively. After obtaining the differential cross section for Compton scattering-which is specified by the Klein-Nishina formula (Klein & Nishina, 1929)-it can be shown that integration over all angles gives the cross section for the Compton interaction: where effectively compares the incident photon energy to the rest mass-energy of the electron, or ∕ e 2 .

Pair Production
Pair production is essentially the reverse process of bremsstrahlung. When in close proximity to a heavy nucleus, a photon can be converted into an electron/positron pair. Conservation of energy for this process dictates where e and p are the energy of the electron and positron, respectively. Because the rest mass of the electron and the positron must come from the energy of the photon, there is a threshold energy of 2 e 2 = 1.022 MeV below which this process cannot occur.
The cross section for pair production originates from quantum electrodynamics, and is given by (Gould & Schréder, 1967): where ] 2 from special relativity. In this case is actually the relativistic energy of the electron/positron in the center-of-mass frame of reference. This cross section is applied for the transitions 40 K → 40 Ar in 40 K (1.461 MeV), 208 Tl → 208 Pb in 232 Th (3.376 MeV), and 214 Bi → 214 Po in 238 U (1.476 MeV).

Activity of Bulk Rock or Sediment
The activity per unit volume vol of bulk rock/sediment in units of Bq/cm 3 is given by where is the mass density of the rock/sediment constituents in g/cm 3 , A is the Avogadro constant, is the molar mass of the radioisotope, g∕g is the concentration of the radioisotope in micrograms per gram (μg/g), and 1∕2 is the half-life of the radioisotope in seconds. Equation 16 is used to convert the concentration of a given radioisotope component (uranium, thorium, or potassium) to the volume-normalized activity in bulk sediment/rock.

Radiolysis
Hydrogen production is evaluated according to the water radiolysis model of Blair and coworkers (Blair et al., 2007): where , , and are the H 2 yields per unit kinetic energy produced in pure water by radiolysis, and , , and are the kinetic energies absorbed by water in bulk rock/sediment due to alpha, beta, and gamma decay, respectively. The values of , , and are 1.57E4, 5.30E3, and 4.50E3 molecules/MeV, respectively (Blair et al., 2007;Spinks & Woods, 1990); the values of , , and (in MeV) are determined by the simulation. These values are affected very little by changes in temperature, pressure, and pH (Draganic & Draganic, 1971). They are, however, enhanced to some extent by salinity due to the reaction of inorganic solutes with hydroxide (Lin et al., 2005;Spinks & Woods, 1990).

Porosity
A porosity correction is needed to account for the fact that only the kinetic energy deposited in the intergranular water (rather than other grains) will contribute to H 2 production. The porosity of a compressible rock or sediment column decreases with increasing depth as the grains are compacted together. The simulations presented in this paper assume that all the porosity is occupied by water; we use Athy's law to quantify porosity as a function of sediment (or rock) depth (Fowler & Yang, 1998): where 0 is the porosity near the surface and is the compaction coefficient and has units of m −1 . In sandstones on Earth, takes an average value close to 5.3E−4 m −1 (McMahon & Parnell, 2014); for geological features on Mars, this study uses a slightly lower compaction coefficient of 3.5E−4 m −1 (Clifford et al., 2010). Equation 18 represents a steady state leading-order solution to equilibrium compaction.
A volume correction factor is required in order to properly scale H 2 production as a function of depth. This factor, 1 − , accounts for the fact that in a given volume of rock/sediment, a large porosity will translate to a smaller proportion being occupied by radioactive mineral grains rather than water. The kinetic energy yield from a given volume of wet rock/sediment is equal to the energy yield of each grain (supplied by the simulation), multiplied by the total number of grains contained in the volume. This total grain number is equal to the volume of the rock/ sediment multiplied by a factor of 1 − , normalized by the volume of each grain (Blair et al., 2007). In other words, one must multiply Equation 17 by Equation 18 (Athy's law), corrected by 1 − = 1 − 0 − . The reader should note that the complete porosity correction used to scale H 2 production as a function of depth is equal to (1 − ) , which contains a global maximum at = 0.5. Figure 3 compares the predicted H 2 production given by the MC method with those calculated by Blair et al. (2007) from marine sediment column measurements as a function of depth. This figure demonstrates the compatibility of our simulations with the deep marine sediment column in which empirically determined radionuclide concentrations, grain sizes, and porosities (measured by the Ocean Drilling Program in the Peru Basin) were provided by Blair and coworkers. Figure 3, in conjunction with Figure 2, further validates our MC methodology. The left panel in Figure 3 displays a comparison after 200 MC trials, while the right panel shows the same after 200,000 trials. Our simulations scale radiolytic hydrogen production using the measured porosity and implemented with (1 − ) described above. In general, this figure shows that the MC method accurately reproduces the calculations of Blair et al. where confidence is high (i.e., standard error is low); any disagreement above ∼55 m is be attributed to (a) the error in measured 232 Th concentrations, and (b) the differing decay energy sums of the two approaches. (Whereas Blair et al. used the endpoint, or maximum, decay energies within a given chain, this study used instead the average to better represent the spectrum of possible outcomes.) Finally, it should be noted that the grain sizes, with an average diameter of 6.2 μm over 120 m, are too small in this case to attenuate the modeled flux. Figure 4 shows modeled H 2 production, normalized to maximum in clay, as a function of depth for the test case of 7.4 μg/g U, 2.8 μg/g Th, and 1.3% K (arbitrary but fairly typical values for geological materials) with 0 = 0.5 and = 5.3E−4 m −1 ; grain size is randomly sampled over intervals corresponding to the grain sizes of clay (0.06-3.9 μm), silt (3.9-63 μm), and sand (63-2,000 μm). Radiolytic hydrogen production in this figure is scaled by the factor (1 − ) described in Section 2.4, which for simplicity is treated as independent of grain size. Figure 4 incorporates the understanding that the porosity of a rock or sediment column-and therefore overall water content-increases nearer the surface. This result shows that the effect of grain size on the H 2 flux is substantial (purely as a result of mineral attenuation of radiation, and ignoring any effect of grain size on porosity).

Effects of Grain Size
Assuming a mass density of 2.6 g/cm 3 results in A vol values of 0.24, 0.030, and 1.1 Bq/cm 3 , respectively, using Equation 16 and the data given in Table 1. Recognizing also that 235 U contributes to the overall activity (0.011 Bq/cm 3 ), the total activity per unit volume of the sediment is simply the sum of the individual activities, or ∼1.4 Bq/cm 3 . From an MC perspective, it is helpful to convert the values above to algorithm-friendly percent fractions of about 18%, 1%, 2%, and 79% for 238 U, 235 U, 232 Th, and 40 K, respectively. Through this example we can see that a random number generator will select potassium predominantly over the course of a given simulation.
Overall, Figure 4 demonstrates the extent to which radiolytic hydrogen production is larger in smaller grains, and so largest in clay when compared with silt or sand. If alpha radiation is used as a baseline with its comparatively large conversion coefficient ( , Equation 17), then this can be understood in terms of an average stopping distance that is at least an entire order of magnitude larger than a characteristic clay grain. By contrast, H 2 flux is reduced in silt-sized grains, where characteristic grain sizes are comparable to the average stopping distance of alpha radiation (∼23 μm). In sand, the production of H 2 is reduced by an entire order of magnitude compared to clay and silt for the porosity-depth relationship considered here (where porosity is itself not dependent on grain size). As before, this phenomenon can be understood in terms of characteristic grain size versus stopping distance and/or attenuation length (in the case of gamma emission). If the average stopping distance of beta radiation from natural radioactive decay is ∼452 μm, then there is the possibility that this decay mode-in conjunction with alpha-is also strongly attenuated, thereby allowing the H 2 contribution due to gamma emission to dominate (attenuation length: ∼3.2 cm or 32,000 μm). It is clear from Figure 4 that the decrease in radiolytic hydrogen production as a function of depth is due in part to a higher order term in the solution to equilibrium compaction, or − 2 0 −2 . Quantitatively, normalized H 2 fluxes near the surface of 100%, 61%, and 10% in clay, silt, and sand, respectively, are reduced to 83%, 51%, and 9% at the maximum depth of 1 km. These values differ according to rock type such that, on average, porosity tends to be greater in coarser rocks; we do not explore this effect here (or the effect of compositional differences between clay-sized and coarser mineral grains typically encountered in siliciclastic rocks). Figure 5 compares the effect of grain size on H 2 production due to alpha, beta, and gamma radiation for the test case presented in Figure 4. (Neither the porosity correction nor the volume correction is applied here since neither would qualitatively alter the comparison.) The top panel shows H 2 production in units of pico (p, E−12) moles per [Earth] year cm 3 as a function of SiO 2 grain radius in cm. Note that H 2 production due to alpha decay is strongly attenuated as a function of grain size. This phenomenon is explained by the relatively short average stopping distance of alphas (∼23 μm or 2.3E−3 cm) compared to beta particles (electrons) and gamma rays (photons). Furthermore, this figure also shows that H 2 production due to alpha decay mirrors the total amount over a logarithmic distribution of grain sizes from clay to sand. This is due to the preponderance of alpha particles in the radioactive emissions of the simulated material, and to the high radiolytic H 2 yield per unit kinetic energy for alpha particles, which dominates that of beta particles and gamma rays by an entire order of magnitude (Spinks & Woods, 1990). Figure 5 (bottom) shows H 2 production in terms of percent fraction as a function of SiO 2 grain radius in cm. This figure (which would have to be adapted for simulations of materials with different U, Th, and K concentrations) is useful in detailing the relative contributions to H 2 production as a result of alpha, beta, and gamma decay in clay, silt, and sand. In general, Figure 5 suggests that the total H 2 production in this test water-SiO 2 matrix is dominated by the influence of alpha decay. Beta and gamma decay, by contrast, make relatively modest contributions, where little or no attenuation is observed. Figure 4. Radiolytic hydrogen production, normalized to maximum in clay, as a function of sediment (or rock) depth for the test case of 7.4 μg/g U, 2.8 μg/g Th, and 1.3% K, with 0 = 0.5 and = 5.3E−4 m −1 . The SiO 2 grain size is sampled over clay (0.06-3.9 μm), silt (3.9-63 μm), and sand (63-2,000 μm In clay, H 2 production from alphas in this simulation is approximately constant at 10.0 ± 0.1 pM/yr cm 3 where the grain size is smaller than the characteristic, or average, stopping distance as a result of natural radioactive decay. By contrast, H 2 production due to beta and gamma decay is comparatively small, with estimates of 0.3377 ± 0.0003 and 0.5220 ± 0.0002 pM/yr cm 3 , respectively. This is explained by the differing radiolytic H 2 yields per unit kinetic energy of these radiation types, that is, , , and factors. The percent contributions in clay due to alpha, beta, and gamma decay are 92%, 3%, and 5%, respectively. The contribution from alpha decay is thus compatible with Blair et al. (2007), who report an estimate of 90 ± 2% in the case of clay-like grain sizes and similar radioactivity levels. These simulations also indicate that the main radioisotope contributor in this particular test material is 238 U, with the remaining components of 235 U, 232 Th, and 40 K producing only trace amounts of H 2 irrespective of the decay mode. This figure shows that these percent contributions are roughly constant due to the comparatively weak attenuation of kinetic energy in the small grains characteristic of clay.
The difference in grain size between clay and silt translates into a marked decrease in H 2 yield by the latter (at a given porosity). This is not surprising given that a typical alpha particle undergoing radioactive decay in SiO 2 will range on the order of E−3 cm. Figure 5 also shows that the degree of attenuation is substantial: The lower grain size bound of silt due to alpha decay produces 9.6 pM/yr cm 3 , while the upper grain size bound produces 3.2 pM/ yr cm 3 . Similar to clay, H 2 production in silt due to beta and gamma radiation are virtually unchanged due to their comparatively large stopping distance and attenuation length, respectively. The effect of increased grain size is more readily apparent in the percent fractions of H 2 production due to alpha, beta, and gamma decay. At the lower grain size bound of silt, alpha radiation is responsible for roughly 92% of H 2 production. This reduces to 79% at the upper bound, with the percent fractions of beta and gamma increasing to 8% and 13%, respectively, as a result.
In sand, the attenuation of alpha particles is more pronounced than that observed in silt. (It should be noted that the reduction of H 2 production in this test water-sand matrix qualitatively resembles the escape efficiency of alpha particles seen in Figure 2.) The lower bound of sand due to alpha decay produces 3.2 pM/yr cm 3 , while the upper bound produces 0.1 pM/yr cm 3 ; between 1 and 2 mm grain size, the production of H 2 due to alpha and Figure 5. Radiolytic hydrogen production via alpha, beta, and gamma decay as a function of SiO 2 grain size for the test case of 7.4 μg/g U, 2.8 μg/g Th, and 1.3% K. Top: total H 2 production per unit volume of water-saturated sediment or rock. Bottom: percentage contribution of each type of radiation to the total H 2 yield.
beta decay is comparable. The effect of alpha particle attenuation is readily apparent in sand, where the percent fraction contributions of beta and gamma radiation compensate proportionally. At the upper bound (∼0.2 cm), alpha, beta, and gamma emissions contribute 14%, 18%, and 68% to H 2 production, respectively. Furthermore, a local maximum in the curve corresponding to beta decay is observed at ∼0.08 cm, suggesting that beta particles "range-out" in grain sizes characteristic of sand. This is consistent with the average beta-particle stopping distance of ∼0.05 cm in natural radioactive decay (National Institute of Standards and Technology, 2021b). No such maximum is observed in the curve corresponding to gamma radiation, however, owing to the comparatively large average attenuation length (∼3.2 cm) of this decay mode in rock/sediment, not to mention the differing physics of photon propagation in matter. Given the grain sizes explored in this study, we have naturally found that all gamma radiation escapes from the originating grain.

Application to Natural Systems
Having validated our method for predicting the radiolytic hydrogen flux in clay, silt, and sand, in this section we consider how grain size could affect radiolytic yields from a variety of geologically interesting natural materials of known (or estimable) porosity and radionuclide concentration. Table 2 summarizes the values of key parameters (other than grain size) used here to simulate porewater radiolysis at Jezero Crater on Mars, the average upper continental crust of the Earth, monazite beach sediment found in Ullal (India), the upper half of the Peru Basin sediment column examined in Figure 3, and Upper Vredefort Crust of South Africa.
The factor tot (decay/s cm 3 ) listed in Table 2 is used in part to convert the H 2 production of Equation 17 (molecules H 2 /decay) to units of M/yr cm 3 . The other activities listed in this table serve as "probability selectors" for the MC simulations. (For example, the radioisotope 40 K, being relatively abundant in minerals, is commonly selected by the random number generator.) In general, these selectors act as proxies for the prevalence of alpha, beta, or gamma emission within the decay chains used by the simulations. The parameters 0 and (surface porosity and compaction coefficient, respectively) are used in Athy's law. For the upper half of the Peru Basin sediment column, we used the average of the porosity measurements above ∼55 m to extrapolate this geological formation to greater depths.
Similar to Figure 4, Figure 6 shows modeled H 2 production, normalized to maximum in clay, as a function of sediment/rock depth for Jezero Crater on Mars, monazite beach sediment found in Ullal (India), and the upper half of the Peru Basin sediment column examined in Figure 3. As before, the grain size is randomly sampled over intervals corresponding to clay (0.06-3.9 μm), silt (3.9-63 μm), and sand (63-2,000 μm). Data corresponding to the upper continental crust of the Earth and Upper Vredefort Crust of South Africa is not shown in this figure due to the null compaction coefficients assumed for this study, which yield constant radiolytic hydrogen fluxes as a function of depth. Relative to maximum production in clay, these fluxes drop to 67% and 68% in silt for the upper continental crust and Upper Vredefort Crust, respectively. In sand, these fluxes are reduced further to 17% and 19%. The total activity per unit volume of rock/sediment tot is the sum of all radioisotope contributions. We assume that hard, crystalline rocks have depth-invariant porosity since they are effectively incompressible over the depth range sampled, i.e., the uppermost one kilometre of crust (Gleeson et al., 2016); for an alternative approach appropriate to larger depth intervals, see Sherwood Lollar et al. (2014). For the other geological settings, we set the compaction coefficient to 5.3E−4 m −1 for Earth and 3.5E−4 m −1 for Mars (Clifford et al., 2010).

Table 2 A Summary of the Input Parameters Used to Generate Figures 6-10
For clay, H 2 production is reduced to 81% relative to maximum in Jezero Crater at a depth of 1 km; in monazite beach sediment-reduced in size to clay-a reduction to 75% is observed. The relatively high surface porosity of the Peru Basin sediment column sees an H 2 production near the surface of 42% relative to maximum, which in this case occurs at approximately 1.06 km. Though not shown in Figure 6, this production will decrease at larger depths in accordance with (1 − ) .
For silt and relative to maximum in clay, the production of molecular hydrogen is 72% and 65% near the surface in Jezero Crater and monazite beach sediment (size-reduced to silt), respectively. The Peru Basin sediment column displays a flux 28% near the surface by virtue of its high porosity. At maximum depth, H 2 production is reduced to 58% and 49% in Jezero Crater and monazite beach sediment, respectively; the Peru Basin sediment column increases to 65% at 1 km.
For sand and again relative to maximum in clay, H 2 production is 27%, 8%, and 6% near the surface in Jezero Crater, monazite beach sediment, and the Peru Basin sediment column, respectively. Whereas these values are reduced at maximum depth to 22% and 6% in Jezero Crater and monazite beach sediment, respectively, by contrast H 2 production increases to 14% in the Peru Basin sediment column. Figure 6. Radiolytic hydrogen production, normalized to maximum in clay, as a function of depth in a variety of geological settings. The SiO 2 grain size is sampled over clay (0.06-3.9 μm), silt (3.9-63 μm), and sand (63-2,000 μm).
The reader should note the difference in H 2 flux between Jezero Crater and monazite beach sediment when uranium enrichment is considered (2% versus 15%, respectively). Since naturally occurring uranium decays by alpha emission, Figure 6 confirms the production of radiolytic hydrogen is comparatively large in clay versus sand. This observation could be relevant to uranium ore environments, in particular those found in the Witwatersrand Basin in South Africa. Figure 7 shows total H 2 production in units of pico (p, E-12) moles per [Earth] year cm 3 as a function of SiO 2 grain radius in cm for the geological formations presented in Figure 6. This figure shows that the monazite beach sediment of Ullal, India produces an H 2 flux that is three orders of magnitude larger than Jezero Crater on Mars. Consistent with the observation in Figure 6, the typical Earth group predicts comparable fluxes of radiolytic hydrogen owing to their similar overall radioactivity. This observation is largely due to their comparable concentrations of thorium and uranium, which decay primarily via alpha emission. If monazite beach sediment were reduced to a grain size in the range of clay, H 2 production from all sources of radiation would be approximately 216 ± 2 pM/yr cm 3 (compared to 86 and 10 pM/yr cm 3 at the upper bounds of silt and sand, respectively). This comparatively large H 2 flux stands in contrast to the typical Earth group, which displays an average of 10.5 ± 1.8 pM/yr cm 3 . At the lower end, Jezero Crater on Mars gives an estimate of 0.57 ± 0.01 pM/yr cm 3 . At the grain size of silt, the typical Earth group gives a reduced H 2 flux of 4.7 ± 0.8 pM/yr cm 3 where grain sizes approach those characterizing sand (i.e., the upper bound of silt). Similarly, Jezero Crater reduces to 0.3 pM/yr cm 3 . This phenomenon is explained by the tendency of alpha emissions to attenuate in silt where the characteristic range of these particles (∼23 μm) is comparable to the grain size. The slope of the decline in H 2 flux with increasing grain size becomes shallower in sand. Specifically, at ∼0.2 cm, radiolytic hydrogen production in the typical Earth group and Jezero Crater is reduced to 1.3 ± 0.3 and 0.1 pM/yr cm 3 , respectively. Finally, since neither the porosity correction nor the volume correction is applied here, the reader should note that the ordering from highest to lowest in Figure 7 naturally differs from that of Figure 6.  When only alpha radiation is considered, our simulations predict that radiolytic hydrogen production would be ∼207 ± 2 pM/yr cm 3 if monazite beach sediment were reduced to a grain size in the range of clay. The typical Earth group displays a middle ground of H 2 flux due to alpha decay, with an average of 8.9 ± 1.7 pM/yr cm 3 in clay. Jezero Crater is expected to produce a flux that is an order of magnitude below the typical Earth group, with an estimate of 0.42 ± 0.01 pM/yr cm 3 . From the standpoint of alpha decay, our simulations show (in descending order) 95.81 ± 0.04% (if reduced to clay), 88.5 ± 0.1%, 85.1 ± 0.2%, and 82.5 ± 0.2% for monazite beach sediment, the Peru Basin sediment column, the average of the upper continental crust of the Earth, and Upper Vredefort Crust, respectively, with an average of 85.4 ± 3.0% for the typical Earth group. This is explained by the 232 Th and 238 U parameters of Table 2, which are derived from the source thorium and uranium measurements of previous studies. (The result for Jezero Crater is slightly lower at 73.9 ± 0.3% for the same reason.) All geological settings studied here attenuate in a predictable manner in silt-sized grains, with estimates of 76.9, 3.2 ± 0.6, and 0.1 pM/yr cm 3 for monazite beach sediment (if appropriately reduced in grain size), the typical Earth group, and Jezero Crater, respectively. Similar to Figure 5 (bottom), in silt the effect of increased grain size is visible on the percent fraction of H 2 production due to alpha decay. The similar shapes of the curves observed in the typical Earth group are due to their analogous radioactive composition; at the upper bound of silt, the percent contribution from alpha decay is reduced to 67.7 ± 5.0%. Bookending the typical Earth group are size-adjusted monazite beach sediment and Jezero Crater, which are reduced to 89.6% and 50.1%, respectively.
In grain sizes characteristic of sand, Figure 8 shows that H 2 production due to alpha decay is reduced substantially: At ∼0.2 cm, this figure gives estimates of 2.59, 0.11 ± 0.02, and 0.005 pM/yr cm 3 for monazite beach sediment, the typical Earth group, and Jezero Crater, respectively. Similar to Figure 5 (bottom), alpha emissions from natural radioactive decay are further attenuated in sand. At the upper bound (∼0.2 cm), alpha emission contributes 26.5%, 8.2 ± 1.8%, and 4.2% to H 2 production in monazite beach sediment, the typical Earth group, and Jezero Crater, respectively. As noted in Figure 5, the upper bound of sand in Figure 8 displays a tapering effect that is consistent with the escape efficiency of alpha particles seen in Figure 2.
Figure 9 (top) shows H 2 production due to beta decay in units of pico (p, E−12) moles per [Earth] year cm 3 as a function of SiO 2 grain radius in cm for a variety of geological formations. In contrast to Figure 8, where the attenuation of alpha particles is observed in silt, in Figure 9 this same phenomenon occurs in grain sizes that char acterize sand. This is explained by the substantially larger stopping distance of beta particles that are routinely emitted via natural radioactive decay (∼0.05 cm). Similar to Figure 8, this figure also shows that-across all grain sizesthe monazite beach sediment of Ullal, India produces the most H 2 overall, a finding that is one order of magnitude above the typical Earth group. Figure 9 (bottom) shows H 2 production in terms of percent fraction as a function of SiO 2 grain radius in cm. Here a global maximum in sand is observed for all geological units.
If monazite beach sediment were reduced to a grain size in the range of clay, our simulations predict that radiolytic hydrogen production would be approximately constant at 3.122 ± 0.002 pM/yr cm 3 . In clay the typical Earth group displays an H 2 flux due to beta decay that is one order of magnitude smaller: 0.5 ± 0.1 pM/yr cm 3 . Jezero Crater on Mars is expected to produce a flux that is an order of magnitude below the typical Earth group, with an estimate of 0.04776 ± 0.00003 pM/yr cm 3 . In descending order, our simulations show percent fraction contributions due to beta decay of 8.37 ± 0.08%, 4.84 ± 0.9%, and 1.45 ± 0.01% for Jezero Crater, the typical Earth group, and monazite beach sediment (in the form of clay), respectively. Similar to Figure 8, these results are explained by the 40 K parameters of Table 2, where potassium decays the most frequently via beta emission. Figure 9 shows that monazite beach sediment, if reduced to grain sizes characteristic of silt, is predicted to produce an H 2 flux due to beta emission that is only slightly diminished compared to that in clay, or 3.0 pM/ yr cm 3 at the upper bound. What is more, the typical Earth group shows almost no variation both qualitatively and quantitatively (at ∼0.006 cm: 0.5 ± 0.1 pM/yr cm 3 ). This is also observed in the case of Jezero Crater, or 0.04685 pM/yr cm 3 at the upper bound of silt. Examining Figure 9 (bottom), if monazite beach sediment were reduced to silt, the percent contribution due to beta decay would increase marginally to 3.5% at the upper bound. Furthermore, Jezero Crater and the typical Earth group also increase to 15.8% and 10.5 ± 1.4%, respectively. Comparing these results to Figure 8 (bottom), we can understand this phenomenon in terms of the tendency of alpha particles to range-out in silt (i.e., beta emission compensates proportionally).
At the upper bound (∼0.2 cm), Figure 9 (top) gives estimates of 1.27, 0.20 ± 0.04, and 0.02 pM/yr cm 3 for monazite beach sediment, the typical Earth group, and Jezero Crater, respectively. Beta emissions are substantially attenuated in sand where the characteristic range of these particles (∼0.05 cm) is comparable to the grain size. That is, beta particles range-out in sand in much the same way as alpha particles in silt. Moreover, it can be seen from Figure 9 that the peak contribution from beta emission occurs at about 0.04, 0.06, and 0.1 cm for Jezero Crater, the typical Earth group, and monazite beach sediment, respectively. Figure 10 shows H 2 production due to gamma decay in units of pico (p, E−12) moles per [Earth] year cm 3 as a function of SiO 2 grain radius in cm for a variety of geological formations. In contrast to alpha and beta emission, no attenuation is observed in radiolytic hydrogen production as a function of grain size. This is explained by (a) the differing physics of photon interaction with matter, and (b) the large attenuation length (∼3.2 cm) characteristic of gamma decay in a natural setting. Similar to Figure 9 (bottom), this figure shows that Jezero Crater produces the most H 2 across all grain sizes due to gamma decay, while monazite beach sediment produces the least, with the typical Earth group presenting a middle ground. These results are sensible given the 40 K parameters of Table 2, which suggest the frequent emission of a residual gamma ray as radioactive potassium decays to 40 Ar.
Gamma rays are photons that lack charge and mass. Consequently they are not subject to electromagnetic slowing as in the case of alpha and beta particles. Thus, the H 2 contribution of gamma rays is constant where the grain size is much less than the characteristic attenuation length found in natural radioactive decay. Our simulations confirm this phenomenon with a summary in Table 3 Figure 10. Radiolytic hydrogen production due to gamma decay as a function of SiO 2 grain size in a variety of geological settings.
In clay (0.06-4 μm), Figure 10 shows percent fraction contributions due to gamma decay of 17.7 ± 0.2%, 9.8 ± 2.2%, and 2.74 ± 0.03% for Jezero Crater, the typical Earth group, and monazite beach sediment (size-reduced), respectively, in descending order. As noted in Figure 9 (bottom), the sequencing observed here is attributed to the relative radioactivity of these geological formations.
In silt (5-60 μm), this figure shows that Jezero Crater, the typical Earth group, and monazite beach sediment (reduced in size to silt) increase to 34.1%, 21.8 ± 3.6%, and 6.9%, respectively. As noted previously, this is a further consequence of the tendency of alpha particles to range-out in silt, where beta emission (Figure 9, bottom) and gamma emission ( Figure 10) compensate proportionally.
Insand (70-2,000 μm), alpha and beta particles are shown to have a diminished effect on the percent contribution to H 2 production. At this point gamma emission dominates radiolytic hydrogen production to the extent that Jezero Crater, the typical Earth group, and monazite beach sediment contribute 80.9%, 76.8 ± 2.0%, and 60.6%, respectively, at the upper bound. Gamma rays characteristic of natural radioactive decay are expected to attenuate in grain sizes of between one and two orders of magnitude beyond those of sand

Discussion and Conclusions
This study has developed and validated a Monte Carlo method that simulates the physics of mineral radiation in order to predict the production of H 2 as a function of mineral grain size and radioisotope composition. The results given here confirm that grain size is a major control on overall radioactivity, and consequently, on radiolytic H 2 yield. The dominant mode of decay (i.e., as controlled by the radioisotope composition) also tends to influence the degree to which H 2 flux is attenuated as a function of grain size; materials whose radioactivity is dominated by their potassium content will show a less severe attenuation of H 2 flux with increasing grain size since 40 K radiation is dominated by beta rather than alpha emissions. When our method is applied to a generic test case, the standard geological classification of grain sizes seems to suggest that clay-sized particles can produce up to an order of magnitude more H 2 per unit time than sand-sized particles (other things being equal). Furthermore, we have used the physics encoded in our simulation method to confirm that alpha particles range-out in silt, while beta particles range-out in grain sizes characteristic of sand. Additionally, all gamma rays, by virtue of their attenuation length, were found to escape to the surrounding water layer in our simulations, with the result that gamma radiation dominates H 2 production for grain sizes larger than sand. Consequently, we do not expect gamma radiation to attenuate until grain sizes exceed those of sand by at least one order of magnitude.
We have applied our numerical method to a variety of interesting natural systems as a means of exploring the magnitude of porewater radiolysis. This exploration was conducted using compositional data from real geological units in order to demonstrate the dependence of natural radionuclide concentration and bulk porosity on radiolytic hydrogen flux. A survey of our simulation results seems to suggest that H 2 production in these geological units is modulated by alpha decay in rock/sediment grain sizes characteristic of clay and silt, while beta and gamma emissions take over as a secondary control in sand-sized grains.
As these results were obtained for homogeneous spherical grains, it is important to note that this work has not addressed the relationships between size and shape obtaining in real geological materials. The importance of these relationships for radiolytic yields is not yet clear; in real rocks and sediments, grains of a size likely to attenuate alpha radiation considerably (e.g., sand) are also likely to be fairly spherical and equant, whereas highly nonspherical grains (e.g., platy clay crystals) are just as likely to be small and thus minimally attenuating. Equally, this work has not considered the broad correlations between grain size and porosity that exist in natural materials; these avenues could be explored in the future. Moreover, future work might also address the effect of  Table 3 Average H 2 Production in Units of pM/yr cm 3 Due to Gamma Decay for a Variety of Geological Settings Predicted Using the Parameters Given in Table 2 surface-area-to-volume ratio rather than grain size on the modeled H 2 production, and the highly heterogeneous distribution of porosity and mineral-water interfaces occurring in some rock types of particular astrobiological interest (e.g., finely fractured crystalline rocks or vesicular basalts; McMahon et al., 2013). It would also be interesting to examine the possible release of radiolytic H 2 from hydrated or hydroxylated mineral grains themselves (which may tend to be more abundant in the clay size-fraction) rather than the porewater. Ultimately, our Monte Carlo method can be readily adapted for materials with constituent grains and pore spaces of arbitrary size and shape, of any composition, and on any planetary body. As a consequence, it should now be possible to estimate more accurately how radiolytic hydrogen might contribute to the habitability of microbial life in porous materials on Earth, Mars, and beyond

Data Availability Statement
The radiolytic hydrogen production data used for the grain size analysis in this study are available via CC BY-SA 4.0 . The Fortran program used for the simulation of porewater radiolysis is available via CC BY-NC-ND 4.0 and developed openly at East Carolina University (DeWitt, 2022).