Setup

Load code libraries

# load libraries
library(tidyverse) # dplyr, tidyr, ggplot
library(readxl) # reading excel files
library(knitr) # generating reports
library(ggpmisc) # add equations of best fit to ggplot
library(scales) # making log scales with exponents
library(lemon) # repeat axis lines for paneled figures in ggplot

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

Load data

# read NSHQ14 borehole log from March 17th, 2018
NSHQ14_profile_raw <- read_xlsx("data_raw/Summary_metadata/NSHQ14_QL40IDRO_20180317_down_5cm.xlsx", na = c("-99999.00", "-99999"))

# read BA3A borehole log from March 18th, 2018
BA3A_profile_raw <- read_xlsx("data_raw/Summary_metadata/BA3A_QL40IDRO_20180318_down_5cm.xlsx", na = c("-99999.00", "-99999"))

# Read in conversion from Ag/AgCl electrode to SHE as function of T
Ag_AgCl_E_T <- read_xlsx("data_raw/Summary_metadata/Ag_AgCl_E_T.xlsx", na = "NA")

NSHQ14 well log

clean up data

# select relevant columns and rename them
NSHQ14_profile <- NSHQ14_profile_raw %>% select(depth_sampled_mbgl = `DEPT[M]`, P_dbar = PRESSURE, T_C = TEMPERATURE, pH = PH, ORP_mV_Ag_AgCl = REDOX, cond_uS_per_cm = `COND(FW)`)

# Filter out nonsense values (i.e. when probe is above water table)
NSHQ14_profile <- NSHQ14_profile %>% filter(!is.na(ORP_mV_Ag_AgCl))

# print some summary statistics
NSHQ14_profile %>% summary()
##  depth_sampled_mbgl     P_dbar             T_C              pH       
##  Min.   : 11.21     Min.   :  4.967   Min.   :35.13   Min.   :10.36  
##  1st Qu.: 81.70     1st Qu.: 74.997   1st Qu.:35.98   1st Qu.:11.05  
##  Median :152.19     Median :144.314   Median :37.62   Median :11.12  
##  Mean   :152.19     Mean   :143.740   Mean   :37.75   Mean   :11.10  
##  3rd Qu.:222.69     3rd Qu.:212.746   3rd Qu.:39.39   3rd Qu.:11.19  
##  Max.   :293.18     Max.   :280.155   Max.   :41.17   Max.   :11.26  
##  ORP_mV_Ag_AgCl    cond_uS_per_cm
##  Min.   :-902.85   Min.   :1595  
##  1st Qu.:-899.70   1st Qu.:3638  
##  Median :-889.79   Median :3733  
##  Mean   :-850.94   Mean   :4116  
##  3rd Qu.:-870.75   3rd Qu.:4905  
##  Max.   : -44.29   Max.   :5027
# subset the data for readability and print
NSHQ14_profile %>% slice(seq(from = 1, to = 30000, by = 1000)) %>% kable(digits=1)
depth_sampled_mbgl P_dbar T_C pH ORP_mV_Ag_AgCl cond_uS_per_cm
11.2 5.1 35.2 10.4 -44.3 1594.8
21.2 15.0 35.1 10.9 -604.0 2541.9
31.2 25.0 35.2 11.2 -768.4 3518.2
41.2 34.8 35.3 11.2 -786.0 3598.0
51.2 44.8 35.4 11.2 -827.5 3620.0
61.2 54.8 35.6 11.2 -847.5 3633.2
71.2 64.6 35.8 11.2 -857.4 3640.1
81.2 74.5 36.0 11.2 -870.3 3634.2
91.2 84.3 36.2 11.2 -876.2 3639.8
101.2 94.3 36.4 11.2 -881.6 3648.9
111.2 104.1 36.6 11.2 -883.0 3640.8
121.2 114.0 36.9 11.2 -884.6 3643.2
131.2 123.7 37.1 11.2 -886.3 3667.8
141.2 133.5 37.3 11.1 -888.2 3701.7
151.2 143.4 37.6 11.2 -890.4 3724.2
161.2 153.1 37.8 11.1 -890.6 3785.9
171.2 162.8 38.1 11.1 -892.7 3961.8
181.2 172.6 38.3 11.1 -898.3 4782.5
191.2 182.4 38.6 11.1 -900.2 4839.6
201.2 191.8 38.8 11.1 -900.0 4866.1
211.2 201.7 39.1 11.1 -900.2 4882.6
221.2 211.4 39.3 11.1 -900.7 4901.8
231.2 220.9 39.6 11.1 -900.0 4921.4
241.2 230.6 39.9 11.0 -899.9 4943.5
251.2 240.2 40.1 11.0 -899.9 4956.6
261.2 249.7 40.4 11.0 -899.4 4969.3
271.2 259.3 40.6 11.0 -899.1 4985.4
281.2 268.8 40.9 11.0 -898.8 5003.5
291.2 278.1 41.1 11.0 -898.5 5024.1

