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.

**Important note:**

The P-model predicts how photosynthesis acclimates to a changing
environment, coordinating stomatal conductance, Vcmax and Jmax. This
yields a model that has the form of a light use efficiency model, where
gross primary production scales linearly with absorbed light, as
described in Stocker
et al. 2020. It is important to note that this implies that the
P-model is valid only for simulating responses to the environment that
evolve over the time scale at which the photosynthetic machinery (e.g.,
Rubisco) can be assumed to acclimate. Sensible choices are on the order
of a couple of weeks to a month. In other words, the arguments (climatic
forcing), provided to `rpmodel()`

should represent typical
daytime mean values, averaged across a couple of weeks. The output is
then representative also for average values across the same time
scale.

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 = 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,
c4 = FALSE,
method_jmaxlim = "wang17",
do_ftemp_kphio = FALSE, # corresponding to setup ORG
do_soilmstress = FALSE, # corresponding to setup ORG
verbose = TRUE
)
```

```
## Warning in rpmodel(tc = 20, 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).
```

```
## $gpp
## [1] 7.119192
##
## $ca
## [1] 40.53
##
## $gammastar
## [1] 3.339251
##
## $kmm
## [1] 46.09928
##
## $ns_star
## [1] 1.125361
##
## $chi
## [1] 0.694352
##
## $xi
## [1] 63.3145
##
## $mj
## [1] 0.7123038
##
## $mc
## [1] 0.3340838
##
## $ci
## [1] 28.14209
##
## $iwue
## [1] 7.742446
##
## $gs
## [1] 0.04784805
##
## $vcmax
## [1] 1.774218
##
## $vcmax25
## [1] 2.78494
##
## $jmax
## [1] 4.001452
##
## $jmax25
## [1] 5.464979
##
## $rd
## [1] 0.02818463
```

Above, we specified the model paramters (arguments `beta`

and `kphio`

). This overrides the defaults, where
`rpmodel()`

uses the parameters as calibrated by Stocker et
al. (2019), depending on the choices for arguments
`do_ftemp_kphio`

and `do_soilmstress`

:

```
kphio = ifelse(do_ftemp_kphio, ifelse(do_soilmstress, 0.087182, 0.081785), 0.049977)
beta = 146.0
apar_soilm = 0.0
bpar_soilm = 0.73300
```

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 <- 30 # 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`

`## [1] 28.14209`

`## [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) \]

`## [1] 0.5927375`

`print( out_pmodel$vcmax * (out_pmodel$ci - out_pmodel$gammastar) / (out_pmodel$ci + out_pmodel$kmm ))`

`## [1] 0.5927375`

`## [1] 0.5927375`

`print( kphio * ppfd * fapar * (out_pmodel$ci - out_pmodel$gammastar) / (out_pmodel$ci + 2 * out_pmodel$gammastar ))`

`## [1] 1.068456`

Yes.

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`

.

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`

`## [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.

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:

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
```

`## [1] 7.642545`

`## [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.

A number of auxiliary functions, which are used within
`rpmodel()`

, are available (public) through the package.

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).
```

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

`## [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}\):

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

`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.):

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

`## [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.):

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

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

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"`

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

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.

Bernacchi, C. J., Pimentel, C., and Long, S. P.: In vivo temperature response func-tions of parameters required to model RuBP-limited photosynthesis, Plant Cell Environ., 26, 1419–1430, 2003

Cai, W., and Prentice, I. C.: Recent trends in gross primary production and their drivers: analysis and modelling at flux-site and global scales, Environ. Res. Lett. 15 124050 https://doi.org/10.1088/1748-9326/abc64e, 2020

Stocker, B. D., Wang, H., Smith, N. G., Harrison, S. P., Keenan, T. F., Sandoval, D., Davis, T., and Prentice, I. C.: P-model v1.0: An optimality-based light use efficiency model for simulating ecosystem gross primary production, Geosci. Model Dev. Discuss., https://doi.org/10.5194/gmd-2019-200, in review, 2019.