The photosynthesis package simulates photosynthetic rate given a set of leaf traits and environmental conditions by solving the Farquhar-von Caemmerer-Berry C3 biochemical model. There are two main steps to using photosynthesis:
photo and photosynthesis
for single and multiple parameter sets, respectively).In this vignette, I’ll show you how to:
You can use the default parameter settings and simulate
photosynthetic rate in a single leaf using the make_*()
functions and photo().
library(dplyr)
library(magrittr)
library(photosynthesis)
# Leaving the make_* functions empty will automatically default to defaults
# parameters.
bake_par = make_bakepar() # temperature response parameters
constants = make_constants(use_tealeaves = FALSE) # physical constants
leaf_par = make_leafpar(use_tealeaves = FALSE) # leaf parameters
enviro_par = make_enviropar(use_tealeaves = FALSE) # environmental parameters
photo(leaf_par, enviro_par, bake_par, constants, quiet = TRUE,
use_tealeaves = FALSE)
#> C_chl value convergence g_tc
#> 1 258.1218 [umol/mol] -1.363227e-06 0 0.1726157 [mol/m^2/s]
#> A J_max J_max25 K_C
#> 1 27.94273 [umol/m^2/s] 200 [umol/m^2/s] 200 [umol/m^2/s] 268.3 [umol/mol]
#> K_C25 K_O K_O25 R_d
#> 1 268.3 [umol/mol] 165084.2 [umol/mol] 165084.2 [umol/mol] 2 [umol/m^2/s]
#> R_d25 T_leaf V_cmax V_cmax25 V_tpu
#> 1 2 [umol/m^2/s] 298.15 [K] 150 [umol/m^2/s] 150 [umol/m^2/s] 200 [umol/m^2/s]
#> V_tpu25 g_mc g_mc25 g_sc
#> 1 200 [umol/m^2/s] 0.4 [mol/m^2/s] 0.4 [mol/m^2/s] 0.4 [mol/m^2/s]
#> g_uc gamma_star gamma_star25 k_mc k_sc k_uc
#> 1 0.01 [mol/m^2/s] 37.9258 [umol/mol] 37.9258 [umol/mol] 1 [1] 1 [1] 1 [1]
#> leafsize phi_J theta_J C_air O P
#> 1 0.1 [m] 0.331 [1] 0.825 [1] 420 [umol/mol] 0.21 [mol/mol] 101.3246 [kPa]
#> PPFD RH wind Ds_Jmax Ds_gmc
#> 1 1500 [umol/m^2/s] 0.5 [1] 2 [m/s] 388.04 [J/K/mol] 487.29 [J/K/mol]
#> Ea_Jmax Ea_KC Ea_KO Ea_Rd
#> 1 56095.18 [J/mol] 80989.78 [J/mol] 23719.97 [J/mol] 40446.75 [J/mol]
#> Ea_Vcmax Ea_Vtpu Ea_gammastar Ea_gmc
#> 1 52245.78 [J/mol] 52245.78 [J/mol] 24459.97 [J/mol] 68901.56 [J/mol]
#> Ed_Jmax Ed_gmc D_c0 D_h0
#> 1 121244.8 [J/mol] 148788.6 [J/mol] 1.29e-05 [m^2/s] 1.9e-05 [m^2/s]
#> D_m0 D_w0 G R eT
#> 1 1.33e-05 [m^2/s] 2.12e-05 [m^2/s] 9.8 [m/s^2] 8.31446 [J/K/mol] 1.75 [1]
#> epsilon sigma C_i
#> 1 0.622 [1] 5.67e-08 [W/K^4/m^2] 350.1432 [umol/mol]You can look at default parameters settings in the manual (run
?make_parameters). These defaults are reasonable, but of
course you will probably want to use different choices and allow some
parameters to vary. Here, I’ll demonstrate how to replace a default. In
the next section, I’ll show you how to set up a gradient of parameter
values over which to solve for leaf temperature.
# Use the `replace` argument to replace defaults. This must be a named list, and
# each named element must have the proper units specified. See `?make_parameters`
# for all parameter names and proper units.
# Temperature response parameters can be updated (but we won't do that here)
bake_par = make_bakepar()
# Physical constants probably do not need to be replaced in most cases,
# that's why we call them 'constants'!
constants = make_constants(use_tealeaves = FALSE)
# First, we'll change photosynthetic photon flux density to 1000 umol / (m^2 s)
enviro_par = make_enviropar(
replace = list(
PPFD = set_units(1000, "umol/m^2/s")
), use_tealeaves = FALSE
)
# Next, we'll change stomatal conductance to 0.3 mol / m^2 / s.
leaf_par = make_leafpar(
replace = list(
g_sc = set_units(0.3, mol / m^2 / s)
), use_tealeaves = FALSE
)
photo <- photo(leaf_par, enviro_par, bake_par, constants, quiet = TRUE,
use_tealeaves = FALSE)
photo |>
select(PPFD, C_chl, A) |>
knitr::kable()| PPFD | C_chl | A |
|---|---|---|
| 1000 [umol/m^2/s] | 252.6879 [umol/mol] | 25.63373 [umol/m^2/s] |
In the previous two examples, I used the photo function
to solve for a single parameter set. In most cases, you’ll want to solve
for many parameter sets. The function photosynthesis
generalizes photo and makes it easy to solve for multiple
parameter sets using the same argument structure. All you need to do is
specify multiple values for one or more leaf or environmental parameters
and photosynthesis uses the purrr::cross
function to fit all combinations1.
# As before, use the `replace` argument to replace defaults, but this time we
# enter multiple values
bake_par = make_bakepar()
constants = make_constants(use_tealeaves = FALSE)
# First, we'll change the PPFD to 1000 and 1500 umol / (m^2 s)
enviro_par = make_enviropar(
replace = list(
PPFD = set_units(c(1000, 1500), umol / m^2 / s)
), use_tealeaves = FALSE
)
# Next, we'll change stomatal conductance to to 0.2 and 0.4 mol / m^2 / s
leaf_par = make_leafpar(
replace = list(
g_sc = set_units(c(0.2, 0.4), mol / m^2 / s)
), use_tealeaves = FALSE
)
# Now there should be 4 combinations (high and low g_sc crossed with high and low PPFD)
ph = photosynthesis(leaf_par, enviro_par, bake_par, constants,
use_tealeaves = FALSE, progress = FALSE, quiet = TRUE)
ph |>
select(g_sc, PPFD, A) |>
knitr::kable()| g_sc | PPFD | A |
|---|---|---|
| 0.2 [mol/m^2/s] | 1000 [umol/m^2/s] | 24.34817 [umol/m^2/s] |
| 0.4 [mol/m^2/s] | 1000 [umol/m^2/s] | 26.27649 [umol/m^2/s] |
| 0.2 [mol/m^2/s] | 1500 [umol/m^2/s] | 25.68222 [umol/m^2/s] |
| 0.4 [mol/m^2/s] | 1500 [umol/m^2/s] | 27.94273 [umol/m^2/s] |
It can take a little while to simulate many different parameter sets.
If you have multiple processors available, you can speed things up by
running simulations in parallel. In the photosynthesis
function, simply use the parallel = TRUE argument to
simulate in parallel. You’ll need to set up a future
plan(). See ?future::plan for more detail.
Here I’ll provide an example simulating an A − Cc
curve.
# NOTE: parallel example is not evaluated because it was causing an issue with CRAN, but you can copy-and-paste the code to run on your own machine.
library(future)
plan("multisession") # Set up plan
# We'll use the `replace` argument to enter multiple atmospheric CO2 concentrations
bake_par = make_bakepar()
constants = make_constants(use_tealeaves = FALSE)
enviro_par = make_enviropar(
replace = list(
C_air = set_units(seq(10, 2000, length.out = 20), umol / mol)
), use_tealeaves = FALSE
)
leaf_par = make_leafpar(use_tealeaves = FALSE)
ph = photosynthesis(leaf_par, enviro_par, bake_par, constants,
use_tealeaves = FALSE, progress = FALSE,
quiet = TRUE, parallel = TRUE)
# Plot C_c versus A
library(ggplot2)
## Drop units for plotting
ph %<>% mutate_if(~ is(.x, "units"), drop_units)
ggplot(ph, aes(C_chl, A)) +
geom_line(size = 2) +
xlab(expression(paste(C[chl], " [ppm]"))) +
ylab(expression(paste("A [", mu, "mol ", m^-2~s^-1, "]"))) +
theme_bw() +
NULLIn experiments, leaf temperature can be kept close to air
temperature, but in nature, leaf temperature can be quite a bit
different than air temperature in the shade depending on environmental
and leaf parameters. If use_tealeaves = TRUE,
photo() and photosynthesis() will call on the
tealeaves
package to calculate leaf temperature using an energy balance model.
# NOTE: parallel example is not evaluated because it was causing an issue with CRAN, but you can copy-and-paste the code to run on your own machine.
# You will need to set use_tealeaves = TRUE when making parameters because additional parameters are needed for tealeaves.
bake_par = make_bakepar()
constants = make_constants(use_tealeaves = TRUE)
enviro_par = make_enviropar(
replace = list(
T_air = set_units(seq(288.15, 313.15, 1), K)
), use_tealeaves = TRUE
)
leaf_par = make_leafpar(replace = list(
g_sc = set_units(c(0.2, 0.4), mol / m^2 / s)
), use_tealeaves = TRUE
)
ph = photosynthesis(leaf_par, enviro_par, bake_par, constants,
use_tealeaves = TRUE, progress = FALSE,
quiet = TRUE, parallel = TRUE)
# Plot temperature and photosynthesis
library(ggplot2)
## Drop units for plotting
ph %<>%
mutate_if(~ is(.x, "units"), drop_units) %>%
mutate(`g[s]` = ifelse(g_sc == 0.2, "low", "high"))
ggplot(ph, aes(T_air, T_leaf, color = `g[s]`)) +
geom_line(size = 2, lineend = "round") +
geom_abline(slope = 1, intercept = 0, linetype = "dotted") +
scale_color_discrete(name = expression(g[s])) +
xlab(expression(paste(T[air], " [K]"))) +
ylab(expression(paste(T[leaf], " [K]"))) +
theme_bw() +
NULL
ggplot(ph, aes(T_air, A, color = `g[s]`)) +
geom_line(size = 2, lineend = "round") +
scale_color_discrete(name = expression(g[s])) +
xlab(expression(paste(T[leaf], " [K]"))) +
ylab(expression(paste("A [", mu, "mol ", m^-2~s^-1, "]"))) +
theme_bw() +
NULLSince optimization is somewhat time-consuming, be
careful about crossing too many combinations. Use
progress = TRUE to show progress bar with estimated time
remaining.↩︎