These code blocks were used to calculate the aragonite sea intensity (ASI) time series

Aragonite Sea Intensity (ASI) calculation

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

Temperature Data

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 ratio

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

Calculate ASI from temperature and Mg/Ca ratios

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