####################################################################################
#SEDIMENTATION FUNCTION
####################################################################################
#burial - vector of average burial increments per time step
#variance of the normal distribution for increments
#tmax = maximum time 
#exp.erosion - multiplication factor that defines erosion in exponential model 
#that is multiplied by burial
#when exp.erosion=0.9, then erosion is 0.9*burial

sedimentation=function(burial, variance, tmax, method, exp.erosion) {
	if (missing(exp.erosion)) exp.erosion <- NULL
	if (missing(variance)) variance <- NULL

average=burial
#keep track of all erosional and depositional events
store.elevation=numeric()
#elevation - final elevation, it is overprinted by erosion events
elevation=numeric()
#starting elevation
elevation[1]=100
#keep track of thickness of each increment (positive or negative)
store.r.burial=numeric()
store.r.exp.burial=numeric()
store.r.exp.erosion=numeric()

#keep track of deposit age
age=numeric()
age[1]=1

for (i in 2:tmax) {
if (method=="normal") store.r.burial[i]=rnorm(1,mean=average[i], sd=sqrt(variance[i]))
if (method=="exponential") {
store.r.exp.burial[i]=rexp(1,rate=1/average[i])
store.r.exp.erosion[i]=rexp(1,rate=1/exp.erosion[i])
store.r.burial[i]=store.r.exp.burial[i]-store.r.exp.erosion[i]
}

#if positive sedimentation
if (store.r.burial[i]>=0) {
	elevation[i]=store.r.burial[i]+elevation[i-1]
	store.elevation[i]=elevation[i]
	age[i]=i
	}

#if erosion
if (store.r.burial[i]<0) {
new.elevation=store.r.burial[i]+elevation[i-1]
store.elevation[i]=new.elevation
#find the uppermost elevation level after erosion
check.erosion=max(which(elevation<new.elevation))
#if erosion would be deeper than initial elevation
if (is.infinite(check.erosion)){check.erosion=1}
if (length(check.erosion)==0) {check.erosion=1}
elevation[(check.erosion):i]=store.r.burial[i]+elevation[i-1]; 
#ages of eroded increments are lost, the recentmost age is assigned to the elevation
age[check.erosion:i]=i
i=check.erosion+1
}

print(elevation[i])
} 

list(final.elevation=elevation, initial.elevation=store.elevation, increment=store.r.burial, age=age)
}


####################################################################################
#RUN SEDIMENTATION, THE METHOD USED IN THE MANUSCRIPT USES NORMAL DISTRIBUTION OF INCREMENTS
####################################################################################
tmax<-11000 #in years, duration of model approximating the duration of sediment cores at Brijuni
METHOD="normal"
#mean long-term sedimentation rate (0.015 cm/y) 
mean.normal=rep(0.015,tmax)
#variance is its standard deviation for three stratigraphic units, with TST, MFZ, and HST values per year
var.normal=c(rep(0.75,4000),rep(0.1,2000),rep(0.005,5000))
if (METHOD=="normal") out=sedimentation(burial=mean.normal,variance=var.normal, tmax=tmax,method=METHOD)

#EXPONENTIAL MODEL NOT USED
#deposition and erosion thickness params for exponential model, not used
#higher burial must be compensated by higher exhumation to keep the long-term sedimentation rate constant
#mean.exp.burial=c(rep(0.015*50,5000),rep(0.015*20,2500),rep(0.015*2,4500))
#mean.exp.erosion=c(rep(0.015*49,5000),rep(0.015*19,2500),rep(0.015*1,4500))
#if (METHOD=="exponential") {out=sedimentation(burial=mean.exp.burial,tmax=tmax, exp.erosion=mean.exp.erosion, method=METHOD)

par(mfrow=c(1,2))
par(mar=c(20,2,1,1))
if (METHOD=="normal") label=paste("burial=", mean.normal[1], "var=",var.normal[1])
if (METHOD=="exponential") label=paste("burial=", mean.exp.burial[1], "erosion=",mean.exp.erosion[1])

