Running the P-model in R

Elevation and pressure

Above, atmospheric pressure (patm) was not provided as an argument, but elevation (elv) was. Hence the warning was printed (only when verbose = TRUE), saying: Atmospheric pressure (patm) not provided. Calculating it as a function of elevation (elv), Assuming standard atmosphere (101325 Pa at sea level).. Alternatively, we can provide atmospheric pressure (patm) as input, which overrides the argument elv.

P-model for time series

The rpmodel() function can also be invoked for time series, where tc, vpd, co2, fapar, patm, and ppfd are vectors.

set.seed(1982)
out_ts_pmodel <- rpmodel( 
  tc             = 20 + rnorm(5, mean = 0, sd = 5),
  vpd            = 1000 + rnorm(5, mean = 0, sd = 50),
  co2            = rep(400, 5),
  fapar          = rep(1, 5),
  ppfd           = 30 + rnorm(5, mean = 0, sd = 3),
  elv            = 0,         
  kphio          = 0.049977,
  beta           = 146,
  c4             = FALSE,
  method_jmaxlim = "none",
  do_ftemp_kphio = TRUE,
  do_soilmstress = FALSE,
  verbose        = FALSE
  )

## Warning in sqrt((1/fact_jmaxlim)^2 - 1): NaNs produced

print(out_ts_pmodel$gpp)

## [1] 8.162522 7.967237 8.148933 7.468270 8.968635

Note that gpp (as well as all other returned variables) are now vectors of the same length as the vectors provided as inputs.

P-model in the tidyverse

We can create a data frame (in tidyverse this is a tibble) and apply the rpmodel() function to each row.

library(dplyr)
library(purrr)

set.seed(1982)

df <- tibble(
  tc             = 20 + rnorm(5, mean = 0, sd = 5),
  vpd            = 1000 + rnorm(5, mean = 0, sd = 50),
  co2            = rep(400, 5),
  fapar          = rep(1, 5),
  ppfd           = 30 + rnorm(5, mean = 0, sd = 3)
  ) %>%
  mutate( out_pmodel = purrr::pmap(., rpmodel, 
    elv            = 0,         
    kphio          = 0.049977,
    beta           = 146,
    c4             = FALSE,
    method_jmaxlim = "none",
    do_ftemp_kphio = FALSE
    ) )

print(df)

Note that the new column out_pmodel now contains the list returned as output of the rpmodel() function applied to each row separately. Additional (constant) arguments are just passed to purrr::pmap as arguments.

If you prefer the elements of these lists to be in separate columns of df, use tidyr to do:

library(tidyr)
df <- df %>% 
  mutate( out_pmodel = purrr::map(out_pmodel, ~as_tibble(.))) %>% 
  unnest(out_pmodel)
print(df)

C4 photosynthesis

Photosynthesis by C4 plants is simulated by the P-model based on the assumption that the “CO2 limitation term” in the FvCB model (Eq. 11 in Stocker et al., 2020) is 1 (no limitation), and based on a distinct parametrisation of the quantum yield efficiency and its temperature dependence \(\varphi_0\), following Cai & Prentice, 2020. Under identical conditions, GPP of C3 and C4 photosynthesis are different:

out_c3 <- rpmodel( 
  tc             = 20,           # temperature, deg C
  vpd            = 1000,         # Pa,
  co2            = 400,          # ppm,
  fapar          = 1,            # fraction  ,
  ppfd           = 30,           # mol/m2/d,
  elv            = 0,            # m.a.s.l.,
  c4             = FALSE
  )
out_c4 <- rpmodel( 
  tc             = 20,           # temperature, deg C
  vpd            = 1000,         # Pa,
  co2            = 400,          # ppm,
  fapar          = 1,            # fraction  ,
  ppfd           = 30,           # mol/m2/d,
  elv            = 0,            # m.a.s.l.,
  c4             = TRUE
  )

## Warning in rpmodel(tc = 20, vpd = 1000, co2 = 400, fapar = 1, ppfd = 30, : rpmodel(): light and Rubisco-limited assimilation rates are not identical.
## Mean relative difference: 1.582923

## Warning in rpmodel(tc = 20, vpd = 1000, co2 = 400, fapar = 1, ppfd = 30, : rpmodel(): Assimilation and GPP are not identical.
## Mean relative difference: 1.582923

print(out_c3$gpp)

## [1] 7.642545

print(out_c4$gpp)

## [1] 126.2565

To accurately simulate C3 photosynthesis, a constant scalar was calibrated by Stocker et al., 2020 to GPP from FLUXNET2015 data and scaled in their ‘BRC’ and ‘FULL’ setup with a temperature dependence factor (their Eq. 20). The implementation of C3 photosynthesis in rpmodel uses kphio as a constant scalar, provided to rpmodel() as an argument, and representing \((a_L b_L)/4\) in their Eq. 20. The temprature dependence is calculated based on Bernacchi et al., 2003 using the function ftemp_kphio().

