These code blocks were used to calculate the aragonite sea intensity (ASI) time series
Mg/Ca ratio and temperature influence on aragonite and calcite precipitation based on experiments from: Balthasar, U. & Cusack, M. Aragonite-calcite seas - Quantifying the gray area. Geology 43, 99-102 (2015). A table with extracted values from the original experimental results was used (“balthasar2015table.csv”)
setwd("C://Users/keichenseer/OneDrive - University of Plymouth/Manuscripts/Nature_Geoscience/Figshare_Uploads/CSVs")
#setwd("D://OneDrive/OneDrive - University of Plymouth/Manuscripts/Nature_Geoscience/Figshare_Uploads/CSVs")# set your appropriate working directory, e.g. ("C://Users/keichenseer/OneDrive - University of Plymouth/Manuscripts/Nature_Geoscience/Figshare_Uploads/CSVs") which contains all the CSV files from this Figshare project
# table with aragonite proportion in each experiment
# load libraries
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
#library(ggplot2)
stagetab <- readRDS("stagetable.rds") # read geol stages
bct <- read.csv(file="balthasar2015table.csv", header=TRUE, sep=",") # read the experimental data
# adjust for the different molar weights of calcite and aragonite
molea <- function(a) (2.95*a)/(2.95*a+2.71*(100-a))
lm(formula =100*molea(aragonite) ~ mgca + temperature, data = bct)
##
## Call:
## lm(formula = 100 * molea(aragonite) ~ mgca + temperature, data = bct)
##
## Coefficients:
## (Intercept) mgca temperature
## -119.610 46.565 4.297
The percentage of aragonite as a function of Mg/Ca ratio and temperature (T) thus is -119.61 + 46.57xMg/Ca + 4.30xT
# a function to calculate aragonite proportion from Mg/Ca and T:
get_aragonite <- function(t,m) -119.61 + 46.57*m + 4.30*t
Palaeotemperatures were calculated from the shallow water dO18 record from the supplementary material of: Veizer, J. & Prokoph, A. Temperatures and oxygen isotopic composition of Phanerozoic oceans. Earth-Science Reviews 146, 92-104 (2015).
# read in the Veizer dO18 data
library(readxl) # create a function that reads in multiple excel spread sheets
read_excel_allsheets <- function(filename) {
sheets <- readxl::excel_sheets(filename)
x <- lapply(sheets, function(X) readxl::read_excel(filename, sheet = X))
names(x) <- sheets
x
}
surf <- read_excel_allsheets("surface_t_veizer_2015.xls") # read in the data
# View(surf[[2]][,2]) # see if it worked
coln <- unique(unlist(sapply(surf,function(x)(colnames(x)))))[-1]
tsurf <- NULL
for(i in 2: length(surf)) {
temp <- matrix(NA, nrow = nrow(surf[[i]]), ncol = length(coln))
for(j in 1:length(coln)) if (coln[j] %in% colnames(surf[[i]])) temp[,j] <- unlist(surf[[i]][,which(colnames(surf[[i]])==coln[j])])
tsurf <- rbind(tsurf,temp)
}
colnames(tsurf) <- coln
tsurf <- as.data.frame(tsurf)
climate_def <- read.csv("veizer2015_surface_climates.csv") # assign latitudinal categories to the region description of the original data
latitude <- sapply(as.character(tsurf$climate), function(x) climate_def$latitude[which(climate_def$climate == x)])
tsurf <- cbind(tsurf, latitude)
tsurf$gts2012 <- as.numeric(as.character(tsurf$gts2012))
#hist(tsurf$gts2012)
colnames(tsurf)[4] <- "d18O"
tsurf$d18O <- as.numeric(as.character(tsurf$d18O))
tsurf <- tsurf[which(!(is.na(tsurf$gts2012))),]
tsurf$fossil <- as.character(tsurf$fossil) # rename the fossils to sensible groups
tsurf$fossil[which(tsurf$fossil %in% c("belemnite", "Belemnite", "belemnites", "Belemnites"))] <- "belemnite"
tsurf$fossil[which(tsurf$fossil %in% c("brachiopod", "Brachiopod", "Brachiopods"))] <- "brachiopod"
tsurf$fossil[which(tsurf$fossil %in% c("PlankticF", "Plankticf"))] <- "plankticF"
tsurf$fossil[which(tsurf$fossil %in% c("bivalve", "Bivalve", "Bivalve-Bent", "Bivalve-Trigonia-", "Gryphea",
"Inoc", "inoceramid", "Inoceramid", "Ostrea", "oyster", "Oyster", "Oysters", "PL-bivalve"))] <- "bivalve"
# Add data from Alberti et al. (2017)
alball <- read.csv("alberti_2017_isotopes.csv")
alb <- alball %>% filter(reference == "present study")
alb$stage <- as.character(alb$stage)
alb <- alb %>% mutate (stagemid = stagetab$stage_midp[sapply(stage, function(x) which(x == stagetab$sstages2))])
alberti <- data.frame(matrix(nrow = nrow(alb), ncol = ncol(tsurf)))
names(alberti) <- names(tsurf)
alberti$d18O <- alb$d18O
alberti$gts2012 <- alb$stagemid
alberti$latitude <- rep("low",nrow(alberti))
alberti$fossil[which(alb$taxonomy == "oyster")] <- "bivalve"
alberti$fossil[which(alb$taxonomy == "brachiopod")] <- "brachiopod"
# View(alberti)
tsurf <- rbind(tsurf, alberti)
# Temperature calculation based on detrending the Phanerozoic d18O trend, Veizer & Prokoph 2015
T_calculate <- function(d,t) {
if(!(identical(d,numeric(0))) & !(identical(t, numeric(0)))) 16.9 - 4*(d- (-0.00003*t^2 + 0.0046*t) -0.27)
else NA
}
temperature <- rep(NA,nrow(tsurf))
temperature <- T_calculate(d=tsurf$d18O, t = tsurf$gts2012)
tsurf <- cbind(tsurf, temperature)
tsurf$temperature <- as.numeric(tsurf$temperature)
Mg/Ca ratios were read in steps of 2 million years from Fig 2A in Demicco, R. V., Lowenstein, T. K., Hardie, L. A. & Spencer, R. J. Model of seawater composition for the Phanerozoic. Geology 33, 877-880 (2005).
demicco2005 <- read.csv("Mg_Ca_Demicco_2005.csv")
mcd <- as.data.frame(matrix(NA,nrow(demicco2005),10))
mcd[,c(1,4)] <- demicco2005[,c(1,2)]
# calculate mgca for the slices:
mgcastage <- as.data.frame(matrix(NA,85,4))
#stage_midp <- apply(cbind(stage_base,stage_top),1,mean)
stagetab <- readRDS("stagetable.rds")
sstages2 <- stagetab$sstages2
stage_base <- stagetab$stage_base
stage_top <- stagetab$stage_top
stage_midp <- (stage_base+stage_top)/2
allstages <- as.data.frame(cbind(sstages2[1:85],stage_midp,stage_base,stage_top))
allstages[,2:4] <- apply(allstages[,2:4],2,as.numeric)
colnames(mgcastage) <- c("time", "slice_base" , "slice_top", "MgCa")
mgcastage[,1] <- allstages[,2]
mgcastage[,2] <- allstages[,3]
mgcastage[,3] <- allstages[,4]
for (j in 1:85){
mgcastage[j,4] <- mean(mcd[which(mcd[,1] < mgcastage[j,2] & mcd[,1] >= mgcastage[j,3]),4], na.rm = T)
if (j == 84) mgcastage[j,4] <- mean(mcd[which(mcd[,1] < mgcastage[j-1,2] & mcd[,1] >= mgcastage[j+1,3]),4], na.rm = T)
}
plot(-stage_midp,mgcastage$MgCa,type = "l", ylab = "Mg/Ca", xlab = "Time (Ma)")
tpara <- array(NA,dim = c(85,13,4)) # first dimension for the values, second for SE, third for SD, fourth for N observations
#tsurf_all <- tsurf
tlow_list <- list(NULL)
ASI_list <- list(NULL)
tsurf <- tsurf[which(tsurf$fossil %in% c("brachiopod", "plankticF", "bivalve")),] # use only brachiopods, planktic Foraminifera and bivalve d18O values
tsurf <- tsurf%>% filter(!(is.na(gts2012)) & !(is.na(d18O)))
colnames(allstages) <- c("name", "midp", "base", "top")
for (i in 1:nrow(allstages)) {
tpara[i,1,1] <- mean(T_calculate(tsurf$d18O[which(tsurf$gts2012 < allstages$base[i] & tsurf$gts2012 >= allstages$top[i])],
tsurf$gts2012[which(tsurf$gts2012 < allstages$base[i] & tsurf$gts2012 >= allstages$top[i])]), na.rm = T)
tpara[i,2,1] <- mean(T_calculate(tsurf$d18O[which(tsurf$latitude == "low" & tsurf$gts2012 < allstages$base[i] & tsurf$gts2012 >= allstages$top[i])],
tsurf$gts2012[which(tsurf$latitude == "low" & tsurf$gts2012 < allstages$base[i] & tsurf$gts2012 >= allstages$top[i])]), na.rm = T)
tpara[i,3,1] <- mean(T_calculate(tsurf$d18O[which(tsurf$latitude == "medium" & tsurf$gts2012 < allstages$base[i] & tsurf$gts2012 >= allstages$top[i])],
tsurf$gts2012[which(tsurf$latitude == "medium" & tsurf$gts2012 < allstages$base[i] & tsurf$gts2012 >= allstages$top[i])]), na.rm = T)
tpara[i,4,1] <- mean(T_calculate(tsurf$d18O[which(tsurf$latitude == "high" & tsurf$gts2012 < allstages$base[i] & tsurf$gts2012 >= allstages$top[i])],
tsurf$gts2012[which(tsurf$latitude == "high" & tsurf$gts2012 < allstages$base[i] & tsurf$gts2012 >= allstages$top[i])]), na.rm = T)
tpara[i,10,1] <- length(tsurf$temperature[which(tsurf$latitude == "low" & tsurf$gts2012 < allstages$base[i] & tsurf$gts2012 >= allstages$top[i])])
tpara[i,11,1] <- length(tsurf$temperature[which(tsurf$latitude == "medium" & tsurf$gts2012 < allstages$base[i] & tsurf$gts2012 >= allstages$top[i])])
tpara[i,1,2] <- sqrt(var(T_calculate(tsurf$d18O[which(tsurf$gts2012 < allstages$base[i] & tsurf$gts2012 >= allstages$top[i])],
tsurf$gts2012[which(tsurf$gts2012 < allstages$base[i] & tsurf$gts2012 >= allstages$top[i])]), na.rm = T)/sqrt(length(tsurf$temperature[which(tsurf$gts2012 < allstages$base[i] & tsurf$gts2012 >= allstages$top[i])]))
)
tpara[i,2,2] <- sqrt(var(T_calculate(tsurf$d18O[which(tsurf$latitude == "low" & tsurf$gts2012 < allstages$base[i] & tsurf$gts2012 >= allstages$top[i])],
tsurf$gts2012[which(tsurf$latitude == "low" & tsurf$gts2012 < allstages$base[i] & tsurf$gts2012 >= allstages$top[i])]), na.rm = T)/sqrt(length(tsurf$temperature[which(tsurf$latitude == "low" & tsurf$gts2012 < allstages$base[i] & tsurf$gts2012 >= allstages$top[i])]))
)
tpara[i,3,2] <- sqrt(var(T_calculate(tsurf$d18O[which(tsurf$latitude == "medium" & tsurf$gts2012 < allstages$base[i] & tsurf$gts2012 >= allstages$top[i])],
tsurf$gts2012[which(tsurf$latitude == "medium" & tsurf$gts2012 < allstages$base[i] & tsurf$gts2012 >= allstages$top[i])]), na.rm = T)/sqrt(length(tsurf$temperature[which(tsurf$latitude == "medium" & tsurf$gts2012 < allstages$base[i] & tsurf$gts2012 >= allstages$top[i])]))
)
tpara[i,4,2] <- sqrt(var(T_calculate(tsurf$d18O[which(tsurf$latitude == "high" & tsurf$gts2012 < allstages$base[i] & tsurf$gts2012 >= allstages$top[i])],
tsurf$gts2012[which(tsurf$latitude == "high" & tsurf$gts2012 < allstages$base[i] & tsurf$gts2012 >= allstages$top[i])]), na.rm = T)/sqrt(length(tsurf$temperature[which(tsurf$latitude == "high" & tsurf$gts2012 < allstages$base[i] & tsurf$gts2012 >= allstages$top[i])]))
)
### SD lowlat T
tpara[i,2,3] <- sqrt(var(T_calculate(tsurf$d18O[which(tsurf$latitude == "low" & tsurf$gts2012 < allstages$base[i] & tsurf$gts2012 >= allstages$top[i])], tsurf$gts2012[which(tsurf$latitude == "low" & tsurf$gts2012 < allstages$base[i] & tsurf$gts2012 >= allstages$top[i])]), na.rm = T)
)
#### N for lowlat T
tpara[i,2,4] <-length(tsurf$temperature[which(tsurf$latitude == "low" & tsurf$gts2012 < allstages$base[i] & tsurf$gts2012 >= allstages$top[i])])
tpara[i,5,1] <- mean(get_aragonite(t = tpara[i,1,1], m = mgcastage[i,4]), na.rm = T)
tpara[i,6,1] <- mean(get_aragonite(t = tpara[i,2,1], m = mgcastage[i,4]), na.rm = T) # lowlat
tpara[i,7,1] <- mean(get_aragonite(t = tpara[i,3,1], m = mgcastage[i,4]), na.rm = T)
tpara[i,8,1] <- mean(get_aragonite(t = tpara[i,4,1], m = mgcastage[i,4]), na.rm = T)
tpara[i,9,1] <- mgcastage[i,4]
tpara[i,12,1] <- get_aragonite(t = tpara[i,2,1]-tpara[i,2,2], m = mgcastage[i,4]) # lowlat - standard error
tpara[i,13,1] <- get_aragonite(t = tpara[i,2,1]+tpara[i,2,2], m = mgcastage[i,4]) # lowlat + standard error
tpara[i,5,2] <- sqrt(var(get_aragonite(t = tsurf$temperature[which(tsurf$gts2012 < allstages$base[i] & tsurf$gts2012 >= allstages$top[i])], m = mgcastage[i,4]), na.rm = T)) / sqrt(length(tsurf$temperature[which(tsurf$gts2012 < allstages$base[i] & tsurf$gts2012 >= allstages$top[i])]))
tpara[i,6,2] <- sqrt(var(get_aragonite(t = tsurf$temperature[which(tsurf$latitude == "low" & tsurf$gts2012 < allstages$base[i] & tsurf$gts2012 >= allstages$top[i])], m = mgcastage[i,4]), na.rm = T)) / sqrt(length(tsurf$temperature[which(tsurf$latitude == "low" & tsurf$gts2012 < allstages$base[i] & tsurf$gts2012 >= allstages$top[i])]))
tpara[i,7,2] <- sqrt(var(get_aragonite(t = tsurf$temperature[which(tsurf$latitude == "medium" & tsurf$gts2012 < allstages$base[i] & tsurf$gts2012 >= allstages$top[i])], m = mgcastage[i,4]), na.rm = T)) / sqrt(length(tsurf$temperature[which(tsurf$latitude == "medium" & tsurf$gts2012 < allstages$base[i] & tsurf$gts2012 >= allstages$top[i])]))
tpara[i,8,2] <- sqrt(var(get_aragonite(t = tsurf$temperature[which(tsurf$latitude == "high" & tsurf$gts2012 < allstages$base[i] & tsurf$gts2012 >= allstages$top[i])], m = mgcastage[i,4]), na.rm = T)) / sqrt(length(tsurf$temperature[which(tsurf$latitude == "high" & tsurf$gts2012 < allstages$base[i] & tsurf$gts2012 >= allstages$top[i])]))
### SD ASI lowlat
tpara[i,6,3] <- sqrt(var(get_aragonite(t = tsurf$temperature[which(tsurf$latitude == "low" & tsurf$gts2012 < allstages$base[i] & tsurf$gts2012 >= allstages$top[i])], m = mgcastage[i,4]), na.rm = T))
# List with actual temperatures
if(nrow(tsurf[which(tsurf$gts2012 < allstages$base[i] & tsurf$gts2012 >= allstages$top[i]),])>= 1) {
tlow_list[[i]] <- T_calculate(tsurf$d18O[which(tsurf$latitude == "low" & tsurf$gts2012 < allstages$base[i] & tsurf$gts2012 >= allstages$top[i])], tsurf$gts2012[which(tsurf$latitude == "low" & tsurf$gts2012 < allstages$base[i] & tsurf$gts2012 >= allstages$top[i])])
} else tlow_list[[i]] <- NA
# List with the actual ASIs for all the temperatures
if(nrow(tsurf[which(tsurf$gts2012 < allstages$base[i] & tsurf$gts2012 >= allstages$top[i]),])>= 1) {
ASI_list[[i]] <- get_aragonite(t = tlow_list[[i]], m = rep(mgcastage[i,4],length(tlow_list[[i]])))
} else ASI_list[[i]] <- NA
}
for (i in 1:nrow(allstages)) { # interpolate T for stages in which no observation fell initially
if(i %in% c(45,58)) tpara[i,2,1] <- mean(c(tpara[i-1,2,1], tpara[i+1,2,1]))
# get ASI for those stages
if(i %in% c(45,58)) tpara[i,6,1] <- get_aragonite(t = tpara[i,2,1], m = mgcastage[i,4])
}
# use low latitude aragonite sea intensity:
asi_l <- tpara[,6,1]
# normalise to values between 0 and 1
library(scales)
asi <- rescale(asi_l,c(0,1))
Save temperature, Mg/Ca ratio and ASI for further use
t_low <- data.frame(matrix(nrow = 85, ncol = 4))
colnames(t_low) <- c("stagemid","Mean","SD","N")
t_low <- t_low %>% mutate(stagemid = -stage_midp, Mean = tpara[,2,1], SD = tpara[,2,3], N = tpara[,2,4])
# assign maximum SD to stages with 0 or 1 value only
t_low <- t_low %>% mutate(SD = case_when(is.na(SD) ~ 5, TRUE ~ SD))
t_low <- t_low %>% mutate(SE = case_when(N > 1 ~ SD/sqrt(N-1)))
t_low <- t_low %>% mutate(SE = case_when(is.na(SE) ~ max(SE, na.rm = T), TRUE ~SE))
# seperate out data with 0 or 1 observation for the plot
na_indices <- which(t_low$N %in% c(0,1))
T_NA <- t_low
T_NA$Mean[na_indices] <- NA
T_NA$SD[na_indices] <- NA
T_NA$SE[na_indices] <- NA
# plot Temperature
plot(T_NA$stagemid,T_NA$Mean, type = "l", col = "red", xlab = "Time (Ma)", ylab = "T (low latitudes)", pch = NA, ylim = c(10,40))
points(c(-stage_base[57]+2,-stage_top[57]-2), rep(T_NA$Mean[57],2), type = "l", lwd = 1, col = rgb(1,0,0,1))
points(-stage_midp[which(t_low$N == 1)], t_low$Mean[which(t_low$N == 1)], type = "p", cex = 0.5, col = rgb(1,0,0,1), pch = 19)
points(-stage_midp[which(t_low$N == 0)], t_low$Mean[which(t_low$N == 0)], type = "p", cex = 0.5, col = rgb(1,0,0,1), pch = 1)
# function to plot error envelopes
error_polygon2 <- function(ep,en,tstart,tend,tmid,color) {
nacheck <- is.na(ep)
segm <- rle(nacheck)
starts <- c(1,cumsum(segm$lengths[1:(length(segm$lengths)-1)])+1)
ends <- cumsum(segm$lengths)
for(i in 1:length(segm$values)) {
if(!(is.na(segm$values[i]))) {
ind <- starts[i]:ends[i]
tmidi <- tmid[ind]
tstarti <- tstart[ind]
tendi <- tend[ind]
polygon( c(tstarti[1], tmidi, tendi[length(tendi)], tendi[length(tendi)], rev(tmidi), tstarti[1]),
c((ep[ind])[1],ep[ind], (ep[ind])[length(ep[ind])], (en[ind])[length(en[ind])], rev(en[ind]), (en[ind])[1]),
border = NA, col = color)
}
}
}
error_polygon2(T_NA$Mean + 2*T_NA$SE, T_NA$Mean - 2*T_NA$SE,
-stage_base,-stage_top, -stage_midp,color = rgb(1,0,0,0.25))
legend("top", legend = c("1 observation", "interpolated"), col = rgb(1,0,0,1), pch = c(19,1), cex = 0.7)
#
ASI <- data.frame(matrix(nrow = 85, ncol = 4))
colnames(ASI) <- c("stagemid","Mean","SD","N")
ASI <- ASI %>% mutate(stagemid = -stage_midp, Mean = tpara[,6,1], SD = tpara[,6,3], N = tpara[,2,4])
# assign maximum SD to stages with 0 or 1 value only
ASI <- ASI %>% mutate(SD = case_when(is.na(SD) ~ max(ASI$SD, na.rm = T), TRUE ~ SD))
ASI <- ASI %>% mutate(SE = case_when(N > 1 ~ SD/sqrt(N-1)))
ASI <- ASI %>% mutate(SE = case_when(is.na(SE) ~ max(ASI$SE, na.rm = T), TRUE ~SE))
# seperate out data with 0 or 1 observation for the plot
na_indices <- which(ASI$N %in% c(0,1))
ASI_NA <- ASI
ASI_NA$Mean[na_indices] <- NA
ASI_NA$SD[na_indices] <- NA
ASI_NA$SE[na_indices] <- NA
# a function to rescale ASI to values between 0 and 1, based on the original curve
asi_rescale_function <- function(x) (x-min(ASI$Mean))/(max(ASI$Mean)-min(min(ASI$Mean)))
# plot ASI
plot(ASI_NA$stagemid,asi_rescale_function(ASI_NA$Mean), type = "l", col = rgb(0,0.3,1,1), xlab = "Time (Ma)", ylab = "T (low latitudes)", pch = NA, ylim = c(-0.075,1.02))
points(c(-stage_base[57]+2,-stage_top[57]-2), rep(asi_rescale_function(ASI_NA$Mean[57]),2), type = "l", lwd = 1, col = rgb(0,0.3,1,1))
points(-stage_midp[which(t_low$N == 1)],asi_rescale_function( ASI$Mean[which(ASI$N == 1)]), type = "p", cex = 0.5, col = rgb(0,0.3,1,1), pch = 19)
points(-stage_midp[which(t_low$N == 0)], asi_rescale_function(ASI$Mean[which(ASI$N == 0)]), type = "p", cex = 0.5, col = rgb(0,0.3,1,1), pch = 1)
error_polygon2(asi_rescale_function(ASI_NA$Mean + 2*ASI_NA$SE), asi_rescale_function(ASI_NA$Mean - 2*ASI_NA$SE),
-stage_base,-stage_top, -stage_midp,color = rgb(0,0.3,1,0.25))
legend("topleft", legend = c("1 observation", "interpolated"), col = rgb(0,0.3,1,1), pch = c(19,1), cex = 0.7)
#
#
Save a table with parameters and SD, SE
# write.csv(ASI,"asi.csv", row.names = F)
# write.csv(t_low,"temperature.csv", row.names = F)
#
# mgca_table <- as.data.frame(cbind(-stage_midp,tpara[,9,1]))
# colnames(mgca_table) <- c("stagemid", "Mean")
# write.csv(mgca_table, "mgca.csv", row.names = F)
#
# saveRDS(ASI_list,"ASI_all_values.rds")
# saveRDS(tlow_list, "Temperature_all_values.rds")
# test <- readRDS("ASI_all_values.rds")
#plot(1:85,log10(unlist(lapply(test,function(x) length(x[which(!(is.na(x)))])))+0.5))