Convert to ORP to SHE

Plot conversion between Ag/AgCl electrode and SHE as function of temperature

formula <- Ag_AgCl_E_T$E_Ag_AgCl_KCl_sat_mV ~ Ag_AgCl_E_T$T_C

plot_redox_offset_T <- Ag_AgCl_E_T %>% ggplot(aes(x = T_C, y = E_Ag_AgCl_KCl_sat_mV))+
  geom_point()+
  geom_smooth(method="lm")+
    stat_poly_eq(aes(label =  paste(stat(eq.label), stat(rr.label), sep = "~~~~")), formula = y ~ x, parse = TRUE, rr.digits = 2, label.x = .2, label.y = .9)+
  theme_bw(base_size = 12)+
  theme(
    panel.grid = element_blank())

plot_redox_offset_T
## Warning: `stat(eq.label)` was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(eq.label)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## `geom_smooth()` using formula = 'y ~ x'

Save a linear model of the relationship between Ag/AgCl electrode and SHE as function of temperature

Ag_AgCl_E_T_lm <- lm(formula)

Ag_AgCl_E_T_lm_summ <- Ag_AgCl_E_T_lm %>% summary()

Ag_AgCl_E_T_lm_summ
## 
## Call:
## lm(formula = formula)
## 
## Residuals:
##     1     2     3     4     5 
## -0.08  0.05  0.08  0.01 -0.06 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     223.230000   0.179722  1242.1 1.15e-09 ***
## Ag_AgCl_E_T$T_C  -1.046000   0.005033  -207.8 2.46e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07958 on 3 degrees of freedom
## Multiple R-squared:  0.9999, Adjusted R-squared:  0.9999 
## F-statistic: 4.319e+04 on 1 and 3 DF,  p-value: 2.457e-07

Save slope

Ag_AgCl_E_T_lm_slope <- coef(Ag_AgCl_E_T_lm_summ)["Ag_AgCl_E_T$T_C",  "Estimate"]

Ag_AgCl_E_T_lm_slope
## [1] -1.046

Save intercept

Ag_AgCl_E_T_lm_int <- coef(Ag_AgCl_E_T_lm_summ)["(Intercept)",  "Estimate"]

Ag_AgCl_E_T_lm_int
## [1] 223.23

Apply conversion

NSHQ14_profile <- NSHQ14_profile %>% mutate(Eh_mV = ORP_mV_Ag_AgCl + Ag_AgCl_E_T_lm_int + Ag_AgCl_E_T_lm_slope * T_C) # convert probe ORP_mV_Ag_AgCl reading from Ag/AgCl to SHE

# print some summary statistics
NSHQ14_profile %>% select(Eh_mV) %>% summary()
##      Eh_mV       
##  Min.   :-721.6  
##  1st Qu.:-718.2  
##  Median :-705.9  
##  Mean   :-667.2  
##  3rd Qu.:-685.2  
##  Max.   : 142.1

Checking values around 20 m and 25 m

NSHQ14_profile %>% filter(depth_sampled_mbgl > 19.99 & depth_sampled_mbgl < 20.01) %>% select(depth_sampled_mbgl, pH, Eh_mV, cond_uS_per_cm)
## # A tibble: 2 × 4
##   depth_sampled_mbgl    pH Eh_mV cond_uS_per_cm
##                <dbl> <dbl> <dbl>          <dbl>
## 1               20.0  10.8 -278.          2409.
## 2               20.0  10.8 -280.          2411.
NSHQ14_profile %>% filter(depth_sampled_mbgl > 24.99 & depth_sampled_mbgl < 25.01) %>% select(depth_sampled_mbgl, pH, Eh_mV, cond_uS_per_cm)
## # A tibble: 2 × 4
##   depth_sampled_mbgl    pH Eh_mV cond_uS_per_cm
##                <dbl> <dbl> <dbl>          <dbl>
## 1               25.0  11.0 -536.          3024.
## 2               25.0  11.0 -537.          3026.

