LithoGas is an R package for Monte Carlo modelling of geologic hydrogen (H₂) and helium (He) production rates via radiolysis and serpentinization. It implements the methods of Warr et al. (2023) and Ardakani et al. (in review), extended to economic source area estimations and deep time production modelling.

Geologic hydrogen and helium are generated in the deep subsurface through two primary mechanisms:

• Radiolysis — the splitting of water molecules into oxygen ions and H2 molecules by alpha, beta, and gamma radiation from the radioactive decay of U, Th, and K in surrounding minerals
• Serpentinization — oxidation of Fe²⁺ to Fe³⁺ during hydration of iron-rich rocks, releasing H₂

LithoGas propagates measurement uncertainties through both models using Monte Carlo simulation, producing distributions of production rates rather than single-point estimates. Results can be scaled to arbitrary source rock volumes for economic feasibility assessment.

LithoGas is not yet on CRAN. Install directly from GitHub:

# install.packages("devtools")
devtools::install_github("yourusername/LithoGas")

library(LithoGas)

# Load the example dataset (example structured dataframe)
data("structuredDF")

# Run Monte Carlo simulation with radiolysis (50 trials per sample)
results <- monteProd(structuredDF, numGen = 50, rad = TRUE)

# Summarise results by sample
sumDF <- monteSum(results, summaryField = "Sample")

# Plot H2 production scaled by source rock volume
plots <- monteH2Plot(sumDF)
plots[[2]]  # view the ggplot2 object

Input dataframe      ← structured dataframe - see examples
      │
      ▼
 monteProd()      ← main wrapper; iterates over samples
      │
      ├──► joinLitProps()      ← function within monteProd() that assigns rock density & porosity from CRPPData
      │                          (when litLith provided, no sample data)
      │
      ├──► monteRad()      ← function within monteProd() that caluclates H2 & He production rates via radiolysis 
      │                      (when rad = TRUE)
      │
      └──► monteSerpFeSpecies()      ← function within monteProd() that calculates H2 production rates via serpentinization from geochemistry data with iron speciation
      │                                (when serp = TRUE)
      │
            └──► monteSerpFeSpecies()      ← function wihtin monteProd() that calculates H2 production rates via serpentinizationfrom geochemistry data with total iron concentration
      │                                      (when serp = TRUE and allowTotalFeSerp = TRUE)
      ▼
 monteSum() ← summarises trials to min/mean/max per sample
      │
            └──► monteH2Plot()        ← log-log source volume plot for H2 from output of monteSum()
            └──► monteHePlot()        ← log-log source volume plot for He from output of monteSum()
      │
      ▼
deepTimeProd() ← Back H2 and He generation rates following average serpentinziation rates, and radioactive decay law

LithoGas accepts a structured dataframe where each row is one sample. All continuous parameters are defined as truncated normal distributions using four columns per variable: Min, Max, Mean, SD. One sample can represent a single geochemical analysis where standard deviation represents the analytical uncertainty, or a one sample can represent multiple geochemical analyses where the range represents the group parameters.

Two approaches are supported for rock properties:

Option 1 — Sample-specific rock properties (use monteDataRockProps as a template)

Provide rockDenMin/Max/Mean/SD and porMin/Max/Mean/SD directly in your input dataframe.

Option 2 — Lithology-based rock properties (use monteDataLithCat as a template)

Provide a litLith column matching a lithology in CRPPData. Rock density and porosity distributions are looked up automatically.

# View available lithologies
data("CRPPData")
unique(CRPPData$Lithology)

Note: monteProd() evaluates each input row/sample for the presence of all required variables prior to performing calculations. It will not run functions where data is not present. As such, multiple methods (radiolysis only, radiolysis + iron speciation serpentinization) can be run from the same input datasheet, as per the example.

monteProd() returns a long-format dataframe with numGen rows per sample, containing all input columns plus computed outputs. Key output columns:

monteSum() returns a wide-format dataframe of H₂/He generation rates summarized by selected field (i.e., by a group_by()). All rates are in mol gas / m³ rock / year Key output columns:

^ RadHeRateMin_molm3yr is equal to HeRateMin_molm3yr. HeRateMin_molm3yr, was included to keep same formatting as H2RateMean_molm3yr

monteH2Plot() and monteHePlot() return a list: [[1]] a long-format dataframe used for plotting, where each summarized entry from monteSum() has multiple rows, and [[2]] being a ggplot. The dataframe includes yearly generation rates each sample scaled from 0.1 km³ to 100 km³. Examples below are for H₂ via monteH2Plot(), but the same columns are generated for He monteHePlot()

[[1]] Dataframe output for plotting

[[2]] ggplot - an example ggplot of log-log scaled sourve volume plot, note that the function can be found using view (e.g., View(monteH2Plot()) so ggplot aesthetics can be manipulated.

If you use LithoGas in your research, please cite:

Coutts, D.S., Warr, O.A., Ardakani, O.A., Sherwood Lollar, B. (in review)

The serpentininzation model is based on:

Ardakani, O.A., Sherwood Lollar, B., Coutts, D.S., Warr, O.A., Delonde, C., Kabanov, P., Lister, C. (in review)

The radiolysis model is based on:

Warr, O., Song, M., Sherwood Lollar, B. (2023). The application of Monte Carlo modelling to quantify in situ hydrogen and associated element production in the deep subsurface. Frontiers in Earth Science, v.11.

Rock property data sourced from:

Enkin, R.J. (2018). The Canadian Rock Physical Property Database: first public release. Geological Survey of Canada, Open File 8460, 68 p. Natural Resources Canada. https://ostrnrcan-dostrncan.canada.ca/entities/publication/c4c0cede-365c-4c87-8077-8e045e874de6

• Daniel Coutts — Geological Survey of Canada, Natural Resources Canada ([email protected])
• Oliver Warr — University of Toronto
• Omid Ardakani — Geological Survey of Canada, Natural Resources Canada
• Barbara Sherwood Lollar — University of Toronto

This package is licensed under CC BY-NC 4.0 — free for non-commercial use with attribution.