#first we want to load in the libraries we will need. We choose seacarb because it is simpler than Phreeqc, because we are not dealing with any ion data. We load readxl to import the data from an excel spreadsheet and use tridyverse as the means for most things for plotting
library(seacarb)
## Loading required package: oce
## Loading required package: gsw
## Loading required package: SolveSAPHE
library(readxl)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.0
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
#this iws what we use to import the data from an excel spreadsheet.
phdata <-read_excel('phexperimentdata.xlsx')
## New names:
## • `` -> `...1`
#now that we loaded the data from the excel spreadsheet we can plot ph from the different conditions over time
ggplot(phdata)+
geom_line(
mapping = aes(x = time_hr, y = ph_HC)
,color = "black"
) +
geom_line(
mapping = aes(x = time_hr, y = ph_HT)
,color = "blue"
) +
geom_line(
mapping = aes(x = time_hr, y = ph_control)
,color = "magenta"
) +
labs(
y= 'pH', x='time (hr)'
) +
#this was just for practive, but here I am adding differnt lines to see what they do.
geom_vline(
xintercept=8.25,
linetype="dashed",
color = "black"
) +
geom_hline(
yintercept=8.25,
linetype="dashed",
color = "coral4"
) +
theme(panel.grid = element_blank()
)
#now we are getting into seacarb. so we are makeing a whole new data sheet based on the parameters that we give seacarb to compute
seacarb_data_HC <- carb(
#the flag tells seacarb which set of carbonate system parameters we're going to give it. refer to page 24 in the seacarb documentation for the full list of options
flag = 21,
#in this case, I chose flag 15, which means I'm going to give seacarb alkalinity and DIC. the units for these are mol/kg, so we need to do a unit conversion:
var1 = 422,
var2 = phdata$ph_HC,
#then we tell it salinity and temperature
S = 35,
T = phdata$temp_HC, # ***add temperature data into excel file, then reload and rerun all code to this poing*** and then repeat for HT and Control.
#surface atmospheric pressure in atm
Patm = 1,
#hydrostatic pressure in bar (surface = 0)
P = 0
)
seacarb_data_HT <- carb(
#the flag tells seacarb which set of carbonate system parameters we're going to give it. refer to page 24 in the seacarb documentation for the full list of options
flag = 21,
#in this case, I chose flag 15, which means I'm going to give seacarb alkalinity and DIC. the units for these are mol/kg, so we need to do a unit conversion:
var1 = 422,
var2 = phdata$ph_HT,
#then we tell it salinity and temperature
S = 35,
T = phdata$temp_HT, # ***add temperature data into excel file, then reload and rerun all code to this poing*** and then repeat for HT and Control.
#surface atmospheric pressure in atm
Patm = 1,
#hydrostatic pressure in bar (surface = 0)
P = 0
)
seacarb_data_control <- carb(
#the flag tells seacarb which set of carbonate system parameters we're going to give it. refer to page 24 in the seacarb documentation for the full list of options
flag = 21,
#in this case, I chose flag 15, which means I'm going to give seacarb alkalinity and DIC. the units for these are mol/kg, so we need to do a unit conversion:
var1 = 422,
var2 = phdata$ph_control,
#then we tell it salinity and temperature
S = 35,
T = phdata$temp_control, # ***add temperature data into excel file, then reload and rerun all code to this poing*** and then repeat for HT and Control.
#surface atmospheric pressure in atm
Patm = 1,
#hydrostatic pressure in bar (surface = 0)
P = 0
)
#now we are plotting the changes in DIC over time form the datafamres that we created in seacarb. We can repeat this for ALk, Ω, or whatever parameter we are interested in.
ggplot(seacarb_data_HC)+
geom_line(
mapping = aes(x = phdata$time_hr, y= DIC*1000)
,color = "purple", ) +
geom_line(
mapping = aes(x = phdata$time_hr, y= seacarb_data_HT$DIC*1000)
,color = "coral1" ) +
geom_line(
mapping = aes(x = phdata$time_hr, y= seacarb_data_control$DIC*1000)
,color = "firebrick1" ) +
labs(y= 'DIC', x='time (hr)'
) +theme(panel.grid = element_blank()
)
ggplot(seacarb_data_HC)+
geom_line(
mapping = aes(x = phdata$time_hr, y= ALK*1000)
, color = "burlywood3", ) +
geom_line(
mapping = aes(x = phdata$time_hr, y= seacarb_data_HT$ALK*1000)
,color = "#96ADFF" ) +
geom_line(
mapping = aes(x = phdata$time_hr, y= seacarb_data_control$ALK*1000)
,color = "gold4" ) +
labs(y= 'Alkalinity', x='time (hr)'
) +
theme(panel.grid = element_blank()
)
ggplot(seacarb_data_HC)+
geom_line(
mapping = aes(x = phdata$time_hr, y=OmegaAragonite)
,color = "burlywood3", ) +
geom_line(
mapping = aes(x = phdata$time_hr, y= seacarb_data_HT$OmegaAragonite)
,color = "#96ADFF" ) +
geom_line(
mapping = aes(x = phdata$time_hr, y= seacarb_data_control$OmegaAragonite)
,color = "gold4" ) +
labs(y= 'Ωarag', x='time (hr)'
)
