The rpmodel package implements the P-model as described in Stocker et al. (2019) Geosci. Mod. Dev.. The main function available through the package is rpmodel(), which returns a list of quantities (see ?rpmodel) for a given set of inputs. An additional set of important functions that are used within rpmodel() are also available through this package. Usage examples are given below.

Example P-model run

Let's run the P-model, without \(J_{\text{max}}\) limitation (argument method_jmaxlim = "none"), for one point. The set of inputs, being temperature (tc), photosynthetic photon flux density (ppfd), vapour pressure deficit (vpd), ambient CO\(_2\) (co2), elevation (elv), and fraction of absorbed photosynthetically active radiation (fapar). The quantum yield efficiency parameter is provided as an argument (kphio) and corresponds to \(\widehat{\varphi_0}\) in Stocker et al. (2019) if the temperature-dependence of this parameter is ignored (argument do_ftemp_kphio = FALSE, corresponding to simulation setup 'ORG' in Stocker et al. (2019)), or to \(\widehat{c_L}\) if the temperature-dependence of the quantum yield efficiency is included (argument do_ftemp_kphio = TRUE, used in simulation setups 'BRC' and 'FULL' in Stocker et al. (2019)). By default the optional argument do_soilmstress is set to FALSE, meaning that the empirical soil moisture stress function is not included. The unit cost ratio (\(\beta\) in Stocker et al. (2019)) is given by argument beta.

To run the rpmodel() function we can do:

library(rpmodel)
out_pmodel <- rpmodel( 
  tc             = 20,           # temperature, deg C
  vpd            = 1000,         # Pa,
  co2            = 400,          # ppm,
  fapar          = 1,            # fraction  ,
  ppfd           = 300,          # mol/m2/d,
  elv            = 0,            # m.a.s.l.,
  kphio          = 0.05,         # quantum yield efficiency,
  beta           = 146,          # unit cost ratio a/b,
  method_optci   = "prentice14",
  method_jmaxlim = "none",
  do_ftemp_kphio = FALSE,
  verbose        = TRUE
  )
## Warning: Atmospheric pressure (patm) not provided. Calculating it as a
## function of elevation (elv), assuming standard atmosphere (101325 Pa at sea
## level).
print(out_pmodel)
## $ca
## [1] 40.53
## 
## $gammastar
## [1] 3.339251
## 
## $kmm
## [1] 46.09928
## 
## $ns_star
## [1] 1.125361
## 
## $chi
## [1] 0.694352
## 
## $mj
## [1] 0.7123038
## 
## $mc
## [1] 0.3340838
## 
## $ci
## [1] 28.14209
## 
## $lue
## [1] 0.4277633
## 
## $gpp
## [1] 128.329
## 
## $iwue
## [1] 7.742446
## 
## $gs
## [1] 0.8624985
## 
## $vcmax
## [1] 31.98167
## 
## $vcmax25
## [1] 50.20073
## 
## $rd
## [1] 0.50805

The function returns a list of variables (see also man page by ?rpmodel), including \(V_{\mathrm{cmax}}\), \(g_s\), and all the parameters of the photosynthesis model (\(K\), \(\Gamma^{\ast}\)), which are all internally consistent, as can be verified for… \[ c_i = c_a - A / g_s = \chi c_a \]

c_molmass <- 12.0107  # molecular mass of carbon
kphio <- 0.05         # quantum yield efficiency, value as used in the function call to rpmodel()
ppfd <- 300           # mol/m2/d, value as used in the function call to rpmodel()
fapar <- 1            # fraction, value as used in the function call to rpmodel()
print( out_pmodel$ci )
## [1] 28.14209
print( out_pmodel$ca - (out_pmodel$gpp / c_molmass) / out_pmodel$gs )
## [1] 28.14209
print( out_pmodel$ca * out_pmodel$chi )
## [1] 28.14209

Yes.

And for… \[ A = V_{\text{cmax}} \frac{c_i-\Gamma^{\ast}}{c_i + K} = \phi_0 I_{\text{abs}} \frac{c_i-\Gamma^{\ast}}{c_i + 2 \Gamma^{\ast}} = g_s (c_a - c_i) \]

print( out_pmodel$gpp / c_molmass )
## [1] 10.68456
print( out_pmodel$vcmax * (out_pmodel$ci - out_pmodel$gammastar) / (out_pmodel$ci + out_pmodel$kmm ))
## [1] 10.68456
print( out_pmodel$gs * (out_pmodel$ca - out_pmodel$ci) )
## [1] 10.68456
print( kphio * ppfd * fapar * (out_pmodel$ci - out_pmodel$gammastar) / (out_pmodel$ci + 2 * out_pmodel$gammastar ))
## [1] 10.68456

Yes.

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.

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           = 300 + rnorm(5, mean = 0, sd = 30),
  elv            = 0,         
  kphio          = 0.05,         
  beta           = 146,
  method_optci   = "prentice14",
  method_jmaxlim = "none",
  do_ftemp_kphio = FALSE,
  verbose        = FALSE
  )
