library(knitr)
opts_chunk$set(echo=FALSE, warning=FALSE, error=FALSE, message=FALSE)

Exploring the Penman-Monteith equation

Here, a couple of key figures using the Penman-Monteith equation. Leaf temperature is calculated by closing the energy balance (numerically), taking into all temperature effects. For the following figures, wind speed (and hence boundary layer conductance) and stomatal conductance were inputs, and transpiration and leaf temperature are simulated.

library(plantecophys)
library(magicaxis)
library(pander)

windlabel <- expression(Wind~speed~~(m~s^-1))


gslow <- 0.05
gshigh <- 0.5

Winds <- exp(seq(log(0.1), log(10), length=25))

These simulations are with plantecophys package version r packageVersion("plantecophys").

# Calculate Tleaf given known gs, from energy balance
rlow <- sapply(Winds, function(x)plantecophys:::FindTleaf(gs=gslow, Wind=x, Tair=25))
rhigh <- sapply(Winds, function(x)plantecophys:::FindTleaf(gs=gshigh, Wind=x, Tair=25))

plot(log10(Winds), rlow, ylim=c(20,30), type='l', col="red", axes=FALSE,
     xlab=windlabel, ylab=expression(T[leaf]~~(degree*C)))
points(log10(Winds), rhigh, type='l', col="blue")

abline(h=25)
magaxis(side=1, unlog=1)
axis(2)
box()
legend("topright", legend=sapply(c(bquote(g[s] == .(gslow)),
                          bquote(g[s] == .(gshigh)),
                          bquote(T[air])
                          ),as.expression), lty=1, col=c("red","blue","black"),
       cex=0.9)
# Calculate ELEAF given known gs; calculate Tleaf from energy balance.
f <- function(w,gs,...){

  tleaf <- plantecophys:::FindTleaf(gs=gs, Wind=w, Tair=25,VPD=1)
  flux <- plantecophys:::LeafEnergyBalance(Tleaf=tleaf, Wind=w, Tair=25, gs=gs, 
                                           VPD=1,
                                           returnwhat="fluxes")
  flux$Tleaf <- tleaf

return(flux)
}

# Wind speed and E/gs
gss <- seq(0.02, 0.5, length=25)
windlow <- 0.1
windhigh <- 10

wlow <- do.call(rbind, lapply(gss, function(x)f(windlow,x)))
whigh <- do.call(rbind, lapply(gss, function(x)f(windhigh,x)))

plot(gss, wlow$ELEAFeb, type='l', col="red", ylim=c(0,5),xlim=c(0,0.6),
     xlab=expression(g[s]~~(mol~m^-2~s^-1)),
     ylab=expression(E[leaf]~~(mmol~m^-2~s^-1)))

points(gss, whigh$ELEAFeb, type='l', col="blue")

legend("topleft", legend=sapply(c(bquote(wind == .(windlow)),
                                   bquote(wind == .(windhigh))),as.expression), 
       lty=1, col=c("red","blue"), cex=0.9)
rlow <- do.call(rbind,lapply(Winds, function(x)f(x, gs=gslow)))
rhigh <- do.call(rbind,lapply(Winds, function(x)f(x, gs=gshigh)))

plot(log10(Winds), rlow$ELEAFeb/gslow/1000,  type='l', col="red", ylim=c(0,0.02),
     axes=FALSE,  xlab=windlabel, ylab=expression(E[leaf]/g[s]~~(mol~mol^-1)))
points(log10(Winds), rhigh$ELEAFeb/gshigh/1000, type='l', col="blue")
abline(h=0.01)

magaxis(side=1, unlog=1)
axis(2)
box()
legend("topright", legend=sapply(c(bquote(g[s] == .(gslow)),
                                   bquote(g[s] == .(gshigh)),
                                   bquote(VPD)), as.expression),
       lty=1, col=c("red","blue","black"))

f <- lapply(Winds, function(x)FARAO2(Wind=x, energybalance=TRUE))
f <- do.call(rbind,f)

FARAO with energy balance

Here, I show some simulations with FARAO (FARquhar And Optimization). A new implementation follows Buckley et al. (2014). Using the leaf gas exchange model (that takes Ci, Tair - and other drivers of course - and calculates A, E, gs and Tleaf), lambda was calculated numerically with,

$$\frac{dA}{dE} == \frac{dA / dC_i}{dE / dC_i}$$

This was done by calculating A and E at a given Ci, then adding a very small number, and calculating $dA = A(C_i+d) - A(C_i)$. Then, the Ci at which the calculated lambda was equal to a preset value was found by optimization.

