Leaf temperature simulations with NicheMapR

Michael Kearney

2024-05-25

Overview

This document provides some examples of how to simulate leaf temperatures with NicheMapR. The NicheMapR package has a general ectotherm heat budget model which can be run through the functions ‘ectotherm’ and ‘ectoR_devel’.

The ‘ectotherm’ function calls a Fortran program that computes the heat and water budget using output from the NicheMapR microclimate model. It’s very fast but it can be a little more involved to customise. The function is fully described in Kearney and Porter (2020).

The ‘ectoR_devel’ function is fully written in R and so is slower, but may be easier to implement and customise. It is of course easier to understand ‘under the hood’ given that it’s written in a high-level, and generally more familiar, language.

Both functions can be used to simulate a leaf and example code is provided below to do just that.

Example code is also given for the leaf temperature calculation described in Campbell and Norman (1998), via the NicheMapR function ‘T_leaf_Campbell’, the ‘FindTLeaf’ function of the plantecophys package (Duursma 2015), and the ‘tleaf’ function of the tealeaves package (Muir, 2019), and the output of all five algorithms models is compared. The ‘ectotherm’ and ‘ectoR_devel’ functions are more flexible, comprehensive (especially in regard to geometry and treatment of infrared radiatoin) and faster to run (esp. the ‘ectotherm’ function).

The NicheMapR package can be integrated with the Photosyn function of the plantecophys package so that stomatal dynamics is driven by photosynthesis-linked models (and allowing more accurate leaf temperatures to drive photosynthesis). An example of how to do this is stepped through in detail.

Finally an example simulation is provided where stomatal closure is additionally driven by water stress due to soil moisture.

Comparing leaf temperature models in R

Microclimate simulation

The first step is to install and load the three relevant packages.

Here’s the code to install them.

library(devtools)
install_github("ropensci/rnoaa")  # this package is no longer on CRAN
install_github("ilyamaclean/microclima")  # install the microclima package
install_github("mrke/NicheMapR")  # install NicheMapR
climate_data_path <- "put your path here"  # place where you want to put the global climate dataset
get.global.climate(folder = climate_data_path)  # download and unpack the global climate data
install.packages("tealeaves")
install.packages("plantecophys")
install.packages("futile.logger")  # needed for micro_usa

Now load the libraries.

library(NicheMapR)
library(tealeaves)
library(plantecophys)

Next we simulate a microclimate with the ‘micro_global’ function of NicheMapR. This function uses the global monthly climatology downloaded above in the installation instructions, and by default simulates the middle day of each month.

Here we run it for a site in the Australian arid zone, ‘Olive Downs’. The model calculates microclimate at a user-specified height ‘Usrhyt’ which in this example is 1 cm. It produces air temperature, wind speed and humidity estimates at this height and at the reference height for the forcing data, in this case 1.2 m.

# simulate a microclimate
loc <- c(141.861086, -29.050343)  # longitude, latitude, decimal degrees
Usrhyt <- 0.01  # leaf height, m
micro <- micro_global(loc = loc, Usrhyt = Usrhyt, microclima = 1)  # set microclima = 1 to get partitioned diffuse and direct solar - needs internet connection for the DEM download
## using microclima and elevatr to adjust solar for topographic and vegetation effects 
## Downloading digital elevation data

The results are returned as a list and are put in an object called ‘micro’. The ‘ectotherm’ function by default expects that ‘micro’ exists and passes items of that list to the microclimate model.

The aboveground microclimatic conditions are in a matrix called ‘metout’. The soil temperatures, including surface temperature, are in the ‘soil’ matrix. Note that the argument ‘microclima’ is set to 1, which means it uses code from the microclima package (Maclean et al., 2019) to partition solar radiation into direct and diffuse, and the diffuse fraction is then part of the microclimate model output. Pressure is obtained for the elevation used in the simulation. Some mock dates are also produced (given that this is based long long-term average monthly means and not particular dates).

metout <- as.data.frame(micro$metout)  # microclimate aboveground conditions
soil <- as.data.frame(micro$soil)  # microclimate soil conditions
PDIFs <- micro$diffuse_frac  # use variable diffuse fraction
# PDIFs <- rep(0.15, nrow(metout)) # assume uniform diffuse
# solar
P_a <- get_pressure(micro$elev)  # air pressure, Pa
P_a <- rep(P_a, nrow(metout))  # create hourly vector
dates <- micro$dates  # mock dates

Define ectotherm model settings

To simulate a leaf we need to turn off the behaviour and metabolism used for simulating animals. This can be achieved most easily by setting ‘live’ to 0. Setting ‘leaf’ to 1 means the model uses stomatal conductance to compute evaporation.

# model settings
live <- 0  # don't simulate behaviour or respiration
leaf <- 1  # use the stomatal conductance formulation for evaporation

Define leaf functional traits

Next the traits of the leaf must be defined. An ellipsoid shape is used in this simulation. Specification of the mass and the ratios of the linear dimensions of the shape (here, semi-minor axes relative to the semi-major axis) together with the density (default value 1 g cm\({^3}\)) allows the geometric functions in the ectotherm model to compute the absolute dimensions and surface areas.

Evaporative exchange is driven by the stomatal conductance parameters, which can be specified for the abaxial and adaxial sides of the leaf as well as the overall base leaf conductance through the cuticle. The ‘pct_wet’ parameter is set to 0.

Radiative exchange parameters can be set including the emissivities of the leaf and the substrate, the solar absorptivitives of the leaf, the configuration factors to the sky and ground, the percentage of the leaf conducting to the ground and the orientation to the sun.