Hydrostatic pressure

Calculate from theory

P_NSHQ14_atmospheric_Pa <- 95490 # measured ambient pressure with Garmin GPSMAP 64S by DBN in 2019
P_NSHQ14_atmospheric_bar <-  P_NSHQ14_atmospheric_Pa * 1e-5

water_level_NSHQ14_m <- as.numeric(NSHQ14_profile %>% select(depth_sampled_mbgl) %>% min())

NSHQ14_profile <- NSHQ14_profile %>% mutate(water_column_above_m = depth_sampled_mbgl - water_level_NSHQ14_m)


# Hydrostatic pressure in a liquid can be calculated as
# 
# p = ρ g h                         (1)
# 
# where
# 
# p = pressure in liquid (N/m2, Pa, lbf/ft2, psf)
# 
# ρ = density of liquid (kg/m3, slugs/ft3)
# 
# g = acceleration of gravity (9.81 m/s2, 32.17405 ft/s2)
# 
# h = height of fluid column - or depth in the fluid where pressure is measured (m, ft)

# https://www.engineeringtoolbox.com/hydrostatic-pressure-water-d_1632.html

density_water <- 1000 # kg/m3
g = 9.81 # acceleration of gravity, m/s2


NSHQ14_profile <- NSHQ14_profile %>% mutate(P_hydrostatic_calc_Pa = P_NSHQ14_atmospheric_Pa + density_water * g * water_column_above_m)

NSHQ14_profile <- NSHQ14_profile %>% mutate(P_hydrostatic_calc_bar = P_hydrostatic_calc_Pa * 1e-5)

# print some summary statistics
NSHQ14_profile %>% select(P_hydrostatic_calc_bar) %>% summary()
##  P_hydrostatic_calc_bar
##  Min.   : 0.9549       
##  1st Qu.: 7.8702       
##  Median :14.7855       
##  Mean   :14.7855       
##  3rd Qu.:21.7008       
##  Max.   :28.6162
# convert measured pressure to bar
NSHQ14_profile <- NSHQ14_profile %>% mutate(P_bar = P_dbar / 10)
# prepare for plotting
NSHQ14_profile_longer_P <- NSHQ14_profile %>% select(depth_sampled_mbgl, P_hydrostatic_calc_bar, P_bar) %>% gather(-depth_sampled_mbgl, key = parameter, value = P_bar) %>% mutate(type = if_else(parameter == "P_hydrostatic_calc_bar", "calculated", "measured"))

Compare theoretical hydrostatic pressure versus measured pressure

plot_NSHQ14_P <- NSHQ14_profile_longer_P %>% ggplot(aes(x = depth_sampled_mbgl))+
  geom_line(aes(y = P_bar, color = type))+
  scale_x_reverse(name = "Depth / [mbgl]", expand = c(0.03,0.02), breaks = c(0, 50, 100, 150, 200, 250, 300))+
    scale_y_continuous(name = latex2exp::TeX("$\\textit{P}_{hydrostatic}\\,/\\,\\lbrack bar \\rbrack$"))+
  coord_flip()+
  theme_bw(base_size = 12)+
  theme(
    panel.grid = element_blank())

plot_NSHQ14_P

They are similar.

BA3A well log

clean up data

# select relevant columns and rename them
BA3A_profile <- BA3A_profile_raw %>% select(depth_sampled_mbgl = `DEPT[M]`, P_dbar = PRESSURE, T_C = TEMPERATURE, pH = PH, ORP_mV_Ag_AgCl = REDOX, cond_uS_per_cm = `COND(FW)`)

# Filter out nonsense values (i.e. when probe is above water table)
BA3A_profile <- BA3A_profile %>% filter(!is.na(ORP_mV_Ag_AgCl))

