########################################################################## # Ecosystem model derived from flux tower data Harvard Forest Ma USA # Pincipal Investigators Munger & Wofsy, Harvard University # Data from http://public.ornl.gov/ameriflux/level2data.html # Joseph Wheatley, Biospherica ########################################################################## # get Data loc=file.path("ftp://cdiac.ornl.gov/pub/ameriflux/data/Level2/Sites_ByID/US-Ha1/with_gaps/usmaharv_1991_L2.csv"); download.file(loc,"temp.dat"); names=read.csv("temp.dat", header=TRUE,skip=17,nrows=1); loc=file.path("ftp://cdiac.ornl.gov/pub/ameriflux/data/Level2/Sites_ByID/US-Ha1/with_gaps/usmaharv_1991_L2.csv"); download.file(loc,"temp.dat"); ha_1991=read.csv("temp.dat", header=FALSE,skip=20); loc=file.path("ftp://cdiac.ornl.gov/pub/ameriflux/data/Level2/Sites_ByID/US-Ha1/with_gaps/usmaharv_1992_L2.csv"); download.file(loc,"temp.dat"); ha_1992=read.csv("temp.dat", header=FALSE,skip=20); loc=file.path("ftp://cdiac.ornl.gov/pub/ameriflux/data/Level2/Sites_ByID/US-Ha1/with_gaps/usmaharv_1993_L2.csv"); download.file(loc,"temp.dat"); ha_1993=read.csv("temp.dat", header=FALSE,skip=20); loc=file.path("ftp://cdiac.ornl.gov/pub/ameriflux/data/Level2/Sites_ByID/US-Ha1/with_gaps/usmaharv_1994_L2.csv"); download.file(loc,"temp.dat"); ha_1994=read.csv("temp.dat", header=FALSE,skip=20); loc=file.path("ftp://cdiac.ornl.gov/pub/ameriflux/data/Level2/Sites_ByID/US-Ha1/with_gaps/usmaharv_1995_L2.csv"); download.file(loc,"temp.dat"); ha_1995=read.csv("temp.dat", header=FALSE,skip=20); loc=file.path("ftp://cdiac.ornl.gov/pub/ameriflux/data/Level2/Sites_ByID/US-Ha1/with_gaps/usmaharv_1996_L2.csv"); download.file(loc,"temp.dat"); ha_1996=read.csv("temp.dat", header=FALSE,skip=20); loc=file.path("ftp://cdiac.ornl.gov/pub/ameriflux/data/Level2/Sites_ByID/US-Ha1/with_gaps/usmaharv_1997_L2.csv"); download.file(loc,"temp.dat"); ha_1997=read.csv("temp.dat", header=FALSE,skip=20); loc=file.path("ftp://cdiac.ornl.gov/pub/ameriflux/data/Level2/Sites_ByID/US-Ha1/with_gaps/usmaharv_1998_L2.csv"); download.file(loc,"temp.dat"); ha_1998=read.csv("temp.dat", header=FALSE,skip=20); loc=file.path("ftp://cdiac.ornl.gov/pub/ameriflux/data/Level2/Sites_ByID/US-Ha1/with_gaps/usmaharv_1999_L2.csv"); download.file(loc,"temp.dat"); ha_1999=read.csv("temp.dat", header=FALSE,skip=20); loc=file.path("ftp://cdiac.ornl.gov/pub/ameriflux/data/Level2/Sites_ByID/US-Ha1/with_gaps/usmaharv_2000_L2.csv"); download.file(loc,"temp.dat"); ha_2000=read.csv("temp.dat", header=FALSE,skip=20); loc=file.path("ftp://cdiac.ornl.gov/pub/ameriflux/data/Level2/Sites_ByID/US-Ha1/with_gaps/usmaharv_2001_L2.csv"); download.file(loc,"temp.dat"); ha_2001=read.csv("temp.dat", header=FALSE,skip=20); loc=file.path("ftp://cdiac.ornl.gov/pub/ameriflux/data/Level2/Sites_ByID/US-Ha1/with_gaps/usmaharv_2002_L2.csv"); download.file(loc,"temp.dat"); ha_2002 =read.csv("temp.dat", header=FALSE,skip=20); loc=file.path("ftp://cdiac.ornl.gov/pub/ameriflux/data/Level2/Sites_ByID/US-Ha1/with_gaps/usmaharv_2003_L2.csv"); download.file(loc,"temp.dat"); ha_2003=read.csv("temp.dat", header=FALSE,skip=20); loc=file.path("ftp://cdiac.ornl.gov/pub/ameriflux/data/Level2/Sites_ByID/US-Ha1/with_gaps/usmaharv_2004_L2.csv"); download.file(loc,"temp.dat"); ha_2004=read.csv("temp.dat", header=FALSE,skip=20); loc=file.path("ftp://cdiac.ornl.gov/pub/ameriflux/data/Level2/Sites_ByID/US-Ha1/with_gaps/usmaharv_2005_L2.csv"); download.file(loc,"temp.dat"); ha_2005=read.csv("temp.dat", header=FALSE,skip=20); loc=file.path("ftp://cdiac.ornl.gov/pub/ameriflux/data/Level2/Sites_ByID/US-Ha1/with_gaps/usmaharv_2006_L2.csv"); download.file(loc,"temp.dat"); ha_2006=read.csv("temp.dat", header=FALSE,skip=20); loc=file.path("ftp://cdiac.ornl.gov/pub/ameriflux/data/Level2/Sites_ByID/US-Ha1/with_gaps/usmaharv_2007_L2.csv"); download.file(loc,"temp.dat"); ha_2007=read.csv("temp.dat", header=FALSE,skip=20); ##################### # pre-processing # ##################### Names <- names(names); harvard <- rbind(ha_1991,ha_1992,ha_1993,ha_1994,ha_1995,ha_1996,ha_1997,ha_1998,ha_1999,ha_2000,ha_2001,ha_2002,ha_2003,ha_2004,ha_2005,ha_2006,ha_2007); # merged dataframes names(harvard) <- Names # available data temp <- (harvard[,] == -9999 | harvard[,] == -6999); #replace -9999 or -6999 with NAs harvard[temp]=NA; ## fortnightly data lims <- seq(from=1,to=365,by=14); lims[27] <- 367; # bound adjust, last fortnight has extra days fortnight <- function(d) { sum(ifelse(d 0.2 m/s #hh <- function(h) {h[w]} #harvardS <- lapply( harvardS, hh) harvardS <- data.frame("YEAR"=harvard$YEAR[w],"FN"=harvard$FN[w],"TA"=harvard$TA[w],"NEE" = harvard$NEE[w], "TS1" = harvard$TS1[w], "VPD" = harvard$VPD[w],"PAR" = harvard$PAR_in[w]) getF <- function(i) {subset(harvardS,FN==i)} # list of dataframes according to fortnight value harvardF <- lapply(seq(1:26),getF) # generate list # plot(harvardF[[20]],harvardF[[20]]$NEE) #plot response to PAR ########################## # Night Only ########################## harvardS.dark <- subset(harvardS,PAR<30) # dark of photonflux < 30mu mol/m2 s-1. #d <- (harvardS$PAR == 0); # dark #harvardS.dark <- data.frame("YEAR"=harvardS$YEAR[d],"FN"=harvardS$FN[d],"TA"=harvardS$TA[d],"NEE" = harvardS$NEE[d], "TS1" = harvardS$TS1[d], "VPD" = harvardS$VPD[d],"PAR" = #harvardS$PAR_in[d]) getF.dark <- function(i) {subset(harvardS.dark,FN==i)} # list of dataframes according to fortnight value harvardF.dark <- lapply(seq(1:26),getF.dark) # generate list ############################# # Day only # ############################# harvardS.light <- subset(harvardS,PAR>30) getF.light <- function(i) {subset(harvardS.light,FN==i)} # list of dataframes according to fortnight value harvardF.light <- lapply(seq(1:26),getF.light) # generate list ############################### # END OF PRE-PROCESSING # ############################### ############################### ############################################# # dark model: 2-stage Arhennius model fit # ############################################# N_params = 1; start_period = 1 end_period = 26 N_periods = end_period-start_period+1; coefs <- rep(0, N_periods*N_params); dim(coefs) <- c(N_periods,N_params); #nls.control(maxiter=100,minFactor=1/4096,tol=0.1) for(i in 1:N_periods) { fit_lowess <- lowess(harvardF.dark[[i]]$TA,harvardF.dark[[i]]$NEE,f=0.7); fit_nls<-nls(y~a, data=as.data.frame(fit_lowess), start=list(a=1)) for (j in 1:N_params){ coefs[i,j] <- as.numeric( coefficients(fit_nls))[j]}} R_E = coefs[,1]; ############################## # mean nocturnal respiration # ############################## getM = function(i) {median(harvardF.dark[[i]]$NEE)}; R_E.dark <- sapply(seq(1:26),getM); ######################################## # plots of response # ######################################## par(mfrow=c(2,2)) plot(harvardF.light[[8]]$PAR, harvardF.light[[8]]$NEE,main="9-22 Apr",pch=16,xlim=c(0,2000),ylim=c(-35,20),font.axis=2,font.lab=2,xlab="PAR",ylab="Daytime NEE") plot(harvardF.light[[14]]$PAR, harvardF.light[[14]]$NEE,main="2-15 Jul",pch=16,xlim=c(0,2000),ylim=c(-35,20),font.axis=2,font.lab=2,xlab="PAR",ylab="Daytime NEE") plot(harvardF.light[[18]]$PAR, harvardF.light[[18]]$NEE,main="27 Aug-09 Sep",pch=16,xlim=c(0,2000),ylim=c(-35,20),font.axis=2,font.lab=2,xlab="PAR",ylab="Daytime NEE") plot(harvardF.light[[22]]$PAR, harvardF.light[[22]]$NEE,main="22 Oct-04 Nov",pch=16,xlim=c(0,2000),ylim=c(-35,20),font.axis=2,font.lab=2,xlab="PAR",ylab="Daytime NEE") ################################################ # daylight model: without T and VPD dependence # ################################################ N_params = 3; start_period = 11 end_period = 20 N_periods = end_period-start_period+1; params.init <- rep(NA, 26*N_params); dim(params.init) <- c(26,N_params); fits.init <- as.list(rep(NA,26)); for(i in start_period:end_period) { fits.init[[i]] <- nls(NEE~ gamma - alpha*beta*PAR/(beta+alpha*PAR), data = harvardF.light[[i]], start=list(gamma=1,alpha=0.02,beta=40),trace=F); # model and VPD & TA for (j in 1:N_params){ params.init[i,j] <- as.numeric( coefficients(fits.init[[i]]))[j]}} fits.init.summary <- lapply(fits.init,summary) #windows() #plot(harvardF.light[[16]]$PAR,harvardF.light[[16]]$NEE,pch=16,main="NEE FN=16",xlab="PAR",ylab="NEE",font.axis=2,font.lab=2) #points(harvardF.light[[16]]$PAR, fitted(fits.init[[16]]),col=2,pch=16) ############################################################################ # T and VPD dependence of efficiency : Parameter determination # ############################################################################ N_params = 6; start_period = 11 end_period = 20 N_periods = end_period-start_period+1; params.mid <- rep(NA, 26*N_params); dim(params.mid) <- c(26,N_params); #nls.control(maxiter=4000,minFactor=1/4096,tol=0.1) fits.mid <- as.list(rep(NA,26)); for(i in start_period:end_period) { fits.mid[[i]] <- nls(NEE~ gamma - alpha*beta*(lambda^2/((TA-25)^2+lambda^2))*mu^2/(VPD^2+mu^2)*PAR/(beta+alpha*PAR), data = harvardF.light[[i]], start=list(gamma=1,alpha=0.02,beta=40,lambda=30,mu=2),trace=F); # model and VPD & TA for (j in 1:N_params){ params.mid[i,j] <- as.numeric( coefficients(fits.mid[[i]]))[j]}} fits.mid.summary <- lapply(fits.mid,summary) # lambda_av <- mean(na.omit(params.mid[,4])) mu_av <- mean(na.omit(params.mid[,5])) windows() plot(harvardF.light[[16]]$PAR,harvardF.light[[16]]$NEE,pch=16,main="NEE FN=16",xlab="PAR",ylab="NEE",font.axis=2,font.lab=2) points(harvardF.light[[16]]$PAR, fitted(fits.mid[[16]]),col=2,pch=16) ################################################ # Final Model # ################################################ windows() N_params = 3; start_period = 11 # 21 May - 3 -Jun end_period = 20 #24 sep - 7 oct N_periods = end_period-start_period+1; params.fin <- rep(NA, 26*N_params); dim(params.fin) <- c(26,N_params); #nls.control(maxiter=4000,minFactor=1/4096,tol=0.1) fits.fin <- as.list(rep(NA,26)); for(i in start_period:end_period) { fits.fin[[i]] <- nls(NEE~ gamma - alpha*beta*(lambda_av^2/((TA-25)^2+lambda_av^2))*mu_av^2/(VPD^2+mu_av^2)*PAR/(beta+alpha*PAR), data = harvardF.light[[i]], start=list(gamma=1,alpha=0.02,beta=40),trace=F); # model and VPD & TA for (j in 1:N_params){ params.fin[i,j] <- as.numeric( coefficients(fits.fin[[i]]))[j]}} fits.summary <- lapply(fits.fin,summary) # plots windows() x<-seq(from=0,to=2000,by=50) par(mfrow=c(1,2)) par(mar=c(5,5,5,0)) ii <- 14 plot(harvardF.light[[ii]]$PAR,harvardF.light[[ii]]$NEE,pch=16,font.axis=2,font.lab=2,xlim=c(0,2000),ylim=c(-40,20),xlab="",ylab=""); title(main="2 Jul-15 Jul",xlab= expression(paste("PAR (",mu," mol ", s^{-1}, m^{-2},")")),ylab=expression(paste("NEE (",mu," mol ", s^{-1}, m^{-2},")"))) points(harvardF.light[[ii]]$PAR, fitted(fits.fin[[ii]]),col=2,pch=16) lines(x, params.fin[ii,1]-params.fin[ii,2]*params.fin[ii,3]*x/(params.fin[ii,3]+params.fin[ii,2]*x),col=3,lwd=2); # optimal NEE lines(x, params.init[ii,1]-params.init[ii,2]*params.init[ii,3]*x/(params.init[ii,3]+params.init[ii,2]*x),col=1,lwd=2,lty=2); text(1000,-40,"\uA9 Biospherica",cex=0.7) legend("top",legend=c("observed NEE","model NEE","model optimal NEE","simple saturation model"),pch=c(16,16,NA,NA),lty=c(NA,NA,1,2),lwd=c(NA,NA,2,2),col=c(1,2,3,1),bty="n") ii <- 18 par(mar=c(5,0,5,5)) plot(harvardF.light[[ii]]$PAR,harvardF.light[[ii]]$NEE,pch=16,xlab="",ylab="",font.axis=2,font.lab=2,yaxt="n",xlim=c(0,2000),ylim=c(-40,20)) title(main="27 Aug-9 Sep",xlab= expression(paste("PAR (",mu," mol ", s^{-1}, m^{-2},")")),ylab="") points(harvardF.light[[ii]]$PAR, fitted(fits.fin[[ii]]),col=2,pch=16) lines(x, params.fin[ii,1]-params.fin[ii,2]*params.fin[ii,3]*x/(params.fin[ii,3]+params.fin[ii,2]*x),col=3,lwd=2); # optimal NEE lines(x, params.init[ii,1]-params.init[ii,2]*params.init[ii,3]*x/(params.init[ii,3]+params.init[ii,2]*x),col=1,lwd=2,lty=2); a_inv <- round(1/params.fin[[ii,2]],1); b<-round(params.fin[[ii,3]],1); g<-round(params.fin[[ii,1]],1); text(1000,-40,"flux tower data from http://public.ornl.gov/ameriflux/",cex=0.7) ###################### # parameter plots # ###################### windows() par( mfrow=c(1,3)) par( mar=c(6, 6, 6, 6) + 0.1) plot(params.fin[,1],ylim=c(0,10),xlim=c(9,22),xlab="Fortnight Index",ylab= expression(paste(mu," mol ", s^{-1}, m^{-2})),pch=16,main=expression(paste("Respiration (R_E)")),font.main=2,font.lab=1,font.axis=2,cex.lab=1.5,cex.main=2) #points(params.fin[,1],pch=16,xlim=c(5,22)) plot(params.fin[,3],ylim=c(0,80),xlim=c(9,22),xlab="Fortnight Index",ylab= expression(paste(mu," mol ", s^{-1}, m^{-2})),pch=16,main=expression(paste("GEP Saturation (",beta,")")),font.main=2,font.lab=1,font.axis=2,cex.lab=1.5,cex.main=2) #points(params.fin[,3],pch=16,xlim=c(9,22)) plot(params.fin[,2],ylim=c(0,0.25),xlim=c(9,22),xlab="Fortnight Index",main=expression(paste("Efficiency (",alpha,")")),pch=16,ylab="",font.main=2,font.lab=1,font.axis=2,cex.lab=1.5,cex.main=2) #points(1/params.fin[,2],pch=16,xlim=c(5,22)) #legend("bottomleft",legend=c("data","model","model optimum epsilon=1","initial fit"),col=c("black","red","green","blue"))