# leaf functional traits
Ww_g <- 1  # wet weight, g
shape <- 2  # 0=plate, 1=cylinder, 2=ellipsoid
shape_b <- 0.0025  # ratio of b axis:a axis for ellipsoid 
shape_c <- 0.1176  # ratio of c axis:a axis for ellipsoid 
g_vs_ab <- 0.2  # leaf vapour conductance, abaxial (bottom of leaf), mol/m2/s
g_vs_ad <- 0  # leaf vapour conductance, adaxial (top of leaf), mol/m2/s
g_vs_ab_max <- 0.3  # leaf vapour conductance, abaxial (bottom of leaf), mol/m2/s
g_vs_ad_max <- 0  # leaf vapour conductance, adaxial (top of leaf), mol/m2/s
epsilon_sub <- 0.95  # emissivity of the substrate (-)
epsilon_sky <- 1  # emissivity of the sky (-), accounted for already in 'TSKY' output from microclimate model so make it 1
epsilon_L <- 0.97  # emissivity of the leaf (-)
alpha_L <- 0.5  # solar absorptivity of the leaf (-)
fatosk <- 0.5  # radiation configuration factor to sky (-)
fatosb <- 0.5  # radiation configuration factor to substrate (-)
conv_enhance <- 1.4  # convective enhancement factor (-)
pct_cond <- 0  # percent of leaf conducting to the ground (%)
postur <- 3  # leaf oriented horizontally, not tracking sun

The ectotherm model needs a vector of hourly stomatal conductance values. I.e., there is no ‘on-the-fly’ response of the stomata other than the stomata opening up to the maximum to avoid leaf temperature reaching the value specified by the argument ‘CT_max’ (default is 40 deg C if the model is run with ‘leaf = 2’). Here we simply close the stomata at night, based on the zenith angle of the sun as calculated by the microclimate model.

# set up vapour conductance vectors and simulate stomatal
# closure at night
f_open <- rep(1, nrow(metout))
f_open[metout$ZEN == 90] <- 0  # close stomata when the sun is set
g_vs_abs <- g_vs_ab * f_open
g_vs_ads <- g_vs_ad * f_open

Run the NicheMapR ectotherm model

Now we can run the model for each hour of the middle day of each month (specified in the microclimate which, as mentioned above, is by default coming from the object ‘micro’ which is assumed to be present in the R environment).

ptm <- proc.time()  # Start timing
ecto <- ectotherm(leaf = leaf, live = live, Ww_g = Ww_g, shape = shape,
    shape_b = shape_b, shape_c = shape_c, g_vs_ab = g_vs_abs,
    g_vs_ad = g_vs_ads, g_vs_ab_max = g_vs_ab_max, g_vs_ad_max = g_vs_ad_max,
    epsilon_sub = epsilon_sub, epsilon_sky = epsilon_sky, epsilon = epsilon_L,
    alpha_max = alpha_L, alpha_min = alpha_L, M_1 = 0, fatosk = fatosk,
    fatosb = fatosb, conv_enhance = conv_enhance, pct_cond = pct_cond,
    postur = postur, preshr = P_a, PDIF = PDIFs)
print(proc.time() - ptm)  # Stop the clock
##    user  system elapsed 
##    0.01    0.00    0.02
environ <- as.data.frame(ecto$environ)
T_leaf_NMR <- environ$TC

Now plot the results against air temperature at leaf height and at the reference height of 1.2 m

month <- 1  # choose a month to plot
sub <- which(floor(dates) + 1 == month)
T_air_1.2m <- metout$TAREF
T_air_1cm <- metout$TALOC
plot(seq(0, 23), T_air_1cm[sub], type = "l", col = "blue", ylab = "temperature, deg C",
    xlab = "hour of day", ylim = c(15, 50))
points(seq(0, 23), T_air_1.2m[sub], type = "l", col = "blue",
    lty = 2)
points(seq(0, 23), T_leaf_NMR[sub], type = "l", ylab = "T_leaf, deg C")

Run NicheMapR ectoR_devel model

Next we make the same calculation with the ‘ectoR_devel’ function, extracting the relevant microclimate variables from the metout and soil output tables and matching the substrate solar absorptivity with that assumed in the microclimate model run. This is done in a loop and is an order of magnitude slower than using the ‘ectotherm’ function.

TAs <- metout$TALOC
TGRDs <- soil$D0cm
TSKYs <- metout$TSKYC
VELs <- metout$VLOC
RHs <- metout$RHLOC
QSOLRs <- metout$SOLR
Zs <- metout$ZEN
alpha_sub <- (1 - micro$REFL)

n <- nrow(metout)
T_leaf_ectoR <- matrix(nrow = n)
ptm <- proc.time() # Start timing
for(i in 1:n){
  ectoR.out <- ectoR_devel(leaf = leaf,
                           Ww_g = Ww_g,
                           M_1 = 0, # make metabolic rate zero
                           shape = shape,
                           shape_b = shape_b,
                           shape_c = shape_c,
                           g_vs_ab = g_vs_abs[i],
                           g_vs_ad = g_vs_ads[i],
                           epsilon = epsilon_L,
                           alpha = alpha_L,
                           fatosk = fatosk,
                           fatosb = fatosb,
                           conv_enhance = conv_enhance,
                           pct_cond = pct_cond,
                           postur = postur,
                           elev = micro$elev,
                           alpha_sub = alpha_sub,
                           epsilon_sub = epsilon_sub,
                           epsilon_sky = epsilon_sky,
                           TA = TAs[i],
                           TGRD = TGRDs[i],
                           TSKY = TSKYs[i],
                           VEL = VELs[i],
                           RH = RHs[i],
                           QSOLR = QSOLRs[i],
                           Z = Zs[i],
                           pres = P_a[i],
                           PDIF = PDIFs[i])
  T_leaf_ectoR[i] <- ectoR.out$TC
}
print(proc.time() - ptm) # Stop the clock
##    user  system elapsed 
##    0.49    0.00    0.50
# plot
plot(seq(0, 23), T_air_1cm[sub], type = 'l', col = 'blue', ylab = 'temperature, deg C', xlab = 'hour of day', ylim = c(15, 50))
points(seq(0, 23), T_air_1.2m[sub], type = 'l', col = 'blue', lty = 2)
points(seq(0, 23), T_leaf_NMR[sub], type = 'l', ylab = 'T_leaf, deg C')
points(seq(0, 23), T_leaf_ectoR[sub], type = 'l', col = 2, lty = 2)