plot(out$initial.elevation, type="l", frame=F, xlab="Time (years)", 
ylab="Stratigraphic elevation (cm)", main="", 
xaxt="n", yaxt="n",ylim=c(min(out$initial.elevation, na.rm=T),max(out$initial.elevation, na.rm=T)+20))
lines(out$final.elevation, col="blue")
#rescale axis labels
axis(1,at=seq(0,11000,by=2500), labels=seq(-11000,0,by=2500))
axis(2,at=seq(100,300,by=50), labels=seq(0,200,by=50))
rect(xleft=0,xright=5000, ytop=max(out$initial.elevation, na.rm=T), ybottom=-100, border=NA, col="gray31")
rect(xleft=5000,xright=7500, ytop=max(out$initial.elevation, na.rm=T), ybottom=-100, border=NA, col="gray91")
lines(out$initial.elevation)
lines(out$final.elevation, col="gray", lwd=3)
abline(v=c(5000,7500), lty=2)
legend(x="topright", col=c("black","gray"), lwd=c(2,3), bty="n", legend=c("Actual sedimentation","Preserved"))

####################################################################################
#split the section into 5 cm increments and assign them to age
####################################################################################
thickness=5
layers=cut(out$final.elevation, breaks=seq(out$final.elevation[1],max(out$final.elevation+thickness),by=thickness), labels=F)
top.age=tapply(out$age,layers,max)
strat.position=tapply(out$final.elevation,layers,max)
strat.layer=tapply(layers,layers,function (x) x[1])*thickness
delta.time=diff(top.age)
strat.separation=diff(strat.position)
rate=diff(strat.position)/diff(top.age)

#true generating rate
mean(mean.normal)
#actual mean sedimentation rate = thickness divided by time
long.term.rate=(max(out$final.elevation)-100)/tmax

###################################
#RESIDENCE TIME AT A GIVEN STRATIGRAPHIC ELEVATION UNDER SHIFTING VARIANCE 
#BUT THE SAME NET SEDIMENTATION RATE
#AT 5 CM INCREMENT THICKNESS
###################################
aggregate.final.5cm=cut(out$final.elevation, breaks=seq(0, 500, by=thickness), labels=F)*thickness
residence.5cm=tapply(out$final.elevation, aggregate.final.5cm, length)
median(residence.5cm)
thickness/long.term.rate
plot(residence.5cm, as.vector(names(residence.5cm)),  yaxt="n", type="l", lty=2, 
ylim=c(min(out$initial.elevation, na.rm=T),max(out$initial.elevation, na.rm=T)+20),
log="", xlim=c(0,3000), ylab="", frame=F, xlab="Time averaging (y)")
points(residence.5cm, as.vector(names(residence.5cm)),  pch=16, cex=1.4)
axis(2,at=seq(100,300,by=50), labels=seq(0,200,by=50))

###################################################
#INCREMENT THICKNESS TO 30 CM - BIOTURBATION EFFECT
###################################################
thickness.2=30
aggregate.final.20cm=cut(out$final.elevation, breaks=seq(0, 500, by=thickness.2), labels=F)*thickness.2
residence.20cm=tapply(out$final.elevation, aggregate.final.20cm, length)
lines(residence.20cm, as.vector(as.numeric(names(residence.20cm)))-(thickness.2/2),  lty=2, col="gray")
points(residence.20cm, as.vector(as.numeric(names(residence.20cm)))-(thickness.2/2),  pch=21, bg="gray", cex=1.4)
legend(x="topright", pch=c(16, 21), pt.bg=c("black","gray"), legend=c("5 cm-thick","30 cm-thick"))

###################################
#SADLER PLOT
###################################
#plot(delta.time, rate, log="xy", xlab="Temporal duration", ylab="Sedimentation rate", pch=16, frame=F)