The implementation of C4 photosynthesis uses a different temperature dependence of \(\varphi_0\) following Eq. 5 in Cai & Prentice (2020), implemented by ftemp_kphio(c4 = TRUE). The calibratable parameter kphio doesn’t appear explicitly in Eq. 5 by Cai & Prentice (2020), but is implemented the same way in rpmodel as for C3 photosynthesis. It should therefore be regarded as a “correction” factor, where kphio = 1.0 means “no correction” and reflects the parametrisation used by Cai & Prentice (2020).

The following shows \(\varphi(T)\), corresponding to kphio * ftemp_kphio(), using parameters as in Stocker et al., 2020 for C3 vegetation (their ‘BRC’ setup, black line in the plot below), and as in Cai & Prentice (2020) for C4 vegetation (blue solid line).

Addendum: As of issue #19 (raised by David Orme), the values provided in Cai & Prentice (2020) are erroneous. And instead, the temperature response function for C4 should be ftemp = -0.064 + 0.03 * tc - 0.000464 * tc^2. This is now also implemented in rpmodel.

Instantaneous temperature scaling

Different instantaneous temperature scaling functions are applied for \(V_\text{cmax}\) and dark respiration (\(R_d\)).

• ftemp_inst_vcmax() calculates the instantaneous temperature response of \(V_\text{cmax}\). Let’s run the P-model for tc = 10 (degrees C). The ratio of \(V_\text{cmax}/V_\text{cmax25}\) should equal the instantaneous temperature scaling function for \(V_\text{cmax}\) at 10 degrees C (calculated by ftemp_inst_vcmax(10)):

out_pmodel <- rpmodel( 
  tc             = 10,           # temperature, deg C
  vpd            = 1000,         # Pa,
  co2            = 400,          # ppm,
  fapar          = 1,            # fraction  ,
  ppfd           = 30,           # mol/m2/d,
  elv            = 0,            # m.a.s.l.,
  kphio          = 0.049977,     # quantum yield efficiency as calibrated for setup ORG by Stocker et al. 2020 GMD,
  beta           = 146,          # unit cost ratio a/b,
  method_jmaxlim = "none",
  do_ftemp_kphio = FALSE,
  verbose        = TRUE
  )

## Warning in rpmodel(tc = 10, vpd = 1000, co2 = 400, fapar = 1, ppfd = 30, : Atmospheric pressure (patm) not provided. Calculating it as a
##       function of elevation (elv), assuming standard atmosphere 
##       (101325 Pa at sea level).

print(paste("Ratio Vcmax/Vcmax25      :", out_pmodel$vcmax/out_pmodel$vcmax25))

## [1] "Ratio Vcmax/Vcmax25      : 0.260975632963417"

print(paste("ftemp_inst_vcmax(10):", ftemp_inst_vcmax(10)))

## [1] "ftemp_inst_vcmax(10): 0.260975632963417"

• ftemp_arrh() Calculates the Arrhenius-type temperature response and is used inside ftemp_inst_vcmax().

• ftemp_inst_rd() calculates the temperature response of dark respiration (\(R_d\)), which is slightly less steep than that for \(V_\text{cmax}\):

print(paste("ftemp_inst_rd(10):", ftemp_inst_rd(10)))

## [1] "ftemp_inst_rd(10): 0.284933345928884"

Parameters in the FvCB model

• calc_gammastar() calculates the CO\(_2\) compensation point (\(\Gamma^\ast\)) in the Farquhar-von Caemmerer-Berry model as a function of temperature (argument tc) and atmospheric pressure (argument patm). This is returned by the rpmodel() function and by the separate auxiliary function calc_gammastar(). calc_gammastar() requires atmospheric pressure (patm) to be given as an argument (in addition to temperature). Corresponding to the rpmodel() call above, let’s calculate this using the auxiliary function calc_patm() with 0 metres above sea level, and assuming standard atmospheric pressure (101325 Pa at 0 m a.s.l.):

print(paste("From rpmodel call :", out_pmodel$gammastar))

## [1] "From rpmodel call : 1.93016150706341"

print(paste("gammastar(10):", calc_gammastar(10, patm = calc_patm(elv = 0))))

## [1] "gammastar(10): 1.93016150706341"

• calc_kmm() calculates the Michaelis Menten coefficient for Rubisco-limited photosynthesis as a function of temperature (argument tc) and atmospheric pressure (argument patm). As above, calc_kmm() requires atmospheric pressure to be given as an argument (in addition to temperature). Corresponding to the rpmodel() call above, let’s calculate this using the auxiliary function calc_patm() with 0 metres above sea level, and assuming standard atmospheric pressure (101325 Pa at 0 m a.s.l.):

print(paste("From rpmodel call:", out_pmodel$kmm))

## [1] "From rpmodel call: 19.6242174524746"

print(paste("kmm(10)     :", calc_kmm(10, patm = calc_patm(elv = 0))))

## [1] "kmm(10)     : 19.6242174524746"

Temperature dependence of quantum yield efficiency