Run the Campbell and Norman model

For comparison we can run the ‘humid operative temperature’ described in Campbell and Norman (1998), via the NicheMapR function ‘leaf_temperature’.

This model requires geometrical traits of the leaf that were internally computed in ‘ectotherm’ and ‘ectoR_devel’ via the ‘GEOM’ subroutine. Specifically ‘leaf_temperature’ needs the leaf width, the total surface area and the silhouette area (i.e. the area receiving direct solar radiation). We can make obtain these via the ‘GEOM_ecto’ function that is part of the ‘ectoR_devel’ algorithm. Note that ‘leaf_temperature’ assumes that the characteristic dimension of the leaf is 0.7 times the width.

The ‘leaf_temperature’ function also needs the base leaf vapor conductance, ‘g_vs_base’. This is not needed explicitly for the ‘ectotherm’ and ‘ectoR_devel’ function because it is derived as a function of the mass transfer coefficient which, in turn, is computed from the convection (heat transfer) coefficient through the Chilton–Colburn analogy.

# get characteristic dimension and areas using ecto_devel
# function GEOM_ecto
GEOM.out <- GEOM_ecto(AMASS = Ww_g/1000, GEOMETRY = shape, SHP = c(1,
    shape_b, shape_c), PTCOND = pct_cond/100, PMOUTH = 0, SKINW = 0/100)
w <- GEOM.out$AL/0.7  # leaf width, m
A <- GEOM.out$AREA  # total leaf surface area, m2
# A_sil <- (GEOM.out$ASILN + GEOM.out$ASILP) / 2 # leaf
# silhouette area to direct solar radiation, m2
A_sil <- A/2  # half the leaf is facing up to direct solar

g_vs_base <- 0.01  # base leaf vapour conductance, mol/m2/s (equivalent to 0.1 (µmol H2O) / (m^2 s Pa) from tealeaves)

Now run the calculation with ‘lapply’. In this case the prediction is almost identical, but it won’t always be so depending on the environmental conditions, due to the subtle differences in how the models treat convection.

# Campbell and Norman Calculation
tm <- proc.time() # Start timing
T_leaf_Campbell <- unlist(lapply(1:length(TAs), function(x){leaf_temperature(
  w = w, # leaf width, m
  A = A, # leaf surface area, m^2
  A_sil = A_sil, # leaf silhouette area, m^2
  heliotropic = 0, # don't track the sun
  alpha_L = alpha_L, # leaf solar absorptivity, -
  alpha_S = alpha_sub, # ground solar absorptivity, -
  epsilon_L = epsilon_L, # leaf emissivity, -
  epsilon_sub = epsilon_sub, # ground emissivity, -
  epsilon_sky = epsilon_sky, # sky emissivity, -
  g_vs_ab = g_vs_abs[x] + g_vs_base / 2, # vapour conductance, abaxial (bottom of leaf), mol/m2/s
  g_vs_ad = g_vs_ads[x] + g_vs_base / 2, # vapour conductance, adaxial (top of leaf), mol/m2/s
  TA = TAs[x], # air temperature at lizard height from microclimate model, deg C
  TGRD = TGRDs[x], # ground temperature from microclimate model, deg C
  TSKY = TSKYs[x], # sky temperature from microclimate model, deg C
  VEL = VELs[x], # wind speed from microclimate model, m/s
  RH = RHs[x], # relative humidity from microclimate model, %
  QSOLR = QSOLRs[x], # total horizontal plane solar radiation from microclimate model, W/m2
  Z = Zs[x], # solar zenith angle from microclimate model, degrees
  PRESS = P_a[x],
  PDIF = PDIFs[x], # proportion solar radiation that is diffuse, -
  conv_enhance = conv_enhance # convective enhancement factor
)})) # run leaf_temperature across environments
print(proc.time() - ptm) # Stop the clock
##    user  system elapsed 
##    0.63    0.11    0.81
# plot
plot(seq(0, 23), T_air_1cm[sub], type = 'l', col = 'blue', ylab = 'temperature, deg C', xlab = 'hour of day', ylim = c(15, 50))
points(seq(0, 23), T_air_1.2m[sub], type = 'l', col = 'blue', lty = 2)
points(seq(0, 23), T_leaf_NMR[sub], type = 'l', ylab = 'T_leaf, deg C')
points(seq(0, 23), T_leaf_Campbell[sub], type = 'l', col = 2, lty = 2)

Run the plantecophys leaf temperature model

