#This script by B. Shuman produces the lake-level reconstruction for Deep Pond, Massachusetts following the approach described in Appendix I of Pribyl and Shuman (2014) A Computational Approach to Quaternary Lake-Level Reconstruction Applied in the Central Rocky Mountains, Wyoming, USA. Quaternary Research. It relies on data and produces the reconstruction like that shown in Marsicek, J., S. Brewer, D. Foster, W.W. Oswald, and B. Shuman. (2013) Moisture and temperature changes associated with the mid-Holocene Tsuga decline in the northeastern United States. Quaternary Science Reviews 80:129-142. http://dx.doi.org/10.1016/j.quascirev.2013.09.001
 
#This version of the script also produces a figure showing the data and reconstructions similar to Fig. 8 in Pribyl and Shuman (2014) or Fig. 7 in Marsicek et al. (2013). Note that modification of script used by Pribyl and Shuman (2014) to account for times when all cores contain littoral sediments yields some lower reconstructed water levels than in Marsicek et al. (2013), and that this version of the reconstruction does not assume that the base of a core is a maximum constraint on water levels prior to the basal age of the core, which had been assumed by Marsicek et al. (2013).

#Read in core data
Core = read.csv("Deep_cores.csv",header=T)
Thresh = read.csv("FaciesThresholds_Deep.csv",header=T)

#Provide inputs
Num.cores=4
total.iterations=30
#version b (allow to rise only to modern sed limit)
SedLim.modern=200
#version a (allow to rise to modern level)
#SedLim.modern=0
Length=12000
Interp.interval=50

#Create interpolation intervals
Int50=((1:(Length/Interp.interval))*Interp.interval)-Interp.interval

#Interpolate core data

Data.i=matrix(NA,length(Int50),Num.cores)
Elev.i=matrix(NA,length(Int50),Num.cores)

i=1

while(i<Num.cores+1){

Core.i=approx(Core$Age[Core$Core.DepthRank==i],Core$LOI[Core$Core.DepthRank==i],method="linear",Int50)
Core.depth.i=approx(Core$Age[Core$Core.DepthRank==i],Core$Depth[Core$Core.DepthRank==i],method="linear",Int50)

Data.i[,i]=Core.i$y
Elev.i[,i]=Core.depth.i$y

i=i+1

}
#Re-scale elevation relative to modern sediment limit

Elev.i=Elev.i-SedLim.modern

#Create matrix for storing iterative output

rows=length(Int50)-2

Recon.Iterations=matrix(NA,rows,total.iterations)
UpperBounds=matrix(NA,rows,total.iterations)
LowerBounds=matrix(NA,rows,total.iterations)

#Initiate loop to run all threshold iterations

iteration=1

while(iteration<total.iterations+1){

Lit.Threshold=Thresh[iteration,1]
Sublit.Threshold=Thresh[iteration,2]

#Identify facies

Core.lit = cbind(Data.i< Lit.Threshold)
Core.sublit = cbind(Data.i< Sublit.Threshold)
Core.profundal = cbind(Data.i>= Sublit.Threshold)
#Estimate littoral elevations

All.lit.max=apply(Core.lit*Elev.i,1,max,na.rm=T)

#Estimate sublittoral elevation 

All.sublit.max=apply(Core.sublit*Elev.i,1,max,na.rm=T)

#Account for intervals with only littoral sediment

Only.lit=(All.lit.max==Elev.i[,1])*(All.sublit.max==All.lit.max)
All.sublit.max[Only.lit>0]=max(Elev.i,na.rm=T)

#Account for intervals with littoral but no sub-littoral sediment by using next deepest core (shallowest profundal core & thus the minimum depth of non-littoral sediments) as lower bound on position of sublittoral zone

NextDeepestCore=Core.profundal*Elev.i
NextDeepestCore[rowMeans(NextDeepestCore)==0]=max(Elev.i[rowMeans(NextDeepestCore)==0])
NextDeepestCore[NextDeepestCore==0]=NA

DeepestCoreElev=apply(Elev.i,1,max,na.rm=T)
NextDeepestCore[is.na(NextDeepestCore)]=DeepestCoreElev[is.na(NextDeepestCore)]
NextDeepestCore=apply(NextDeepestCore,1,min,na.rm=T)

All.sublit.max[All.sublit.max==All.lit.max]=NextDeepestCore[All.sublit.max==All.lit.max]
All.sublit.max[is.na(All.sublit.max)]=DeepestCoreElev[is.na(All.sublit.max)]

#Account for intervals with no littoral or sub-littoral sediment by using depth of shallowest core as lower bound on position of sublittoral zone

ShallowestCoreElev=apply(Elev.i,1,min,na.rm=T)
All.sublit.max=All.sublit.max+((All.sublit.max==0)*ShallowestCoreElev)

#Estimate sediment limit depths
All.sublit.smooth=(All.sublit.max[1:(length(Int50)-2)]+All.sublit.max[2:(length(Int50)-1)]+All.sublit.max[3:length(Int50)])/3
All.lit.smooth=(All.lit.max[1:(length(Int50)-2)]+All.lit.max[2:(length(Int50)-1)]+All.lit.max[3:length(Int50)])/3

SedLim = (All.sublit.smooth+All.lit.smooth)/2

#Add results to matrix of iterations

Recon.Iterations[,iteration]=SedLim
UpperBounds[,iteration]=All.lit.smooth
LowerBounds[,iteration]=All.sublit.smooth

iteration=iteration+1

}