The temperature dependence of quantum yield efficiency is modelled following Bernacchi et al. (2003), if the argument to the rpmodel() call do_ftemp_kphio = TRUE. This affects several quantities returned by the rpmodel() call (GPP, LUE, Vcmax), and can be calculated direction using ftemp_kphio().

out_pmodel_ftemp_kphio_ON <- rpmodel( 
  tc             = 20,           # temperature, deg C
  vpd            = 1000,         # Pa,
  co2            = 400,          # ppm,
  fapar          = 1,            # fraction  ,
  ppfd           = 30,           # mol/m2/d,
  elv            = 0,            # m.a.s.l.,
  do_ftemp_kphio = TRUE
  )
out_pmodel_ftemp_kphio_OFF <- rpmodel( 
  tc             = 20,           # temperature, deg C
  vpd            = 1000,         # Pa,
  co2            = 400,          # ppm,
  fapar          = 1,            # fraction  ,
  ppfd           = 30,           # mol/m2/d,
  elv            = 0,            # m.a.s.l.,
  do_ftemp_kphio = FALSE
  )
print(paste("LUE ftemp_ON /LUE ftemp_OFF =", out_pmodel_ftemp_kphio_ON$lue / out_pmodel_ftemp_kphio_OFF$lue))

## [1] "LUE ftemp_ON /LUE ftemp_OFF = "

print(paste("GPP ftemp_ON /GPP ftemp_OFF =", out_pmodel_ftemp_kphio_ON$gpp / out_pmodel_ftemp_kphio_OFF$gpp))

## [1] "GPP ftemp_ON /GPP ftemp_OFF = 1.07351301598735"

print(paste("Vcmax ftemp_ON /Vcmax ftemp_OFF =", out_pmodel_ftemp_kphio_ON$vcmax / out_pmodel_ftemp_kphio_OFF$vcmax))

## [1] "Vcmax ftemp_ON /Vcmax ftemp_OFF = 1.07351301598735"

print(paste("ftemp_kphio(20) =", ftemp_kphio(20)))

## [1] "ftemp_kphio(20) = 0.656"

Soil moisture stress

The soil moisture stress function is available as a separate public function.

vec_soilm <- seq(from = 1.0, to = 0.0, by = -0.05)
vec_soilmstress <- calc_soilmstress( vec_soilm, meanalpha = 1.0, apar_soilm = 0.0, bpar_soilm = 0.7330 )
plot(vec_soilm, vec_soilmstress)

Similar to above, the soil moisture dependence of LUE (and hence GPP, and Vcmax) can be calculated directly using the function calc_soilmstress() and affects several quantities returned by the rpmodel() call (GPP, LUE, Vcmax):

out_pmodel_soilmstress_OFF <- rpmodel( 
  tc             = 20,           # temperature, deg C
  vpd            = 1000,         # Pa,
  co2            = 400,          # ppm,
  fapar          = 1,            # fraction  ,
  ppfd           = 30,           # mol/m2/d,
  elv            = 0,            # m.a.s.l.,
  do_ftemp_kphio = FALSE,
  do_soilmstress = FALSE
  )
out_pmodel_soilmstress_ON <- rpmodel( 
  tc             = 20,           # temperature, deg C
  vpd            = 1000,         # Pa,
  co2            = 400,          # ppm,
  fapar          = 1,            # fraction  ,
  ppfd           = 30,           # mol/m2/d,
  elv            = 0,            # m.a.s.l.,
  do_ftemp_kphio = FALSE,
  do_soilmstress = TRUE,
  soilm          = 0.2,
  apar_soilm     = 0.1,
  bpar_soilm     = 0.7,
  meanalpha      = 0.2 
  )
print(paste("LUE soilmstress_ON /LUE soilmstress_OFF =", out_pmodel_soilmstress_ON$lue / out_pmodel_soilmstress_OFF$lue))

## [1] "LUE soilmstress_ON /LUE soilmstress_OFF = "

print(paste("GPP soilmstress_ON /GPP soilmstress_OFF =", out_pmodel_soilmstress_ON$gpp / out_pmodel_soilmstress_OFF$gpp))

## [1] "GPP soilmstress_ON /GPP soilmstress_OFF = 0.662222222222222"

print(paste("Vcmax soilmstress_ON /Vcmax soilmstress_OFF =", out_pmodel_soilmstress_ON$vcmax / out_pmodel_soilmstress_OFF$vcmax))

## [1] "Vcmax soilmstress_ON /Vcmax soilmstress_OFF = 0.662222222222222"

print(paste("soilmstress(0.2, apar_soilm = 0.1, bpar_soilm = 0.7, meanalpha = 0.2) =", calc_soilmstress(0.2, apar_soilm = 0.1, bpar_soilm = 0.7, meanalpha = 0.2)))

## [1] "soilmstress(0.2, apar_soilm = 0.1, bpar_soilm = 0.7, meanalpha = 0.2) = 0.662222222222222"

ftemp_arrh() Calculates the Arrhenius-type temperature response.