The next bit of code shows how to use the ‘FindTleaf’ function in plantecophys and compares its predictions to the ‘ectotherm’ function. Because ‘FindTleaf’ needs vapour pressure deficit, the WETAIR function is used to compute this from relative humidity, air temperature and atmospheric pressure. It also takes photon flux density in \(\mu \space mol \space m^{-2}\) rather than solar radiation in \(W m^{-2}\), so a conversion is needed.

Note that ‘FindTleaf’ doesn’t consider sky and ground longwave radiation fluxes as an explicit function of ‘sky temperature’ and ground temperature, which is the main reason for the discrepancies.

# plant ecophys calculation
UMOLPERJ <- 4.57 # mu_mol photons / J
WETAIR_out <- WETAIR(db = TAs, rh = RHs, bp = P_a)
VPD_Pa <- WETAIR_out$esat - WETAIR_out$e
tm <- proc.time() # Start timing
T_leaf_ecophys <- unlist(lapply(1:length(TAs), function(x){FindTleaf(
  Tair = TAs[x],
  gs = g_vs_abs[x] + g_vs_ads[x] + g_vs_base,
  PPFD = QSOLRs[x] * 0.5 * UMOLPERJ, # assuming PAR is 0.5 of global radiation
  VPD = VPD_Pa[x] / 1000,
  Patm = P_a[x] / 1000,
  Wind = VELs[x],
  Wleaf = w,
  StomatalRatio = 1,
  LeafAbs = alpha_L)})) # run FindTleaf across environments
print(proc.time() - ptm) # Stop the clock
##    user  system elapsed 
##    0.72    0.17    1.00
# plot
plot(seq(0, 23), T_air_1cm[sub], type = 'l', col = 'blue', ylab = 'temperature, deg C', xlab = 'hour of day', ylim = c(15, 50))
points(seq(0, 23), T_air_1.2m[sub], type = 'l', col = 'blue', lty = 2)
points(seq(0, 23), T_leaf_NMR[sub], type = 'l', ylab = 'T_leaf, deg C')
points(seq(0, 23), T_leaf_ecophys[sub], type = 'l', col = 2, lty = 2)

Run the tealeaves model

The last comparison is against the tealeaves package. For this we have to alter the units involved. Note that it is very slow, likely because of the use of the units package to attach units to the inputs.

# Leaving the make_* functions empty will automatically
# default to defaults parameters.
leaf_par <- make_leafpar()  # leaf parameters
enviro_par <- make_enviropar()  # environmental parameters
constants <- make_constants()  # physical constants
T_leaf <- tleaf(leaf_par, enviro_par, constants, quiet = TRUE)
T_leaf_tealeaves <- matrix(data = NA, nrow = length(TAs), ncol = 12)

pctdones <- seq(0, length(TAs), round(0.01 * length(TAs), 0))
ptm <- proc.time()  # Start timing
for (i in 1:length(TAs)) {
    enviro_par <- make_enviropar(replace = list(RH = set_units(RHs[i]/100,
        "1"), S_sw = set_units(QSOLRs[i], "W/m^2"), T_air = set_units(TAs[i] +
        273.15, "K"), T_sky = set_units(TSKYs[i] + 273.15, "K"),
        wind = set_units(VELs[i], "m/s"), P = set_units(P_a[i]/1000,
            "kPa"), r = set_units(micro$REF, "1")))
    leaf_par <- make_leafpar(replace = list(g_sw = set_units((g_vs_abs[i] +
        g_vs_ads[i]) * 1e+06/P_a[i], "umol/m^2/s/Pa"), g_uw = set_units(g_vs_base *
        1e+06/P_a[i], "umol/m^2/s/Pa"), abs_l = set_units(epsilon_L,
        "1"), abs_s = set_units(alpha_L, "1"), leafsize = set_units(GEOM.out$AL,
        "m")))
    T_leaf <- tleaf(leaf_par, enviro_par, constants, quiet = TRUE)
    T_leaf_tealeaves[i, ] <- unlist(T_leaf)
    pctdone <- round(i/length(TAs) * 100, 1)
    if (i %in% pctdones) {
        # cat(round(i/length(TAs)*100,1), 'percent done
        # \n')
    }
}
print(proc.time() - ptm)  # Stop the clock
##    user  system elapsed 
##   40.22    0.99   41.99
colnames(T_leaf_tealeaves) <- colnames(T_leaf)
T_leaf_tealeaves <- as.data.frame(T_leaf_tealeaves)$T_leaf -
    273.15

# plot
plot(seq(0, 23), T_air_1cm[sub], type = "l", col = "blue", ylab = "temperature, deg C",
    xlab = "hour of day", ylim = c(15, 50))
points(seq(0, 23), T_air_1.2m[sub], type = "l", col = "blue",
    lty = 2)
points(seq(0, 23), T_leaf_NMR[sub], type = "l", ylab = "T_leaf, deg C")
points(seq(0, 23), T_leaf_tealeaves[sub], type = "l", col = 2,
    lty = 2)

Compare all the outputs

Create Figure 1 from the manuscript.

par(mfrow = c(1, 1))
month <- 1  # choose a month to plot
sub <- which(floor(dates) + 1 == month)
plot(seq(0, 23), T_air_1cm[sub], type = "l", col = "blue", ylab = "temperature, deg C",
    xlab = "hour of day", ylim = c(15, 50))
points(seq(0, 23), T_air_1.2m[sub], type = "l", col = "blue",
    lty = 2)
points(seq(0, 23), T_leaf_NMR[sub], type = "l")
points(seq(0, 23), T_leaf_Campbell[sub], col = "darkgreen", lty = 2,
    type = "l")
points(seq(0, 23), T_leaf_ecophys[sub], col = "brown", lty = 1,
    type = "l")
