#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 Newby, P., B. N. Shuman, J. P Donnelly, and D. MacDonald (2011). Evidence of centennial-scale droughts during the Younger Dryas and Holocene near the Hudson River watershed, U.S.A. Quaternary Research 75: 523-530. 

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

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

#Provide inputs
Num.cores=6
total.iterations=18
#version b (allow to rise only to modern sed limit)
SedLim.modern=40
#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)