vpds <- seq(1, 3.5, length=25)
windlow <- 0.4
windhigh <- 10
runcon <- FARAO2(energybalance=FALSE, VPD=vpds, Tleaf=25)
run1 <- FARAO2(energybalance=TRUE, VPD=vpds, Wind=windlow, Tair=25)
run2 <- FARAO2(energybalance=TRUE, VPD=vpds, Wind=windhigh, Tair=25)

run1 <- run1[is.finite(run1$ELEAF),]
run2 <- run2[is.finite(run2$ELEAF),]

f2 <- function(x)10^-3*x$ELEAF/(x$VPD/101)
runcon$GSinf <- f2(runcon)
run1$GSinf <- f2(run1)
run2$GSinf <- f2(run2)

f <- function(x, vpdname="VPD")x$ALEAF/(x$Ca*sqrt(x[,vpdname]))
runcon$bb <- f(runcon)
run1$bb <- f(run1)
run2$bb <- f(run2)
run1$bb2 <- f(run1,"VPDleaf")
run2$bb2 <- f(run2, "VPDleaf")

with(runcon, plot(VPD, GS, type='l', ylim=c(0,0.16),
                  xlim=c(1,5),
                  xlab="VPD (air) (kPa)",
                  ylab=expression(g[s]~~(mol~m^-2~s^-1))))
with(run1, points(VPDleaf, GS, type='l', col="red"))
with(run2, points(VPDleaf, GS, type='l', col="blue"))
legend("bottomleft", c("No energy balance","Wind = 0.4", "Wind = 10"),
       lty=1, col=c("black","red","blue"))

with(runcon, plot(VPD, ALEAF/ELEAF, type='l', ylim=c(0,10),
                  xlim=c(1,5),
                  xlab="VPD (air) (kPa)",
                  ylab=expression(g[s]~~(mol~m^-2~s^-1))))
with(run1, points(VPDleaf, ALEAF/ELEAF, type='l', col="red"))
with(run2, points(VPDleaf, ALEAF/ELEAF, type='l', col="blue"))
legend("bottomleft", c("No energy balance","Wind = 0.4", "Wind = 10"),
       lty=1, col=c("black","red","blue"))



with(runcon, plot(bb, GS, xlim=c(0,0.04), pch=19))
with(run1, points(bb, GS, col="red", pch=19))
with(run2, points(bb, GS, col="blue", pch=19))

with(runcon, plot(bb, GS, xlim=c(0,0.04), pch=19))
with(run1, points(bb, GSinf, col="red", pch=19))
with(run2, points(bb, GSinf, col="blue", pch=19))
n <- 101
set.seed(1111)
pars <- runif(n, 200,1000)
vpds <- runif(n,1,3.5)
tairs <- runif(n,15,28)
cas <- runif(n, 380, 500)
windlow <- 0.4
windhigh <- 10


rcon <- FARAO2(energybalance=FALSE, VPD=vpds, PPFD=pars, Tleaf=tairs, Ca=cas)
r1 <- FARAO2(energybalance=TRUE, VPD=vpds, Tair=tairs, PPFD=pars, Ca=cas, Wind=windlow)
r2 <- FARAO2(energybalance=TRUE, VPD=vpds, Tair=tairs, PPFD=pars, Ca=cas, Wind=windhigh)
r1 <- r1[is.finite(r1$ELEAF),]
r2 <- r2[is.finite(r2$ELEAF),]

f <- function(x, vpdname="VPD")x$ALEAF/(x$Ca*sqrt(x[,vpdname]))
rcon$bb <- f(rcon)
r1$bb <- f(r1)
r2$bb <- f(r2)
r1$bb2 <- f(r1,"VPDleaf")
r2$bb2 <- f(r2, "VPDleaf")

# inferred GS using air VPD (this is what we might do for eddy data or when gbl is unknown)
f2 <- function(x)10^-3*x$ELEAF/(x$VPD/101)
rcon$GSinf <- f2(rcon)
r1$GSinf <- f2(r1)
r2$GSinf <- f2(r2)


with(rcon, plot(bb, GS, pch=19, ylim=c(0,0.25), xlim=c(0,0.04),
                  xlab=expression(A/C[a]*sqrt(D[leaf])), 
                  ylab=expression(Optimized~g[s]~~(mol~m^-2~s^-1))))
with(r1, points(bb2, GS, pch=19, col="red"))
with(r2, points(bb2, GS, pch=19, col="blue"))
legend("topleft", c("No energy balance","Wind = 0.4", "Wind = 10"),
       lty=1, col=c("black","red","blue"))

abline(lm(GS ~ bb, data=rcon))
abline(lm(GS ~ bb2, data=r1), col="red")
abline(lm(GS ~ bb2, data=r2), col="blue")
with(rcon, plot(VPD, ALEAF/ ELEAF, xlim= c(0,5), ylim=c(0,15)))
with(r1, points(VPD, ALEAF/ELEAF, col="red"))
with(r2, points(VPD, ALEAF/ELEAF, col="blue"))