points(seq(0, 23), T_leaf_tealeaves[sub], type = "l", col = "darkgreen")
legend(0, 50, legend = c("T_air 1.2m", "T_air 1cm", "T_leaf NMR",
    "T_leaf C&N", "T_leaf ecophys", "T_leaf TL"), lty = c(2,
    1, 1, 2, 1, 1), col = c("blue", "blue", "black", "darkgreen",
    "brown", "darkgreen"), bty = "n")

Now compare the different models for six different months.

months <- c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul",
    "Aug", "Sep", "Oct", "Nov", "Dec")
par(mfrow = c(1, 1))
for (i in seq(1, 12, 2)) {
    month <- i  # choose a month to plot
    sub <- which(floor(dates) + 1 == month)
    plot(seq(0, 23), T_air_1cm[sub], type = "l", col = "blue",
        ylab = "temperature, deg C", xlab = "hour of day", ylim = c(0,
            50), main = months[i])
    points(seq(0, 23), T_air_1.2m[sub], type = "l", col = "blue",
        lty = 2)
    points(seq(0, 23), T_leaf_NMR[sub], type = "l")
    points(seq(0, 23), T_leaf_Campbell[sub], col = "darkgreen",
        lty = 2, type = "l")
    points(seq(0, 23), T_leaf_ecophys[sub], col = "brown", lty = 1,
        type = "l")
    points(seq(0, 23), T_leaf_tealeaves[sub], type = "l", col = "darkgreen")
    legend(0, 50, legend = c("T_air 1.2m", "T_air 1cm", "T_leaf NMR",
        "T_leaf C&N", "T_leaf ecophys", "T_leaf TL"), lty = c(2,
        1, 1, 2, 1, 1), col = c("blue", "blue", "black", "darkgreen",
        "brown", "darkgreen"), bty = "n")
}

Incorporating stomatal dynamics with the plantecophys package

In this section we use the ‘Photosyn’ function of the plantecophys package to obtain more realistic stomatal dynamics. This code runs Photosyn assuming the leaf temperature calculated with the ectotherm model to start with.

# now include simulated stomatal conductance from plantecophys
photo_out <- Photosyn(Tleaf = T_leaf_NMR,
                      g1 = 5.2, # Parameter of Ball-Berry type stomatal conductance models
                      Ca = 400, # Atmospheric CO2 concentration (ppm)
                      VPD = VPD_Pa/1000, # Vapour pressure deficit (kPa)
                      vpdmin = 0.5, # lower threshold on VPD
                      PPFD = QSOLRs * UMOLPERJ * 1, # Photosynthetic photon flux density ('PAR') (mu mol m-2 s-1)
                      Patm = P_a/1000 # Atmospheric pressure (kPa)
                      )
g_vs_abs_plantecophys <- photo_out$GS

Plot the stomatal conductance values for the chosen month.

month <- 1  # choose a month to plot
sub <- which(floor(dates) + 1 == month)
par(mfrow = c(1, 1))
plot(seq(0, 23), g_vs_abs_plantecophys[sub], type = "l", ylim = c(0,
    g_vs_ab_max), ylab = "stomatal conductance, mol/m2/s", xlab = "hour of day")
abline(h = g_vs_ab_max, lty = 2, col = 3)

Now recalculate leaf temperature using these stomatal conductance values from Photosyn.

ptm <- proc.time()  # Start timing
ecto <- ectotherm(leaf = leaf, live = live, Ww_g = Ww_g, shape = shape,
    shape_b = shape_b, shape_c = shape_c, g_vs_ab = g_vs_abs_plantecophys,
    g_vs_ad = g_vs_ads, g_vs_ab_max = g_vs_ab_max, g_vs_ad_max = g_vs_ad_max,
    epsilon_sub = epsilon_sub, epsilon_sky = epsilon_sky, epsilon = epsilon_L,
    alpha_max = alpha_L, alpha_min = alpha_L, M_1 = 0, fatosk = fatosk,
    fatosb = fatosb, conv_enhance = conv_enhance, pct_cond = pct_cond,
    postur = postur, preshr = P_a, PDIF = PDIFs)
print(proc.time() - ptm)  # Stop the clock
##    user  system elapsed 
##    0.02    0.00    0.01
environ <- as.data.frame(ecto$environ)
T_leaf_NMR_plantecophys <- environ$TC

Compare the two calculations of leaf temperature. Stomatal closure in the middle of the day has caused leaf temperature to rise substantially.

plot(seq(0, 23), T_leaf_NMR_plantecophys[sub], type = "l", lty = 2,
    ylab = "temperature, deg C", xlab = "hour of day", ylim = c(15,
        50))
points(seq(0, 23), T_leaf_NMR[sub], type = "l")
abline(h = 45, lty = 2, col = 2)
abline(h = 50, lty = 1, col = 2)

But, the problem is, we need to know the leaf temperature to determine stomatal behaviour and vice versa. So now we do a loop feeding the results back and forth between the two functions until we get convergence.

# now do a loop of n iterations to see when convergence
# occurs

