Load Packages and Soil Data
require(ggplot2);require(SoilR);library(SoilR);library(FME);library(MASS);library(forecast);library(lattice);
library(RColorBrewer);source("multiplot.R");library(tidyverse);library(dplyr);library(tidyr);
library(stringr);library(lubridate);source("systemAge.R");source("transitTime.R");library(expm);
library(Metrics)
files<-("ERV_SoilR_Inputs_vFinal.csv")
Enter Model Parameters and Initial Guesses
#Select Site (options: "BM", "RV")
plot<-""
#Select Depth Interval(options:"10","20","30","40","60","80","100")
depth<-""
#Select Model ("P1","P2","P3")
mdl<-""
#Select # MCMC Iterations
nit=
#Enter Partitioning of C Inputs Between Pools
p1_in= #LF
p2_in= #OF
p3_in= #HF (calculated in model,value not used)
#Enter Initial Guess of Pool Turnover Constants
k1= #LF
k2= #OF
k3= #HF
#Enter Initial Guess of Pool Radiocarbon Contents
F1= #LF
F2= #OF
F3= #HF
#Pool Carbon Initialization (fraction of present day)
if(depth<){
pI=
}else{ if (depth>){
pI=
}}
Build Soil Carbon Dataframe
data<-tibble()
data<-read_csv(files)
names(data)<-c("sample","site","surface_depth","base_depth","sample_type","elevation","latitude", "longitude","OC_400","ROC","TIC","TC","TOC","sample_collection_year","14C_measured","bulkdensity_sample","ROCf","ROCfshale","SOC_stock","OCpetro","f_inert","f_petro","f_passive","mass_percent_LF","mass_percent_OF","mass_percent_HF","14C_bio","14C_effective")
if(mdl=="P1"){
plot_data<-data%>%filter(str_detect(sample,plot))
depth_exact<-(paste(depth,"\\b",sep=""))
depth_data<-plot_data%>%filter(str_detect(base_depth,depth_exact))
#Soil Bulk Density
blk_name<-(paste(plot,depth,"\\b",sep=""))
blk_data<-plot_data%>%filter(str_detect(sample,blk_name))
blk_bd<-blk_data$bulkdensity_sample
#Extract Carbon and Radiocarbon Content
blk_toc<-blk_data$TOC# bulk soil TOC
data_cal<-depth_data%>%as_tibble()%>%mutate(blk_cstock=blk_toc/100*blk_bd*100*100)# bulk soil C stock
obs_t<-as.numeric((slice(dplyr::select(data_cal,time=sample_collection_year),1))) #observation year
blk_cstock<-as.numeric(slice(dplyr::select(data_cal,blk_cstock),1))
blk14<-cbind(time=obs_t,blk14=as.numeric(slice(dplyr::select(blk_data,blk14=14C_measured),1))) # bulk soil 14C
lf_data<-data_cal%>%filter(str_detect(sample,"LF")); lf_data<-lf_data%>%as_tibble()%>%mutate(cstock=blk_cstock*mass_percent_LF/100)#pool C stock
lf<-cbind(time=obs_t,lf=as.numeric(slice(dplyr::select(lf_data,lf=cstock),1))) # LF C stock
lf14<-cbind(time=obs_t,lf14=as.numeric(slice(dplyr::select(lf_data,lf14=14C_measured),1))) # LF 14C
of_data<-data_cal%>%filter(str_detect(sample,"OF")); of_data<-of_data%>%as_tibble()%>%mutate(cstock=blk_cstock*mass_percent_OF/100)#pool C stock
of<-cbind(time=obs_t,of=as.numeric(slice(dplyr::select(of_data,of=cstock),1))) # OF C stock
of14<-cbind(time=obs_t,of14=as.numeric(slice(dplyr::select(of_data,of14=14C_measured),1))) # OF 14C
hf_data<-data_cal%>%filter(str_detect(sample,"HF")); hf_data<-hf_data%>%as_tibble()%>%mutate(cstock=blk_cstock*mass_percent_HF/100)#pool C stock
hf<-cbind(time=obs_t,hf=as.numeric(slice(dplyr::select(hf_data,hf=cstock),1))) # HF C stock
hf14<-cbind(time=obs_t,hf14=as.numeric(slice(dplyr::select(hf_data,hf14=14C_measured),1))) # HF 14C
}else{if(mdl=="P2"){
#shale OC calc.
data<-data%>%as_tibble()%>%mutate(P2_SOCbio=TOC-OCpetro)
plot_data<-data%>%filter(str_detect(sample,plot))
depth_exact<-(paste(depth,"\\b",sep=""))
depth_data<-plot_data%>%filter(str_detect(base_depth,depth_exact))
#Soil Bulk Density
blk_name<-(paste(plot,depth,"\\b",sep=""))
blk_data<-plot_data%>%filter(str_detect(sample,blk_name))
blk_bd<-blk_data$bulkdensity_sample
#Extract Carbon and Radiocarbon Content
blk_toc<-blk_data$P2_SOCbio# bulk soil TOC
data_cal<-depth_data%>%as_tibble()%>%mutate(blk_cstock=blk_toc/100*blk_bd*100*100)# bulk soil C stock
orig_blk_toc<-blk_data$TOC# bulk soil TOC
data_cal<-data_cal%>%as_tibble()%>%mutate(orig_blk_cstock=blk_toc/100*blk_bd*100*100)# bulk soil C stock
obs_t<-as.numeric((slice(dplyr::select(data_cal,time=sample_collection_year),1))) #observation year
blk_cstock<-as.numeric(slice(dplyr::select(data_cal,blk_cstock),1))
orig_blk_cstock<-as.numeric(slice(dplyr::select(data_cal,orig_blk_cstock),1))
blk14<-cbind(time=obs_t,blk14=as.numeric(slice(dplyr::select(blk_data,blk14=14C_bio),1))) # bulk soil 14C
lf_data<-data_cal%>%filter(str_detect(sample,"LF")); lf_data<-lf_data%>%as_tibble()%>%mutate(cstock=orig_blk_cstock*mass_percent_LF/100)#pool C stock
lf<-cbind(time=obs_t,lf=as.numeric(slice(dplyr::select(lf_data,lf=cstock),1))) # LF C stock
lf14<-cbind(time=obs_t,lf14=as.numeric(slice(dplyr::select(lf_data,lf14=14C_measured),1))) # LF 14C
of_data<-data_cal%>%filter(str_detect(sample,"OF")); of_data<-of_data%>%as_tibble()%>%mutate(cstock=orig_blk_cstock*mass_percent_OF/100)#pool C stock
of<-cbind(time=obs_t,of=as.numeric(slice(dplyr::select(of_data,of=cstock),1))) # OF C stock
of14<-cbind(time=obs_t,of14=as.numeric(slice(dplyr::select(of_data,of14=14C_measured),1))) # OF 14C
hf_data<-data_cal%>%filter(str_detect(sample,"HF")); hf_data<-hf_data%>%as_tibble()%>%mutate(cstock=blk_cstock*mass_percent_HF/100)#pool C stock
hf<-cbind(time=obs_t,hf=as.numeric(slice(dplyr::select(hf_data,hf=cstock),1))) # HF C stock
hf14<-cbind(time=obs_t,hf14=as.numeric(slice(dplyr::select(hf_data,hf14=14C_bio),1))) # HF 14CpI=0.2
}else{ if (mdl=="P3"){
#shale OC calc.
data<-data%>%as_tibble()%>%mutate(P3_SOCeff=TOC-ROC)
plot_data<-data%>%filter(str_detect(sample,plot))
depth_exact<-(paste(depth,"\\b",sep=""))
depth_data<-plot_data%>%filter(str_detect(base_depth,depth_exact))
#Soil Bulk Density
blk_name<-(paste(plot,depth,"\\b",sep=""))
blk_data<-plot_data%>%filter(str_detect(sample,blk_name))
blk_bd<-blk_data$bulkdensity_sample
#Extract Carbon and Radiocarbon Content
blk_toc<-blk_data$P3_SOCeff# bulk soil TOC
data_cal<-depth_data%>%as_tibble()%>%mutate(blk_cstock=blk_toc/100*blk_bd*100*100)# bulk soil C stock
orig_blk_toc<-blk_data$TOC# bulk soil TOC
data_cal<-data_cal%>%as_tibble()%>%mutate(orig_blk_cstock=blk_toc/100*blk_bd*100*100)# bulk soil C stock
obs_t<-as.numeric((slice(dplyr::select(data_cal,time=sample_collection_year),1))) #observation year
blk_cstock<-as.numeric(slice(dplyr::select(data_cal,blk_cstock),1))
orig_blk_cstock<-as.numeric(slice(dplyr::select(data_cal,orig_blk_cstock),1))
blk14<-cbind(time=obs_t,blk14=as.numeric(slice(dplyr::select(blk_data,blk14=14C_effective),1))) # bulk soil 14C
lf_data<-data_cal%>%filter(str_detect(sample,"mass_percent_LF")); lf_data<-lf_data%>%as_tibble()%>%mutate(cstock=orig_blk_cstock*mass_percent_LF/100)#pool C stock
lf<-cbind(time=obs_t,lf=as.numeric(slice(dplyr::select(lf_data,lf=cstock),1))) # LF C stock
lf14<-cbind(time=obs_t,lf14=as.numeric(slice(dplyr::select(lf_data,lf14=14C_measured),1))) # LF 14C
of_data<-data_cal%>%filter(str_detect(sample,"OF")); of_data<-of_data%>%as_tibble()%>%mutate(cstock=orig_blk_cstock*mass_percent_OF/100)#pool C stock
of<-cbind(time=obs_t,of=as.numeric(slice(dplyr::select(of_data,of=cstock),1))) # OF C stock
of14<-cbind(time=obs_t,of14=as.numeric(slice(dplyr::select(of_data,of14=14C_measured),1))) # OF 14C
hf_data<-data_cal%>%filter(str_detect(sample,"HF")); hf_data<-hf_data%>%as_tibble()%>%mutate(cstock=blk_cstock*mass_percent_HF/100)#pool C stock
hf<-cbind(time=obs_t,hf=as.numeric(slice(dplyr::select(hf_data,hf=cstock),1))) # HF C stock
hf14<-cbind(time=obs_t,hf14=as.numeric(slice(dplyr::select(hf_data,hf14=14C_effective),1))) # HF 14C
}}
R_data<-data_cal%>%filter(str_detect(sample,"SI"))
R14<-cbind(time=obs_t,R14=as.numeric(slice(dplyr::select(R_data,R14=14C_measured),1))) # R 14C
R<-cbind(time=obs_t,R=NA)# Respiration C
Build Atmospheric 14CO2 Curve
AtmBomb=function(){
preandpost=bind.C14curves(prebomb=IntCal13,postbomb=Hua2013$NHZone1,time.scale="AD")
bombcurve=preandpost[preandpost[,1]>=500,1:2]
yrs=seq(1966,2009.5,by=1/4) # A series of years by quarters
nz2=spline(Hua2013$NHZone2[,c(1,4)],xout=yrs) #Spline interpolation of the NH_Zone 2 dataset at a quaterly basis
nhz2=ts((nz2$y-1)*1000,start=1966,freq=4) #Transformation into a time-series object
m=ets(nhz2) #Fits an exponential smoothing state space model to the time series
f2=forecast(m,h=11*4) #Uses the fitted model to forecast 11 years into the future
bc=data.frame(Year=c(bombcurve[-dim(bombcurve)[1],1],
seq(tsp(f2$mean)[1],tsp(f2$mean)[2], by=1/tsp(f2$mean)[3])),
Delta14C=c(bombcurve[-dim(bombcurve)[1],2],as.numeric(f2$mean)))
return(bc)
}
Run SoilR with Parameter Optimization
#Soil R Model: 3-pools in parallel
ER14C=function(pars){
mod=ThreepParallelModel14(
t=years,
ks=pars,
C0=C0_pI,
F0_Delta14C=F0,
In=c_in,
gam1=p1_in,
gam2=p2_in,
#xi=, (environmental scaler)
inputFc=Fatm,
#lag, (radiocarbon lag)
pass=TRUE
)
n.R14m=getF14R(mod)
n.C14t=getF14(mod)
n.R=getReleaseFlux(mod)
n.C=getC(mod)
m.R=data.frame(time=years, R=rowSums(n.R))
m.R14=data.frame(time=years, R14=n.R14m)
m.lf=data.frame(time=years, lf=n.C[,1])
m.of=data.frame(time=years, of=n.C[,2])
m.hf=data.frame(time=years, hf=n.C[,3])
m.lf14=data.frame(time=years, lf14=n.C14t[,1])
m.of14=data.frame(time=years, of14=n.C14t[,2])
m.hf14=data.frame(time=years, hf14=n.C14t[,3])
modelOutput<-list("R14"=m.R14, "lf14"=m.lf14, "of14"=m.of14, "hf14"=m.hf14, "R"=m.R, "lf"=m.lf, "of"=m.of, "hf"=m.hf)
return(modelOutput)
}
#Cost Function
ERcost=function(pars){
modelOutput=ER14C(pars)
if (!is.na(lf14[2])) {
cost=modCost(model=modelOutput$lf14, obs=lf14, x="time", err=NULL)
if (!is.na(of14[2])) {cost=modCost(model=modelOutput$of14, obs=of14, x="time", err=NULL,cost=cost)
}
if (!is.na(hf14[2])) {cost=modCost(model=modelOutput$hf14, obs=hf14, x="time", err=NULL,cost=cost)
}
if (!is.na(R14[2])) {cost=modCost(model=modelOutput$R14, obs=R14, x="time", err=NULL,cost=cost)
}
if (!is.na(lf[2])) {cost=modCost(model=modelOutput$lf, obs=lf, x="time", err=NULL,cost=cost)
}
if (!is.na(of[2])) {cost=modCost(model=modelOutput$of, obs=of, x="time", err=NULL,cost=cost)
}
if (!is.na(hf[2])) {cost=modCost(model=modelOutput$hf, obs=hf, x="time", err=NULL,cost=cost)
}
}else{
if (!is.na(of14[2])) {
cost=modCost(model=modelOutput$of14, obs=of14, x="time", err=NULL)
if (!is.na(hf14[2])) {cost=modCost(model=modelOutput$hf14, obs=hf14, x="time", err=NULL,cost=cost)
}
if (!is.na(R14[2])) {cost=modCost(model=modelOutput$R14, obs=R14, x="time", err=NULL,cost=cost)
}
if (!is.na(lf[2])) {cost=modCost(model=modelOutput$lf, obs=lf, x="time", err=NULL,cost=cost)
}
if (!is.na(of[2])) {cost=modCost(model=modelOutput$of, obs=of, x="time", err=NULL,cost=cost)
}
if (!is.na(hf[2])) {cost=modCost(model=modelOutput$hf, obs=hf, x="time", err=NULL,cost=cost)
}
}else{
if (!is.na(hf14[2])) {
cost=modCost(model=modelOutput$hf14, obs=hf14, x="time", err=NULL)
if (!is.na(R14[2])) {cost=modCost(model=modelOutput$R14, obs=R14, x="time", err=NULL,cost=cost)
}
if (!is.na(lf[2])) {cost=modCost(model=modelOutput$lf, obs=lf, x="time", err=NULL,cost=cost)
}
if (!is.na(of[2])) {cost=modCost(model=modelOutput$of, obs=of, x="time", err=NULL,cost=cost)
}
if (!is.na(hf[2])) {cost=modCost(model=modelOutput$hf, obs=hf, x="time", err=NULL,cost=cost)
}
}else{
if (!is.na(R14[2])) {
cost=modCost(model=modelOutput$R14, obs=R14, x="time", err=NULL)
if (!is.na(lf[2])) {cost=modCost(model=modelOutput$lf, obs=lf, x="time", err=NULL,cost=cost)
}
if (!is.na(of[2])) {cost=modCost(model=modelOutput$of, obs=of, x="time", err=NULL,cost=cost)
}
if (!is.na(hf[2])) {cost=modCost(model=modelOutput$hf, obs=hf, x="time", err=NULL,cost=cost)
}
}else{
if (!is.na(lf[2])) {
cost=modCost(model=modelOutput$lf, obs=lf, x="time", err=NULL)
if (!is.na(of[2])) {cost=modCost(model=modelOutput$of, obs=of, x="time", err=NULL,cost=cost)
}
if (!is.na(hf[2])) {cost=modCost(model=modelOutput$hf, obs=hf, x="time", err=NULL,cost=cost)
}
}else{
if (!is.na(of[2])) {
cost=modCost(model=modelOutput$of, obs=of, x="time", err=NULL)
if (!is.na(hf[2])) {cost=modCost(model=modelOutput$hf, obs=hf, x="time", err=NULL,cost=cost)
}
}else{
if (!is.na(hf[2])) {cost=modCost(model=modelOutput$hf, obs=hf, x="time", err=NULL)
}
}
}
}
}
}
}
return(cost)
}
#Turnover Rate Parameter Fitting
upper=c(0.99,0.99,0.99) #upper turnover rate boundaries
lower=c(0.00001,0.00001,0.00001) #lower turnover rate boundaries
#Fitted Parameter Matrix to Obtain Best Parameter Set (Optimization algorithm)
#Model fitting using FME with Pseudo-random search method
ERfit1<-modFit(f=ERcost, p=k0, method="Pseudo", upper=upper, lower=lower)
kfit<-ERfit1$par
k_opt="Pseudo"
#Run Model with Best Fit Parameters
FullFit=ThreepParallelModel14(
t=years,
ks=kfit,
C0=C0_pI,
F0_Delta14C=F0,
In=c_in,
gam1=p1_in,
gam2=p2_in,
#xi=, (environmental scaler)
inputFc=Fatm,
#lag=0
pass=TRUE
)
#Get Model Output Variables
ff.R14m=getF14R(FullFit) #respired 14CO2
ff.C14m=getF14C(FullFit) #bulk Soil 14C
ff.C14t=getF14(FullFit) #14C of indiviudal pools
AccR=getAccumulatedRelease(FullFit)
ff.C=getC(FullFit) #pool Sizes for lf, of, and hf respectively
ff.rC=getReleaseFlux(FullFit)
ff.R=data.frame(time=years, R=rowSums(ff.rC))
ff.bC<-rowSums(ff.C) #calculate total bulk c
# Run a Bayesian optimization to estimate uncertainty ranges
set.seed(3); rnorm(3)
res=ERfit1$residuals
var0=ERfit1$var_ms_unweighted
kfit_p=ERfit1$par
ERmcmc=modMCMC(f=ERcost, p=kfit_p,niter=nit,jump=NULL,var0=var0,lower=lower,upper=upper,burninlength=nit/10)
#Look at the output
mcmc.summary=summary(ERmcmc)
#Rerun model using MCMC derived best fit parameters in order to plot results
#Run Model with Best Fit Turnover Parameters
FullFit2=ThreepParallelModel14(
t=years,
#ks=ERmcmc$bestpar,#the parameter set that gave the highest probability
ks=mcmc.summary[1,1:3],#the average parameter set
C0=C0_pI,
F0_Delta14C=F0,
In=c_in,
gam1=p1_in,
gam2=p2_in,
#xi=, (environmental scaler)
inputFc=Fatm,
#lag=0
pass=TRUE
)
#Get Model Output Variables
ff.R14m2=getF14R(FullFit2) #respired 14CO2
ff.C14m2=getF14C(FullFit2) #bulk Soil 14C
ff.C14t2=getF14(FullFit2) #14C of indiviudal pools
AccR2=getAccumulatedRelease(FullFit2)
ff.C2=getC(FullFit2) #pool Sizes for lf, of, and hf respectively
ff.rC2=getReleaseFlux(FullFit2)
ff.R2=data.frame(time=years, R=rowSums(ff.rC2))
ff.bC2<-rowSums(ff.C2) #calculate total bulk c
#Calculate Transit Times and System Ages
k=mcmc.summary[1,1:3]
ages=seq(0,100)
A2=diag(-k)
u2=matrix(c(p1_in*c_in,p2_in*c_in,p3_in*c_in),ncol=1) #vector of inputs
#System age
SA2=systemAge(A=A2, u=u2, a=ages)
SA2_mean=SA2$meanSystemAge
#Transit time
TT2=transitTime(A=A2, u=u2, a=ages)
TT2_mean=TT2$meanTransitTime