with(rcon, plot(VPD, ALEAF/ ELEAF, xlim= c(0,5), ylim=c(0,15)))
with(r1, points(VPDleaf, ALEAF/ELEAF, col="red"))
with(r2, points(VPDleaf, ALEAF/ELEAF, col="blue"))
with(rcon, plot(bb, GS, pch=19, ylim=c(0,0.25), xlim=c(0,0.04),
                  xlab=expression(A/C[a]*sqrt(D[air])), 
                  ylab=expression(g[s]~~(mol~m^-2~s^-1))))
with(r1, points(bb, GS, pch=19, col="red"))
with(r2, points(bb, GS, pch=19, col="blue"))
legend("topleft", c("No energy balance","Wind = 0.4", "Wind = 10"),
       lty=1, col=c("black","red","blue"))

abline(lm(GS ~ bb, data=rcon))
abline(lm(GS ~ bb, data=r1), col="red")
abline(lm(GS ~ bb, data=r2), col="blue")
with(rcon, plot(bb, GSinf, pch=19, xlim=c(0,0.04),  ylim=c(0,0.2), 
                  xlab=expression(A/C[a]*sqrt(D[air])), 
                  ylab=expression(Inferred~g[s]~~(mol~m^-2~s^-1))))
with(r1, points(bb, GSinf, pch=19, col="red"))
with(r2, points(bb, GSinf, pch=19, col="blue"))
abline(lm(GSinf ~ bb, data=rcon))
abline(lm(GSinf ~ bb, data=r1), col="red")
abline(lm(GSinf ~ bb, data=r2), col="blue")

legend("topleft", c("No energy balance","Wind = 0.4", "Wind = 10"),
       lty=1, col=c("black","red","blue"))

I also fit g1 with the Medlyn et al. 2011 model using non-linear regression (and assuming g0=0). The table below shows two estimates: the first is using leaf-to-air VPD, and fitting actual optimized gs. The second estimate is using the air VPD, and using the inferred gs from Eleaf (by gs = E/VPD).

fitg1 <- function(x, vpdname, gsname="GS", g0=TRUE ){
  x$VPD <- x[,vpdname]
  x$gs <- x[,gsname]
  if(g0){
    fit <- nls(gs ~ g0 + 1.6*(1+g1/sqrt(VPD))*ALEAF/Ca, data=x,
             start=list(g0=0, g1=4))
    return(coef(fit))
  } else {
    fit <- nls(gs ~ 1.6*(1+g1/sqrt(VPD))*ALEAF/Ca, data=x,
             start=list(g1=4))
    return(coef(fit))
  }

}


df <- data.frame(g1_estimate1=c(fitg1(rcon, "VPD", g0=F),
                      fitg1(r2, "VPDleaf", g0=F),
                      fitg1(r1, "VPDleaf", g0=F)),
                 g1_estimate2=c(fitg1(rcon, "VPD", g0=F),
                      fitg1(r2, "VPD", "GSinf", g0=F),
                      fitg1(r1, "VPD", "GSinf", g0=F)))
rownames(df) <- c("Infinite gbl","Wind = 10", "Wind = 0.4")

pander(df)

Caveats

# This is to check it works; here we found a discnontinuity before.
windlow <- 0.5
windhigh <- 10

Cis <- seq(100,375, length=101)

rlow <- PhotosynEB(Ci=Cis, Wind=windlow)
rhigh <- PhotosynEB(Ci=Cis, Wind=windhigh)
rcon <- Photosyn(Ci=Cis)

with(rcon, plot(Ci, ALEAF, type='l'))
with(rlow, points(Ci, ALEAF, type='l', col="red"))
with(rhigh, points(Ci, ALEAF, type='l', col="blue"))

with(rcon, plot(Ci, ELEAF, type='l'))
with(rlow, points(Ci, ELEAF, type='l', col="red"))
with(rhigh, points(Ci, ELEAF, type='l', col="blue"))

with(rcon, plot(Ci, Tleaf, type='l'))
with(rlow, points(Ci, Tleaf, type='l', col="red"))
with(rhigh, points(Ci, Tleaf, type='l', col="blue"))
w <- seq(0.2, 10, length=25)

r <- FARAO2(energybalance=TRUE, Wind=w, Ca=400, Tair=25, VPD=2, Wleaf=0.05)

f2 <- function(x)10^-3*x$ELEAF/(x$VPD/101)

g1 <- function(x, vpdname="VPD", gsname="GS"){
  with(x, sqrt(x[,vpdname])*((x[,gsname]/1.6)*(Ca/ALEAF) -1))
}