n <- 6  # iterations to do to find steady state stomatal conductance
par(mfrow = c(2, 3))
g_vs_abs_iterate <- g_vs_abs
T_leaf_NMR_iterate <- NULL
for (i in 1:n) {
    ecto <- ectotherm(leaf = leaf, live = live, Ww_g = Ww_g,
        shape = shape, shape_b = shape_b, shape_c = shape_c,
        g_vs_ab = g_vs_abs_iterate, g_vs_ad = g_vs_ads, g_vs_ab_max = g_vs_ab_max,
        g_vs_ad_max = g_vs_ad_max, epsilon_sub = epsilon_sub,
        epsilon_sky = epsilon_sky, epsilon = epsilon_L, alpha_max = alpha_L,
        alpha_min = alpha_L, M_1 = 0, fatosk = fatosk, fatosb = fatosb,
        conv_enhance = conv_enhance, pct_cond = pct_cond, postur = postur,
        preshr = P_a, PDIF = PDIFs)
    environ <- as.data.frame(ecto$environ)
    if (i > 1) {
        T_leaf_prev <- T_leaf_NMR_iterate
    }
    T_leaf_NMR_iterate <- environ$TC
    photo_out <- Photosyn(Tleaf = T_leaf_NMR_iterate, g1 = 5.2,
        Ca = 400, VPD = VPD_Pa/1000, vpdmin = 0.5, PPFD = QSOLRs *
            UMOLPERJ * 1, Patm = P_a/1000)
    g_vs_abs_iterate <- photo_out$GS
    plot(seq(0, 23), T_leaf_NMR_iterate[sub], type = "l", col = 1,
        ylim = c(15, 55), main = i, ylab = "leaf temperature, deg C",
        xlab = "hour of day")
    points(seq(0, 23), T_air_1cm[sub], type = "l", col = "blue")
    if (i > 1) {
        points(seq(0, 23), T_leaf_prev[sub], type = "l", lty = 2,
            col = "grey")
    }
    abline(h = 45, lty = 2, col = 2)
    abline(h = 50, lty = 1, col = 2)
}

Photosynthesis consequences

We can look at the consequences for photosynthesis as predicted by the ‘Photosyn’ function. First we calculate and plot different aspects of metabolism with our finalised leaf temperature; the rubisco-limited photosynthesis, the RUBP-limited photosynthesis, the respiration and the net photosynthesis.

# plot photosynthesis outputs
par(mfrow = c(1, 1))
A_j <- photo_out$Aj
A_c <- photo_out$Ac
R_d <- photo_out$Rd
A_n <- photo_out$ALEAF  # A_n = pmin(A_c, A_j) - R_d
plot(seq(0, 23), A_j[sub], type = "l", ylab = "mu mol m-2 s-1",
    main = "computed leaf temperature", xlab = "hour of day",
    ylim = c(0, 25), col = "grey")
points(seq(0, 23), A_c[sub], type = "l", lty = 2)
points(seq(0, 23), R_d[sub], type = "l", ylab = "mu mol m-2 s-1",
    col = "brown", lty = 2)
points(seq(0, 23), A_n[sub], type = "l", col = "red")
legend(2, 25, c("Rubisco-limited", "RUBP-limited", "Net", "respiration"),
    col = c("grey", "black", "red", "brown"), lty = c(1, 2, 1,
        2), bty = "n")

Now we can see what we’d have predicted if we assumed leaf temperature was equal to air temperature.

# redo with T_leaf = T_air
photo_out_T_air <- Photosyn(Tleaf = T_air_1cm, g1 = 5.2, Ca = 400,
    VPD = VPD_Pa/1000, vpdmin = 0.5, PPFD = QSOLRs * UMOLPERJ *
        1, Patm = P_a/1000)

# plot photosynthesis outputs
A_j2 <- photo_out_T_air$Aj
A_c2 <- photo_out_T_air$Ac
R_d2 <- photo_out_T_air$Rd
A_n2 <- photo_out_T_air$ALEAF
plot(seq(0, 23), A_j2[sub], type = "l", ylab = "mu mol m-2 s-1",
    main = "leaf temperature = air temperature", xlab = "hour of day",
    ylim = c(0, 25), col = "grey")
points(seq(0, 23), A_c2[sub], type = "l", lty = 2)
points(seq(0, 23), R_d2[sub], type = "l", ylab = "mu mol m-2 s-1",
    col = "brown", lty = 2)
points(seq(0, 23), A_n2[sub], type = "l", col = "red")
legend(2, 25, c("Rubisco-limited", "RUBP-limited", "Net", "respiration"),
    col = c("grey", "black", "red", "brown"), lty = c(1, 2, 1,
        2), bty = "n")

Finally, we can plot the difference in net photosynthesis.

# compare net photosynthesis using calculated leaf
# temperature or air temperature
plot(seq(0, 23), A_n[sub], type = "l", ylab = "mu mol m-2 s-1",
    main = "leaf vs air temperature", xlab = "hour of day", ylim = c(-2,
        15))
points(seq(0, 23), A_n2[sub], type = "l", col = "red", lty = 2)
legend(2, 12, c("calculated T_leaf", "T_leaf = T_air"), col = c("black",
    "red"), lty = c(1, 2), bty = "n")

Soil moisture stress and stomatal closure

In this last example application, we use the SPAC model embedded in the microclimate model to assess when stomata would be simulated to close due to reaching a critical leaf water potential (due to the risk of cavitation, i.e. bubbles in the xylem).

The SPAC model is in the microclimate model to properly account for transpiration when simulating soil moisture. It includes the parameters ‘PC’, the critical leaf water potential for stomatal closure, and ‘SP’, the stability parameter for stomatal closure. The idea is to drive our simulated leaf’s stomatal closure by what the SPAC model did in response to the calculated soil moisture during the microclimate simulation.

Simulate microclimate in the Mojave Desert

For this example we’ll simulate the same leaf but in the Mojave Desert, Palm Springs, using the ‘micro_usa’ function. The ‘micro_global’ function uses long-term average monthly climate whereas the ‘micro_usa’ function uses daily historical data via the gridMET database (Abatzoglou 2011). The ‘micro_usa’ function draws data from the web via opendap. We specifically choose the year 2017 (with 2016 as a ‘spin-up’ year) when soil moisture becomes sufficiently low in the simulation to drive stomatal closure.