# print some summary statistics
BA3A_profile %>% summary()
##  depth_sampled_mbgl     P_dbar            T_C              pH       
##  Min.   :  8.293    Min.   :  4.38   Min.   :35.00   Min.   : 9.37  
##  1st Qu.: 79.243    1st Qu.: 75.20   1st Qu.:35.96   1st Qu.:10.62  
##  Median :150.193    Median :145.13   Median :37.65   Median :10.74  
##  Mean   :150.193    Mean   :144.62   Mean   :37.75   Mean   :10.67  
##  3rd Qu.:221.143    3rd Qu.:214.26   3rd Qu.:39.43   3rd Qu.:10.81  
##  Max.   :292.093    Max.   :282.88   Max.   :41.23   Max.   :10.93  
##  ORP_mV_Ag_AgCl   cond_uS_per_cm  
##  Min.   :-894.7   Min.   : 985.9  
##  1st Qu.:-870.5   1st Qu.:1573.4  
##  Median :-857.2   Median :1622.2  
##  Mean   :-855.3   Mean   :1998.7  
##  3rd Qu.:-850.4   3rd Qu.:2382.0  
##  Max.   :-728.3   Max.   :3563.7
# subset the data for readability and print
BA3A_profile %>% slice(seq(from = 1, to = 30000, by = 1000)) %>% kable(digits=1)
depth_sampled_mbgl P_dbar T_C pH ORP_mV_Ag_AgCl cond_uS_per_cm
8.3 4.4 35.0 9.5 -729.5 987.2
18.3 14.6 35.1 9.9 -777.0 1012.9
28.3 24.5 35.1 10.6 -814.0 1237.3
38.3 34.5 35.2 10.8 -833.8 1425.8
48.3 44.4 35.4 10.9 -840.9 1483.1
58.3 54.4 35.5 10.9 -844.1 1538.0
68.3 64.3 35.7 10.9 -848.1 1580.5
78.3 74.2 35.9 10.9 -850.0 1612.0
88.3 84.1 36.2 10.9 -853.7 1592.4
98.3 94.0 36.4 10.8 -852.3 1575.5
108.3 103.9 36.6 10.8 -851.9 1571.2
118.3 113.8 36.9 10.8 -852.7 1579.2
128.3 123.6 37.1 10.8 -854.8 1600.9
138.3 133.5 37.4 10.8 -856.4 1634.5
148.3 143.2 37.6 10.8 -857.1 1616.6
158.3 153.1 37.8 10.7 -856.0 1620.6
168.3 162.8 38.1 10.7 -858.7 1656.1
178.3 172.6 38.4 10.6 -862.9 1815.3
188.3 182.4 38.6 10.6 -863.4 2022.8
198.3 192.1 38.8 10.6 -868.9 2382.3
208.3 201.8 39.1 10.7 -884.4 3205.1
218.3 211.7 39.4 10.7 -886.9 3321.6
228.3 221.3 39.6 10.7 -889.1 3446.7
238.3 231.0 39.9 10.8 -890.8 3496.6
248.3 240.6 40.1 10.7 -883.9 3055.5
258.3 250.4 40.4 10.6 -876.7 2705.6
268.3 259.9 40.6 10.6 -872.1 2402.5
278.3 269.5 40.9 10.5 -866.8 2193.3
288.3 279.2 41.1 10.5 -863.7 2078.1

Convert to ORP to SHE

Apply conversion

BA3A_profile <- BA3A_profile %>% mutate(Eh_mV = ORP_mV_Ag_AgCl + Ag_AgCl_E_T_lm_int + Ag_AgCl_E_T_lm_slope * T_C) # convert probe ORP_mV_Ag_AgCl reading from Ag/AgCl to SHE

# print some summary statistics
BA3A_profile %>% select(Eh_mV) %>% summary()
##      Eh_mV       
##  Min.   :-713.3  
##  1st Qu.:-689.7  
##  Median :-673.3  
##  Mean   :-671.6  
##  3rd Qu.:-664.8  
##  Max.   :-541.7

Hydrostatic pressure

Calculate from theory

P_BA3A_atmospheric_Pa <- 95490 # measured ambient pressure with Garmin GPSMAP 64S by DBN in 2019
P_BA3A_atmospheric_bar <-  P_BA3A_atmospheric_Pa * 1e-5

water_level_BA3A_m <- as.numeric(BA3A_profile %>% select(depth_sampled_mbgl) %>% min())

BA3A_profile <- BA3A_profile %>% mutate(water_column_above_m = depth_sampled_mbgl - water_level_BA3A_m)


# Hydrostatic pressure in a liquid can be calculated as
# 
# p = ρ g h                         (1)
# 
# where
# 
# p = pressure in liquid (N/m2, Pa, lbf/ft2, psf)
# 
# ρ = density of liquid (kg/m3, slugs/ft3)
# 
# g = acceleration of gravity (9.81 m/s2, 32.17405 ft/s2)
# 
# h = height of fluid column - or depth in the fluid where pressure is measured (m, ft)

# https://www.engineeringtoolbox.com/hydrostatic-pressure-water-d_1632.html

density_water <- 1000 # kg/m3
g = 9.81 # acceleration of gravity, m/s2