print(out_ts_pmodel$gpp)
## [1] 121.6522 137.6573 130.6448 126.0823 148.5923

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 rpmode() function to each row.

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(purrr)
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           = 300 + rnorm(5, mean = 0, sd = 30)
  ) %>%
  mutate( out_pmodel = purrr::pmap(., rpmodel, 
    elv            = 0,         
    kphio          = 0.05,         
    beta           = 146,
    method_optci   = "prentice14",
    method_jmaxlim = "none",
    do_ftemp_kphio = FALSE
    ) )
print(df)
## # A tibble: 5 x 6
##      tc   vpd   co2 fapar  ppfd out_pmodel       
##   <dbl> <dbl> <dbl> <dbl> <dbl> <list>           
## 1  33.2 1008.   400     1  256. <named list [15]>
## 2  15.6 1014.   400     1  346. <named list [15]>
## 3  21.6 1081.   400     1  341. <named list [15]>
## 4  21.9  953.   400     1  275. <named list [15]>
## 5  19.9 1064.   400     1  363. <named list [15]>

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)
## # A tibble: 5 x 20
##      tc   vpd   co2 fapar  ppfd    ca gammastar   kmm ns_star   chi    mj
##   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>     <dbl> <dbl>   <dbl> <dbl> <dbl>
## 1  33.2 1008.   400     1  256.  40.5      6.52 143.    0.838 0.833 0.582
## 2  15.6 1014.   400     1  346.  40.5      2.64  31.6   1.26  0.636 0.745
## 3  21.6 1081.   400     1  341.  40.5      3.62  52.7   1.08  0.706 0.697
## 4  21.9  953.   400     1  275.  40.5      3.69  54.2   1.08  0.722 0.698
## 5  19.9 1064.   400     1  363.  40.5      3.33  45.8   1.13  0.687 0.711
## # … with 9 more variables: mc <dbl>, ci <dbl>, lue <dbl>, gpp <dbl>,
## #   iwue <dbl>, gs <dbl>, vcmax <dbl>, vcmax25 <dbl>, rd <dbl>

Auxiliary functions

A number of auxiliary functions, which are used within rpmodel(), are available (public) through the package.

Instantaneous temperature scaling

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

out_pmodel <- rpmodel( 
  tc             = 10,           # temperature, deg C
  vpd            = 1000,         # Pa,
  co2            = 400,          # ppm,
  fapar          = 1,            # fraction  ,
  ppfd           = 300,          # mol/m2/d,
  elv            = 0,            # m.a.s.l.,
  kphio          = 0.05,         # quantum yield efficiency,
  beta           = 146,          # unit cost ratio a/b,
  method_optci   = "prentice14",
  method_jmaxlim = "none",
  do_ftemp_kphio = FALSE,
  verbose        = TRUE
  )
## Warning: 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("calc_ftemp_inst_vcmax(10):", calc_ftemp_inst_vcmax(10)))
## [1] "calc_ftemp_inst_vcmax(10): 0.260975632963417"
print(paste("calc_ftemp_inst_rd(10):", calc_ftemp_inst_rd(10)))
## [1] "calc_ftemp_inst_rd(10): 0.284933345928884"

Parameters in the FvCB model

print(paste("From rpmodel call :", out_pmodel$gammastar))
## [1] "From rpmodel call : 1.93016150706341"
print(paste("calc_gammastar(10):", calc_gammastar(10, patm = calc_patm(elv = 0))))
## [1] "calc_gammastar(10): 1.93016150706341"
print(paste("From rpmodel call:", out_pmodel$kmm))
## [1] "From rpmodel call: 19.6242174524746"
print(paste("calc_kmm(10)     :", calc_kmm(10, patm = calc_patm(elv = 0))))
## [1] "calc_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 calc_ftemp_kphio().

out_pmodel_ftemp_kphio_ON <- rpmodel( 
  tc             = 20,           # temperature, deg C
  vpd            = 1000,         # Pa,
  co2            = 400,          # ppm,
  fapar          = 1,            # fraction  ,
  ppfd           = 300,          # 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           = 300,          # 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 = 1.08933333333333"
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.08933333333333"
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.08933333333333"
print(paste("calc_ftemp_kphio(20) =", calc_ftemp_kphio(20)))
## [1] "calc_ftemp_kphio(20) = 0.656"

Soil moisture stress

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           = 300,          # mol/m2/d,
  elv            = 0,            # m.a.s.l.,
  do_ftemp_kphio = FALSE,
  do_soilmstress = TRUE
  )
out_pmodel_soilmstress_ON <- rpmodel( 
  tc             = 20,           # temperature, deg C
  vpd            = 1000,         # Pa,
  co2            = 400,          # ppm,
  fapar          = 1,            # fraction  ,
  ppfd           = 300,          # 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 = 0.662222222222222"
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("calc_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] "calc_soilmstress(0.2, apar_soilm = 0.1, bpar_soilm = 0.7, meanalpha = 0.2) = 0.662222222222222"

calc_ftemp_arrh() Calculates the Arrhenius-type temperature response.