Note that in this example we are explicitly stating the values of ‘PC’ and ‘SP’ (these are what are used by default in the microclimate model but you can vary them to simulate plants with varying resistance to low soil water potential).

loc <- c(-116.53, 33.83)
Usrhyt <- 0.3  # leaf height
PC <- -1500  # critical leaf water potential for stomatal closure, J/kg
SP <- 10  # stability parameter for stomatal closure, -

micro <- micro_usa(loc = loc, Usrhyt = Usrhyt, PC = PC, SP = SP,
    dstart = "01/01/2016", dfinish = "31/12/2017")
## extracting elevation via opendaps 
## extracting weather data 
## reading weather input for 01/01/2016 to 31/12/2017 
## tmin weather input 
## tmax weather input 
## rhmin weather input 
## rhmax weather input 
## rain weather input 
## solar weather input 
## wind weather input 
## running micro_global to get clear sky solar 
## running microclimate model for 731 days from 01/01/2016  to  31/12/2017  at site  long -116.53 lat 33.83 
##    user  system elapsed 
##   36.57    0.20   36.88
# define time zone
tz <- paste0("Etc/GMT-", floor(micro$longlat[1]/15 * -1))
dates <- micro$dates
attr(dates, "tzone") <- tz
dates2 <- micro$dates2

Now extract the aboveground microclimate, the soil water potential predictions and the plant predictions (from the SPAC model) from the ‘micro’ object.

metout <- as.data.frame(micro$metout)
soilpot <- as.data.frame(micro$soilpot)
plant <- as.data.frame(micro$plant)

Compute stomatal dynamics with Photosyn

The ‘plant’ object contains the hourly predictions of transpiration rate, leaf water potential and root water potentials (for each depth in the soil). We need the leaf water potential to reconstruct the scaling factor (0 to 1) used in the SPAC model (see eq. 15 in Appendix 1). We will use this scaling factor as a metric of stomatal openness.

# calculate stomata opennes as a function of leaf water
# potential
f_open <- 1/(1 + (plant$LPT/PC)^SP)
plot(dates, f_open, type = "l", ylab = "fraction open", xlab = "",
    main = "stomatal closure due to water stress")

Next we run the Photosyn function to obtain the photosynthesis-driven stomatal dynamics, as in the previous section (note we’ve also redefined g_vs_ab and g_vs_ab_max to have a higher base level than before).

g_vs_ab <- 0.4  # leaf vapour conductance, abaxial (bottom of leaf), mol/m2/s
g_vs_ab_max <- 0.4  # leaf vapour conductance, abaxial (bottom of leaf), mol/m2/s

# use Photosyn to get stomatal dynamics
P_a <- get_pressure(micro$elev)  # air pressure, Pa
P_a <- rep(P_a, nrow(metout))  # create hourly vector
TAs <- metout$TALOC
RHs <- metout$RHLOC
QSOLRs <- metout$SOLR
WETAIR_out <- WETAIR(db = TAs, rh = RHs, bp = P_a)
VPD_Pa <- WETAIR_out$esat - WETAIR_out$e
UMOLPERJ <- 4.57  # mu_mol photons / J

photo_out <- Photosyn(Tleaf = TAs, g1 = 5.2, Ca = 400, VPD = VPD_Pa/1000,
    vpdmin = 0.5, PPFD = QSOLRs * UMOLPERJ * 1, Patm = P_a/1000)
if (g_vs_ab > 0) {
    g_vs_abs <- photo_out$GS
} else {
    g_vs_abs <- photo_out$GS * 0
}
if (g_vs_ad > 0) {
    g_vs_ads <- photo_out$GS
} else {
    g_vs_ads <- photo_out$GS * 0
}
g_vs_abs_photosyn <- g_vs_abs  # save for comparison plot
sub <- which(as.numeric(format(dates, "%Y")) >= 2017)
plot(dates[sub], g_vs_abs[sub], type = "l", ylab = "stomatal conductance, mol/m2/s",
    xlab = "time")

Here’s what the leave temperatures would be if the plant wasn’t responding to soil moisture stress.

ptm <- proc.time()  # Start timing
ecto <- ectotherm(leaf = leaf, live = live, Ww_g = Ww_g, shape = shape,
    shape_b = shape_b, shape_c = shape_c, g_vs_ab = g_vs_abs,
    g_vs_ad = g_vs_ads, g_vs_ab_max = g_vs_ab_max, g_vs_ad_max = g_vs_ad_max,
    epsilon_sub = epsilon_sub, epsilon_sky = epsilon_sky, epsilon = epsilon_L,
    alpha_max = alpha_L, alpha_min = alpha_L, M_1 = 0, fatosk = fatosk,
    fatosb = fatosb, conv_enhance = conv_enhance, pct_cond = pct_cond,
    postur = postur, preshr = P_a, PDIF = PDIFs)
print(proc.time() - ptm)  # Stop the clock
##    user  system elapsed 
##    1.38    0.08    1.46
environ <- as.data.frame(ecto$environ)
T_leaf_photosyn <- environ$TC
plot(dates[sub], T_leaf_photosyn[sub], type = "l", ylim = c(-5,
    70), ylab = "leaf temperature, deg C", xlab = "time", col = "blue")

Impose effects of water stress on stomatal dynamics

Now we can adjust this so that the photosynthesis dynamics are overridden by the dynamics in response to leaf water stress.