r$g1_VPD <- g1(r, "VPD")
r$g1_VPDleaf <- g1(r, "VPDleaf")
r$GSinf <- f2(r)
r$g1_VPD_gsinf <- g1(r, "VPD", "GSinf")

with(r, plot(Gbh, g1_VPD, type='l', xlim=c(0,6), ylim=c(0,5)))
with(r, points(Gbh, g1_VPDleaf, type='l', lty=5))

with(r, points(Gbh, g1_VPD_gsinf, type='l', lwd=2))





getg1 <- function(w, dfr){

  rcontrol <- FARAO2(energybalance=FALSE, Ca=dfr$Ca, Tleaf=dfr$Tair, 
              VPD=dfr$VPD, PPFD=dfr$PPFD)
  r <- FARAO2(energybalance=TRUE, Wind=w, Ca=dfr$Ca, Tair=dfr$Tair, 
              VPD=dfr$VPD, Wleaf=0.05, PPFD=dfr$PPFD)

  fitg1 <- function(x, vpdname, gsname="GS"){
    x$VPD <- x[,vpdname]
    x$gs <- x[,gsname]
    fit <- nls(gs ~ 1.6*(1+g1/sqrt(VPD))*ALEAF/Ca, data=x,
             start=list(g1=4))
    return(coef(fit))
  }
  r$GSinf <- with(r, 10^-3*ELEAF/(VPD/101))

  l <- list(data=r, datacontrol=rcontrol, g1=data.frame(g1_0=fitg1(rcontrol,"VPD"),
                              g1_VPDair=fitg1(r, "VPD"),
                              g1_VPDleaf=fitg1(r, "VPDleaf"),
                              g1_VPDair_GSinf=fitg1(r, "VPD", "GSinf")
                              ))
  class(l) <- "gblrun"
  return(l)

}

plot.gblrun <- function(x,...){


  with(x$data, plot(ALEAF/(sqrt(VPD)*Ca), GS, pch=19,
                    ylim=c(0, max(GS)), xlim=c(0,0.05)))
  with(x$datacontrol, points(ALEAF/(sqrt(VPD)*Ca), GS, pch=19, col="darkgrey"))

  with(x$data, plot(ALEAF/(sqrt(VPDleaf)*Ca), GS, pch=19,col="blue",
                    ylim=c(0, max(GS)), xlim=c(0,0.05)))
  with(x$datacontrol, points(ALEAF/(sqrt(VPD)*Ca), GS, pch=19, col="darkgrey"))

  with(x$data, plot(ALEAF/(sqrt(VPD)*Ca), GSinf, pch=19,col="forestgreen",
                    ylim=c(0, max(GS)), xlim=c(0,0.05)))
  with(x$datacontrol, points(ALEAF/(sqrt(VPD)*Ca), GS, pch=19, col="darkgrey"))

}

vpds <- c(seq(1, 3.5, by=0.5))
rhs <- seq(20,60,by=5)
cas <- c(380,600)
tairs <- 22
ppfds <- c(400,800,1200)
dfr <- expand.grid(RH=rhs, Tair=tairs, Ca=cas, PPFD=ppfds)
dfr$VPD <- RHtoVPD(dfr$RH,dfr$Tair)

g <- getg1(0.8, dfr)


w <- seq(0.5, 5, length=10)
res <- lapply(w, function(x)getg1(x,dfr))

g1df <- do.call(rbind, lapply(res, "[[", "g1"))
g1df$Wind <- w

with(g1df, plot(Wind, g1_0, type='l', lwd=2, ylim=c(1.5,4.2)))
with(g1df, points(Wind, g1_VPDair, type='l', col="red"))
with(g1df, points(Wind, g1_VPDleaf, type='l', col="blue"))
with(g1df, points(Wind, g1_VPDair_GSinf, type='l', col="forestgreen"))

d <- do.call(rbind, lapply(res, "[[", "data"))
d$bb_VPD <- with(d, ALEAF / (sqrt(VPD)*Ca))
d$bb_VPDleaf <- with(d, ALEAF / (sqrt(VPDleaf)*Ca))

palette(rainbow(10))
with(d, plot(bb_VPD, GS, pch=15, col=as.factor(Wind), ylim=c(0,0.2), xlim=c(0,0.05)))
with(d, plot(bb_VPDleaf, GS, pch=15, col=as.factor(Wind), ylim=c(0,0.2), xlim=c(0,0.05)))
with(d, plot(bb_VPD, GSinf, pch=15, col=as.factor(Wind), ylim=c(0,0.2), xlim=c(0,0.05)))


RemkoDuursma/plantecophys documentation built on April 6, 2021, 2:45 a.m.