Oman methane isotope plotting
Source file: 200813_Oman_14CH4_paper_figs.Rmd
Daniel Nothaft
08 Jun 2025
1 Setup
1.1 Load code libraries
# load libraries
# install if needed
# install.packages("tidyverse")
library(tidyverse) # dplyr, tidyr, ggplot
# install.packages("readxl")
library(readxl) # reading excel files
# install.packages("knitr")
library(knitr) # generating reports
# install.packages("scales")
library(scales) # making log scales with exponents
# install.packages("ggrepel")
library(ggrepel) # algorithm to avoid overlapping labels
# install.packages("cowplot")
library(cowplot) # extracting legends from plots
# install.packages("lemon")
library(lemon) # repeat axis lines for paneled figures in ggplot
# install.packages("ggpmisc")
library(ggpmisc) # adding linear regression equations to ggplot
# install.packages("latex2exp")
library(latex2exp) # plot labels with latex
# global knitting options for automatic saving of all plots as .png and .pdf
knitr::opts_chunk$set(
dev = c("png", "pdf"), fig.keep = "all",
dev.args = list(pdf = list(encoding = "WinAnsi", useDingbats = FALSE)),
fig.path = file.path("fig_output/", paste0(gsub("\\.[Rr]md", "/", knitr::current_input()))),
cache.path = file.path("cache/", paste0(gsub("\\.[Rr]md", "/", knitr::current_input())))
)# source all relevant scripting files
source(file.path("scripts", "plotting_functions.R")) # includes LaTeX labelling function1.2 Load data
Load in Oman well data
# Load in Oman aqueous geochemistry etc.
data <- read_xlsx("data_raw/Summary_metadata/Oman_Geochem_2012-2019_14CH4.xlsx", na = "NA") # read in data from excel file
# add site column for plotting against other literature data
data <- data %>% mutate(site = "Oman/UAE")
# print data
dataLoad in comparison data from the literature
# read comparison data from literature
data_literature <- read_excel("data_raw/Summary_metadata/CH4_isotope_literature_data.xlsx", na = "NA")
# Filter for only samples mentioning chromite or chromitite for Greek rock crushings
data_literature <- data_literature %>% filter(if_else(site == "Greece", str_detect(Lithology, "hromit") == TRUE, !is.na(site)))
# print
data_literatureLoad in Henry’s constants from Sander, 2015
# read Henry's constants from literature
Henrys <- read_csv("data_raw/Henrys_constants_Sander.csv")## Rows: 8 Columns: 3
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): compound
## dbl (2): Hcp_mol_per_m3_Pa_25C, dlnHcp_over_d_inv_T
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.1.3 Set plotting aesthetics
# colors
microbial_color <- "#009E73" # forest/mint green
thermogenic_color <- "#F0E442" # yellow
abiotic_color <- "#999999" # grey
microbial_opacity <- 0.05
thermo_opacity <- 0.1
abiotic_opacity <- 0.08
gabbro_color <- "#CC79A7" # magenta
Mg_HCO3_color <- "#E69F00" # gold
Ca_OH_color <- "#0072B2" # dark blue
cbPalette_water_types <- c(gabbro_color, Mg_HCO3_color, Ca_OH_color)
names(cbPalette_water_types) <- c("gabbro", "Mg_HCO3", "Ca_OH")
cbPalette_MWL <- c("#D55E00", "#0072B2", "#000000", "#999999")
color_rock_crush_named <- c("rock crushing" = "#D55E00", "vent/well fluid" = "black")
color_rock_crush_named_new_line <- c("rock\ncrushing" = "#D55E00", "vent/well\nfluid" = "black")
fill_4_wells_named <- c("WAB71" = "#999999", "NSHQ14" = "#000000", "NSHQ04" = "#56B4E9", "CM2A" = "#0072B2")
fill_5_wells_named <- c("WAB71" = "#999999", "NSHQ14" = "#000000", "NSHQ04" = "#56B4E9", "CM2A" = "#0072B2", "WAB56" = NA)
# levels (for order of appearance in legend)
data_literature$site <- factor(data_literature$site, levels = c("Oman/UAE", "Philippines", "Mid-Cayman Rise", "Ashadze II","Greece", "Kidd Creek"))
# shapes
shapes_all_sites_named <- c("Oman/UAE" = 21, "Philippines" = 0, "Mid-Cayman Rise" = 2, "Ashadze II" = 5, "Kidd Creek" = 4, "Greece" = 3)
shapes_fills <- c(22, 23, 24, 21)2 C\(_1\) / (C\(_2\) + C\(_3\)) vs. \(\delta^{13}\text{C}_{\text{CH}_4}\) (“Bernard”) Plot
2.1 Prepare data
Prepare Oman well data
# select samples with quantifiable C2H6
data_C1_C2_C3 <- data %>% filter(!is.na(C2H6_uM)) %>% select(sampling_site, notes, year_sampled, depth_fluid_intake_mbct, water_type, CH4_uM, CH4_was_measured, CH4_LQ_uM, CH4_unc_uM, C2H6_uM, C2H6_was_measured, C2H6_LQ_uM, C2H6_unc_uM, C3H8_uM, C3H8_was_measured, C3H8_LQ_uM, C3H8_unc_uM,
CH4_d13C_permil_VPDB_CUB, CH4_d13C_unc_permil_CUB) # select relevant data for CH4, C2H6, C3H8
data_C1_C2_C3 <- data_C1_C2_C3 %>% filter(CH4_uM > 0.5) # filter for samples for with greater than 0.5 µmol/L CH4
# calculate C1/(C2+C3)
data_C1_C2_C3 <- data_C1_C2_C3 %>% mutate(C1_over_C2_plus_C3 = CH4_uM / (C2H6_uM + if_else(!is.na(C3H8_uM), C3H8_uM, 0)))
# # summarize wells with multiple replicates in a given year (NSHQ14 2017)
# data_C1_C2_C3 <- data_C1_C2_C3 %>% group_by(sampling_site, year_sampled) %>% summarise(water_type = first(water_type), C1_over_C2_plus_C3 = mean(C1_over_C2_plus_C3), CH4_d13C_permil_VPDB_CUB = mean(CH4_d13C_permil_VPDB_CUB, na.rm = TRUE), CH4_uM = mean(CH4_uM), C2H6_uM = mean(C2H6_uM), C3H8_uM = mean(C3H8_uM)) %>% ungroup()
data_C1_C2_C3 <- data_C1_C2_C3 %>% mutate(site = "Oman/UAE")
# print
data_C1_C2_C3Get literature data
# select samples with quantifiable C2H6
data_C1_C2_C3_lit <- data_literature %>% filter(!is.na(C2H6_conc)) %>% select(site, sample_type, CH4_conc, C2H6_conc, C3H8_conc,
CH4_d13C_permil_VPDB) # select relevant data for CH4, C2H6, C3H8
# calculate C1/(C2+C3)
data_C1_C2_C3_lit <- data_C1_C2_C3_lit %>% mutate(C1_over_C2_plus_C3 = CH4_conc / (C2H6_conc + if_else(!is.na(C3H8_conc), C3H8_conc, 0)))
# print
data_C1_C2_C3_lit2.2 Plot
# Set coordinates of fields of origins (Milkov and Etiope 2018)
primary_microbial_bernard <- data.frame(
x = c(-90,-70,-50, -50, -90),
y = c(2e2,2e2,7e2, 1e5, 1e5)
)
secondary_microbial_bernard <- data.frame(
x = c(-55, -50, -45, -35, -40, -55),
y = c(4e3, 4e0, 2e0, 7e0, 1e5, 1e5)
)
thermogenic_bernard <- data.frame(
x = c(-50, -65, -70, -50, -47, -42, -40, -40, -15, -15, -30, -50),
y = c(1e-1, 1e0, 2e2, 7e2, 2e2, 1e2, 2e2, 1e5, 1e5, 5e2, 1e-1, 1e-1)
)
abiotic_bernard <- data.frame(
x = c(-50, -30, -10, 10, -30),
y = c(8e1, 1e-1, 1e-1, 1e5, 1e5)
)
# create plot
C1_C2_C3_plot <-
ggplot() +
# add origin fields
geom_polygon(data = primary_microbial_bernard, aes(x=x, y=y), fill = microbial_color, alpha = microbial_opacity, color = microbial_color, size = 0.2) +
geom_polygon(data = secondary_microbial_bernard, aes(x=x, y=y), fill = NA, color = "black", size = 0.2, linetype = "dashed") +
geom_polygon(data = thermogenic_bernard, aes(x=x, y=y), alpha = thermo_opacity, fill = thermogenic_color, color = thermogenic_color, size = 0.5) +
geom_polygon(data = abiotic_bernard, aes(x=x, y=y), alpha= abiotic_opacity, fill = abiotic_color, color = abiotic_color, size = .2) +
# add text annotations
annotate("text", y = 3e4, x = -80, size = 3, label = "PM", color = microbial_color) +
annotate("text", y = 1e4, x = -45, size = 3, label = "SM", color = "grey20") +
annotate("text", y = 1e0, x = -56, size = 3, label = "T", color = "grey20") +
annotate("text", y = 3e0, x = -15, size = 3, label = "A", color = "grey20") +
# plot literature data
geom_point(data = data_C1_C2_C3_lit,
aes(
x = CH4_d13C_permil_VPDB,
y = C1_over_C2_plus_C3,
shape = site,
color = sample_type,
), alpha = 1) +
# plot Oman well data
geom_point(data = data_C1_C2_C3,
aes(
x = CH4_d13C_permil_VPDB_CUB ,
y = C1_over_C2_plus_C3,
shape = site,
fill = sampling_site
), alpha = 1) +
# axis styling
scale_y_continuous(trans = 'log10',
breaks = trans_breaks('log10', function(x) 10 ^ x, n=5),
labels = trans_format('log10', math_format(10^.x)),
limits = c(.1, 1e5), expand = c(0,0),
name = latex2exp::TeX("$C_1$/($C_{2}$+$C_{3}$)"))+
scale_x_continuous(name = latex2exp::TeX("$\\delta ^{13}C_{CH_4}$$\\,$/$\\,$$\\lbrack$‰ VPDB$\\rbrack$"), limits = c(-90, 15), expand = c(0,0)) +
# symbol styling
scale_shape_manual(values = shapes_all_sites_named) +
scale_color_manual(values = color_rock_crush_named) +
scale_fill_manual(values = fill_4_wells_named, name = "water type", guide = guide_legend(order = 2)) +
# design
theme_classic(base_size = 12)+
theme(
# axis.title.x = element_blank(),
# axis.text.x.bottom = element_blank(),
# axis.ticks.x.bottom = element_blank(),
# axis.line.x.bottom = element_blank(),
plot.margin = margin(5, 1, 5, 1, "pt"),
legend.position = "none"
)
# print plot
C1_C2_C3_plot
3 Alkanes \(\delta^{13}\)C plot
3.1 Prepare data
# create data frame of CH4, C2H6, and C3H8 d13C data
data_alkanes_d13C <- data %>% filter(!is.na(C2H6_d13C_permil_VPDB)) %>% select(sampling_site, year_sampled, depth_fluid_intake_mbct, CH4_d13C_permil_VPDB_CUB, CH4_d13C_unc_permil_CUB, CH4_d13C_permil_VPDB_MIT, CH4_d13C_unc_permil_MIT, CH4_d13C_permil_VPDB_ISOTECH, CH4_d13C_unc_permil_ISOTECH, C2H6_d13C_permil_VPDB, C2H6_d13C_unc_permil, C3H8_d13C_permil_VPDB, C3H8_d13C_unc_permil)
# print
data_alkanes_d13CGather alkane \(\delta^{13}\text{C}\) data so it can be plotted
# gather data making a column for the compound
data_alkanes_d13C_gathered <- data_alkanes_d13C %>% select(-depth_fluid_intake_mbct) %>% gather("compound_lab", "d13C_permil_VPDB", -sampling_site, -year_sampled, -CH4_d13C_unc_permil_CUB, -CH4_d13C_unc_permil_MIT, -CH4_d13C_unc_permil_ISOTECH, -C2H6_d13C_unc_permil, -C3H8_d13C_unc_permil)
# make a column that gives the # of C atoms in the compound based on the compounds chemcal formula
data_alkanes_d13C_gathered <- data_alkanes_d13C_gathered %>% mutate(C_atoms = case_when(
str_detect(compound_lab, "CH") == TRUE ~ 1,
str_detect(compound_lab, "C2") == TRUE ~ 2,
str_detect(compound_lab, "C3") == TRUE ~ 3
))
# make another column for the uncertainty of each d13C analysis
data_alkanes_d13C_gathered <- data_alkanes_d13C_gathered %>% mutate(uncertainty_d13C_permil = case_when(
compound_lab == "CH4_d13C_permil_VPDB_CUB" ~ CH4_d13C_unc_permil_CUB,
compound_lab == "CH4_d13C_permil_VPDB_MIT" ~ CH4_d13C_unc_permil_MIT,
compound_lab == "CH4_d13C_permil_VPDB_ISOTECH" ~ CH4_d13C_unc_permil_ISOTECH,
compound_lab == "C2H6_d13C_permil_VPDB" ~ C2H6_d13C_unc_permil,
compound_lab == "C3H8_d13C_permil_VPDB" ~ C3H8_d13C_unc_permil
))
# select only the most relevant columns
data_alkanes_d13C_gathered <- data_alkanes_d13C_gathered %>% select(sampling_site, year_sampled, C_atoms, compound_lab, d13C_permil_VPDB, uncertainty_d13C_permil) %>% filter(!is.na(d13C_permil_VPDB))
# print
data_alkanes_d13C_gatheredAdd a site and sample type to Oman well data so that it is compatible for plotting with literature data
data_alkanes_d13C_gathered <- data_alkanes_d13C_gathered %>% mutate(site = "Oman/UAE", sample_type = "vent/well\nfluid")
data_alkanes_d13C_gatheredGet literature data
# get relevant data for C1-C5 alkane analyses
data_alkanes_d13C_lit <- data_literature %>% filter(!is.na(C2H6_d13C_permil_VPDB)) %>% filter(!is.na(C2H6_d13C_permil_VPDB)) %>% select(site, Sample, sample_type, CH4_d13C_permil_VPDB, C2H6_d13C_permil_VPDB, C3H8_d13C_permil_VPDB, n_C4H10_d13C_permil_VPDB, n_C5H12_d13C_permil_VPDB)
# print
data_alkanes_d13C_litAs previously done for Oman well data, gather literature data by compound for plotting purposes
# select most relevant data from literature
data_alkanes_d13C_lit <- data_alkanes_d13C_lit %>% select(site, Sample, sample_type, CH4_d13C_permil_VPDB, C2H6_d13C_permil_VPDB, C3H8_d13C_permil_VPDB, n_C4H10_d13C_permil_VPDB, n_C5H12_d13C_permil_VPDB)
# gather by compound
data_alkanes_d13C_gathered_lit <- data_alkanes_d13C_lit %>% gather("compound", "d13C_permil_VPDB", -site, -Sample, -sample_type)
# assign number of C atoms per molecule based on the compound formula
data_alkanes_d13C_gathered_lit <- data_alkanes_d13C_gathered_lit %>% mutate(C_atoms = case_when(
str_detect(compound, "CH4") == TRUE ~ 1,
str_detect(compound, "C2H6") == TRUE ~ 2,
str_detect(compound, "C3H8") == TRUE ~ 3,
str_detect(compound, "C4H10") == TRUE ~ 4,
str_detect(compound, "C5H12") == TRUE ~ 5
))
# change sample type to have new line character for plotting
data_alkanes_d13C_gathered_lit <- data_alkanes_d13C_gathered_lit %>% mutate(sample_type = if_else(sample_type == "vent/well fluid", "vent/well\nfluid", "rock\ncrushing"))
# print
data_alkanes_d13C_gathered_lit3.2 Plot
3.2.1 Plot data
# create plot
alkanes_d13C_oman_wells_and_lit_plot <- ggplot() +
# Add NSHQ14 label
annotate("text", label = "2017", x = 1.15, y = 1.2, label.size = 0.22) +
annotate("text", label = "2019", x = 1.15, y = 8.2, label.size = 0.22) +
# Plot literature data
# lines
geom_line(data = data_alkanes_d13C_gathered_lit %>% filter(C_atoms < 4),
mapping = aes(
x = C_atoms,
y = d13C_permil_VPDB,
color = sample_type,
group = Sample), alpha = 0.1) +
# points
geom_point(data = data_alkanes_d13C_gathered_lit %>% filter(C_atoms < 4),
mapping = aes(
x = C_atoms,
y = d13C_permil_VPDB,
shape = site,
color = sample_type), alpha = 0.7) +
# Plot Oman well data
# lines
geom_line(data = data_alkanes_d13C_gathered %>% filter(C_atoms < 4),
mapping = aes(
x = C_atoms,
y = d13C_permil_VPDB,
group = year_sampled,
color = sample_type), alpha = 1, size = .6) +
# error bars
geom_linerange(data = data_alkanes_d13C_gathered,
mapping = aes(
x = C_atoms,
ymax = d13C_permil_VPDB + uncertainty_d13C_permil,
ymin = d13C_permil_VPDB - uncertainty_d13C_permil,
color = sample_type), show.legend = FALSE) +
# points
geom_point(data = data_alkanes_d13C_gathered,
mapping = aes(
x = C_atoms,
y = d13C_permil_VPDB,
shape = site,
color = sample_type,
fill = sampling_site)) +
# symbol styling
scale_shape_manual(name = "sample site", values = shapes_all_sites_named, guide = guide_legend(order = 1)) +
scale_fill_manual(values = fill_4_wells_named, guide = FALSE) +
scale_color_manual(values = color_rock_crush_named_new_line, name = "sample type", guide = guide_legend(nrow = 1, order = 3)) +
guides(color = guide_legend(override.aes=list(shape=21))) +
# axis styling
scale_x_continuous(name = "C atoms in alkane", breaks = c(1,2,3), expand = expand_scale(add = c(0.04)))+
scale_y_continuous(name = latex2exp::TeX("$\\delta ^{13}C$$\\,$/$\\,$\\lbrack$$‰ VPDB$\\rbrack$"))+
# design
theme_classic(base_size = 12) +
theme(
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank(),
legend.direction = "vertical",
legend.position = "bottom",
legend.text = element_text(margin = margin(t=3, b=3)),
plot.margin = margin(1, 3, 0, 1, "pt")
)
# print plot
alkanes_d13C_oman_wells_and_lit_plot
# create plot
alkanes_d13C_oman_only_plot <- ggplot() +
# Add NSHQ14 label
annotate("text", label = "NSHQ14\n2017", x = 1.2, y = -4, label.size = 0.22, hjust = 0.5) +
annotate("text", label = "NSHQ14\n2019", x = 1.2, y = 7.5, label.size = 0.22, hjust = 0.5) +
# Add Nizwa label
annotate("text", label = "Nizwa", x = 1.2, y = -11, label.size = 0.22, hjust = 0.5, color = "#999999") +
# Plot literature data
# lines
geom_line(data = data_alkanes_d13C_gathered_lit %>% filter(C_atoms < 4) %>% filter(site == "Oman/UAE"),
mapping = aes(
x = C_atoms,
y = d13C_permil_VPDB,
color = sample_type,
group = Sample), alpha = 0.1) +
# points
geom_point(data = data_alkanes_d13C_gathered_lit %>% filter(C_atoms < 4) %>% filter(site == "Oman/UAE"),
mapping = aes(
x = C_atoms,
y = d13C_permil_VPDB,
shape = site,
color = sample_type), alpha = 0.7) +
# Plot Oman well data
# lines
geom_line(data = data_alkanes_d13C_gathered %>% filter(C_atoms < 4),
mapping = aes(
x = C_atoms,
y = d13C_permil_VPDB,
group = year_sampled,
color = sample_type), alpha = 1, size = .6) +
# error bars
geom_linerange(data = data_alkanes_d13C_gathered,
mapping = aes(
x = C_atoms,
ymax = d13C_permil_VPDB + uncertainty_d13C_permil,
ymin = d13C_permil_VPDB - uncertainty_d13C_permil,
color = sample_type), show.legend = FALSE) +
# points
geom_point(data = data_alkanes_d13C_gathered,
mapping = aes(
x = C_atoms,
y = d13C_permil_VPDB,
shape = site,
color = sample_type,
fill = sampling_site)) +
# symbol styling
scale_shape_manual(name = "sample site", values = shapes_all_sites_named, guide = guide_legend(order = 1)) +
scale_fill_manual(values = fill_4_wells_named, guide = FALSE) +
scale_color_manual(values = color_rock_crush_named_new_line, name = "sample type", guide = guide_legend(nrow = 1, order = 3)) +
guides(color = guide_legend(override.aes=list(shape=21))) +
# axis styling
scale_x_continuous(name = "C atoms in alkane", breaks = c(1,2,3), expand = expand_scale(add = c(0.04)))+
scale_y_continuous(name = latex2exp::TeX("$\\delta ^{13}C$$\\,$/$\\,$\\lbrack$$‰ VPDB$\\rbrack$"))+
# design
theme_classic(base_size = 12) +
theme(
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank(),
legend.direction = "vertical",
legend.position = "none",
legend.text = element_text(margin = margin(t=3, b=3)),
plot.margin = margin(1, 3, 0, 1, "pt")
)
# print plot
alkanes_d13C_oman_only_plot
4 CH\(_4\) Fraction modern facet grid plot
4.1 Prepare data
# create data frame of CH4, C2H6, and C3H8 d13C data
data_14CH4 <- data %>% filter(!is.na(CH4_Fm_UCIAMS)) %>% select(sampling_site, year_sampled, notes, CH4_Fm_UCIAMS, CH4_Fm_unc_UCIAMS, CH4_Fm_NOSAMS, CH4_Fm_unc_NOSAMS, CH4_d13C_permil_VPDB_CUB, CH4_d13C_unc_permil_CUB, CH4_d13C_permil_VPDB_ISOTECH, CH4_d13C_unc_permil_ISOTECH, CH4_dD_permil_VSMOW_CUB, CH4_dD_unc_permil_CUB)
# print
data_14CH4(data_14CH4_NSHQ14 <- data_14CH4 %>% filter(sampling_site == "NSHQ14"))data_14CH4_NSHQ14_longer <- data_14CH4_NSHQ14 %>%
# combine F14C data from different labs: UCIAMS and NOSAMS
pivot_longer(c(CH4_Fm_UCIAMS, CH4_Fm_NOSAMS), names_to = "parameter", values_to = "CH4_Fm", values_drop_na = TRUE) %>%
# rename "parameter" column (e.g. "CH4_Fm_UCIAMS") to the laboratory name
mutate(lab = case_when(
str_detect(parameter, "UCIAMS") == TRUE ~ "CUB/UCI",
str_detect(parameter, "NOSAMS") == TRUE ~ "Isotech/\nNOSAMS"
)) %>%
# drop "parameter column
select(-parameter) %>%
# assign Fm uncertainties based on lab
mutate(CH4_Fm_unc = if_else(lab == "Isotech/\nNOSAMS", CH4_Fm_unc_NOSAMS, CH4_Fm_unc_UCIAMS)) %>%
select(-CH4_Fm_unc_NOSAMS, -CH4_Fm_unc_UCIAMS) %>%
# relocate lab column after CH4_Fm_unc_UCIAMS
relocate(lab, .before = CH4_Fm) %>%
# assign CH4 d13C based on lab
mutate(CH4_d13C_permil_VPDB = if_else(lab == "Isotech/\nNOSAMS", CH4_d13C_permil_VPDB_ISOTECH, CH4_d13C_permil_VPDB_CUB)) %>%
select(-CH4_d13C_permil_VPDB_ISOTECH, -CH4_d13C_permil_VPDB_CUB) %>%
# assign CH4 d13C uncertainty based on lab
mutate(CH4_d13C_unc_permil = if_else(lab == "Isotech/\nNOSAMS", CH4_d13C_unc_permil_ISOTECH, CH4_d13C_unc_permil_CUB)) %>%
select(-CH4_d13C_unc_permil_ISOTECH, -CH4_d13C_unc_permil_CUB) %>%
# assign CH4 dD based on lab (CH4 dD was not measured by Isotech)
mutate(CH4_dD_permil_VSMOW = if_else(lab == "Isotech/\nNOSAMS", NA_real_, CH4_dD_permil_VSMOW_CUB)) %>%
select(-CH4_dD_permil_VSMOW_CUB) %>%
# assign CH4 dD uncertainty based on lab
mutate(CH4_dD_unc_permil = if_else(lab == "Isotech/\nNOSAMS", NA_real_, CH4_dD_unc_permil_CUB)) %>%
select(-CH4_dD_unc_permil_CUB)
data_14CH4_NSHQ14_longerdata_14CH4_NSHQ14_longer_longer <- data_14CH4_NSHQ14_longer %>%
pivot_longer(c(CH4_d13C_permil_VPDB, CH4_dD_permil_VSMOW), names_to = "parameter", values_to = "value") %>%
# assign uncertainty based on parameter
mutate(uncertainty_permil = if_else(parameter == "CH4_d13C_permil_VPDB", CH4_d13C_unc_permil, CH4_dD_unc_permil)) %>%
select(-CH4_d13C_unc_permil, -CH4_dD_unc_permil)
data_14CH4_NSHQ14_longer_longer# add tex formatting
data_14CH4_NSHQ14_longer_longer <- data_14CH4_NSHQ14_longer_longer %>%
mutate(
tex_params = as_factor(parameter)%>%
fct_recode(
"$\\delta ^{13}C_{CH_4}\\,/\\,\\[$‰$ \\, VPDB \\]$" = "CH4_d13C_permil_VPDB",
"$\\delta D_{CH_4}\\,/\\,\\[$‰$ \\, VSMOW \\]$" = "CH4_dD_permil_VSMOW"
)) %>%
# make shorter note for NSHQ14
mutate(notes_short = if_else(notes == "no purge (NP)", "NP", NA_character_))4.2 Plot
CH4_F14C_d13C_dD_plot <- data_14CH4_NSHQ14_longer_longer %>% ggplot(aes(
x = CH4_Fm,
y = value,
# x error bars are smaller than data points
# xmin = CH4_Fm - CH4_Fm_unc,
# xmax = CH4_Fm + CH4_Fm_unc,
ymin = value - uncertainty_permil,
ymax = value + uncertainty_permil,
shape = as.character(year_sampled),
color = lab,
label = notes_short
)) +
# error bars
geom_linerange(orientation = "x") +
# geom_linerange(orientation = "y") +
# points
geom_point() +
# text
geom_text(nudge_x = -0.026, size = 2.75, show.legend = FALSE) +
# grid
facet_grid(rows = vars(tex_params), scales = "free_y", switch = "y", labeller = latex_labeller) +
# scale adjustments
scale_x_continuous(name = latex2exp::TeX("CH$_4$ fraction modern"), limits = c(0, 0.32), expand = c(0,0)) +
scale_y_continuous(name = NULL) +
scale_shape_discrete(name = "Year sampled") +
scale_color_manual(name = "Laboratory", values = c("CUB/UCI" = "black", "Isotech/\nNOSAMS" = "#0072B2")) +
# styling
theme_bw() +
theme(
panel.grid = element_blank(),
legend.position = "right",
legend.direction = "vertical",
legend.box.margin = margin(r=2),
plot.margin = margin(l=2, t=2, b=1, r=2),
legend.box.spacing = unit(0.1, "cm")
)
CH4_F14C_d13C_dD_plot
CH4_F14C_d13C_dD_plot_year_color_shape <- data_14CH4_NSHQ14_longer_longer %>% ggplot(aes(
x = CH4_Fm,
y = value,
# xmin = CH4_Fm - CH4_Fm_unc,
# xmax = CH4_Fm + CH4_Fm_unc,
ymin = value - uncertainty_permil,
ymax = value + uncertainty_permil,
shape = as.character(year_sampled),
color = as.character(year_sampled),
label = notes_short
)) +
# error bars
geom_linerange(orientation = "x") +
# geom_linerange(orientation = "y") +
# points
geom_point() +
# text
geom_text(nudge_x = -0.015, size = 3, show.legend = FALSE) +
# grid
facet_grid(rows = vars(tex_params), scales = "free_y", switch = "y", labeller = latex_labeller) +
# scale adjustments
scale_x_continuous(name = latex2exp::TeX("CH$_4$ fraction modern"), limits = c(0, 0.32), expand = c(0,0)) +
scale_y_continuous(name = NULL) +
scale_shape_discrete(name = "Year sampled") +
scale_color_manual(name = "Year sampled", values = cbPalette_MWL) +
# styling
theme_bw() +
theme(
panel.grid = element_blank(),
legend.position = "bottom",
legend.box.margin = margin(r=2),
plot.margin = margin(l=2, t=2, b=1, r=2),
legend.box.spacing = unit(0.1, "cm"),
legend.text = element_text(size = 8.2),
legend.title = element_text(size = 10)
)
CH4_F14C_d13C_dD_plot_year_color_shape
CH4_F14C_d13C_dD_plot_year_color_shape_clean <- CH4_F14C_d13C_dD_plot_year_color_shape +
facet_rep_grid(rows = vars(tex_params), scales = "free_y", switch = "y", labeller = latex_labeller) +
theme_classic() +
theme(
strip.placement = "outside",
strip.background = element_blank(),
strip.text = element_text(size = 11, vjust = 0),
legend.position = "bottom",
legend.box.margin = margin(r=14),
plot.margin = margin(l=0, t=0, b=1, r=0),
legend.box.spacing = unit(0.1, "cm"),
legend.text = element_text(size = 8.2),
legend.title = element_text(size = 10)
)
CH4_F14C_d13C_dD_plot_year_color_shape_clean
CH4_F14C_d13C_dD_plot_year_color_shape_clean %>% tag_facet_keep_strip(
open = "",
close = "",
tag_pool = LETTERS
)
5 Correlations and regressions
5.1 \(\delta^{13}\text{C}_{\text{CH}_4}\) vs. F\(^{14}\)C\(_{\text{CH}_4}\)
cor.test( ~ CH4_Fm + CH4_d13C_permil_VPDB,
data=data_14CH4_NSHQ14_longer,
method = "pearson",
conf.level = 0.95)##
## Pearson's product-moment correlation
##
## data: CH4_Fm and CH4_d13C_permil_VPDB
## t = -5.4521, df = 3, p-value = 0.01212
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.9969981 -0.4448515
## sample estimates:
## cor
## -0.9530624# generate linear model
model_d13C_F14C = lm(CH4_d13C_permil_VPDB ~ CH4_Fm,
data = data_14CH4_NSHQ14_longer)
# print model summary
summary(model_d13C_F14C)##
## Call:
## lm(formula = CH4_d13C_permil_VPDB ~ CH4_Fm, data = data_14CH4_NSHQ14_longer)
##
## Residuals:
## 1 2 3 4 5
## 1.3368 0.7942 -1.6385 0.4595 -0.9519
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.850 1.066 7.363 0.00518 **
## CH4_Fm -34.923 6.405 -5.452 0.01212 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.44 on 3 degrees of freedom
## Multiple R-squared: 0.9083, Adjusted R-squared: 0.8778
## F-statistic: 29.73 on 1 and 3 DF, p-value: 0.01212Predict \(\delta^{13}\text{C}_{\text{CH}_4}\) (and 95% confidence interval) for F\(^{14}\)C\(_{\text{CH}_4}\) of 0 and 1
predict(model_d13C_F14C, newdata = data.frame(CH4_Fm=c(0,1)), interval = c("confidence"))## fit lwr upr
## 1 7.85011 4.457081 11.243138
## 2 -27.07277 -44.871755 -9.273792data_14CH4_NSHQ14_longer %>% ggplot(aes(x = CH4_Fm, y = CH4_d13C_permil_VPDB))+
geom_smooth(method = "lm", formula = y ~ x, fullrange = TRUE) +
stat_poly_eq(formula = y ~ x, parse = TRUE, aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~~")), label.x = "right") +
geom_point()+
scale_x_continuous(limits = c(0,1), expand = c(0,0), name = TeX("CH$_4$ fraction modern")) +
scale_y_continuous(name = TeX("$\\delta ^{13}C_{CH_4}\\,/\\,\\[$‰$ \\, VPDB \\]$")) +
theme_classic() +
theme(
plot.margin = margin(r = 12)
)
5.2 \(\delta\text{D}_{\text{CH}_4}\) vs. F\(^{14}\)C\(_{\text{CH}_4}\)
cor.test( ~ CH4_Fm + CH4_dD_permil_VSMOW,
data=data_14CH4_NSHQ14_longer,
method = "pearson",
conf.level = 0.95)##
## Pearson's product-moment correlation
##
## data: CH4_Fm and CH4_dD_permil_VSMOW
## t = -3.7882, df = 2, p-value = 0.06316
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.9987068 0.2433714
## sample estimates:
## cor
## -0.9368447# generate linear model
model_dD_F14C = lm(CH4_dD_permil_VSMOW ~ CH4_Fm,
data = data_14CH4_NSHQ14_longer)
# print model summary
summary(model_dD_F14C)##
## Call:
## lm(formula = CH4_dD_permil_VSMOW ~ CH4_Fm, data = data_14CH4_NSHQ14_longer)
##
## Residuals:
## 1 3 4 5
## -8.381 2.661 -6.741 12.461
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -242.71 10.30 -23.575 0.00179 **
## CH4_Fm -212.48 56.09 -3.788 0.06316 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11.79 on 2 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.8777, Adjusted R-squared: 0.8165
## F-statistic: 14.35 on 1 and 2 DF, p-value: 0.06316Predict \(\delta^{13}\text{C}_{\text{CH}_4}\) (and 95% confidence interval) for F\(^{14}\)C\(_{\text{CH}_4}\) of 0 and 1
predict(model_dD_F14C, newdata = data.frame(CH4_Fm=c(0,1)), interval = c("confidence"))## fit lwr upr
## 1 -242.7067 -287.0036 -198.4098
## 2 -455.1897 -661.7777 -248.6018data_14CH4_NSHQ14_longer %>% ggplot(aes(x = CH4_Fm, y = CH4_dD_permil_VSMOW))+
geom_smooth(method = "lm", formula = y ~ x, fullrange = TRUE) +
stat_poly_eq(formula = y ~ x, parse = TRUE, aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~~")), label.x = "right") +
geom_point()+
scale_x_continuous(limits = c(0,1), expand = c(0,0), name = TeX("CH$_4$ fraction modern")) +
scale_y_continuous(name = TeX("$\\delta D_{CH_4}\\,/\\,\\[$‰$ \\, VSMOW \\]$")) +
theme_classic() +
theme(
plot.margin = margin(r = 12)
)
5.3 \(\delta\text{D}_{\text{CH}_4}\) vs. \(\delta^{13}\text{C}_{\text{CH}_4}\)
5.3.1 Only data from NSHQ14 with paired \(\text{F}^{14}\text{C}_{\text{CH}_4}\) measurements
cor.test( ~ CH4_d13C_permil_VPDB + CH4_dD_permil_VSMOW,
data=data_14CH4_NSHQ14_longer,
method = "pearson",
conf.level = 0.95)##
## Pearson's product-moment correlation
##
## data: CH4_d13C_permil_VPDB and CH4_dD_permil_VSMOW
## t = 1.9057, df = 2, p-value = 0.197
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.6925588 0.9956742
## sample estimates:
## cor
## 0.8030376# generate linear model
model_dD_d13C = lm(CH4_dD_permil_VSMOW ~ CH4_d13C_permil_VPDB,
data = data_14CH4_NSHQ14_longer)
# print model summary
summary(model_dD_d13C)##
## Call:
## lm(formula = CH4_dD_permil_VSMOW ~ CH4_d13C_permil_VPDB, data = data_14CH4_NSHQ14_longer)
##
## Residuals:
## 1 3 4 5
## -13.43 13.28 -14.94 15.08
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -287.122 11.980 -23.966 0.00174 **
## CH4_d13C_permil_VPDB 5.192 2.724 1.906 0.19696
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 20.09 on 2 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.6449, Adjusted R-squared: 0.4673
## F-statistic: 3.632 on 1 and 2 DF, p-value: 0.197data_14CH4_NSHQ14_longer %>% ggplot(aes(x = CH4_d13C_permil_VPDB, y = CH4_dD_permil_VSMOW))+
geom_smooth(method = "lm", formula = y ~ x, fullrange = TRUE) +
stat_poly_eq(formula = y ~ x, parse = TRUE, aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~~")), label.x = "right") +
geom_point()+
scale_x_continuous(name = TeX("$\\delta ^{13}C_{CH_4}\\,/\\,\\[$‰$ \\, VSMOW \\]$")) +
scale_y_continuous(name = TeX("$\\delta D_{CH_4}\\,/\\,\\[$‰$ \\, VSMOW \\]$")) +
theme_classic() +
theme(
plot.margin = margin(r = 12)
)
5.3.2 All NSHQ14 data 2014-2019 data
Prepare Oman well data
# make data frame of CH4 d13C and dD data for all wells
data_CD <- data %>% select(sampling_site, notes, year_sampled, depth_fluid_intake_mbct, water_type, CH4_d13C_permil_VPDB_LBNL, CH4_dD_permil_VSMOW_LBNL,
CH4_d13C_permil_VPDB_CUB, CH4_dD_permil_VSMOW_CUB, CH4_d13C_permil_VPDB_MIT, CH4_dD_permil_VSMOW_MIT, CH4_d13C_permil_VPDB_UCLA, CH4_dD_permil_VSMOW_UCLA, site)
# compute means between analyses in multiple labs
data_CD <- data_CD %>%
rowwise() %>%
mutate(CH4_d13C_permil_VPDB_interlab_mean = mean(c(CH4_d13C_permil_VPDB_LBNL, CH4_d13C_permil_VPDB_CUB, CH4_d13C_permil_VPDB_MIT, CH4_d13C_permil_VPDB_UCLA), na.rm = TRUE))
data_CD <- data_CD %>%
rowwise() %>%
mutate(CH4_dD_permil_VSMOW_interlab_mean = mean(c(CH4_dD_permil_VSMOW_LBNL, CH4_dD_permil_VSMOW_CUB, CH4_dD_permil_VSMOW_MIT, CH4_dD_permil_VSMOW_UCLA), na.rm = TRUE))
# remove samples without C and D isotope data
data_CD <- data_CD %>% filter(!is.nan(CH4_d13C_permil_VPDB_interlab_mean) & !is.nan(CH4_dD_permil_VSMOW_interlab_mean))
# print
data_CD %>% select(sampling_site, notes, year_sampled, CH4_d13C_permil_VPDB_interlab_mean, CH4_dD_permil_VSMOW_interlab_mean)# make data frame of CH4 d13C data
data_C_gathered <- data_CD %>% select(sampling_site, notes, year_sampled, water_type, CH4_d13C_permil_VPDB_CUB, CH4_d13C_permil_VPDB_MIT, CH4_d13C_permil_VPDB_UCLA, CH4_d13C_permil_VPDB_LBNL, site, depth_fluid_intake_mbct) %>%
pivot_longer(cols = c(-sampling_site, -notes, -water_type, -year_sampled, -site, -depth_fluid_intake_mbct), names_to = "parameter", values_to = "CH4_d13C_permil_VPDB")
# make a column for analytical laboratory
data_C_gathered <- data_C_gathered %>% mutate(lab = case_when(
str_detect(parameter, "LBNL") == TRUE ~ "LBNL",
str_detect(parameter, "CUB") == TRUE ~ "CUB",
str_detect(parameter, "MIT") == TRUE ~ "MIT",
str_detect(parameter, "UCLA") == TRUE ~ "UCLA"
)) %>%
select(-parameter)
# data_C_gathered
# repeat the above for dD
data_D_gathered <- data_CD %>% select(sampling_site, notes, water_type, year_sampled, CH4_dD_permil_VSMOW_CUB, CH4_dD_permil_VSMOW_MIT, CH4_dD_permil_VSMOW_UCLA, CH4_dD_permil_VSMOW_LBNL) %>%
pivot_longer(cols = c(-sampling_site, -notes, -water_type, -year_sampled), names_to = "parameter", values_to = "CH4_dD_permil_VSMOW")
data_D_gathered <- data_D_gathered %>% mutate(lab = case_when(
str_detect(parameter, "LBNL") == TRUE ~ "LBNL",
str_detect(parameter, "CUB") == TRUE ~ "CUB",
str_detect(parameter, "MIT") == TRUE ~ "MIT",
str_detect(parameter, "UCLA") == TRUE ~ "UCLA"
)) %>% select(-parameter)
# combine 13C and D data
data_CD_gathered <- full_join(data_C_gathered, data_D_gathered) %>%
filter(!is.na(CH4_d13C_permil_VPDB) | !is.na(CH4_dD_permil_VSMOW))## Joining with `by = join_by(sampling_site, notes, year_sampled, water_type,
## lab)`data_CD_gathered_NSHQ14 <- data_CD_gathered %>%
# only samples from NSHQ14
filter(sampling_site == "NSHQ14") %>%
# only samples with CH4 d13C and CH4 dD measurements
filter(!is.na(CH4_d13C_permil_VPDB) & !is.na(CH4_dD_permil_VSMOW))
# print
data_CD_gathered_NSHQ14# get sample size
data_CD_gathered_NSHQ14 %>% summarise(n = n())cor.test( ~ CH4_d13C_permil_VPDB + CH4_dD_permil_VSMOW,
data = data_CD_gathered_NSHQ14,
method = "pearson",
conf.level = 0.95)##
## Pearson's product-moment correlation
##
## data: CH4_d13C_permil_VPDB and CH4_dD_permil_VSMOW
## t = 4.0233, df = 8, p-value = 0.003823
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.3886158 0.9555266
## sample estimates:
## cor
## 0.8180736# generate linear model
model_dD_d13C = lm(CH4_dD_permil_VSMOW ~ CH4_d13C_permil_VPDB,
data = data_CD_gathered_NSHQ14)
# print model summary
summary(model_dD_d13C)##
## Call:
## lm(formula = CH4_dD_permil_VSMOW ~ CH4_d13C_permil_VPDB, data = data_CD_gathered_NSHQ14)
##
## Residuals:
## Min 1Q Median 3Q Max
## -24.029 -6.857 3.163 5.851 30.721
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -278.146 5.323 -52.252 1.99e-11 ***
## CH4_d13C_permil_VPDB 5.142 1.278 4.023 0.00382 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 16.83 on 8 degrees of freedom
## Multiple R-squared: 0.6692, Adjusted R-squared: 0.6279
## F-statistic: 16.19 on 1 and 8 DF, p-value: 0.003823data_CD_gathered_NSHQ14 %>% ggplot(aes(x = CH4_d13C_permil_VPDB, y = CH4_dD_permil_VSMOW))+
geom_smooth(method = "lm", formula = y ~ x, fullrange = TRUE) +
stat_poly_eq(formula = y ~ x, parse = TRUE, aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~~"))) +
geom_point()+
scale_x_continuous(name = TeX("$\\delta ^{13}C_{CH_4}\\,/\\,\\[$‰$ \\, VSMOW \\]$")) +
scale_y_continuous(name = TeX("$\\delta D_{CH_4}\\,/\\,\\[$‰$ \\, VSMOW \\]$")) +
theme_classic()
6 \(\alpha_{\text{methane/water}}\)
# solve CH4-H2O equilibrium (Horibe and Craig, 1995)
eq_CH4_H2O <- c(35, 270) # input temperatures, 35˚C for Samail Ophiolite groundwater, 270˚C for CH4-H2 closure temperature (Wang et al., 2018)
eq_CH4_H2O <- as.data.frame(eq_CH4_H2O) # convert above vector into data frame
eq_CH4_H2O <- eq_CH4_H2O %>% rename(temp_C = eq_CH4_H2O) # rename column
eq_CH4_H2O <- eq_CH4_H2O %>% mutate(temp_K = temp_C + 273.15) # add column for temperature in Kelvin
eq_CH4_H2O <- eq_CH4_H2O %>% mutate(alph_H2O_CH4 = 1.0997 + 8456 / temp_K^2 + 0.9611e9 / temp_K^4 - 27.82e12 / temp_K^6) # calculate alpha H2O/CH4 as function of temperature using equation from Horibe and Craig (1995)
eq_CH4_H2O <- eq_CH4_H2O %>% mutate(delta_H2O_VSMOW = 0) # use VSMOW as H2O dD
eq_CH4_H2O <- eq_CH4_H2O %>% mutate(delta_CH4_VSMOW = (delta_H2O_VSMOW + 1)/(alph_H2O_CH4) - 1) # calculate dD CH4
eq_CH4_H2O <- eq_CH4_H2O %>% mutate(delta_CH4_VSMOW_permil = delta_CH4_VSMOW * 1000) # convert to permil
eq_CH4_H2O7 C\(_1\) / (C\(_2\) + C\(_3\)) vs. \(\text{F}^{14}\text{C}_{\text{CH}_4}\) at NSHQ14
data_C1_C2_C3_NSHQ14_2017_thru_2019 <- data_C1_C2_C3 %>% filter(sampling_site == "NSHQ14") %>% filter(year_sampled >= 2017 & year_sampled <= 2019) %>% filter(depth_fluid_intake_mbct == 85) %>% select(year_sampled, C1_over_C2_plus_C3, notes)
data_C1_C2_C3_NSHQ14_2017_thru_2019data_C1_C2_C3_NSHQ14_2017_thru_2019_14C <- data_C1_C2_C3_NSHQ14_2017_thru_2019 %>% left_join(data_14CH4_NSHQ14 %>% select(year_sampled, notes, CH4_Fm_UCIAMS), by = c("year_sampled", "notes"))
data_C1_C2_C3_NSHQ14_2017_thru_2019_14Cdata_C1_C2_C3_NSHQ14_2017_thru_2019_14C %>% ggplot(aes(x = CH4_Fm_UCIAMS, y = C1_over_C2_plus_C3, color = as.character(year_sampled), shape = notes))+
geom_point()+
scale_shape_discrete(na.value = 17) +
scale_color_manual(name = "year sampled", values = cbPalette_MWL) +
scale_x_continuous(name = TeX("CH$_4$ fraction modern"), limits = c(0,0.31), expand = c(0,0)) +
scale_y_continuous(name = TeX("$C_1$/($C_{2}$+$C_{3}$)"), limits = c(0,1200), expand = c(0,0)) +
theme_classic()
8 Rate comparisons
Methane concentration in 2017 in \(\mu\text{mol}\,\text{L}^{-1}\)
conc_CH4_NSHQ14_2017_umol_per_L <- data %>% filter(sampling_site == "NSHQ14") %>%
filter(year_sampled == 2017) %>%
filter(depth_fluid_intake_mbct == 85) %>%
select(CH4_uM) %>% as.numeric()
conc_CH4_NSHQ14_2017_umol_per_L## [1] 105.5295Methane fraction modern in 2017
Fm_CH4_NSHQ14_2017 <- data %>% filter(sampling_site == "NSHQ14") %>%
filter(year_sampled == 2017) %>%
filter(depth_fluid_intake_mbct == 85) %>%
select(CH4_Fm_UCIAMS) %>% as.numeric()
Fm_CH4_NSHQ14_2017## [1] 0.1918Concentration of modern methane in 2017
conc_modern_CH4_NSHQ14_2017_umol_per_L <- conc_CH4_NSHQ14_2017_umol_per_L * Fm_CH4_NSHQ14_2017
conc_modern_CH4_NSHQ14_2017_umol_per_L## [1] 20.24056Rate comparison with Fones et al. 2019 “Physiological adaptations to serpentinization in the Samail Ophiolite, Oman” ISME J Calculate minimum time to form modern methane by microbial methanogenesis, in years
NSHQ14_2017_max_methanogenesis_rate_pmol_per_mL_per_day <- 18 # Fones et al. 2019
NSHQ14_2017_max_methanogenesis_rate_pmol_per_L_per_day <- NSHQ14_2017_max_methanogenesis_rate_pmol_per_mL_per_day * 1000
NSHQ14_2017_max_methanogenesis_rate_umol_per_L_per_day <- NSHQ14_2017_max_methanogenesis_rate_pmol_per_L_per_day / 1e6
NSHQ14_2017_max_methanogenesis_rate_umol_per_L_per_year <- NSHQ14_2017_max_methanogenesis_rate_umol_per_L_per_day * 365
NSHQ14_2017_max_methanogenesis_rate_umol_per_L_per_year## [1] 6.57min_time_to_form_modern_CH4_by_methanogenesis_NSHQ14_2017_year <- conc_modern_CH4_NSHQ14_2017_umol_per_L / NSHQ14_2017_max_methanogenesis_rate_umol_per_L_per_year
min_time_to_form_modern_CH4_by_methanogenesis_NSHQ14_2017_year## [1] 3.0807559 Sulfate and hydrogen comparisons
data %>% filter(!is.na(CH4_Fm_UCIAMS)) %>% select(
sampling_site, year_sampled, month_sampled, day_sampled, depth_fluid_intake_mbct, CH4_Fm_UCIAMS, Sulfate_uM, Sulfate_was_measured, H2_uM, H2_was_measured, H2_LQ_uM)10 Henrys constant calculations
R_m3_Pa_per_K_mol <- 8.31446261815324 # Gas constant
# Calculate Henry's constants at 22 ˚C and 37 ˚C
Henrys <- Henrys %>% mutate(Hcp_mol_per_m3_Pa_22C = Hcp_mol_per_m3_Pa_25C * exp(dlnHcp_over_d_inv_T * (1/(22+273.15) - 1/(25+273.15))))
Henrys <- Henrys %>% mutate(Hcp_mol_per_m3_Pa_37C = Hcp_mol_per_m3_Pa_25C * exp(dlnHcp_over_d_inv_T * (1/(37+273.15) - 1/(25+273.15))))
# Calculate Henry's constants at 22 ˚C and 37 ˚C
Henrys <- Henrys %>% mutate(Hcc_22C = Hcp_mol_per_m3_Pa_22C * R_m3_Pa_per_K_mol * (22+273.15))
Henrys <- Henrys %>% mutate(Hcc_37C = Hcp_mol_per_m3_Pa_37C * R_m3_Pa_per_K_mol * (37+273.15))
# Calculate % difference at 37 ˚C vs 22 ˚C
Henrys <- Henrys %>% mutate(Hcp_perc_diff_37C_22C = (Hcp_mol_per_m3_Pa_37C - Hcp_mol_per_m3_Pa_22C) / Hcp_mol_per_m3_Pa_22C * 100)
Henrys <- Henrys %>% mutate(Hcc_perc_diff_37C_22C = (Hcc_37C - Hcc_22C) / Hcc_22C * 100)
HenrysGet relevant gas data
data_14CH4_2019_gas_concs <- data %>% filter(year_sampled == 2019) %>% filter(sampling_site == "NSHQ14" | sampling_site == "WAB71") %>% select(sampling_site, notes, CH4_uM, CH4_LQ_uM, C2H6_uM, C2H6_LQ_uM, C3H8_uM, C3H8_LQ_uM, i_C4H10_uM, i_C4H10_LQ_uM, n_C4H10_uM, n_C4H10_LQ_uM, i_C5H12_uM, i_C5H12_LQ_uM, n_C5H12_uM, n_C5H12_LQ_uM, C6H14_uM, C6H12_LQ_uM)
data_14CH4_2019_gas_concscalculate equilibrium gas phase concentration
data_14CH4_2019_gas_concs_gas_phase <- data_14CH4_2019_gas_concs %>%
mutate(CH4_uM = CH4_uM / Henrys %>% filter(compound == "methane") %>% select(Hcc_22C) %>% as.numeric(),
CH4_LQ_uM = CH4_LQ_uM / Henrys %>% filter(compound == "methane") %>% select(Hcc_22C) %>% as.numeric(),
C2H6_uM = C2H6_uM / Henrys %>% filter(compound == "ethane") %>% select(Hcc_22C) %>% as.numeric(),
C2H6_LQ_uM = C2H6_LQ_uM / Henrys %>% filter(compound == "ethane") %>% select(Hcc_22C) %>% as.numeric(),
C3H8_uM = C3H8_uM / Henrys %>% filter(compound == "propane") %>% select(Hcc_22C) %>% as.numeric(),
C3H8_LQ_uM = C3H8_LQ_uM / Henrys %>% filter(compound == "propane") %>% select(Hcc_22C) %>% as.numeric(),
i_C4H10_uM = i_C4H10_uM / Henrys %>% filter(compound == "iso-butane") %>% select(Hcc_22C) %>% as.numeric(),
i_C4H10_LQ_uM = i_C4H10_LQ_uM / Henrys %>% filter(compound == "iso-butane") %>% select(Hcc_22C) %>% as.numeric(),
n_C4H10_uM = n_C4H10_uM / Henrys %>% filter(compound == "n-butane") %>% select(Hcc_22C) %>% as.numeric(),
n_C4H10_LQ_uM = n_C4H10_LQ_uM / Henrys %>% filter(compound == "n-butane") %>% select(Hcc_22C) %>% as.numeric(),
i_C5H12_uM = i_C5H12_uM / Henrys %>% filter(compound == "iso-pentane") %>% select(Hcc_22C) %>% as.numeric(),
i_C5H12_LQ_uM = i_C5H12_LQ_uM / Henrys %>% filter(compound == "iso-pentane") %>% select(Hcc_22C) %>% as.numeric(),
n_C5H12_uM = n_C5H12_uM / Henrys %>% filter(compound == "n-pentane") %>% select(Hcc_22C) %>% as.numeric(),
n_C5H12_LQ_uM = n_C5H12_LQ_uM / Henrys %>% filter(compound == "n-pentane") %>% select(Hcc_22C) %>% as.numeric(),
C6H14_uM = C6H14_uM / Henrys %>% filter(compound == "hexane") %>% select(Hcc_22C) %>% as.numeric(),
C6H12_LQ_uM = C6H12_LQ_uM / Henrys %>% filter(compound == "hexane") %>% select(Hcc_22C) %>% as.numeric()
)
data_14CH4_2019_gas_concs_gas_phase# calculate C1/(C2+C3)
data_14CH4_2019_gas_concs_gas_phase <- data_14CH4_2019_gas_concs_gas_phase %>% mutate(C1_over_C2_plus_C3 = CH4_uM / (C2H6_uM + if_else(!is.na(C3H8_uM), C3H8_uM, 0)))
data_14CH4_2019_gas_concs_gas_phase %>% select(sampling_site, C1_over_C2_plus_C3, notes)# calculate C1/(C moles in total C1-C6)
data_14CH4_2019_gas_concs_gas_phase <- data_14CH4_2019_gas_concs_gas_phase %>% mutate(mol_C_in_C1_over_mol_C_C1_thru_C6 = CH4_uM / (CH4_uM + if_else(is.na(C2H6_uM) == FALSE, C2H6_uM, 0) * 2 + if_else(is.na(C3H8_uM) == FALSE, C3H8_uM, 0) * 3 + if_else(is.na(i_C4H10_uM) == FALSE, i_C4H10_uM, 0) * 4 + if_else(is.na(n_C4H10_uM) == FALSE, n_C4H10_uM, 0) * 4 + if_else(is.na(i_C5H12_uM) == FALSE, i_C4H10_uM, 0) * 5 + if_else(is.na(n_C5H12_uM) == FALSE, n_C5H12_uM, 0) * 5 + if_else(is.na(C6H14_uM) == FALSE, C6H14_uM, 0) * 6))
data_14CH4_2019_gas_concs_gas_phase %>% select(sampling_site, mol_C_in_C1_over_mol_C_C1_thru_C6, notes)Calculate what the F\(^{14}\)C of the total gaseous hydrocarbon pool would be if CH\(_4\) had F\(^{14}\)C = 0 and non-methane hydrocarbon gases (NMHC) had F\(^{14}\)C = 1 assuming gas concentrations ratios from samples in prior chunk
F14C_NMHC_hypo <- 1 # hypothetical assumption of fraction modern of total non-methane hydrocarbon pool
F14C_CH4_hypo <- 0 # hypothetical assumption of fraction modern of total non-methane hydrocarbon pool
data_14CH4_2019_gas_concs_gas_phase_hypo_dead_CH4 <- data_14CH4_2019_gas_concs_gas_phase %>%
mutate(F14C_tot_HC_hypo = mol_C_in_C1_over_mol_C_C1_thru_C6 * F14C_CH4_hypo + (1 - mol_C_in_C1_over_mol_C_C1_thru_C6) * F14C_NMHC_hypo)
data_14CH4_2019_gas_concs_gas_phase_hypo_dead_CH4 %>% select(sampling_site, F14C_tot_HC_hypo, mol_C_in_C1_over_mol_C_C1_thru_C6, notes)Calculate what the F\(^{14}\)C of the total gaseous hydrocarbon pool would be if CH\(_4\) had F\(^{14}\)C = 1 and non-methane hydrocarbon gases (NMHC) had F\(^{14}\)C = 0 assuming gas concentrations ratios from samples in prior chunk
F14C_NMHC_hypo <- 0 # hypothetical assumption of fraction modern of total non-methane hydrocarbon pool
F14C_CH4_hypo <- 1 # hypothetical assumption of fraction modern of total non-methane hydrocarbon pool
data_14CH4_2019_gas_concs_gas_phase_hypo_dead_NMHC <- data_14CH4_2019_gas_concs_gas_phase %>%
mutate(F14C_tot_HC_hypo = mol_C_in_C1_over_mol_C_C1_thru_C6 * F14C_CH4_hypo + (1 - mol_C_in_C1_over_mol_C_C1_thru_C6) * F14C_NMHC_hypo)
data_14CH4_2019_gas_concs_gas_phase_hypo_dead_NMHC %>% select(sampling_site, F14C_tot_HC_hypo, mol_C_in_C1_over_mol_C_C1_thru_C6, notes)Calculate what the F\(^{14}\)C of the total gaseous hydrocarbon pool would be if CH\(_4\) had F\(^{14}\)C of measured value of NSHQ14 2019 prepared at Isotech and measured at NOSAMS and non-methane hydrocarbon gases (NMHC) had F\(^{14}\)C = 0 assuming gas concentrations ratios from samples in prior chunk
F14C_NMHC_hypo <- 0 # hypothetical assumption of fraction modern of total non-methane hydrocarbon pool
(F14C_CH4_measured_NOSAMS <- data %>% filter(sampling_site == "NSHQ14" & year_sampled == 2019 & is.na(notes) == T) %>% select(CH4_Fm_NOSAMS) %>% as.character() %>% as.numeric()) # hypothetical assumption of fraction modern of total non-methane hydrocarbon pool## [1] 0.0614data_14CH4_2019_gas_concs_gas_phase_hypo_dead_NMHC_measured_NSHQ14_f14C_CH4 <- data_14CH4_2019_gas_concs_gas_phase %>% filter(sampling_site == "NSHQ14" & is.na(notes) == T) %>%
mutate(F14C_tot_HC_hypo = mol_C_in_C1_over_mol_C_C1_thru_C6 * F14C_CH4_measured_NOSAMS + (1 - mol_C_in_C1_over_mol_C_C1_thru_C6) * F14C_NMHC_hypo)
data_14CH4_2019_gas_concs_gas_phase_hypo_dead_NMHC_measured_NSHQ14_f14C_CH4 %>% select(sampling_site, F14C_tot_HC_hypo, mol_C_in_C1_over_mol_C_C1_thru_C6, notes)F14C_CH4_hypo <- data_14CH4_2019_gas_concs_gas_phase_hypo_dead_NMHC_measured_NSHQ14_f14C_CH4 %>% select(F14C_tot_HC_hypo) %>% as.character %>% as.numeric()
# calculate difference between measured and hypothetical F14C CH4 at NSHQ14 in 2019
(F14C_CH4_diff_measured_hypo <- F14C_CH4_measured_NOSAMS - F14C_CH4_hypo)## [1] 0.0002654512F\(^{14}\)C of total hydrocarbon gases in NSHQ14 in 2019 measured at UC-Irvine and prepared at CU-Boulder
(F14C_CH4_measured_UCI <- data %>% filter(sampling_site == "NSHQ14" & year_sampled == 2019 & is.na(notes) == T) %>% select(CH4_Fm_UCIAMS) %>% as.character() %>% as.numeric())## [1] 0.0566# calculate difference between measured F14C CH4 at NSHQ14 in 2019 at NOSAMS and UCI
(F14C_CH4_diff_NOSAMS_UCI <- F14C_CH4_measured_NOSAMS - F14C_CH4_measured_UCI)## [1] 0.0048% of difference between NOSAMS and UCI F\(^{14}\)C attributable to dilution of NMHCs
F14C_CH4_diff_measured_hypo / F14C_CH4_diff_NOSAMS_UCI * 100## [1] 5.530233Calculate relative difference (%) of the UCI and NOSAMS samples
(F14C_CH4_measured_NOSAMS - F14C_CH4_measured_UCI) / mean(F14C_CH4_measured_NOSAMS, F14C_CH4_measured_UCI) * 100## [1] 7.81759