# impose water stress effects
g_vs_abs[g_vs_ab * f_open < g_vs_abs] <- g_vs_ab * f_open[g_vs_ab *
    f_open < g_vs_abs]
g_vs_ads[g_vs_ad * f_open < g_vs_ads] <- g_vs_ad * f_open[g_vs_ad *
    f_open < g_vs_ads]

plot(dates[sub], g_vs_abs_photosyn[sub], type = "l", ylab = "stomatal conductance, mol/m2/s",
    xlab = "time")
points(dates[sub], g_vs_abs[sub], type = "l", ylab = "stomatal conductance, mol/m2/s",
    xlab = "time", col = "grey")

Now we can run the ectotherm model with these adjusted stomatal dynamics.

ptm <- proc.time()  # Start timing
ecto <- ectotherm(leaf = leaf, live = live, Ww_g = Ww_g, shape = shape,
    shape_b = shape_b, shape_c = shape_c, g_vs_ab = g_vs_abs,
    g_vs_ad = g_vs_ads, g_vs_ab_max = g_vs_ab_max, g_vs_ad_max = g_vs_ad_max,
    epsilon_sub = epsilon_sub, epsilon_sky = epsilon_sky, epsilon = epsilon_L,
    alpha_max = alpha_L, alpha_min = alpha_L, M_1 = 0, fatosk = fatosk,
    fatosb = fatosb, conv_enhance = conv_enhance, pct_cond = pct_cond,
    postur = postur, preshr = P_a, PDIF = PDIFs)
print(proc.time() - ptm)  # Stop the clock
##    user  system elapsed 
##    1.34    0.05    1.39
environ <- as.data.frame(ecto$environ)
T_leaf <- environ$TC
# plot(dates[sub], T_leaf[sub], type = 'l', ylim = c(-5,
# 70), ylab = 'leaf #temperature, deg C', xlab = 'time',
# col = 'red') points(dates[sub], T_leaf_photosyn[sub],
# type = 'l', col = 'blue')

Now recalculate, feeding leaf temperature output into Photosyn (ideally this would be iterated a few more times).

photo_out <- Photosyn(Tleaf = T_leaf, g1 = 5.2, Ca = 400, VPD = VPD_Pa/1000,
    vpdmin = 0.5, PPFD = QSOLRs * UMOLPERJ * 1, Patm = P_a/1000)
if (g_vs_ab > 0) {
    g_vs_abs <- photo_out$GS
} else {
    g_vs_abs <- photo_out$GS * 0
}
if (g_vs_ad > 0) {
    g_vs_ads <- photo_out$GS
} else {
    g_vs_ads <- photo_out$GS * 0
}
g_vs_abs[g_vs_abs * f_open < g_vs_abs] <- g_vs_abs[g_vs_abs *
    f_open < g_vs_abs] * f_open[g_vs_abs * f_open < g_vs_abs]
g_vs_ads[g_vs_ads * f_open < g_vs_ads] <- g_vs_ads[g_vs_ads *
    f_open < g_vs_ads] * f_open[g_vs_ads * f_open < g_vs_ads]

ptm <- proc.time()  # Start timing
ecto <- ectotherm(leaf = leaf, live = live, Ww_g = Ww_g, shape = shape,
    shape_b = shape_b, shape_c = shape_c, g_vs_ab = g_vs_abs,
    g_vs_ad = g_vs_ads, g_vs_ab_max = g_vs_ab_max, g_vs_ad_max = g_vs_ad_max,
    epsilon_sub = epsilon_sub, epsilon_sky = epsilon_sky, epsilon = epsilon_L,
    alpha_max = alpha_L, alpha_min = alpha_L, M_1 = 0, fatosk = fatosk,
    fatosb = fatosb, conv_enhance = conv_enhance, pct_cond = pct_cond,
    postur = postur, preshr = P_a, PDIF = PDIFs)
print(proc.time() - ptm)  # Stop the clock
##    user  system elapsed 
##    1.66    0.10    1.75
environ <- as.data.frame(ecto$environ)
T_leaf <- environ$TC
plot(dates[sub], T_leaf[sub], type = "l", ylim = c(-5, 70), ylab = "leaf temperature, deg C",
    xlab = "time", col = "red")
points(dates[sub], T_leaf_photosyn[sub], type = "l", col = "blue")

You can see that leaves are getting hotter when accounting for stomatal closure due to water stress. Plotting the differences between them reveals that the magnitude of this effect is around 1 or 2 \(^\circ\)C.

plot(dates[sub], T_leaf[sub] - T_leaf_photosyn[sub], type = "l",
    ylab = "difference, deg C", xlab = "time")

References

Abatzoglou, J. T. 2011. Development of gridded surface meteorological data for ecological applications and modelling. International Journal of Climatology 33:121–131.

Campbell, G. S., and J. M. Norman. 1998. Environmental Biophysics. Springer, New York.

Duursma, R. A. 2015. Plantecophys - An R Package for Analysing and Modelling Leaf Gas Exchange Data. PLOS ONE 10:e0143346.

Kearney, M. R., and W. P. Porter. 2020. NicheMapR – an R package for biophysical modelling: the ectotherm and Dynamic Energy Budget models. Ecography 43:85–96.

Maclean, I. M. D., J. R. Mosedale, and J. J. Bennie. 2019. Microclima: An r package for modelling meso- and microclimate. Methods in Ecology and Evolution 10:280–290.

Muir, C. D. 2019. tealeaves: an R package for modelling leaf temperature using energy budgets. AoB PLANTS 11:plz054.