BA3A_profile <- BA3A_profile %>% mutate(P_hydrostatic_calc_Pa = P_BA3A_atmospheric_Pa + density_water * g * water_column_above_m)

BA3A_profile <- BA3A_profile %>% mutate(P_hydrostatic_calc_bar = P_hydrostatic_calc_Pa * 1e-5)

# print some summary statistics
BA3A_profile %>% select(P_hydrostatic_calc_bar) %>% summary()
##  P_hydrostatic_calc_bar
##  Min.   : 0.9549       
##  1st Qu.: 7.9151       
##  Median :14.8753       
##  Mean   :14.8753       
##  3rd Qu.:21.8355       
##  Max.   :28.7957
# convert measured pressure to bar
BA3A_profile <- BA3A_profile %>% mutate(P_bar = P_dbar / 10)
# prepare for plotting
BA3A_profile_longer_P <- BA3A_profile %>% select(depth_sampled_mbgl, P_hydrostatic_calc_bar, P_bar) %>% gather(-depth_sampled_mbgl, key = parameter, value = P_bar) %>% mutate(type = if_else(parameter == "P_hydrostatic_calc_bar", "calculated", "measured"))

Compare theoretical hydrostatic pressure versus measured pressure

plot_BA3A_P <- BA3A_profile_longer_P %>% ggplot(aes(x = depth_sampled_mbgl))+
  geom_line(aes(y = P_bar, color = type))+
  scale_x_reverse(name = "Depth / [mbgl]", expand = c(0.03,0.02), breaks = c(0, 50, 100, 150, 200, 250, 300))+
    scale_y_continuous(name = latex2exp::TeX("$\\textit{P}_{hydrostatic}\\,/\\,\\lbrack bar \\rbrack$"))+
  coord_flip()+
  theme_bw(base_size = 12)+
  theme(
    panel.grid = element_blank())

plot_BA3A_P

They are similar.

Geotherm

NSHQ14

plot_T_NSHQ14 <- NSHQ14_profile %>% ggplot( mapping = aes(
    x = T_C,
    y = depth_sampled_mbgl)) +
  geom_point() +
    geom_smooth(method="lm", color = "blue") +
  stat_poly_eq(aes(label =  paste(stat(eq.label), stat(rr.label), sep = "~~~~")),
               formula = x ~ y, parse = TRUE, rr.digits = 2, color = "blue", label.x = 0.2)+
  scale_y_reverse(name = latex2exp::TeX("depth$\\,$/$\\,$$\\left[$m$\\right]$"))+
    scale_x_continuous(name = latex2exp::TeX("temperature$\\,$/$\\,$$\\left[$˚C$\\right]$"))+
  theme_bw()+
  theme(
    panel.grid = element_blank()
  )

plot_T_NSHQ14 
## `geom_smooth()` using formula = 'y ~ x'

NSHQ14_profile_below_100m <- NSHQ14_profile %>% filter(depth_sampled_mbgl > 100)

plot_T_NSHQ14_below_100m <- NSHQ14_profile_below_100m %>% ggplot( mapping = aes(
    x = T_C,
    y = depth_sampled_mbgl)) +
  geom_point() +
    geom_smooth(method="lm", color = "blue") +
  stat_poly_eq(aes(label =  paste(stat(eq.label), stat(rr.label), sep = "~~~~")),
               formula = x ~ y, parse = TRUE, rr.digits = 2, color = "blue", label.x = 0.2)+
  scale_y_reverse(name = latex2exp::TeX("depth$\\,$/$\\,$$\\left[$m$\\right]$"))+
    scale_x_continuous(name = latex2exp::TeX("temperature$\\,$/$\\,$$\\left[$˚C$\\right]$"))+
    theme_bw()+
  theme(
    panel.grid = element_blank()
  )

plot_T_NSHQ14_below_100m 
## `geom_smooth()` using formula = 'y ~ x'

# generate linear model
model_T_d_NSHQ14_below_100m = lm(T_C ~ depth_sampled_mbgl,
           data = NSHQ14_profile_below_100m)

# print model summary
summary(model_T_d_NSHQ14_below_100m)
## 
## Call:
## lm(formula = T_C ~ depth_sampled_mbgl, data = NSHQ14_profile_below_100m)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.053155 -0.014659 -0.002923  0.013311  0.060838 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        3.381e+01  5.165e-04   65463   <2e-16 ***
## depth_sampled_mbgl 2.508e-02  2.528e-06    9923   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01959 on 19316 degrees of freedom
## Multiple R-squared:  0.9998, Adjusted R-squared:  0.9998 
## F-statistic: 9.846e+07 on 1 and 19316 DF,  p-value: < 2.2e-16