SedLim = rowMeans(Recon.Iterations)
All.sublit.smooth=rowMeans(LowerBounds)
All.lit.smooth=rowMeans(UpperBounds)

#All.sublit.smooth=apply(LowerBounds,1,quantile,prob=0.95,na.rm=T)
#All.lit.smooth=apply(UpperBounds,1,quantile,prob=0.05,na.rm=T)

SedLim.smooth=cbind(Int50[2:(length(Int50)-1)],SedLim,All.lit.smooth, All.sublit.smooth)
colnames(SedLim.smooth)=c("Age","LakeElev.Anom","UpperBound","LowerBound")
SedLim.smooth=data.frame(SedLim.smooth)

#Write output file
write.csv(Recon.Iterations, "Iterations.csv", row.names = TRUE)
write.csv(SedLim.smooth, "LLrecon.csv", row.names = FALSE)


################ Plot #####################

#Make two panel plot

quartz(width=(4*0.66667)+1,height=6)

t=layout(matrix(c(1,2,2), 3, 1, byrow = TRUE))
Length=12000
par(mai=c(0,0.75,0.5,0.25))

plot(Core$Age[Core$Core.DepthRank==1],Core$LOI[Core$Core.DepthRank==1],type="l",xlim=c(Length,0),xaxp=c(Length,0,Length/3000),ylim=c(0,100),xaxt="n",xlab="",ylab="% LOI",frame.plot=F,lwd=1.5,col="dark gray")
abline(h=Thresh[1,2],lwd=1,col="gray")
abline(h=Thresh[1,1],lwd=1,col="gray")
title("A. Deep Pond core data")

i=1

while(i<Num.cores+1){

lines(Core$Age[Core$Core.DepthRank==i+1],Core$LOI[Core$Core.DepthRank==i+1],lwd=2/(i+1))

i=i+1

}

par(mai=c(1,0.75,0.5,0.25))

plot(Core$Age[Core$Core.DepthRank==1],Core$Depth[Core$Core.DepthRank==1]-SedLim.modern,type="l",xlim=c(Length,0),xaxp=c(Length,0,Length/3000),ylim=c(max(Core$Depth[Core$Age<Length])-SedLim.modern,0),xlab="Cal yr before AD 1950",ylab="Depth below modern (cm)",frame.plot=F)
title("B. Reconstruction and Age-Depth-Facies Relationships")

SedLim.smooth$LowerBound[1]=max(SedLim.smooth$LowerBound,na.rm=TRUE)
SedLim.smooth$UpperBound[1]=max(SedLim.smooth$UpperBound,na.rm=TRUE)
SedLim.smooth$LowerBound[length(SedLim.smooth$LowerBound)]=max(SedLim.smooth$LowerBound,na.rm=TRUE)
SedLim.smooth$UpperBound[length(SedLim.smooth$UpperBound)]=max(SedLim.smooth$UpperBound,na.rm=TRUE)
polygon(SedLim.smooth$Age,SedLim.smooth$UpperBound,col="light blue",border=FALSE)
polygon(SedLim.smooth$Age,SedLim.smooth$LowerBound,col="white",border=FALSE)

iteration=1

while(iteration<total.iterations+1){

lines(Int50[2:length(Int50)-2],Recon.Iterations[,iteration],col="dark gray",lwd=0.5)

iteration=iteration+1

}
i=1

while(i<Num.cores+1){

lines(Core$Age[Core$Core.DepthRank==i],Core$Depth[Core$Core.DepthRank==i]-SedLim.modern)
points(Int50,Elev.i[,i]*(Data.i[,i]<=Thresh[1,2]),pch=19,col="tan")
points(Int50,Elev.i[,i]*(Data.i[,i]<=Thresh[1,1]),pch=19,col="orange")
points(Int50,matrix(0,length(Int50),1),pch=19,cex=1.5,col="white")

i=i+1

}
lines(SedLim.smooth$Age,SedLim.smooth$LakeElev.Anom,col="blue",lwd=2)
points(Length*0.25,(max(Core$Depth[Core$Age<Length])-SedLim.modern)*0.8,pch=19,col="tan")
text(Length*0.23,(max(Core$Depth[Core$Age<Length])-SedLim.modern)*0.81,"Sub-littoral",pos=4)
points(Length*0.25,(max(Core$Depth[Core$Age<Length])-SedLim.modern)*0.9,pch=19,col="orange")
text(Length*0.23,(max(Core$Depth[Core$Age<Length])-SedLim.modern)*0.91,"Littoral",pos=4)