Nothing
Running the P-model in R
In rpmodel: P-Model
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.
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 = 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
)
print(out_pmodel)
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 )
print( out_pmodel$ca - (out_pmodel$gpp / c_molmass) / out_pmodel$gs )
print( out_pmodel$ca * out_pmodel$chi )
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 )
print( out_pmodel$vcmax * (out_pmodel$ci - out_pmodel$gammastar) / (out_pmodel$ci + out_pmodel$kmm ))
print( out_pmodel$gs * (out_pmodel$ca - out_pmodel$ci) )
print( kphio * ppfd * fapar * (out_pmodel$ci - out_pmodel$gammastar) / (out_pmodel$ci + 2 * out_pmodel$gammastar ))
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.
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
)
print(out_ts_pmodel$gpp)
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
)
print(out_c3$gpp)
print(out_c4$gpp)
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.
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$).
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
)
print(paste("Ratio Vcmax/Vcmax25 :", out_pmodel$vcmax/out_pmodel$vcmax25))
print(paste("ftemp_inst_vcmax(10):", ftemp_inst_vcmax(10)))
-
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)))
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))
print(paste("gammastar(10):", calc_gammastar(10, patm = calc_patm(elv = 0))))
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))
print(paste("kmm(10) :", calc_kmm(10, patm = calc_patm(elv = 0))))
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))
print(paste("GPP ftemp_ON /GPP ftemp_OFF =", out_pmodel_ftemp_kphio_ON$gpp / out_pmodel_ftemp_kphio_OFF$gpp))
print(paste("Vcmax ftemp_ON /Vcmax ftemp_OFF =", out_pmodel_ftemp_kphio_ON$vcmax / out_pmodel_ftemp_kphio_OFF$vcmax))
print(paste("ftemp_kphio(20) =", ftemp_kphio(20)))
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))
print(paste("GPP soilmstress_ON /GPP soilmstress_OFF =", out_pmodel_soilmstress_ON$gpp / out_pmodel_soilmstress_OFF$gpp))
print(paste("Vcmax soilmstress_ON /Vcmax soilmstress_OFF =", out_pmodel_soilmstress_ON$vcmax / out_pmodel_soilmstress_OFF$vcmax))
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)))
ftemp_arrh() Calculates the Arrhenius-type temperature response.
References
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.
Try the rpmodel package in your browser
Any scripts or data that you put into this service are public.
rpmodel documentation built on June 28, 2024, 1:06 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
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.
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 = 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 ) print(out_pmodel)
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 ) print( out_pmodel$ca - (out_pmodel$gpp / c_molmass) / out_pmodel$gs ) print( out_pmodel$ca * out_pmodel$chi )
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 ) print( out_pmodel$vcmax * (out_pmodel$ci - out_pmodel$gammastar) / (out_pmodel$ci + out_pmodel$kmm )) print( out_pmodel$gs * (out_pmodel$ca - out_pmodel$ci) ) print( kphio * ppfd * fapar * (out_pmodel$ci - out_pmodel$gammastar) / (out_pmodel$ci + 2 * out_pmodel$gammastar ))
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.
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 ) print(out_ts_pmodel$gpp)
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 ) print(out_c3$gpp) print(out_c4$gpp)
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.
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$).
ftemp_inst_vcmax()calculates the instantaneous temperature response of $V_\text{cmax}$. Let's run the P-model fortc = 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 byftemp_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 ) print(paste("Ratio Vcmax/Vcmax25 :", out_pmodel$vcmax/out_pmodel$vcmax25)) print(paste("ftemp_inst_vcmax(10):", ftemp_inst_vcmax(10)))
-
ftemp_arrh()Calculates the Arrhenius-type temperature response and is used insideftemp_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)))
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 (argumenttc) and atmospheric pressure (argumentpatm). This is returned by therpmodel()function and by the separate auxiliary functioncalc_gammastar().calc_gammastar()requires atmospheric pressure (patm) to be given as an argument (in addition to temperature). Corresponding to therpmodel()call above, let's calculate this using the auxiliary functioncalc_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)) print(paste("gammastar(10):", calc_gammastar(10, patm = calc_patm(elv = 0))))
calc_kmm()calculates the Michaelis Menten coefficient for Rubisco-limited photosynthesis as a function of temperature (argumenttc) and atmospheric pressure (argumentpatm). As above,calc_kmm()requires atmospheric pressure to be given as an argument (in addition to temperature). Corresponding to therpmodel()call above, let's calculate this using the auxiliary functioncalc_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)) print(paste("kmm(10) :", calc_kmm(10, patm = calc_patm(elv = 0))))
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)) print(paste("GPP ftemp_ON /GPP ftemp_OFF =", out_pmodel_ftemp_kphio_ON$gpp / out_pmodel_ftemp_kphio_OFF$gpp)) print(paste("Vcmax ftemp_ON /Vcmax ftemp_OFF =", out_pmodel_ftemp_kphio_ON$vcmax / out_pmodel_ftemp_kphio_OFF$vcmax)) print(paste("ftemp_kphio(20) =", ftemp_kphio(20)))
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)) print(paste("GPP soilmstress_ON /GPP soilmstress_OFF =", out_pmodel_soilmstress_ON$gpp / out_pmodel_soilmstress_OFF$gpp)) print(paste("Vcmax soilmstress_ON /Vcmax soilmstress_OFF =", out_pmodel_soilmstress_ON$vcmax / out_pmodel_soilmstress_OFF$vcmax)) 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)))
ftemp_arrh() Calculates the Arrhenius-type temperature response.
References
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.
Try the rpmodel package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.