Predict depth of temperature limit of life

# generate linear model (inverse geothemal gradient - easier to work with for predict function)
model_T_d_NSHQ14_below_100m_inv_lm = lm(depth_sampled_mbgl ~ T_C,
           data = NSHQ14_profile_below_100m)

# print model summary
summary(model_T_d_NSHQ14_below_100m_inv_lm)
## 
## Call:
## lm(formula = depth_sampled_mbgl ~ T_C, data = NSHQ14_profile_below_100m)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.4440 -0.5251  0.1156  0.5814  2.1143 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.348e+03  1.557e-01   -8654   <2e-16 ***
## T_C          3.986e+01  4.017e-03    9923   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.781 on 19316 degrees of freedom
## Multiple R-squared:  0.9998, Adjusted R-squared:  0.9998 
## F-statistic: 9.846e+07 on 1 and 19316 DF,  p-value: < 2.2e-16
# make a data frame for use with predict() function to predict the depth at which the temperature limit of life is reached based on the NSHQ14 geothermal gradient below 100 m
# T limit of life = 122 ˚C (Takai et al., 2008 "Cell proliferation at 122°C and isotopically heavy CH4 production by a hyperthermophilic methanogen under high-pressure cultivation")
T_limit_of_life_C <- data.frame(T_C = c(122))

# predict based on modeled geothermal and print
(model_T_d_NSHQ14_below_100m_inv_lm %>% predict.lm(T_limit_of_life_C))
##        1 
## 3515.643

BA3A

plot_T_BA3A <- BA3A_profile %>% ggplot( mapping = aes(
    x = T_C,
    y = depth_sampled_mbgl)) +
  geom_point() +
    geom_smooth(method="lm", color = "blue") +
  stat_poly_eq(aes(label =  paste(stat(eq.label), stat(rr.label), sep = "~~~~")),
               formula = x ~ y, parse = TRUE, rr.digits = 2, color = "blue", label.x = 0.2)+
  scale_y_reverse(name = latex2exp::TeX("depth$\\,$/$\\,$$\\left[$m$\\right]$"))+
    scale_x_continuous(name = latex2exp::TeX("temperature$\\,$/$\\,$$\\left[$˚C$\\right]$"))+
  theme_bw()+
  theme(
    panel.grid = element_blank()
  )

plot_T_BA3A 
## `geom_smooth()` using formula = 'y ~ x'

BA3A_profile_below_100m <- BA3A_profile %>% filter(depth_sampled_mbgl > 100)

plot_T_BA3A_below_100m <- BA3A_profile_below_100m %>% ggplot( mapping = aes(
    x = T_C,
    y = depth_sampled_mbgl)) +
  geom_point() +
    geom_smooth(method="lm", color = "blue") +
  stat_poly_eq(aes(label =  paste(stat(eq.label), stat(rr.label), sep = "~~~~")),
               formula = x ~ y, parse = TRUE, rr.digits = 2, color = "blue", label.x = 0.2)+
  scale_y_reverse(name = latex2exp::TeX("depth$\\,$/$\\,$$\\left[$m$\\right]$"))+
    scale_x_continuous(name = latex2exp::TeX("temperature$\\,$/$\\,$$\\left[$˚C$\\right]$"))+
    theme_bw()+
  theme(
    panel.grid = element_blank()
  )

plot_T_BA3A_below_100m 
## `geom_smooth()` using formula = 'y ~ x'

# generate linear model
model_T_d_BA3A_below_100m = lm(T_C ~ depth_sampled_mbgl,
           data = BA3A_profile_below_100m)

# print model summary
summary(model_T_d_BA3A_below_100m)
## 
## Call:
## lm(formula = T_C ~ depth_sampled_mbgl, data = BA3A_profile_below_100m)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.031482 -0.007958 -0.003657  0.005497  0.039192 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        3.388e+01  3.134e-04  108111   <2e-16 ***
## depth_sampled_mbgl 2.510e-02  1.538e-06   16320   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01182 on 19208 degrees of freedom
## Multiple R-squared:  0.9999, Adjusted R-squared:  0.9999 
## F-statistic: 2.663e+08 on 1 and 19208 DF,  p-value: < 2.2e-16