The conductivity of an infinitesimally short xylem element can be described as a function of the water potential (\(\Psi\)), declining to half at a level \(\Psi_{50}\).
\[
k(\Psi) = 0.5 ^ {\left( \Psi / \Psi_{50} \right) ^ b}
\] Plotted with different shape parameter values \(b\).
par_conductivity_std <- list(psi50 = -2, b = 2)
calc_conductivity <- function(psi, par){
0.5^((psi/par$psi50)^par$b)
}
ggplot( data = tibble(x = 0), aes(x = x) ) +
stat_function(fun = calc_conductivity, args = list(par = list(psi50 = -1, b = 1)), col = 'black') +
stat_function(fun = calc_conductivity, args = list(par = list(psi50 = -1, b = 2)), col = 'green') +
stat_function(fun = calc_conductivity, args = list(par = list(psi50 = -1, b = 5)), col = 'red') +
xlim( -3, 0 ) +
labs( title = "Conductivity as a function of water potential" )
The total conductance along the soil-to-leaf pathway is then given by the integral \[ \int_{\Psi_S}^{\Psi_L}k(\Psi) d \Psi \]
Plotted here for different soil water potentials.
calc_conductance = function(dpsi, psi_soil, ...){
-integrate(calc_conductivity, lower = psi_soil, upper = (psi_soil - dpsi), ...)$value
}
tibble( dpsi = seq(0, 5, length.out = 30)) %>%
mutate(total_conductance_soil1 = purrr::map_dbl(dpsi, ~calc_conductance(., psi_soil = -1, par = par_conductivity_std))) %>%
mutate(total_conductance_soil2 = purrr::map_dbl(dpsi, ~calc_conductance(., psi_soil = -2, par = par_conductivity_std))) %>%
ggplot() +
geom_line(aes( x = dpsi, y = total_conductance_soil1)) +
geom_line(aes( x = dpsi, y = total_conductance_soil2), col = 'red') +
labs( title = "Total conductance as a function of leaf water potential" )
Short term regulation of \(g_s\) at time scales of minutes to days.
This yields \(g_s^\ast\) as a function of hydraulic parameters, the “isohydricity” constant \(\Delta \Psi^\ast\), and the environment: \[ g_s^\ast = \frac{\nu_H \; k_{s0}}{1.6 H \eta D} \int_{\Psi_s}^{\Psi_s - \Delta\Psi^\ast} k(\Psi) \; d\Psi \]
calc_gs_star <- function(psi_soil, vpd, viscosity, dpsi_star, v_huber, par, par_conductivity){
(v_huber * par$conductivity_base / (1.6 * par$height * viscosity * vpd)) * calc_conductance(dpsi = dpsi_star, psi_soil = psi_soil, par = par_conductivity_std)
}
par_conductance_std <- list(conductivity_base = 1, height = 1 )
The regulated stomatal conductance \(g_s^\ast\) responds to a drying soil as shown below (plotted as a function of soil water potential).
tibble(psi_soil = seq(0, -5, length.out = 30)) %>%
mutate(gs_star_1 = purrr::map_dbl(psi_soil, ~calc_gs_star(., vpd = 100, viscosity = 1, dpsi_star = 1, v_huber = 1, par = par_conductance_std, par_conductivity = par_conductivity_std))) %>%
mutate(gs_star_2 = purrr::map_dbl(psi_soil, ~calc_gs_star(., vpd = 100, viscosity = 1, dpsi_star = 2, v_huber = 1, par = par_conductance_std, par_conductivity = par_conductivity_std))) %>%
ggplot() +
geom_line(aes(x = psi_soil, y = gs_star_1)) +
geom_line(aes(x = psi_soil, y = gs_star_2, col = 'red')) +
labs(title = "Regulated stomatal conductance", subtitle = "At two different dpsi_star levels")
With soil water potential, expressed as a function of volumetric soil water content, using (from BiomeE): \[ \Psi_s = \Psi_\text{sat} \left( \frac{W_\text{FC}}{W}\right) ^b \]
## functional form and values from BiomeE model
calc_psi_soil <- function(w_vol, par){
par$psi_soil_sat * (par$w_vol_fc / w_vol)^par$b
}
par_psi_soil_std = list(psi_soil_sat = -0.6, w_vol_fc = 0.4, b = 1) # b is 2.2 in BiomeE
… it can be plotted as a function of soil volumetric water content.
df_w_vol <- tibble( w_vol = seq(0.05, 0.6, length.out = 100)) %>%
mutate(psi_soil = calc_psi_soil(w_vol, par_psi_soil_std)) %>%
mutate(gs_star_1 = purrr::map_dbl(psi_soil, ~calc_gs_star(., vpd = 100, viscosity = 1, dpsi_star = 1, v_huber = 1, par = par_conductance_std, par_conductivity = par_conductivity_std)))
df_w_vol %>%
ggplot() +
geom_line(aes(x = w_vol, y = gs_star_1)) +
labs(title = "Regulated stomatal conductance", subtitle = "As a funtion of volumetric soil water content")
This yields a GPP response to drying soil.
calc_assim <- function(gs, vcmax, ca, par){
## Rubisco is limiting
## Solve Eq. system
## A = gs (ca- ci)
## A = Vcmax * (ci - gammastar)/(ci + Kmm)
## This leads to a quadratic equation:
## A * ci^2 + B * ci + C = 0
## 0 = a + b*x + c*x^2
## with
A <- -1.0 * gs
B <- gs * ca - gs * par$kmm - vcmax
C <- gs * ca * par$kmm + vcmax * par$gammastar
ci <- QUADM(A, B, C)
a_c <- vcmax * (ci - par$gammastar) / (ci + par$kmm)
return(a_c)
}
## Set P-model parameters
beta <- 146 # unit cost ratio a/b
c_cost <- 0.41
gamma <- 0.105 # unit cost ratio c/b
kphio <- 0.05 # quantum yield efficiency
c_molmass <- 12.0107 # molar mass, g / mol
## Define environmental conditions
tc <- 20 # temperature, deg C
ppfd <- 300 # mol/m2/d
vpd <- 1000 # Pa
co2 <- 400 # ppm
elv <- 0 # m.a.s.l.
fapar <- 1 # fraction
out_analytical <- rpmodel::rpmodel(
tc = tc,
vpd = vpd,
co2 = co2,
elv = elv,
kphio = kphio,
beta = beta,
fapar = fapar,
ppfd = ppfd,
method_optci = "prentice14",
method_jmaxlim = "none",
do_ftemp_kphio = FALSE
)
par_photosynth_std <- list(
kmm = out_analytical$kmm,
gammastar = out_analytical$gammastar
)
df_w_vol <- df_w_vol %>%
mutate(a_c = purrr::map_dbl(gs_star_1, ~calc_assim(., vcmax = out_analytical$vcmax, ca = out_analytical$ca, par = par_photosynth_std)))
df_w_vol %>%
ggplot() +
geom_line(aes(x = w_vol, y = a_c)) +
labs(title = "Assimilation", subtitle = "As a funtion of volumetric soil water content")
df_vpd <- tibble( vpd = seq(0, 1000, length.out = 100)) %>%
mutate(gs_star = purrr::map_dbl(vpd, ~calc_gs_star(-1, vpd = ., viscosity = 1, dpsi_star = 1, v_huber = 1, par = par_conductance_std, par_conductivity = par_conductivity_std))) %>%
mutate(a_c_inst = purrr::map_dbl(gs_star, ~calc_assim(., vcmax = out_analytical$vcmax, ca = out_analytical$ca, par = par_photosynth_std)))
df_vpd %>%
ggplot() +
geom_line(aes(x = vpd, y = a_c_inst)) +
labs(title = "Assimilation A", subtitle = "As a funtion of VPD")
## Warning: Removed 1 rows containing missing values (geom_path).
xxx
On time scales of weeks to months, \(V_{\mathrm{cmax}}\) and \(\nu_H\) are acclimated/regulated to satisfy some optimality criterion.
At this time scale, other hydraulic parameters are regarded as constant, while stomatal conductance is instantaneously regulated with fast-varying environmental conditions (\(g_s^\ast\)).
A least-cost (LC) principle could be formulated as \[ b V_{\mathrm{cmax}} / A + (cA_s /\nu_H)/A = min. \] The second term is the leaf construction cost which is a function of the Huber value \(\nu_H\) and the sapwood cross-sectional area \(A_s\). Note that the water cost is implicitly included in \(A\) because it’s a function of the regulated \(g_s^\ast\).
Treat \(\nu_H\) as a constant and optimise \(V_{\mathrm{cmax}}\) so that net assimilation is maximised: \[ A_C - b V_{\mathrm{cmax}} = max. \] Does it even have a maximum?
fn_target <- function(vcmax, gs_star, ca, par_cost, par_photosynth, do_optim = FALSE){
## Required parameter sets:
## par_cost : b
## par_photosynth : kmm, gammastar
## Rubisco-limited assimilation
a_c <- calc_assim(gs_star, vcmax, ca, par = par_photosynth)
## Profit Maximisation
out <- a_c - par_cost$b * vcmax
if (do_optim){
return(-out)
} else {
return(out)
}
}
optim_vcmax_midterm_constvhuber <- function(fn_target, gs_star, ca, par_cost, par_photosynth, vcmax_init, return_all = FALSE){
out_optim <- optimr::optimr(
par = vcmax_init,
lower = 0.00001 * vcmax_init,
upper = 100000 * vcmax_init,
fn = fn_target,
gs_star = gs_star,
ca = ca,
par_cost = par_cost,
par_photosynth = par_photosynth,
do_optim = TRUE,
method = "L-BFGS-B",
control = list( maxit = 100, maximize = TRUE )
)
out_optim$value <- -out_optim$value
if (return_all){
out_optim
} else {
return(out_optim$par)
}
}
df_test <- tibble( vcmax = seq(out_analytical$vcmax * 0.1, out_analytical$vcmax * 10, length.out = 30) ) %>%
mutate(assim = purrr::map_dbl(vcmax, ~calc_assim(gs = 1, vcmax = ., ca = out_analytical$ca, par = par_photosynth_std))) %>%
mutate(target = purrr::map_dbl(vcmax, ~fn_target(., gs_star = 1, ca = out_analytical$ca, par_cost = list(b = 0.05), par = par_photosynth_std)))
out_midterm_simpl <- optim_vcmax_midterm_constvhuber(fn_target, gs_star = 1, ca = out_analytical$ca, par_cost = list(b = 0.05), par_photosynth = par_photosynth_std, vcmax_init = out_analytical$vcmax, return_all = TRUE)
df_test %>%
ggplot(aes(x = vcmax, y = target)) +
geom_line() +
geom_point(aes(x = out_midterm_simpl$par, y = out_midterm_simpl$value), col = 'red')
Yes. Ok.
Btw: A least-cost criterion (\(b V_{\mathrm{cmax}} / A_C = min.\)) doesn’t work here
This allows us to predict how \(V_{\mathrm{cmax}}\) would acclimate to a drying soil.
df_w_vol <- tibble( w_vol = seq(0.05, 0.6, length.out = 100)) %>%
mutate(psi_soil = calc_psi_soil(w_vol, par_psi_soil_std)) %>%
mutate(gs_star = purrr::map_dbl(psi_soil, ~calc_gs_star(., vpd = 100, viscosity = 1, dpsi_star = 1, v_huber = 1, par = par_conductance_std, par_conductivity = par_conductivity_std))) %>%
mutate(vcmax_opt = purrr::map_dbl(gs_star, ~optim_vcmax_midterm_constvhuber(fn_target, gs_star = ., ca = out_analytical$ca, par_cost = list(b = 0.05), par_photosynth = par_photosynth_std, vcmax_init = out_analytical$vcmax))) %>%
mutate(a_c_accl = purrr::map2_dbl(gs_star, vcmax_opt, ~calc_assim(.x, vcmax = .y, ca = out_analytical$ca, par = par_photosynth_std))) %>%
mutate(a_c_inst = purrr::map_dbl(gs_star, ~calc_assim(., vcmax = 0.886016918, ca = out_analytical$ca, par = par_photosynth_std)))
df_w_vol %>%
ggplot() +
geom_line(aes(x = w_vol, y = vcmax_opt)) +
labs(title = "Acclimated Vcmax", subtitle = "As a funtion of volumetric soil water content")
df_w_vol %>%
tidyr::gather("source", "a_c", c(a_c_accl, a_c_inst)) %>%
ggplot() +
geom_line(aes(x = w_vol, y = a_c, linetype = source)) +
labs(title = "Acclimated vs. instantaneous A", subtitle = "As a funtion of volumetric soil water content")
The instantaneous dependence of \(A\) on VPD goes with \(1/D\). This is probably not reasonable. Compare this to the acclimated response.
df_vpd <- tibble( vpd = seq(0, 1000, length.out = 100)) %>%
mutate(gs_star = purrr::map_dbl(vpd, ~calc_gs_star(-1, vpd = ., viscosity = 1, dpsi_star = 1, v_huber = 1, par = par_conductance_std, par_conductivity = par_conductivity_std))) %>%
mutate(a_c_inst = purrr::map_dbl(gs_star, ~calc_assim(., vcmax = out_analytical$vcmax, ca = out_analytical$ca, par = par_photosynth_std))) %>%
filter(!is.infinite(gs_star)) %>%
mutate(vcmax_opt = purrr::map_dbl(gs_star, ~optim_vcmax_midterm_constvhuber(fn_target, gs_star = ., ca = out_analytical$ca, par_cost = list(b = 0.05), par_photosynth = par_photosynth_std, vcmax_init = out_analytical$vcmax))) %>%
mutate(a_c_accl = purrr::map2_dbl(gs_star, vcmax_opt, ~calc_assim(.x, vcmax = .y, ca = out_analytical$ca, par = par_photosynth_std)))
df_vpd %>%
tidyr::gather("source", "a_c", c(a_c_accl, a_c_inst)) %>%
ggplot() +
geom_line(aes(x = vpd, y = a_c, linetype = source)) +
labs(title = "Assimilation vs. instantaneous A", subtitle = "As a funtion of VPD")
XXX todo: