---
title: "Oman methane isotope plotting"
subtitle: "Source file: 200813_Oman_14CH4_paper_figs.Rmd"
author: "Daniel Nothaft"
date: "`r format(Sys.Date(), '%d %b %Y')`"
output:
  html_document: 
    df_print: paged # omit to disable paged data table output
    css: stylesheet.css # omit if no need for custom stylesheet
    number_sections: yes # change to no for unnumbered sections
    toc: yes # change to no to disable table of contents
    toc_float: true # change to false to keep toc at the top
    toc_depth: 3 # change to specify which headings to include in toc
    code_folding: show # change to hide to hide code by default
editor_options:
  chunk_output_type: inline
---
# Setup
## Load code libraries
```{r setup, message=FALSE, warning=FALSE}
# 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())))
)
```

```{r source}
# source all relevant scripting files
source(file.path("scripts", "plotting_functions.R")) # includes LaTeX labelling function
```

## Load data
Load in Oman well data
```{r load-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
data
```

Load in comparison data from the literature
```{r load-lit-data}
# 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_literature
```

Load in Henry's constants from Sander, 2015
```{r load-henrys}
# read Henry's constants from literature
Henrys <- read_csv("data_raw/Henrys_constants_Sander.csv")
```

## Set plotting aesthetics

```{r set-plot-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)
```



# C$_1$ / (C$_2$ + C$_3$) vs. $\delta^{13}\text{C}_{\text{CH}_4}$ ("Bernard") Plot
## Prepare data
Prepare Oman well data
```{r get-Oman-C1-C2-C3-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_C3
```

Get literature data
```{r get-lit-C1-C2-C3-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_lit
```

## Plot
```{r C1_C2_C3_plot, fig.width=3.4, fig.height=3.4, warning=FALSE}
# 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
```

# Alkanes $\delta^{13}$C plot

## Prepare data
```{r get-alkane-d13C-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_d13C
```

Gather alkane $\delta^{13}\text{C}$ data so it can be plotted
```{r gather-alkane-d13C-data}
# 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_gathered
```

Add a site and sample type to Oman well data so that it is compatible for plotting with literature data
```{r mutate-alkane-d13C-data-site}
data_alkanes_d13C_gathered <- data_alkanes_d13C_gathered %>% mutate(site = "Oman/UAE", sample_type = "vent/well\nfluid")

data_alkanes_d13C_gathered
```

Get literature data
```{r get-alkane-d13C-lit-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_lit
```

As previously done for Oman well data, gather literature data by compound for plotting purposes
```{r gather-alkane-d13C-data-lit}
# 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_lit
```

## Plot
### Plot data
```{r alkanes-d13C-oman-wells-and-lit-plot, fig.height=6.5, fig.width=3.42, warning=FALSE}
# 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
```

```{r alkanes-d13C-oman-only-plot, fig.height=3, fig.width=3.42, warning=FALSE}
# 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
```


# CH$_4$ Fraction modern facet grid plot

## Prepare data
```{r get-14CH4-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
```

```{r filter-14CH4-for-NSHQ14}
(data_14CH4_NSHQ14 <- data_14CH4 %>% filter(sampling_site == "NSHQ14"))
```

```{r pivot-longer-14CH4-data}
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_longer
```

```{r pivot-longer-14CH4-data-again}
data_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
```

```{r 14CH4-tex}
# 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_))
```

## Plot
```{r, "plot-CH4_F14C-d13C-dD-lab-color", fig.width = 3.42, fig.height = 3.5, warning=FALSE}
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
```




```{r, "plot-CH4_F14C-d13C-dD-year-color-shape", fig.width = 3.42, fig.height = 4.5, warning=FALSE}
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
```

```{r, "plot-CH4_F14C-d13C-dD-year-color-shape-clean", fig.width = 3.42, fig.height = 4.5, warning=FALSE}
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
```

```{r, "plot-CH4_F14C-d13C-dD-year-color-shape-clean-letter-tagged", fig.width = 3.42, fig.height = 4.5, warning=FALSE}
CH4_F14C_d13C_dD_plot_year_color_shape_clean %>% tag_facet_keep_strip(
  open = "",
  close = "",
  tag_pool = LETTERS
)
```

# Correlations and regressions
## $\delta^{13}\text{C}_{\text{CH}_4}$ vs. F$^{14}$C$_{\text{CH}_4}$
```{r cor-test-d13C-Fm}
cor.test( ~ CH4_Fm + CH4_d13C_permil_VPDB,
         data=data_14CH4_NSHQ14_longer,
         method = "pearson",
         conf.level = 0.95)
```

```{r lm-d13C-Fm}
# 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)
```

Predict $\delta^{13}\text{C}_{\text{CH}_4}$ (and 95% confidence interval) for F$^{14}$C$_{\text{CH}_4}$ of 0 and 1
```{r predict-d13C}
predict(model_d13C_F14C, newdata = data.frame(CH4_Fm=c(0,1)), interval = c("confidence"))
```

```{r plot-d13C-Fm-w-reg, warning=FALSE}
data_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)
  )
```

## $\delta\text{D}_{\text{CH}_4}$ vs. F$^{14}$C$_{\text{CH}_4}$
```{r cor-test-D-Fm}
cor.test( ~ CH4_Fm + CH4_dD_permil_VSMOW,
         data=data_14CH4_NSHQ14_longer,
         method = "pearson",
         conf.level = 0.95)
```

```{r lm-D-Fm}
# 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)
```
Predict $\delta^{13}\text{C}_{\text{CH}_4}$ (and 95% confidence interval) for F$^{14}$C$_{\text{CH}_4}$ of 0 and 1
```{r predict-dD}
predict(model_dD_F14C, newdata = data.frame(CH4_Fm=c(0,1)), interval = c("confidence"))
```

```{r plot-D-Fm-w-reg, warning=FALSE}
data_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)
  )
```

## $\delta\text{D}_{\text{CH}_4}$ vs. $\delta^{13}\text{C}_{\text{CH}_4}$
### Only data from NSHQ14 with paired $\text{F}^{14}\text{C}_{\text{CH}_4}$ measurements
```{r cor-test-D-13C-w-paired-14CH4}
cor.test( ~ CH4_d13C_permil_VPDB + CH4_dD_permil_VSMOW,
         data=data_14CH4_NSHQ14_longer,
         method = "pearson",
         conf.level = 0.95)
```

```{r lm-D-13C-w-paired-14CH4}
# 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)
```

```{r plot-D-13C-w-reg-w-paired-14CH4, warning=FALSE}
data_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)
  )
```

### All NSHQ14 data 2014-2019 data
Prepare Oman well data

```{r data_CD}
# 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)
```


```{r data_CD_gathered}
# 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))
```


```{r get-CD-isotope-data-for-NSHQ14}
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
```

```{r CD_NSHQ14_n}
# get sample size
data_CD_gathered_NSHQ14 %>% summarise(n = n())
```


```{r cor-test-D-13C-2014-2019}
cor.test( ~ CH4_d13C_permil_VPDB + CH4_dD_permil_VSMOW,
         data = data_CD_gathered_NSHQ14,
         method = "pearson",
         conf.level = 0.95)
```

```{r lm-D-13C-2014-2019}
# 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)
```

```{r plot-D-13C-w-reg-2014-2019, warning=FALSE}
data_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()
```

# $\alpha_{\text{methane/water}}$

```{r}
# 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_H2O
```
# C$_1$ / (C$_2$ + C$_3$) vs. $\text{F}^{14}\text{C}_{\text{CH}_4}$ at NSHQ14

```{r data_C1_C2_C3_NSHQ14_2017_thru_2019}
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_2019
```

```{r join_14C_to_C1_C2_C3_data}
data_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_14C
```

```{r plot-C1_C2_C3-14CH4-NSHQ14, warning=FALSE}
data_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()
```

# Rate comparisons
Methane concentration in 2017 in $\mu\text{mol}\,\text{L}^{-1}$
```{r conc-CH4-NSHQ14-2017}
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
```

Methane fraction modern in 2017
```{r F14C-CH4-NSHQ14-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
```

Concentration of modern methane in 2017
```{r conc-modern-CH4-NSHQ14-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
```

Rate 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
```{r methanogenesis-rates-Fones}
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
```


```{r min-methanogenesis-time}
min_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
```


# Sulfate and hydrogen comparisons

```{r}
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)
```

# Henrys constant calculations
```{r}
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)

Henrys
```

Get relevant gas data
```{r}
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_concs
```
calculate equilibrium gas phase concentration
```{r}
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
```

```{r}
# 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)

```

```{r}
# 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

```{r}
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

```{r}
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
```{r}
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

data_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)
```

F$^{14}$C of total hydrocarbon gases in NSHQ14 in 2019 measured at UC-Irvine and prepared at CU-Boulder

```{r}
(F14C_CH4_measured_UCI <- data %>% filter(sampling_site == "NSHQ14" & year_sampled == 2019 & is.na(notes) == T) %>% select(CH4_Fm_UCIAMS) %>% as.character() %>% as.numeric())

# 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)
```

% of difference between NOSAMS and UCI F$^{14}$C attributable to dilution of NMHCs
```{r}
F14C_CH4_diff_measured_hypo / F14C_CH4_diff_NOSAMS_UCI * 100
```

Calculate relative difference (%) of the UCI and NOSAMS samples
```{r}
(F14C_CH4_measured_NOSAMS - F14C_CH4_measured_UCI) / mean(F14C_CH4_measured_NOSAMS, F14C_CH4_measured_UCI) * 100
```

