# GSA Data Repository 201Xxxx # Appendix S2 # # Nawrot R., Berensmeier M., Gallmetzer I., Haselmair A., Tomašových A. and Zuschin M. 2022. # Multiple phyla, one time resolution? Similar time averaging in benthic foraminifera, mollusk, # echinoid, crustacean and otolith fossil assemblages # # # -------------------------------------------------------------------------------------------- # R script for analyzing the age data # Rafał Nawrot, Department of Palaeontology, University of Vienna # email: rafal.nawrot@univie.ac.at, paleo.nawrot@gmail.com # # Last updated: 14 March 2022 #--------------------------------------------------------------------------------------------- ##### Required R packages #################################################################### library(rcarbon) library(moments) library(boot) ##### Uploading the data ##################################################################### ages<-read.csv("Appendix S1.csv", skip=4, stringsAsFactors=T) str(ages) ages$Taxon<-factor(ages$Taxon, levels=c("Corbula gibba", "Gouldia minima", "Adelosina intricata", "Elphidium crispum", "Echinocyamus pusillus", "Psammechinus/Paracentrotus", "Pilumnus sp.", "Gobius niger")) taxa<-levels(ages$Taxon) ##### Calibrating 14C ages ################################################################## # Calibration based on Marine20 curve (Heaton et al., 2020; doi: 10.1017/RDC.2020.68) # the regional marine reservoir correction (ΔR) of -61 (SD = 50) (Tomašových et al., 2019; # doi: 10.1111/sed.12516) cal.ages<-calibrate(x=ages$X14C.age[!is.na(ages$X14C.age)], errors=ages$X14C.error[!is.na(ages$X14C.age)], ids=ages$Specimen.ID[!is.na(ages$X14C.age)], calCurves="marine20", resOffsets=-61, resErrors=50, timeRange=c(55000,0)) #median probability ages ages$Median.calibrated.age.BP[!is.na(ages$X14C.age)]<-medCal(cal.ages) #calibration age errors (1 and 2 sigma) errors.cal.ages<-summary(cal.ages) ##### Basic parameters of the age distributions ############################################## # Combine AAR ages and calibrated 14C ages ages$Age.BP<-ages$Median.calibrated.age.BP ages$Age.BP[is.na(ages$Age.BP)]<-ages$AAR.age.BP[is.na(ages$Age.BP)] # Basic summary sum.age<-tapply(ages$Age.BP, ages$Taxon, summary) sum.age<-as.data.frame(do.call(rbind, sum.age)) #number of dated specimens per taxon n.dated<-tapply(ages$Age.BP, ages$Taxon, length) #interquartile age range (IQR; measure of time averaging) per taxon iqr.age<-tapply(ages$Age.BP, ages$Taxon, IQR) #median age per taxon med.age<-tapply(ages$Age.BP, ages$Taxon, median) #total age range per taxon sum.age$Range<-sum.age$Max.-sum.age$Min. #skewness sum.age$Skew<-tapply(ages$Age.BP, ages$Taxon, skewness) #Correlation between median age and IQR cor.test(med.age, iqr.age) #Age offsets between the taxa (pairwise differences between median ages) dist(med.age) #Kruskal-Wallis test and Wilcoxon test with Holm's correction comparing median ages kruskal.test(ages$Age.BP, g=ages$Taxon) pwil<-pairwise.wilcox.test(ages$Age.BP, ages$Taxon, p.adjust.method="holm") #save the results pw.table<-round(pwil$p.value, 3) pw.table[pwil$p.value<0.001]<-"<0.001" write.csv(pw.table, file="Table S2 pairwise Wilcoxson test.csv") #Bonferroni correction (for comparison) pairwise.wilcox.test(ages$Age.BP, ages$Taxon, p.adjust.method="bonferroni") ##### Bootstrap confidence intervals ########################################################## # 95% bias-corrected bootstrap percentile confidence intervals (CIs) for median age and IQR med.age.boot<-data.frame(Taxon=taxa, Median.age=med.age, Median.age.lCI=NA, Median.age.uCI=NA) iqr.age.boot<-data.frame(Taxon=taxa, IQR=iqr.age, IQR.lCI=NA, IQR.uCI=NA) boot.fun<-function(x, nrep=10000, stat, seed){ statfun<-function(a, b) stat(a[b]) if(!is.na(seed)) set.seed(seed) res<-boot(x, statistic=statfun, R=nrep, sim="ordinary") cis<-boot.ci(res, conf=0.95, type="bca") attributes(cis$bca)<-NULL cis$bca[4:5] } #random number generator seed set to 1 to allow reproduction of the results med.out<-tapply(ages$Age.BP, ages$Taxon, boot.fun, stat=median, seed=1) med.age.boot[,3:4]<-as.data.frame(do.call(rbind, med.out)) iqr.out<-tapply(ages$Age.BP, ages$Taxon, boot.fun, stat=IQR, seed=1) iqr.age.boot[,3:4]<-as.data.frame(do.call(rbind, iqr.out)) # Summary table: age distribution m<-apply(med.age.boot[,3:4], 1, function(x) paste(round(x[1]), "-", round(x[2]))) q<-apply(iqr.age.boot[,3:4], 1, function(x) paste(round(x[1]), "-", round(x[2]))) summary.table<-data.frame(Taxon=names(med.age), Number.of.specimens=n.dated, Median.age=paste0(round(med.age),"\n", "(",m,")"), Interquartile.age.range=paste0(round(iqr.age),"\n", "(",q,")"), Minimum.age=sum.age$Min., Maximum.age=sum.age$Max., Age.range=sum.age$Range, Skewness=round(sum.age$Skew,2)) write.csv(summary.table, file="Table 1.csv", row.names=F) ##### Size of skeletal elements ############################################################## # Summary table sum.size<-tapply(ages$Size.mm, ages$Taxon, summary) sum.size<-as.data.frame(do.call(rbind, sum.size)) sum.mass<-tapply(ages$Mass.mg, ages$Taxon, summary) sum.mass[1:2]<-NA #no data on mass for the two bivalves sum.mass<-as.data.frame(do.call(rbind, sum.mass)) sum.size1<-round(sum.size, 3) sum.mass1<-round(sum.mass, 3) summary.table.size<-data.frame(Taxon=rownames(sum.size), Number.of.specimens=n.dated, Mean.mass.mg=sum.mass1$Mean, Median.mass.mg=sum.mass1$Median, Mass.range.mg=paste(sum.mass1$Min., "-", sum.mass1$Max.), Mean.size.mm=sum.size1$Mean, Median.size.mm=sum.size1$Median, Size.range.mm=paste(sum.size1$Min., "-", sum.size1$Max.)) write.csv(summary.table.size, "Table S1.csv", row.names=F) ##### Figures ############################################################################### ### Figure 2: age distributions with shading of sea-level phases #colors for higher taxa cols<-c("seashell3","seashell3","skyblue","skyblue", "seagreen3", "seagreen3", "moccasin", "salmon") sealev<-c(10000, 7000, 5500, 2000, 0) #ages of sea-level phases sealevcol<-c("gray75", "gray85", "gray95", "gray99") bin.size<-500 cairo_pdf("Figure 2.pdf", height=4, width=7.125) op<-par(mfcol=c(2,4), mar=c(2.5, 2, 4, 0), oma=c(0.5, 1, 0, 1), mgp=c(1.5, 0.5, 0), las=1, cex.main=0.85, cex.axis=0.8, cex.lab=0.8, font.main=4, tcl=-0.3) for(i in 1:8){ ok<-ages$Taxon==taxa[i] plot(1:1, xlim=c(0,10000), ylim=c(0, 9.8), ylab="", xlab="", type="n", axes=F) for (s in 2:length(sealev)){ rect(sealev[s], 0, sealev[s-1], 9, col=sealevcol[s-1], border=sealevcol[s-1]) } par(new=T) hist(ages$Age.BP[ok], breaks=seq(0,10000, by=bin.size), ylim=c(0, 9.8), ylab="", xlab="Post-mortem age (yr B.P.)", yaxt="n", main=taxa[i], col=cols[i]) axis(2, at=seq(0,8, by=2)) segments(med.age[i], 0, med.age[i], 9, lty=2, lwd=1.5) if(i %in% c(1,2)) mtext("Number of specimens", side=2, line=1.5, cex=0.6, las=0, at=4) if(i==1) mtext("Bivalvia", side=3, line=2.5, cex=0.7) if(i==3) mtext("Foraminifera", side=3, line=2.5, cex=0.7) if(i==5) mtext("Echinoidea", side=3, line=2.5, cex=0.7) if(i==7) mtext("Brachyura", side=3, line=2.5, cex=0.7) if(i==8) mtext("Teleostei", side=3, line=2.5, cex=0.7) mtext(paste0("N = ", n.dated[i], "\nMedian age = ", round(med.age[i]), " yr B.P.", "\nIQR = ", round(iqr.age[i])), line=-1, cex=0.55) } par(op) dev.off() ### Figure 3: median age vs. IQR with 95% CIs #coding for mineralogy: aragonite 21, high-Mg calcite 24, low-Mg calcite 25, pchmin<-c(21, 21, 24, 25, 24, 24, 25, 21) cairo_pdf("Figure 3.pdf", height=4.66, width=4.66) op<-par(pty="s", mgp=c(1.7,0.5,0), tcl=-0.3, mar=c(3, 3, 2, 2)) plot(med.age, iqr.age, xlim=c(0, 7500), ylim=c(0, 7500), type="n", xlab="Median age (yr B.P.)", ylab="Interquartile range (yr)") segments(med.age.boot$Median.age.lCI, iqr.age, med.age.boot$Median.age.uCI, iqr.age) segments(med.age, iqr.age.boot$IQR.lCI, med.age, iqr.age.boot$IQR.uCI) points(med.age, iqr.age, cex=1.2, pch=pchmin, bg=cols) text(med.age-170, iqr.age+170, labels=1:8, cex=0.9) legend("topleft", legend=c("Bivalvia","Foraminifera","Echinoidea","Brachyura","Teleostei (otoliths)"), bty="n", pch=15, col=unique(cols), cex=0.9, pt.cex=1.2) legend("topright", legend=c("Aragonite", "High-Mg calcite", "Low-Mg calcite"), bty="n", pch=c(21,24,25), cex=0.9, pt.cex=1.2) par(op) dev.off() ### Figure S2: mean size vs. median age and IQR cairo_pdf("Figure S2.pdf", height=4, width=7.125) op<-par(mfrow=c(1,2), pty="s", mgp=c(1.7,0.5,0), tcl=-0.3, mar=c(3, 3, 1, 1), cex.axis=0.9) plot(sum.size$Mean, iqr.age, type="n", xlim=c(0.1, 5), ylim=c(0, 5000), xlab="Mean size (mm)", ylab="Interquartile age range (yr)") points(sum.size$Mean, iqr.age, cex=1, pch=pchmin, bg=cols) text(sum.size$Mean, iqr.age, labels=1:8, cex=0.8, pos=3) x<-cor.test(sum.size$Mean, iqr.age, method="pearson") legend('topright', leg=paste0("r = ", round(x$estimate, 2),", p = ", round(x$p.value, 2)), bty='n', cex=0.9) legend("bottomleft", legend=c("Bivalvia","Foraminifera","Echinoidea","Brachyura","Teleostei (otoliths)"), bty="n", pch=15, col=unique(cols), cex=0.8) legend("bottomright", legend=c("Aragonite", "High-Mg calcite", "Low-Mg calcite"), bty="n", pch=c(21,24,25), cex=0.8) mtext("A", side=3, line=0.1, at=-0.9, cex=1.4) plot(sum.size$Mean, med.age, type="n", xlim=c(0.1, 5), ylim=c(0, 8000), xlab="Mean size (mm)", ylab="Median age (yr B.P.)") points(sum.size$Mean, med.age, cex=1, pch=pchmin, bg=cols) text(sum.size$Mean, med.age, labels=1:8, cex=0.8, pos=3) x<-cor.test(sum.size$Mean, med.age, method="pearson") legend('topright', leg=paste0("r = ", round(x$estimate, 2),", p = ", round(x$p.value, 2)), bty='n', cex=0.9) mtext("B", side=3, line=0.1, at=-0.9, cex=1.4) par(op) dev.off() # -------------------------------------------------------------------------------------------- # R script for: 1) estimating IQR uncertainty for AAR ages # 2) estimating IQR expected under a given sedimentation rate and mixing depth # Adam Tomašových, Earth Science Institute, Slovak Academy of Sciences # email: Adam.Tomasovych@savba.sk # # Last updated: 16 February 2022 #--------------------------------------------------------------------------------------------- ############################################################################################# # 1) estimating IQR uncertainty for AAR ages ############################################################################################# bootci.iqr=function(temp) { x.resample<-numeric(); x25.resample<-numeric();x75.resample<-numeric();iqr.resample<-numeric(); for (i in 1:1000) { for (j in 1:length(temp)) {x.resample[j]<-sample(temp,1,replace=T)} x25.resample[i]<-quantile(x.resample, c(0.25), na.rm=T) x75.resample[i]<-quantile(x.resample, c(0.75), na.rm=T) iqr.resample[i]=x75.resample[i]-x25.resample[i] } out=list(resample25=x25.resample, resample75=x75.resample, resampleiqr=iqr.resample) out } ############################################################################################# #variance parameter for lognormal calibration model (TDK1) and for gamma calibration model (SPK0) for Corbula Corbula.log.d=-2.896 Corbula.gamma.d=4.784 #variance parameters for gamma and lognormal calibration models for Gouldia Gouldia.gamma.d=3.935 Gouldia.log.d=-2.52 Gouldia.ages=ages$Age.BP[ages$Taxon=="Gouldia minima"] Corbula.ages=ages$Age.BP[ages$Taxon=="Corbula gibba"] ############################################################################################# #GOULDIA ############################################################################################# par(mfrow=c(2,2)) Gouldia.Piran2.IQR=IQR(Gouldia.ages) boot.Gouldia.IQR=bootci.iqr(Gouldia.ages)$resampleiqr SIM.meaniqr.log.age=numeric(); SIM.meaniqr.gamma.age=numeric() for (j in 1:1000) { meaniqr.x.age=numeric(); x.log.age=numeric() x.gamma.age=numeric() for (x in 1:length(ages)) { #we note that log(ages[x]) is median sampled from the posterior distribution, and thus #equivalent to the mean parameter that defines the lognormally-distributed ages temp=rlnorm(length(Gouldia.ages), meanlog=log(Gouldia.ages[x]),sdlog=sqrt(exp(Gouldia.log.d))) x.log.age[x]=IQR(temp) temp=rgamma(length(Gouldia.ages), shape=Gouldia.ages[x]/exp(Gouldia.gamma.d),scale=exp(Gouldia.gamma.d)) x.gamma.age[x]=IQR(temp) } SIM.meaniqr.log.age[j]=mean(x.log.age); SIM.meaniqr.gamma.age[j]=mean(x.gamma.age); } par(mfrow=c(2,2)) hist(Gouldia.Piran2.IQR-SIM.meaniqr.log.age, breaks=seq(-100,1500, by=50), xlim=c(0,3000), col="gray", cex.main=0.9, xlab="Inter-quartile age range (yr)", main="Gouldia-IQR corrected-lognormal uncertainty") abline(v=Gouldia.Piran2.IQR-mean(SIM.meaniqr.log.age), lwd=2) #mean estimate of IQR corrected for Corbula with lognormal uncertainty Gouldia.Piran2.IQR-mean(SIM.meaniqr.log.age) hist(boot.Gouldia.IQR, breaks=seq(0,5000, by=50), add=T, col="white") hist(Gouldia.Piran2.IQR-SIM.meaniqr.gamma.age, breaks=seq(-100,1500, by=50), xlim=c(0,3000), col="gray", cex.main=0.9, xlab="Inter-quartile age range (yr)", main="Gouldia-IQR corrected-gamma uncertainty") abline(v=Gouldia.Piran2.IQR-mean(SIM.meaniqr.gamma.age), lwd=2) #mean estimate of IQR corrected for Corbula with lognormal uncertainty Gouldia.Piran2.IQR-mean(SIM.meaniqr.gamma.age) hist(boot.Gouldia.IQR, breaks=seq(0,5000, by=50), add=T, col="white") ############################################################################################# #CORBULA ############################################################################################# Corbula.Piran.IQR=IQR(Corbula.ages) boot.Corbula.IQR=bootci.iqr(Gouldia.ages)$resampleiqr SIM.meaniqr.log.age=numeric(); SIM.meaniqr.gamma.age=numeric() for (j in 1:1000) { meaniqr.x.age=numeric(); x.log.age=numeric() x.gamma.age=numeric() for (x in 1:length(Corbula.ages)) { temp=rlnorm(length(Corbula.ages), meanlog=log(Corbula.ages[x]),sdlog=sqrt(exp(Corbula.log.d))) x.log.age[x]=IQR(temp) temp=rgamma(length(Corbula.ages), shape=Corbula.ages[x]/exp(Corbula.gamma.d),scale=exp(Corbula.gamma.d)) x.gamma.age[x]=IQR(temp) } SIM.meaniqr.log.age[j]=mean(x.log.age); SIM.meaniqr.gamma.age[j]=mean(x.gamma.age); } hist(Corbula.Piran.IQR-SIM.meaniqr.log.age, breaks=seq(-100,3000, by=50), xlim=c(0,3000),col="gray", cex.main=0.9, xlab="Inter-quartile age range (years)",main="Corbula-IQR corrected-lognormal uncertainty") abline(v=Corbula.Piran.IQR-mean(SIM.meaniqr.log.age), lwd=2) Corbula.Piran.IQR-mean(SIM.meaniqr.log.age) hist(boot.Corbula.IQR, breaks=seq(0,4000, by=50), add=T, col="white") hist(Corbula.Piran.IQR-SIM.meaniqr.gamma.age, breaks=seq(-100,3000, by=50), xlim=c(0,3000),col="gray", cex.main=0.9, xlab="Inter-quartile age range (years)",main="Corbula-IQR corrected-gamma uncertainty") abline(v=Corbula.Piran.IQR-mean(SIM.meaniqr.gamma.age), lwd=2) Corbula.Piran.IQR-mean(SIM.meaniqr.gamma.age) hist(boot.Corbula.IQR, breaks=seq(0,4000, by=50), add=T, col="white") ######################################################################################## # 2) Compute IQR expected under a given sedimentation rate and mixing depth ######################################################################################## # vectors with Corbula and Gouldia ages (since 2013) and the corresponding 5 cm increments in the core Piran 2 # age data from: Tomašových A., Gallmetzer I., Haselmair A., Kaufman D.S., Mavrič B. and Zuschin M. 2019. # A decline in molluscan carbonate production driven by the loss of vegetated habitats encoded in the Holocene # sedimentary record of the Gulf of Trieste. Sedimentology, 66(3), pp.781-807. Corbula.Piran.ages=c(10695,9867,6048,10096,8138,7013,7244,7735,5739,6795,6051,8127,5212,7989,7326,7865,7616,10104,10758,9473,7558, 9209,9998,10355,7700,9891,5626,8594,8554,7969,6446,8792,6707,8539,9708,5650,8249,9410,9748,7961,10159,10451, 8384,7669,6783,8723,7390,6061,9147,6493,9385,8729,10738,6094,6281,9055,6064,8750,8594,3782,6576,5124,9838, 3942,8794,7275,6119,4257,7757,5973,9905,8482,3969,9171,6821,6349,7754,6038,7751,9578,9459,8297,6230,7305, 6760,6387,8154,9230,7979,7836,7662,6299,7800,6433,5706,7463,6402,7137,8997,6420,8725,6896,8807,8033,8141, 5808,8674,4381,8689,7098,5819,8756,6971,8614,9414,9291,6930,9381,7734,8823,4964,8396,9132,60,1016,16, 2184,1108,42,177,1303,152,7168,5246,2798,110,1141,1891,715,392,3135,4298,1355,1464,1088,980,2145, 140,115,253,7342,15,926,2597,1264,4186,2081,6796,4420,122,961,2169,3005,1144,1215,2000,483,3762, 6854,448,2275,6761,5849,7292,5672,5522,7917,5015,5607,8770,4567,6227,6913,3381,4787,7886,8182,5240,5758, 7469,4376,3560,6187,8107,5257,6235,5403,8980,8149,5455,3246,9079,7092,6080,3092,3970,7114,3841,7639,677, 7299,5965,3217,6832,3207,11024,6674,5575,7570,5706,7272,6868,7123,5989,7816,9651,9097,9213,6212,5023,4914,7977) Corbula.Piran.5depths=c(67.5,67.5,67.5,67.5,67.5,67.5,67.5,67.5,67.5,67.5,67.5,67.5,67.5,67.5,67.5,67.5,67.5,67.5,67.5,67.5,67.5, 67.5,67.5,67.5,67.5,67.5,67.5,87.5,87.5,87.5,87.5,87.5,87.5,87.5,87.5,87.5,87.5,87.5,87.5,87.5,87.5,87.5, 87.5,87.5,87.5,87.5,87.5,87.5,87.5,87.5,87.5,87.5,87.5,87.5,87.5,87.5,102.5,102.5,102.5,102.5,102.5,102.5,102.5, 102.5,102.5,102.5,102.5,102.5,102.5,102.5,102.5,102.5,102.5,102.5,102.5,102.5,102.5,127.5,127.5,127.5,127.5,127.5,127.5,127.5, 127.5,127.5,127.5,127.5,127.5,127.5,127.5,127.5,127.5,127.5,127.5,127.5,127.5,127.5,147.5,147.5,147.5,147.5,147.5,147.5,147.5, 147.5,147.5,147.5,147.5,147.5,147.5,147.5,147.5,147.5,147.5,147.5,147.5,147.5,147.5,147.5,147.5,147.5,147.5,2.0,2.0,2.0, 2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0, 2.0,2.0,2.0,2.0,2.0,2.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0, 10.0,10.0,10.0,10.0,10.0,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5, 27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,47.5,47.5,47.5,47.5,47.5,47.5,47.5,47.5,47.5,47.5,47.5,47.5, 47.5,47.5,47.5,47.5,47.5,47.5,47.5,47.5,47.5,47.5,47.5,47.5,127.5,102.5,102.5,102.5,127.5,127.5,127.5,127.5,127.5,127.5) Gouldia.Piran2.ages=c(3814,3204,2798,2317,537,677,1697,3256,3010,2679,1135,1589,3396,1483,3346,3066,4268,2745,2473,876,2883, 1139,4873,1853,620,2928,1059,475,3062,4577,2379,2529,4283,3587,2618,4966,11006,6675,4921,1255,5473,3741, 4221,5091,3145,3125,3547,3002,1338,2536,516,1416,1106,1180,1185,1279,1730,660,598,757,3815,2781,992, 1038,1770,598,741,1027,403,954,1038,308,345,3335,1090,329,303,3212,1444,162,685,1696,640,819, 1552,1816,1405,1412,3424,229,2934,1359,1359,1742,2126,2681,1467,1446,4183,2027,1800,658,4224,311,7477, 1584,746,473,3725,2828,3425,3973,5828,5428,4714,4981,2418,3066,3220,4975,3123,3713,3351,2362,3559,6000, 3642,5080,4931,2303,4674,6855,6023,6144,2429,4352,3586,7079,5281,3666,3718,3642,5527,2622,3547,3639,4943, 3859,2432,3492,4515,2910,4307,4649,5748,3370,4861,9719,7065,4940,5724,4612,3828,3851,4424,5610,3729,3566, 5384,4233,6622,3088,3260,3113,4036,4873,545,5860,3044,416,3841,6066,4858,5295,5774,4880,7198,5102,5291, 2840,3575,6129,846,4821,3788,6797,3930,3042,4822,8685,6277,10402,6926,6762,3401,6844,5863,6016,3173,7518, 10494,5593,1251,7271,11502,7779,5069,5404,7063,7472,2761,4040,1206,7180,6445,5481,7384,7310,4506,6606,4900, 5293,7744,8015,3472,5615,2846,6322,6129,2531,4327,4479,7221,3514,5286,9471,3383,3246,6612,5613,4844,7645, 4994,4997,6675,7696,3594,4718,5680,3564,4722,6603,6025,5815,7860,7828,2769,978,5967,4352,3196,3662,6736, 5371,5942,4938,4510,7666,4450,6308,4722,5268,1797,7873,4561,1546,6068,4779,2448,1209,2795,4209,4182,5705, 3581,3733,4259,7367,5421,2394,4953,5733,4281,5481,4229,3974,3950,3577,5349,4420,2055,6950,3444,991,3406,2819) Gouldia.Piran2.5depths=c(14.0,14.0,14.0,14.0,14.0,14.0,14.0,14.0,14.0,14.0,14.0,14.0,14.0,14.0,14.0,14.0,14.0,14.0,14.0,14.0,14.0, 14.0,14.0,14.0,14.0,14.0,14.0,14.0,18.0,18.0,18.0,18.0,18.0,18.0,18.0,18.0,18.0,18.0,18.0,18.0,18.0,18.0, 18.0,18.0,18.0,18.0,18.0,18.0,18.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0, 2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,10.0,10.0,10.0,10.0, 10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0, 10.0,10.0,10.0,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5, 27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,47.5,47.5,47.5,47.5,47.5,47.5,47.5,47.5,47.5, 47.5,47.5,47.5,47.5,47.5,47.5,47.5,47.5,47.5,47.5,47.5,47.5,47.5,47.5,47.5,47.5,47.5,47.5,47.5,47.5,47.5, 67.5,67.5,67.5,67.5,67.5,67.5,67.5,67.5,67.5,67.5,67.5,67.5,67.5,67.5,67.5,67.5,67.5,67.5,67.5,67.5,67.5, 67.5,67.5,67.5,67.5,67.5,67.5,67.5,67.5,67.5,87.5,87.5,87.5,87.5,87.5,87.5,87.5,87.5,87.5,87.5,87.5,87.5, 87.5,87.5,87.5,87.5,87.5,87.5,87.5,87.5,87.5,87.5,87.5,87.5,87.5,87.5,87.5,87.5,87.5,87.5,102.5,102.5,102.5, 102.5,102.5,102.5,102.5,102.5,102.5,102.5,102.5,102.5,102.5,102.5,102.5,102.5,102.5,102.5,102.5,102.5,102.5,102.5,102.5,102.5, 102.5,102.5,127.5,127.5,127.5,127.5,127.5,127.5,127.5,127.5,127.5,127.5,127.5,127.5,127.5,127.5,127.5,127.5,127.5,127.5,127.5, 127.5,127.5,127.5,127.5,127.5,127.5,127.5,127.5,127.5,127.5,127.5,147.5,147.5,147.5,147.5,147.5,147.5,147.5,147.5,147.5,147.5, 147.5,147.5,147.5,147.5,147.5,147.5,147.5,147.5,147.5,147.5,147.5,147.5,147.5,147.5,147.5,147.5,147.5,147.5,147.5,18.0,18.0, 18.0) ################################################################################################################################## #AGE MODEL FOR SEDIMENT CORE BASED ON GOULDIA AND CORBULA #TIME AVERAGING EXPECTED ON THE BASIS OF NET SEDIMENTATION RATE ESTIMATED ON THE BASIS OF #DOWNCORE CHANGES IN MEDIAN AGE ################################################################################################################################## out=split(c(Gouldia.Piran2.ages,Corbula.Piran.ages),c(Gouldia.Piran2.5depths,Corbula.Piran.5depths)) par(mfcol=c(6,3)) par(mar=c(3,4,1,1)) hist(c(unlist(out$"2"),unlist(out$"10")),col="gray",breaks=seq(0,12000,by=500),cex.main=0.9,main="0-12,cm",xlab="",ylim=c(0,30)) par(mar=c(4,4,1,1)) hist(c(unlist(out$"14"),unlist(out$"18")),col="gray",breaks=seq(0,12000,by=500),cex.main=0.9,main="12-20,cm",xlab="",ylim=c(0,20)) par(mar=c(4,4,1,1)) hist(unlist(out$"27.5"),col="gray",breaks=seq(0,12000,by=500),cex.main=0.9,main="20-30,cm",xlab="Postmortem,age,(y)",ylim=c(0,20)) par(mfrow=c(1,2)) par(mar=c(20,4,2,1)) out3=boxplot(split(c(Gouldia.Piran2.ages,Corbula.Piran.ages),c(Gouldia.Piran2.5depths,Corbula.Piran.5depths)), at=-as.numeric(names(split(c(Gouldia.Piran2.ages,Corbula.Piran.ages),c(Gouldia.Piran2.5depths,Corbula.Piran.5depths)))), col="gray81",range=1.5,boxwex=3,horizontal=TRUE,ylab="Sediment depth (cm)",xlab="Time (years)",main="",xlim=c(-90,0), ylim=c(0,12000),cex.main=0.9,outline=F,frame=F,yaxt="n") axis(2,at=-c(0,25,50,75,100,125,150), lab=-c(0,25,20,75,100,125,150), las=2) median.ages=tapply(c(Gouldia.Piran2.ages,Corbula.Piran.ages),c(Gouldia.Piran2.5depths,Corbula.Piran.5depths),median) mean.ages=tapply(c(Gouldia.Piran2.ages,Corbula.Piran.ages),c(Gouldia.Piran2.5depths,Corbula.Piran.5depths),mean) cm.levels=as.numeric(names(mean.ages)) lines(median.ages[c(1:5,7,8)],-cm.levels[c(1:5,7,8)]) thicknesses=diff(cm.levels[c(1:5,7,8)]) median.age.diffs=diff(median.ages[c(1:5,7,8)]) sed.rates=round(thicknesses/median.age.diffs,digits=3) text(9000,-5,labels=paste(sed.rates[1],"cm/y"),cex=0.9) text(9000,-11,labels=paste(sed.rates[2],"cm/y"),cex=0.9) text(9000,-15,labels=paste(sed.rates[3],"cm/y"),cex=0.9) text(9000,-20,labels=paste(sed.rates[4],"cm/y"),cex=0.9) text(9000,-35,labels=paste(sed.rates[5],"cm/y"),cex=0.9) text(9000,-75,labels=paste(sed.rates[6],"cm/y"),cex=0.9) mixed.depth=c(10,20,50) exp.residence=array(NA,dim=c(length(mixed.depth),length(sed.rates))) rownames(exp.residence)=c("Mixing depth=10 cm","Mixing depth=20 cm","Mixing depth=50 cm") for (i in 1:length(mixed.depth)) { exp.residence[i,]=round(mixed.depth[i]/sed.rates) } expected.IQR=round(log(3)/(1/exp.residence)) expected.IQR #add expected iqr to the 6th increment at 45-50 cm, extrapolating from sedimentation between 30 and 70 cm expected.IQR.50cm=rep(NA,7) expected.IQR.50cm[c(1:5,7)]=expected.IQR[3,] expected.IQR.50cm[6]=expected.IQR.50cm[5] expected.IQR.20cm=rep(NA,7) expected.IQR.20cm[c(1:5,7)]=expected.IQR[2,] expected.IQR.20cm[6]=expected.IQR.20cm[5] expected.IQR.10cm=rep(NA,7) expected.IQR.10cm[c(1:5,7)]=expected.IQR[1,] expected.IQR.10cm[6]=expected.IQR.10cm[5] multi.IQR=tapply(c(Gouldia.Piran2.ages,Corbula.Piran.ages),c(Gouldia.Piran2.5depths,Corbula.Piran.5depths),IQR) plot(expected.IQR.20cm,multi.IQR[2:8],pch=21,bg="gray",xlim=c(100,14000),ylim=c(100,10000),log="xy",frame=F,cex=1.2, xlab="Expected IQR (sedim. rate and ML)",ylab="Molluscan (Gouldia-Varicorbula) IQR (years)") lines(c(1,12000),c(1,12000),lty=2) points(expected.IQR.10cm,multi.IQR[2:8],pch=21,cex=1.2,bg="white") points(expected.IQR.50cm,multi.IQR[2:8],pch=21,cex=1.2,bg="black") legend(x="bottomright",legend=c("ML-10 cm","ML-20 cm","ML-50 cm"),pch=21,pt.bg=c("white","gray","black"),cex=1.2)