diff --git a/RAnalysis/Output/1_Survival_Growth/Seawater_chemistry/F1_Chem_95CI.pdf b/RAnalysis/Output/1_Survival_Growth/Seawater_chemistry/F1_Chem_95CI.pdf
new file mode 100644
index 0000000..7223a15
Binary files /dev/null and b/RAnalysis/Output/1_Survival_Growth/Seawater_chemistry/F1_Chem_95CI.pdf differ
diff --git a/RAnalysis/Output/1_Survival_Growth/Seawater_chemistry/F2_Chem_95CI.pdf b/RAnalysis/Output/1_Survival_Growth/Seawater_chemistry/F2_Chem_95CI.pdf
new file mode 100644
index 0000000..654d7c9
Binary files /dev/null and b/RAnalysis/Output/1_Survival_Growth/Seawater_chemistry/F2_Chem_95CI.pdf differ
diff --git a/RAnalysis/Scripts/Dstage_Abnormality.Rmd b/RAnalysis/Scripts/Dstage_Abnormality.Rmd
index a0720b7..9529ec7 100644
--- a/RAnalysis/Scripts/Dstage_Abnormality.Rmd
+++ b/RAnalysis/Scripts/Dstage_Abnormality.Rmd
@@ -1068,6 +1068,7 @@ capture.output(aovMOD_Exp2_LowvMod,file="Output/2_Dstage_Shell_Abnormalities/Exp
 
 ```{r  Statistics; F2 Experiment #4} 
 
+# ANOVA first
 aovMOD_Exp4 <- lm(Percent.deformed ~ Larvae_pCO2, data=df_Exp4)
 shapiro.test(resid(aovMOD_Exp4)) # 0.2875 - norm
 leveneTest(aovMOD_Exp4) # 0.5552- pass
@@ -1094,6 +1095,25 @@ summary(Exp4_Dstage_BetaReg, type = "quantile")  # same exact outomce without th
 # Larvae_pCO2Moderate pCO2  0.03702    0.41964   0.088     0.93  
 
 
+# T-test 
+
+# normality Shapiro Wilk
+(df_Exp4 %>% 
+        group_by(as.factor(Larvae_pCO2)) %>% 
+        rstatix::shapiro_test(Percent.deformed))
+# Low pCO2	Percent.deformed	0.7714106		# low pCO2 non normal!
+# Moderate pCO2	Percent.deformed	0.9003712	0.38668421		            
+
+# equal variance 
+(df_Exp4 %>% 
+    rstatix::levene_test(Percent.deformed ~ as.factor(Larvae_pCO2)))$p[1] # 0.5552307 pass variance
+
+# we gotta run the Wilkcox! non normal under low pCO2!  
+wilcox.test(df_Exp4$Percent.deformed ~ (as.factor(df_Exp4$Larvae_pCO2)))
+# Wilcoxon rank sum exact test
+# data:  df_Exp4$Percent.deformed by as.factor(df_Exp4$Larvae_pCO2)
+# W = 4, p-value = 1
+?wilcox.test
 # output 
 
 ## anova - no effect
diff --git a/RAnalysis/Scripts/SeawaterChemistry.Rmd b/RAnalysis/Scripts/SeawaterChemistry.Rmd
new file mode 100644
index 0000000..6a0a142
--- /dev/null
+++ b/RAnalysis/Scripts/SeawaterChemistry.Rmd
@@ -0,0 +1,613 @@
+---
+title: "Seawater Chemistry"
+author: "Samuel Gurr"
+date: "8/7/2024"
+output: html_document
+---
+
+### SET UP
+```{r setup, include=FALSE}
+
+knitr::opts_chunk$set(echo = TRUE)
+# SET WORKING DIRECTORY 
+
+#knitr::opts_knit$set(root.dir = "C:/Users/samjg/Documents/Github_repositories/EAD-ASEB-Airradians_Larvae_crossgen_OA/RAnalysis") # Sam's
+knitr::opts_knit$set(root.dir = 
+                       "C:/Users/samjg/Documents/Github_repositories/EAD-ASEB-Airradians_Larvae_crossgen_OA/RAnalysis/") # Sam's work
+
+
+```
+
+
+
+
+
+```{r load data}
+# install.packages('gcookbook')
+# LOAD PACKAGES :::::::::::::::::::::::::::::::::::::::::::::::::::::::
+library(dplyr)
+library(ggplot2)
+library(reshape2)
+library(knitr)
+library(kableExtra)
+library(car)
+library(gcookbook)
+library(ggpubr)
+
+
+
+# LOAD DATA :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+# 8.16/2025 shannon assembled a means and SD dataset for the DIC+TA to calculate 
+CarbChem_Meseck_menas <- read.csv(file="Data/Seawater_chemistry/raw/Seawater_Chemsitry_means_Meseck_20240816.csv", 
+                                  header=TRUE) # load the chem data
+# This is all the data, we can back calculate to get to what Shannon assembled or just use the one above..
+# if we come toa conclusion that TA+pH or pH+DIC is desired than we can do that
+CarbChem_all <- read.csv(file="Data/Seawater_chemistry/raw/Seawater_Chemsitry_raw.csv", 
+                         header=TRUE) # load the chem data
+
+```
+
+
+## Objective 
+
+* For the F1 F2 Larvae paper we want two 8-panel figures (one per generation) containing the following: 
+
+  - Salinity 
+  - temperature bucket
+  - total alkalinity
+  - total dissolved carbon 
+  - pH measured
+  - pH out calculated
+  - pCO2 out
+  - omega aragonite
+
+
+## Master files 
+
+CarbChem_all - contains the chemistry for all larval periods, note the generation and the run #s 
+the manuscript on F1 and F2 larvae used experiment #s 2 and 4
+
+```{r master files}
+
+# Filter for experiment 2 and 4
+CarbChem_Exp2_Exp4 <- CarbChem_all %>% dplyr::filter(Run == c(2,4))
+
+# call columns with TADIC in the column name - these are the CO2sys calculated values 
+CarbChem_Exp2_Exp4_select <- cbind( 
+                                  (CarbChem_Exp2_Exp4[,c(1,3,7:11)]), # columns, core measured data  and factors
+                                  (CarbChem_Exp2_Exp4[ ,grepl( ("TADIC"), names(CarbChem_Exp2_Exp4) ) ]) # TADIC colmns
+                                )
+
+# long format by group 
+CarbChem_Exp2_Exp4_long <- CarbChem_Exp2_Exp4_select %>% 
+                          melt(id.vars=c('Date', 'Generation', 'Larvae_OA')) %>% 
+                          na.omit() %>% 
+                          dplyr::mutate(Date = 
+                                          as.POSIXct(Date, format = "%m/%d/%Y"),
+                                        Month = months(Date),
+                                        Year  = lubridate::year(Date)) 
+
+CarbChem_Exp2_Exp4_long$value <- as.numeric(CarbChem_Exp2_Exp4_long$value) # make numeric
+
+# MEANS PER SAMPLING DAY (reducing replicates within treatment to one mean) 
+# aquire means by group for mean SE plotting by Date and treamtent
+CarbChem_Exp2_Exp4_means.by.day <- CarbChem_Exp2_Exp4_long %>% # calc means and standard error
+                           dplyr::group_by(Date, Generation, Larvae_OA, variable) %>% 
+                           dplyr::summarise(mean = mean(value),
+                                            sd   = sd(value),
+                                            se   = sd/(sqrt(n())), 
+                                            n = n()) %>% 
+                            dplyr::rename(pCO2_treatment = Larvae_OA) %>% 
+                            dplyr::mutate(pCO2_treatment = 
+                                            case_when(
+                                                      pCO2_treatment %in% 'Mid' ~ 'Moderate',
+                                                      TRUE ~ as.factor(pCO2_treatment)
+                                            ),
+                                          
+                            )
+CarbChem_Exp2_Exp4_means.by.day$pCO2_treatment <-
+  forcats::fct_relevel(CarbChem_Exp2_Exp4_means.by.day$pCO2_treatment,c("Low","Moderate","High"))                           
+# Note there are quite a few dates that have only a signle same for that day...
+
+
+```
+
+```{r SmoothRibbon function with gam}
+# how do we get smooth geom ribbon figures
+# in summary run gam and use a fine scale for x axis 
+# https://stackoverflow.com/questions/61029929/smoothing-geom-ribbon
+# A principled way to achieve what you want is to fit a GAM model to your data using the gam() function in mgcv and then apply the predict() function to that model over a finer grid of values for your predictor variable. The grid can cover the span defined by the range of observed values for your predictor variable.
+# load R packages
+library(ggplot2)
+library(mgcv)
+
+
+SmoothRibbon_twoOA <- function(datafilename, col_name) {
+
+data <- as.data.frame(
+  # carbchem.MEANS.by.day %>% 
+                datafilename %>% 
+                na.omit() %>% 
+                # dplyr::filter(!Type %in% 'Larvae') %>carbchem.MASTER_MEANS.MONTLY% 
+                dplyr::filter(pCO2_treatment %in% c("Low", "Moderate")) %>% 
+                dplyr::filter(variable %in% (
+                                               # 'Temperature_bucket_C',
+                                               # 'Salinity',
+                                               # 'DO_mg_l',
+                                               # 'pH',
+                                               # 'pH_out'
+                                               col_name
+                                               # 'TCO2_mmol_kgSW'
+                                               # 'TA_mmol_kgSW',
+                                               # 'pCO2_out_matm'
+                                               #'War_out_DIC_pH'
+                                               # 'Wca_out'
+                                               ))
+                                          ) %>% 
+                 mutate(Date_asnum = dense_rank(Date))
+
+
+# calc the num days because the Date_asnumin ascending order dot not 
+# represent the elapsed time properly as an x axis
+data$NumDays  <- as.numeric(
+                            difftime(as.Date(data$Date),
+                            as.Date(data$Date[1]),
+                            units="days")
+                            )
+
+data.low      <- data %>% dplyr::filter(pCO2_treatment %in% 'Low')
+data.moderate <- data %>% dplyr::filter(pCO2_treatment %in% 'Moderate')
+data.all.ordered  <- rbind(data.low, data.moderate)
+
+
+# fit GAM model
+m.low      <- gam(mean ~ s(NumDays, k = length(table(data.low$mean)) ), data = data.low) 
+m.moderate <- gam(mean ~ s(NumDays, k = length(table(data.moderate$mean)) ), data = data.moderate) 
+
+# define finer grid of predictor values
+xnew.low      <- seq(min(data.low$NumDays), 
+                     max(data.low$NumDays), by = 0.1) 
+xnew.moderate <- seq(min(data.moderate$NumDays), 
+                     max(data.moderate$NumDays), by = 0.1) 
+
+# apply predict() function to the fitted GAM model
+# using the finer grid of x values
+p.low      <- predict(m.low, newdata = data.frame(NumDays = xnew.low), se = TRUE) 
+p.moderate <- predict(m.moderate, newdata = data.frame(NumDays = xnew.moderate), se = TRUE) 
+
+# plot the estimated mean values of y (fit) at given x values
+# over the finer grid of x values;
+# superimpose approximate 95% confidence band for the true 
+# mean values of y at given x values in the finer grid
+g.low <- data.frame(NumDays = xnew.low, 
+            fit = p.low$fit,
+            lwr = p.low$fit - 1.96*p.low$se.fit, # 1.96 * SE = 95% CI
+            upr = p.low$fit + 1.96*p.low$se.fit  # 1.96 * SE = 95% CI
+            ) %>% 
+         dplyr::mutate(pCO2_treatment = "Low")
+
+g.moderate <- data.frame(NumDays = xnew.moderate, 
+            fit = p.moderate$fit,
+            lwr = p.moderate$fit - 1.96*p.moderate$se.fit, 
+            upr = p.moderate$fit + 1.96*p.moderate$se.fit) %>% 
+         dplyr::mutate(pCO2_treatment = "Moderate")
+
+# r bind the data sets
+# use the first date of data and add the elapsed numberic time, not that 1 = seconds, we 
+# convert to days by mutliplying by 86400 (seconds in a day)
+g.master <- rbind(g.low, g.moderate) %>% 
+                  dplyr::mutate(Date = data$Date[1] + NumDays*86400)
+
+theme_set(theme_classic())
+
+p <- ggplot(data = g.master, aes(Date, fit)) + 
+  geom_ribbon(aes(ymin = lwr, ymax = upr, fill=factor(pCO2_treatment)), alpha = 0.2) + 
+  geom_line(aes(y = lwr,
+                colour=factor(pCO2_treatment)),
+            linetype = "solid") +
+  geom_line(aes(y = upr,
+                colour=factor(pCO2_treatment)),
+            linetype = "solid") +
+  geom_point(data = data.all.ordered, 
+             aes(Date, mean,
+                 colour=factor(pCO2_treatment)), 
+             shape = 3, size =1) +
+  ylab(col_name) +
+  theme(legend.position="none") +
+  scale_fill_manual(values=c("forestgreen", "orange"), name="fill") +
+  scale_colour_manual(values=c("forestgreen", "orange"), name="fill") 
+
+print(p)
+}
+
+
+
+
+SmoothRibbon_threeOA <- function(datafilename, col_name) {
+
+data <- as.data.frame(
+  #carbchem.MEANS.by.day.F1 %>% 
+                datafilename %>% 
+                na.omit() %>% 
+                # dplyr::filter(!Type %in% 'Larvae') %>carbchem.MASTER_MEANS.MONTLY% 
+                # dplyr::filter(pCO2_treatment %in% c("Low", "Moderate", "High")) %>% 
+                dplyr::filter(variable %in% (
+                                               #'TADIC_WAr_out_in_situ'
+                                               # 'Salinity',
+                                               # 'DO_mg_l',
+                                               # 'pH',
+                                               # 'pH_out'
+                                               col_name
+                                               # 'TCO2_mmol_kgSW'
+                                               # 'TA_mmol_kgSW',
+                                               # 'pCO2_out_matm'
+                                               #'War_out_DIC_pH'
+                                               # 'Wca_out'
+                                               ))
+                                          ) %>% 
+                 mutate(Date_asnum = dense_rank(Date))
+
+
+# calc the num days because the Date_asnumin ascending order dot not 
+# represent the elapsed time properly as an x axis
+data$NumDays  <- as.numeric(
+                            difftime(as.Date(data$Date),
+                            as.Date(data$Date[1]),
+                            units="days")
+                            )
+
+data.low      <- data %>% dplyr::filter(pCO2_treatment %in% 'Low')
+data.moderate <- data %>% dplyr::filter(pCO2_treatment %in% 'Moderate')
+data.high     <- data %>% dplyr::filter(pCO2_treatment %in% 'High')
+data.all.ordered  <- rbind(data.low, data.moderate, data.high)
+
+
+# fit GAM model
+m.low      <- gam(mean ~ s(NumDays, k = length(table(data.low$mean)) ), data = data.low) 
+m.moderate <- gam(mean ~ s(NumDays, k = length(table(data.moderate$mean)) ), data = data.moderate) 
+m.high     <- gam(mean ~ s(NumDays, k = length(table(data.high$mean)) ), data = data.high) 
+
+# define finer grid of predictor values
+xnew.low      <- seq(min(data.low$NumDays), 
+                     max(data.low$NumDays), by = 0.1) 
+xnew.moderate <- seq(min(data.moderate$NumDays), 
+                     max(data.moderate$NumDays), by = 0.1) 
+xnew.high     <- seq(min(data.high$NumDays), 
+                     max(data.high$NumDays), by = 0.1) 
+# apply predict() function to the fitted GAM model
+# using the finer grid of x values
+p.low      <- predict(m.low, newdata = data.frame(NumDays = xnew.low), se = TRUE) 
+p.moderate <- predict(m.moderate, newdata = data.frame(NumDays = xnew.moderate), se = TRUE) 
+p.high     <- predict(m.high, newdata = data.frame(NumDays = xnew.high), se = TRUE) 
+
+# plot the estimated mean values of y (fit) at given x values
+# over the finer grid of x values;
+# superimpose approximate 95% confidence band for the true 
+# mean values of y at given x values in the finer grid
+g.low <- data.frame(NumDays = xnew.low, 
+            fit = p.low$fit,
+            lwr = p.low$fit - 1.96*p.low$se.fit, # 1.96 * SE = 95% CI
+            upr = p.low$fit + 1.96*p.low$se.fit  # 1.96 * SE = 95% CI
+            ) %>% 
+         dplyr::mutate(pCO2_treatment = "Low")
+
+g.moderate <- data.frame(NumDays = xnew.moderate, 
+            fit = p.moderate$fit,
+            lwr = p.moderate$fit - 1.96*p.moderate$se.fit, 
+            upr = p.moderate$fit + 1.96*p.moderate$se.fit) %>% 
+         dplyr::mutate(pCO2_treatment = "Moderate")
+
+g.high <- data.frame(NumDays = xnew.high, 
+            fit = p.high$fit,
+            lwr = p.high$fit - 1.96*p.high$se.fit, 
+            upr = p.high$fit + 1.96*p.high$se.fit) %>% 
+         dplyr::mutate(pCO2_treatment = "High")
+# r bind the data sets
+# use the first date of data and add the elapsed numberic time, not that 1 = seconds, we 
+# convert to days by mutliplying by 86400 (seconds in a day)
+g.master <- rbind(g.low, g.moderate, g.high) %>% 
+                  dplyr::mutate(Date = data$Date[1] + NumDays*86400)
+
+theme_set(theme_classic())
+
+levels(g.master$pCO2_treatment) = c("Low", "Moderate", "High")
+
+p <- ggplot(data = g.master, aes(Date, fit)) + 
+  geom_ribbon(aes(ymin = lwr, ymax = upr, 
+                  fill=(forcats::fct_relevel(pCO2_treatment,c("Low", "Moderate", "High")))), 
+                  alpha = 0.2) + 
+  geom_line(aes(y = lwr,
+                colour=(pCO2_treatment)),
+            linetype = "solid") +
+  geom_line(aes(y = upr,
+                colour=(pCO2_treatment)),
+            linetype = "solid") +
+  geom_point(data = data.all.ordered, 
+             aes(Date, mean,
+                 colour=(pCO2_treatment)), 
+             shape = 3, size =1) +
+  ylab(col_name) +
+  theme(legend.position="none") +
+  scale_fill_manual(values=c("forestgreen", "orange","purple"), name="fill") +
+  scale_colour_manual(values=c("forestgreen","orange", "purple"), name="fill") 
+
+print(p)
+}
+```
+
+* USe Smooth Ribbon and GAM to make some narly plots!
+```{r F1 three OA treamtents - Nice plots!}
+library(ggpubr)
+carbchem.MEANS.by.day.F1 <- CarbChem_Exp2_Exp4_means.by.day %>% filter(Generation %in% 'F1') 
+# pdf("C:/Users/samjg/Documents/Github_repositories/Airradians_multigen_OA/RAnalysis/Output/SeawaterChemistry/F1/plots/F1_Chem_95CI.pdf", width=8, height=11)
+pdf("Output/1_Survival_Growth/Seawater_chemistry/F1_Chem_95CI.pdf", width=8, height=11)
+    ggarrange(
+          SmoothRibbon_threeOA(carbchem.MEANS.by.day.F1, "Temperature_bucket_C"),
+          SmoothRibbon_threeOA(carbchem.MEANS.by.day.F1, "Salinity"),
+          SmoothRibbon_threeOA(carbchem.MEANS.by.day.F1, "TA_mmol_kgSW"),
+          SmoothRibbon_threeOA(carbchem.MEANS.by.day.F1, "TCO2_mmol_kgSW"),
+          # SmoothRibbon_threeOA(carbchem.MEANS.by.day.F1, ""),
+          SmoothRibbon_threeOA(carbchem.MEANS.by.day.F1, "TADIC_pH_out_in_situ"),          
+          SmoothRibbon_threeOA(carbchem.MEANS.by.day.F1, "TADIC_pCO2_out_.matm._in_situ"),          
+          SmoothRibbon_threeOA(carbchem.MEANS.by.day.F1, "TADIC_WAr_out_in_situ"),          
+          nrow = 4,
+          ncol = 2
+    )
+graphics.off()
+
+
+```
+
+
+```{r F2 two OA treamtents - Nice plots!}
+
+carbchem.MEANS.by.day.F2 <- CarbChem_Exp2_Exp4_means.by.day %>% filter(Generation %in% 'F2')
+
+# pdf("C:/Users/samjg/Documents/Github_repositories/Airradians_multigen_OA/RAnalysis/Output/SeawaterChemistry/F2/plots/F2_Chem_95CI.pdf", width=8, height=11)
+pdf("Output/1_Survival_Growth/Seawater_chemistry/F2_Chem_95CI.pdf", width=8, height=11)
+    ggarrange(
+          SmoothRibbon_twoOA(carbchem.MEANS.by.day.F1, "Temperature_bucket_C"),
+          SmoothRibbon_twoOA(carbchem.MEANS.by.day.F1, "Salinity"),
+          SmoothRibbon_twoOA(carbchem.MEANS.by.day.F1, "TA_mmol_kgSW"),
+          SmoothRibbon_twoOA(carbchem.MEANS.by.day.F1, "TCO2_mmol_kgSW"),
+          # SmoothRibbon_twoOA(carbchem.MEANS.by.day.F1, ""),
+          SmoothRibbon_twoOA(carbchem.MEANS.by.day.F1, "TADIC_pH_out_in_situ"),          
+          SmoothRibbon_twoOA(carbchem.MEANS.by.day.F1, "TADIC_pCO2_out_.matm._in_situ"),          
+          SmoothRibbon_twoOA(carbchem.MEANS.by.day.F1, "TADIC_WAr_out_in_situ"), 
+          nrow = 4,
+          ncol = 2
+    )
+graphics.off()
+
+```
+
+
+
+# TABLE (all data)
+
+```{r use dplyr to build chemistry table, include=FALSE}
+carbchem.MASTER_MEANS <- carbchem.MASTER_long %>% # calc means and standard error
+   na.omit() %>% 
+   dplyr::group_by(Date, Type, pCO2_treatment, variable) %>% 
+   dplyr::summarise(mean = mean(value),
+                    sd   = sd(value),
+                    se   = sd/(sqrt(n())) )
+
+chem_meanDate <- carbchem.MASTER_MEANS %>% 
+     na.omit() %>% 
+   dplyr::group_by(pCO2_treatment, variable) %>% 
+   dplyr::summarise(mean2 = mean(mean),
+                    sd   = sd(mean),
+                    se   = sd/(sqrt(n())) ,
+                    n = n())
+
+carbchem.MASTER_wideMEANS <- dcast(chem_meanDate, pCO2_treatment ~ variable, value.var="mean2")
+carbchem.MASTER_wideStErr <- dcast(chem_meanDate, pCO2_treatment ~ variable, value.var="se")
+carbchem.MASTER_wideN     <- dcast(chem_meanDate, pCO2_treatment ~ variable, value.var="n")
+
+# final table
+
+colnames(carbchem.MASTER_wideMEANS)
+FINAL_TABLE                <- data.frame(matrix(nrow = nrow(carbchem.MASTER_wideMEANS), ncol = 1))
+
+FINAL_TABLE$pCO2_Treatment <- carbchem.MASTER_wideMEANS$pCO2_treatment
+
+FINAL_TABLE$Salinity       <-  paste( (paste( (signif(carbchem.MASTER_wideMEANS$Salinity, digits=3)),
+                                     (signif(carbchem.MASTER_wideStErr$Salinity, digits=3)), sep=" ± ")),
+                                     " (",
+                                     carbchem.MASTER_wideN$Salinity,
+                                     ")",
+                                     sep = "")
+
+FINAL_TABLE$pCO2           <-  paste( (paste( (signif(carbchem.MASTER_wideMEANS$xCO2_out._dry_at._1_atm_.ppm_.DIC_pH, digits=3)),
+                                     (signif(carbchem.MASTER_wideStErr$xCO2_out._dry_at._1_atm_.ppm_.DIC_pH, digits=3)), sep=" ± ")),
+                                     " (",
+                                     carbchem.MASTER_wideN$xCO2_out._dry_at._1_atm_.ppm_.DIC_pH,
+                                     ")",
+                                     sep = "")
+
+FINAL_TABLE$Temperature    <- paste( (paste( (signif(carbchem.MASTER_wideMEANS$Temperature_bucket_C, digits=3)),
+                                     (signif(carbchem.MASTER_wideStErr$Temperature_bucket_C, digits=3)), sep=" ± ")),
+                                     " (",
+                                     carbchem.MASTER_wideN$Temperature_bucket_C,
+                                     ")",
+                                     sep = "")
+
+FINAL_TABLE$pH             <- paste( (paste( (signif(carbchem.MASTER_wideMEANS$pH_out, digits=3)),
+                                     (signif(carbchem.MASTER_wideStErr$pH_out, digits=3)), sep=" ± ")),
+                                     " (",
+                                     carbchem.MASTER_wideN$pH_out,
+                                     ")",
+                                     sep = "")
+
+FINAL_TABLE$HCO3           <-  paste( (paste( (signif(carbchem.MASTER_wideMEANS$HCO3_out_mmol_kgSW, digits=3)),
+                                     (signif(carbchem.MASTER_wideStErr$HCO3_out_mmol_kgSW, digits=3)), sep=" ± ")),
+                                     " (",
+                                     carbchem.MASTER_wideN$HCO3_out_mmol_kgSW,
+                                     ")",
+                                     sep = "")
+
+FINAL_TABLE$CO3           <-  paste( (paste( (signif(carbchem.MASTER_wideMEANS$CO3_out_mmol_kgSW, digits=3)),
+                                             "±",
+                                     (signif(carbchem.MASTER_wideStErr$CO3_out_mmol_kgSW, digits=3)), sep="  ")),
+                                     " (",
+                                     carbchem.MASTER_wideN$CO3_out_mmol_kgSW,
+                                     ")",
+                                     sep = "")
+
+FINAL_TABLE$TA             <- paste( (paste( (signif(carbchem.MASTER_wideMEANS$TA_mmol_kgSW, digits=3)),
+                                     (signif(carbchem.MASTER_wideStErr$TA_mmol_kgSW, digits=3)), sep=" ± ")),
+                                     " (",
+                                     carbchem.MASTER_wideN$TA_mmol_kgSW,
+                                     ")",
+                                     sep = "")
+
+FINAL_TABLE$Aragonite.Sat  <-  paste( (paste( (signif(carbchem.MASTER_wideMEANS$War_out_DIC_pH, digits=3)),
+                                     (signif(carbchem.MASTER_wideStErr$War_out_DIC_pH, digits=3)), sep=" ± ")),
+                                     " (",
+                                     carbchem.MASTER_wideN$War_out_DIC_pH,
+                                     ")",
+                                     sep = "")
+
+FINAL_TABLE                <- FINAL_TABLE[,-1] # view table
+```
+
+
+## knitR to print the table 
+```{r lable output, echo=TRUE}
+FINAL_TABLE %>%
+  kbl(caption = "Table 1. Seawater chemistry") %>%
+  kable_classic(full_width = T, html_font = "Cambria")
+```
+
+
+#### save data
+```{r save tables, echo=TRUE}
+# save output table
+# write.table(chem_meanDate,"C:/Users/samjg/Documents/Github_repositories/Cvirginica_multistressor/RAnalysis/Output/Water_Chem/Chem.Table_sep.csv",sep=",", row.names=FALSE)  # write table to 
+# write.table(FINAL_TABLE,"C:/Users/samjg/Documents/Github_repositories/Cvirginica_multistressor/RAnalysis/Output/Water_Chem/Chem.Table.csv",sep=",", row.names=FALSE)  # write table to 
+
+```
+
+
+
+
+
+
+
+
+
+
+# TABLE (all data)
+
+```{r use dplyr to build chemistry table, include=FALSE}
+
+# change to long format to properly facet
+carbchem.MASTER_long <- carbchem.MASTER %>% 
+  dplyr::select(!c(Treatment,Replicate)) %>% 
+  melt(id.vars=c('Date', 'Type', 'pCO2_treatment'))
+
+carbchem.MASTER_long$value <- as.numeric(carbchem.MASTER_long$value) # make numeric
+
+carbchem.MASTER_MEANS <- carbchem.MASTER_long %>% # calc means and standard error
+   na.omit() %>% 
+   dplyr::group_by(Date, Type, pCO2_treatment, variable) %>% 
+   dplyr::summarise(mean = mean(value),
+                    sd   = sd(value),
+                    se   = sd/(sqrt(n())) )
+
+chem_meanDate <- carbchem.MASTER_MEANS %>% 
+     na.omit() %>% 
+   dplyr::group_by(pCO2_treatment, variable) %>% 
+   dplyr::summarise(mean2 = mean(mean),
+                    sd   = sd(mean),
+                    se   = sd/(sqrt(n())) ,
+                    n = n())
+
+carbchem.MASTER_wideMEANS <- dcast(chem_meanDate, pCO2_treatment ~ variable, value.var="mean2")
+carbchem.MASTER_wideStErr <- dcast(chem_meanDate, pCO2_treatment ~ variable, value.var="se")
+carbchem.MASTER_wideN     <- dcast(chem_meanDate, pCO2_treatment ~ variable, value.var="n")
+
+# final table
+
+colnames(carbchem.MASTER_wideMEANS)
+FINAL_TABLE                <- data.frame(matrix(nrow = nrow(carbchem.MASTER_wideMEANS), ncol = 1))
+
+FINAL_TABLE$pCO2_Treatment <- carbchem.MASTER_wideMEANS$pCO2_treatment
+
+FINAL_TABLE$Salinity       <-  paste( (paste( (signif(carbchem.MASTER_wideMEANS$Salinity, digits=3)),
+                                     (signif(carbchem.MASTER_wideStErr$Salinity, digits=3)), sep=" ± ")),
+                                     " (",
+                                     carbchem.MASTER_wideN$Salinity,
+                                     ")",
+                                     sep = "")
+
+FINAL_TABLE$pCO2           <-  paste( (paste( (signif(carbchem.MASTER_wideMEANS$xCO2_out._dry_at._1_atm_.ppm_.DIC_pH, digits=3)),
+                                     (signif(carbchem.MASTER_wideStErr$xCO2_out._dry_at._1_atm_.ppm_.DIC_pH, digits=3)), sep=" ± ")),
+                                     " (",
+                                     carbchem.MASTER_wideN$xCO2_out._dry_at._1_atm_.ppm_.DIC_pH,
+                                     ")",
+                                     sep = "")
+
+FINAL_TABLE$Temperature    <- paste( (paste( (signif(carbchem.MASTER_wideMEANS$Temperature_bucket_C, digits=3)),
+                                     (signif(carbchem.MASTER_wideStErr$Temperature_bucket_C, digits=3)), sep=" ± ")),
+                                     " (",
+                                     carbchem.MASTER_wideN$Temperature_bucket_C,
+                                     ")",
+                                     sep = "")
+
+FINAL_TABLE$pH             <- paste( (paste( (signif(carbchem.MASTER_wideMEANS$pH_out, digits=3)),
+                                     (signif(carbchem.MASTER_wideStErr$pH_out, digits=3)), sep=" ± ")),
+                                     " (",
+                                     carbchem.MASTER_wideN$pH_out,
+                                     ")",
+                                     sep = "")
+
+FINAL_TABLE$HCO3           <-  paste( (paste( (signif(carbchem.MASTER_wideMEANS$HCO3_out_mmol_kgSW, digits=3)),
+                                     (signif(carbchem.MASTER_wideStErr$HCO3_out_mmol_kgSW, digits=3)), sep=" ± ")),
+                                     " (",
+                                     carbchem.MASTER_wideN$HCO3_out_mmol_kgSW,
+                                     ")",
+                                     sep = "")
+
+FINAL_TABLE$CO3           <-  paste( (paste( (signif(carbchem.MASTER_wideMEANS$CO3_out_mmol_kgSW, digits=3)),
+                                             "±",
+                                     (signif(carbchem.MASTER_wideStErr$CO3_out_mmol_kgSW, digits=3)), sep="  ")),
+                                     " (",
+                                     carbchem.MASTER_wideN$CO3_out_mmol_kgSW,
+                                     ")",
+                                     sep = "")
+
+FINAL_TABLE$TA             <- paste( (paste( (signif(carbchem.MASTER_wideMEANS$TA_mmol_kgSW, digits=3)),
+                                     (signif(carbchem.MASTER_wideStErr$TA_mmol_kgSW, digits=3)), sep=" ± ")),
+                                     " (",
+                                     carbchem.MASTER_wideN$TA_mmol_kgSW,
+                                     ")",
+                                     sep = "")
+
+FINAL_TABLE$Aragonite.Sat  <-  paste( (paste( (signif(carbchem.MASTER_wideMEANS$War_out_DIC_pH, digits=3)),
+                                     (signif(carbchem.MASTER_wideStErr$War_out_DIC_pH, digits=3)), sep=" ± ")),
+                                     " (",
+                                     carbchem.MASTER_wideN$War_out_DIC_pH,
+                                     ")",
+                                     sep = "")
+
+FINAL_TABLE                <- FINAL_TABLE[,-1] # view table
+```
+
+
+## knitR to print the table 
+```{r lable output, echo=TRUE}
+FINAL_TABLE %>%
+  kbl(caption = "Table 1. Seawater chemistry") %>%
+  kable_classic(full_width = T, html_font = "Cambria")
+```
+
+
+#### save data (food x OA manuscript)
+```{r save tables, echo=TRUE}
+# save output table
+# write.table(chem_meanDate,"C:/Users/samjg/Documents/Github_repositories/Cvirginica_multistressor/RAnalysis/Output/Water_Chem/Chem.Table_sep.csv",sep=",", row.names=FALSE)  # write table to 
+# write.table(FINAL_TABLE,"C:/Users/samjg/Documents/Github_repositories/Cvirginica_multistressor/RAnalysis/Output/Water_Chem/Chem.Table.csv",sep=",", row.names=FALSE)  # write table to 
+
+```