Leaf gas exchange from R

library(phylloptim)
library(dplyr)
library(tidyr)
library(ggplot2)

theme_set(theme_bw(base_size = 11) +
            theme(panel.grid.minor = element_blank(),
                  strip.background = element_rect(fill = "grey95",
                                                  colour = NA)))

What this model does differently

Most leaf gas-exchange models obtain stomatal conductance from an empirical function of assimilation, VPD and CO₂. Ball-Berry, Medlyn and their descendants all work that way. This model has no such function.

Instead it builds an explicit hydraulic path from the soil, through the roots and the root collar, up the stem to the leaf. Then it asks which leaf water potential maximises profit, where profit is assimilation minus a hydraulic cost that rises as the xylem loses conductivity. Conductance is whatever falls out of the answer.

Stomatal closure under drought is therefore not fitted here. It is derived, from two Weibull vulnerability curves, the soil potential, and the plant’s own conductance.

What is on this page

Section	What it covers
One solve to Light and VPD together	driving the model and reading the result
Traits and numerical settings are separate	which arguments are physiology and which are numerics
Water potentials are positive magnitudes	the one convention that will break a ported script
Two supply paths, two routes	the root network `plant` uses, and the collar route
Trait gradients	exact derivatives, and when they are worth using
A faster method	the closed form: what it costs in accuracy, and when it refuses

vignette("the-models") is the reference for what is being solved. vignette("fitting") works through a calibration.

One solve

leaf_solve(psi_soil = 2.0, PPFD = 900) |>
  as_tibble() |>
  glimpse()
#> Rows: 1
#> Columns: 23
#> $ psi_soil       <dbl> 2
#> $ layers         <int> 1
#> $ PPFD           <dbl> 900
#> $ atm_vpd        <dbl> 2
#> $ ca             <dbl> 40
#> $ leaf_temp      <dbl> 25
#> $ atm_kpa        <dbl> 101.3
#> $ psi_leaf       <dbl> 3.864854
#> $ collar         <dbl> NaN
#> $ ci             <dbl> 18.71375
#> $ A              <dbl> 12.8818
#> $ E              <dbl> 3.641312e-05
#> $ gc             <dbl> 0.0613037
#> $ profit         <dbl> 9.188099
#> $ hydraulic_cost <dbl> 3.693699
#> $ E_up           <dbl> NaN
#> $ uptake         <dbl> 0
#> $ lambda         <dbl> 191735.2
#> $ g1_eff         <dbl> 1.243302
#> $ Tleaf          <dbl> 25
#> $ shadow_cost    <dbl> 0
#> $ residual_flux  <dbl> 0
#> $ carbon         <dbl> 9.188099

Drivers in, operating point out. psi_soil is a positive magnitude in MPa, which a later section covers. The defaults are plant’s model, TF24_floor, which at its default price is TF24, maximised over the leaf potential with the potential upstream of the stem held at psi_soil (the leaf route). collar is NaN because that route never solves for a collar potential; the root network plant uses, and the collar route that solves for one, come later.

Most of those columns are the usual gas-exchange quantities. Two are less usual.

lambda is the marginal cost of water, dA/dE, at the operating point. It is how much carbon the leaf gives up per unit of water it declines to transpire, and it is the quantity the optimisation equalises.

g1_eff re-expresses the solved conductance as the g1 of the Medlyn model, which is a convenient common scale if that is the model you are used to. It is computed from the solve rather than supplied, so it varies with conditions.

A drought response

leaf_solve() is vectorised over its drivers, so a response curve is one call.

drought <- leaf_solve(psi_soil = seq(0.2, 5.5, length.out = 40), PPFD = 900) |>
  as_tibble()

drought |>
  select(psi_soil, A, gc) |>
  pivot_longer(c(A, gc), names_to = "variable") |>
  mutate(variable = recode(variable,
                           A = "assimilation~(mu*mol~m^-2~s^-1)",
                           gc = "stomatal~conductance~(mol~m^-2~s^-1)")) |>
  ggplot(aes(psi_soil, value)) +
  geom_line(linewidth = 0.8, colour = "steelblue4") +
  facet_wrap(~variable, scales = "free_y", labeller = label_parsed) +
  labs(x = "soil water potential (MPa)", y = NULL)

Both fall as the soil dries, because each unit of water costs more as the stem nears the end of its vulnerability curve. Neither is held up by a minimum-conductance floor: g_min is zero by default.

The same sweep on a familiar scale

g1_eff across the same drought, over the range where the stomata are open:

drought |>
  filter(gc > 0, is.finite(g1_eff)) |>
  ggplot(aes(psi_soil, g1_eff)) +
  geom_line(linewidth = 0.8, colour = "firebrick") +
  labs(x = "soil water potential (MPa)",
       y = expression(effective~g[1]~(kPa^0.5)))

It declines as the soil dries, because it is computed from the solve rather than being a parameter of it.

Light and VPD together

Because leaf_solve() recycles its drivers elementwise, a factorial design is an expand_grid() and a single call. There is no loop, and nothing to join back afterwards, because the drivers come back as columns of the answer.

expand_grid(PPFD = seq(50, 1800, length.out = 30),
            atm_vpd = c(0.5, 1, 2, 3)) |>
  with(leaf_solve(psi_soil = 1.0, PPFD = PPFD, atm_vpd = atm_vpd)) |>
  ggplot(aes(PPFD, A, colour = factor(atm_vpd))) +
  geom_line(linewidth = 0.8) +
  scale_colour_viridis_d(end = 0.9, direction = -1) +
  labs(x = expression(PPFD~(mu*mol~m^-2~s^-1)),
       y = expression(A~(mu*mol~m^-2~s^-1)),
       colour = "VPD (kPa)")

The curves saturate at different heights because the hydraulic cost of the transpiration needed to support a given assimilation rises with VPD.

Traits and numerical settings are separate

The C++ constructor takes fifteen positional arguments, four of which are tolerances sitting among the physiology. There are also four traits it does not take at all, which have to be assigned afterwards. leaf_traits() and leaf_control() split the first fifteen and supply the missing four, so a calibration loop varying traits never has to know either fact.

vulnerable <- leaf_traits(stem_P50 = 2.2, root_P50 = 2.2)

bind_rows(
  default = leaf_solve(psi_soil = 3.0, PPFD = 900),
  vulnerable = leaf_solve(psi_soil = 3.0, PPFD = 900, traits = vulnerable),
  .id = "traits"
) |>
  select(traits, psi_leaf, A, gc, g1_eff)
#>       traits psi_leaf         A          gc    g1_eff
#> 1    default 4.760361 9.0675958 0.035387735 0.7651323
#> 2 vulnerable 3.278536 0.8185905 0.002452242 0.2586561

stem_P50/stem_c and root_P50/root_c are two different vulnerability curves. The stem pair drives the hydraulic cost; the root pair drives uptake. They are named explicitly rather than left as an unmarked b/c pair, because using one where the other belongs is easy to do and produces a plausible wrong answer rather than an error.

Tolerances go the other way, through leaf_control():

tibble(
  setting = c("default", "tight"),
  psi_leaf = c(
    leaf_solve(psi_soil = 2.0, PPFD = 900)$psi_leaf,
    leaf_solve(psi_soil = 2.0, PPFD = 900,
               control = leaf_control(GSS_tol_abs = 1e-5))$psi_leaf
  )
)
#> # A tibble: 2 × 2
#>   setting psi_leaf
#>   <chr>      <dbl>
#> 1 default     3.86
#> 2 tight       3.86

Those two are identical, because the solve root-finds its own first-order condition rather than searching profit, to about 1e-12. GSS_tol_abs does not reach the leaf route at all.

It is read in two places, both on the collar route: the threshold below which the feasible collar interval is treated as a single point, and the tolerance of a fallback search used when neither bracket endpoint has a usable gradient.

This is not only a precision story. It is what makes the trait gradients below possible. An argmax resolved only to 1e-3 moves in discontinuous jumps as a trait varies, so its derivative comes back as noise, or worse, as something smooth, plausible and sign-inverted.

Water potentials are positive magnitudes

There is one representation throughout: every ψ is a positive number of MPa, including psi_soil, psi_crit and the collar potential. Where an equation needs one potential to oppose another, the minus sign is in the equation.

This is asserted rather than documented, which matters if you are porting a script that passed signed potentials:

leaf_solve(psi_soil = -2.0, PPFD = 900)
#> Error:
#> ! `psi_soil` must be positive magnitudes in MPa, not signed potentials. If this came from an older script, drop the minus sign rather than negating the result.

That is an error rather than a quiet abs(). A model accepting both conventions would let such a script return a plausible wrong number instead.

Two supply paths, two routes

Everything above used the single-potential supply, which collapses the whole soil-to-collar journey into one series resistance. The other path is the multi-layer root network that plant uses. Both are driven through the same root_network argument, because both take their resistances as a per-call driver.

On the leaf route that resistance never reaches the solve: the potential upstream of the stem is pinned at psi_soil. The collar route maximises over the root-collar potential instead, so water has to be pulled through the roots, and the resistance costs the leaf. Every cost curve takes either route; the collar route is the only one on the multi-layer supply.

leaf_solve(psi_soil = 1.5, PPFD = 900, route = "collar",
           root_network = series_resistance(1e3)) |>
  as_tibble() |>
  select(psi_soil, psi_leaf, collar, A, gc, g1_eff)
#> # A tibble: 1 × 6
#>   psi_soil psi_leaf collar     A     gc g1_eff
#>      <dbl>    <dbl>  <dbl> <dbl>  <dbl>  <dbl>
#> 1      1.5     3.49   2.52  8.22 0.0310  0.691

The resistance is per unit leaf area, because nothing in this model is allowed to scale with plant size. On the collar route, raising it makes water harder to get, and the leaf closes down:

tibble(resistance = 10^seq(2.5, 4.5, length.out = 25)) |>
  rowwise() |>
  mutate(A = leaf_solve(psi_soil = 1.5, PPFD = 900, route = "collar",
                        root_network = series_resistance(resistance))$A) |>
  ungroup() |>
  ggplot(aes(resistance, A)) +
  geom_line(linewidth = 0.8, colour = "steelblue4") +
  scale_x_log10() +
  labs(x = expression(soil-to-collar~resistance~(MPa~s~mol^-1~m^2)),
       y = expression(A~(mu*mol~m^-2~s^-1)))

The single potential is what makes comparison against other optimality models meaningful. Medlyn, Prentice least-cost and Cowan-Farquhar are all formulated against a single soil potential, so comparing them against a multi-layer root network would compare two things at once.

The supply is chosen when the leaf is built rather than afterwards. There is no settable leaf$supply_kind:

l <- leaf_model()
l$supply_kind <- "multilayer"
#> Error:
#> ! Leaf$supply_kind is read-only

Flipping the tag alone would leave the other path’s state configured and silently ignored. Both leaf_supply_*() functions therefore reconfigure the object completely, so it can never be in a state where the tag and the supply disagree.

Multiple soil layers

On the multi-layer supply, pass a vector for a layered profile, or a list of vectors for several profiles. Layers default to 1 m thick with the root carbon split evenly.

leaf_solve(psi_soil = list(c(1.0, 1.5, 2.0), c(3.0, 3.5, 4.0)), PPFD = 900,
           supply = leaf_supply_multilayer(), route = "collar") |>
  as_tibble() |>
  select(psi_soil, layers, A, gc, uptake)
#> # A tibble: 2 × 5
#>   psi_soil layers     A      gc   uptake
#>      <dbl>  <int> <dbl>   <dbl>    <dbl>
#> 1        1      3  5.00 0.0168  0.000555
#> 2        3      3  1.42 0.00431 0.000142

psi_soil in the result is the wettest layer, with layers alongside, so the column stays numeric and the frame stays rectangular.

The stateful interface

leaf_solve() is a wrapper over an R6 object you can drive yourself. This is what plant does, and it is what you want when you care about the intermediate state rather than only the operating point.

It is also the cheaper route per row: driving one object by hand, as below, measures at about half of leaf_solve()’s cost per row. Both are dominated by the same thing, which is that a call across the R boundary costs about half as much as the whole C++ solve. If you need an inner loop faster than that, what has to go is the number of R calls per row rather than which R function makes them. The measured figures are in the README’s Performance from R.

l <- leaf_model()
set_drivers(l, psi_soil = 2.0, PPFD = 900)
l$optimise()

tibble(profit = l$profit_, cost = l$hydraulic_cost_, A = l$assim_colimited_)
#> # A tibble: 1 × 3
#>   profit  cost     A
#>    <dbl> <dbl> <dbl>
#> 1   9.19  3.69  12.9

The object exposes the whole model, not just the operating point. The stem vulnerability curve is defined for any potential:

tibble(psi = seq(0, 6, length.out = 200)) |>
  mutate(k = vapply(psi, l$proportion_of_conductivity, numeric(1))) |>
  ggplot(aes(psi, k)) +
  geom_line(linewidth = 0.8, colour = "darkgreen") +
  labs(x = expression(psi~(MPa)), y = "proportion of conductivity",
       title = "stem vulnerability")

The profit function is not defined everywhere. Water flows from the soil to the leaf, so psi_leaf must be at least the upstream potential. Below it the leaf would be wetter than the soil, so transpiration(), stom_cond_CO2() and psi_leaf_to_ci() stop with an error rather than return a negative flux. profit_psi_leaf_TF24() reads that pair as zero flux, as the solve’s evaluator does, so it returns a number, but not one on the curve the solve maximises. psi_leaf equal to the upstream potential is the shut stomata’s state, with no stomatal flux. Plot the domain the solve actually searches:

psi_upstream <- 2.0

tibble(psi = seq(psi_upstream, 6, length.out = 200)) |>
  mutate(profit = vapply(psi, l$profit_psi_leaf_TF24, numeric(1),
                         psi_upstream = psi_upstream)) |>
  ggplot(aes(psi, profit)) +
  geom_line(linewidth = 0.8, colour = "darkgreen") +
  geom_vline(xintercept = l$opt_psi_leaf_, linetype = 2, colour = "grey40") +
  labs(x = expression(psi[leaf]~(MPa)), y = "profit",
       title = "the function being maximised")

The dashed line is where the solve landed, and the maximum it sits on is flat. That flatness has a practical consequence: the profit value is much better determined than its location. An error dp in the value displaces the location by sqrt(dp/k) for curvature k, so locating the peak by comparing heights costs precision that solving dprofit = 0 keeps. It is why the solve does not search this curve. vignette("the-models") Appendix C has the measurements.

Trait gradients, and where they stop working

Because the operating point is now located precisely, it can be differentiated. leaf_gradient() returns the derivative of each solved output with respect to each trait:

g <- leaf_gradient(psi_soil = 2.0, PPFD = 900,
                   pars = c("vcmax_25", "jmax_25", "stem_P50", "TF24_cost_scale"))
round(g$gradient, 5)
#>                        A       gc psi_leaf collar   profit
#> vcmax_25         0.05412  0.00007  0.00374    NaN  0.04565
#> jmax_25          0.01923  0.00008  0.00408    NaN  0.01000
#> stem_P50         3.78753  0.03326  0.57970    NaN  5.04967
#> TF24_cost_scale -0.29384 -0.00258 -0.12968    NaN -0.49249

Choosing which gradients to compute

pars selects them, and by default you get all of them. Omit it and leaf_gradient() differentiates every parameter it can:

dim(leaf_gradient(psi_soil = 2.0, PPFD = 900)$gradient)
#> [1] 18  5
rownames(leaf_gradient(psi_soil = 2.0, PPFD = 900)$gradient)
#>  [1] "vcmax_25"                      "stem_c"                       
#>  [3] "stem_P50"                      "root_c"                       
#>  [5] "root_P50"                      "TF24_beta2"                   
#>  [7] "jmax_25"                       "phi_J"                        
#>  [9] "curv_fact_elec_trans"          "curv_fact_colim"              
#> [11] "TF24_cost_scale"               "R_d_25"                       
#> [13] "JS22_gamma"                    "CMax_a"                       
#> [15] "CMax_b"                        "g_min"                        
#> [17] "leaf_specific_conductance_max" "TF24_floor_lambda_o"

That is 16 traits from leaf_traits(), plus leaf_specific_conductance_max, which is not a trait but is fitted like one, and TF24_floor_lambda_o, the default model’s price of water: 18 in all. The collar route over the single-potential path adds resistance, its series resistance, which the leaf route never reads; model = "TF24" drops the price, and model = "CF77" swaps in CF77_lambda_, that curve’s prescribed marginal value of water. These counts are computed rather than written out, for the reason the hidden chunk below records.

R_d_25, dark respiration at 25 °C, is one of the traits and behaves like any other. So dY/dvcmax_25 is a partial derivative at fixed respiration. A fit that wants respiration to follow Vcmax moves both and adds the two columns.

Selecting fewer parameters costs no less. The work is one solve plus one pass that carries the derivative in every parameter at once, so asking for one parameter costs what asking for all of them does.

request	ms per call
`pars = NULL` (all 18)	2.80
4 parameters	2.49
1 parameter	2.48

So name the parameters you are fitting for the shape of the result, not for speed.

The five outputs are not selectable. Every call returns the derivative of A, gc, psi_leaf, collar and profit with respect to each requested parameter, because they are read off the same solve and cost nothing extra. pars chooses the rows of that matrix, never the columns.

The first four are what a gas-exchange calibration observes. profit is the objective, and it is there for a demographic consumer: plant bills carbon from the leaf’s profit rather than its assimilation. It is also the one output the envelope theorem reaches, so at an interior optimum its gradient is the direct partial alone. ?leaf_gradient explains that, along with the measurement that makes it worth doing rather than merely true.

A collar potential you supply

If your model carries the collar potential rather than solving for it each step, as plant does, pass it and the gradient is taken there. This and the rest of the section are on the collar route over the multi-layer supply:

collar <- list(supply = leaf_supply_multilayer(), route = "collar")
g <- leaf_gradient(psi_soil = 2.0, PPFD = 900, pars = "vcmax_25", psi = 2.7,
                   supply = collar$supply, route = collar$route)
g$status
#> [1] "prescribed"

psi is exogenous, so the answer is the partial derivative at fixed collar, plus dpsi_dtheta if the collar you imposed itself moves with the traits. M and dY_dpsi come back too, so a caller integrating its own sensitivity of psi has the coefficients.

Hand back the collar the solver found, with the response it derived, and the two paths agree to rounding:

a <- leaf_gradient(psi_soil = 2.0, PPFD = 900, pars = "vcmax_25",
                   supply = collar$supply, route = collar$route)
b <- leaf_gradient(psi_soil = 2.0, PPFD = 900, pars = "vcmax_25",
                   psi = a$psi, dpsi_dtheta = -a$M / a$H,
                   supply = collar$supply, route = collar$route)
all.equal(a$gradient, b$gradient, tolerance = 1e-12)
#> [1] TRUE

A psi outside the feasible collar interval is refused: the leaf does not move it onto the interval’s end, where the outputs would follow the bound rather than your dpsi_dtheta. status reports "refused", g$psi is the one you gave, and everything else is NA. A feasible psi whose leaf potential lies in the stem curve’s tail, just inside the dry end, is "no-gradient": it has a value, but every row would be rounding there.

Naming something that cannot be differentiated is an error rather than a silent omission, and the message says what is available instead:

leaf_gradient(psi_soil = 2.0, PPFD = 900, pars = "resistance",
              supply = collar$supply, route = collar$route)
#> Error:
#> ! `pars` names things this cannot differentiate: resistance. Available: the traits from leaf_traits(), plus `leaf_specific_conductance_max`. `resistance` is a parameter of the single-potential path only; on the multi-layer path the soil-to-collar resistances are derived from root carbon.

These are not finite differences of the solve. The outputs are read at the profit-maximising collar potential, so a trait moves them twice over: directly, and by moving the optimum itself. Differentiating the optimality condition dprofit/dψ = 0, rather than the solved output, captures both terms, and one pass of automatic differentiation carries every parameter at once.

For TF24_cost_scale, TF24_beta2, stem_P50 and stem_c the second term is the entire answer. Change any of them and assimilation moves only because the leaf chooses a different operating point.

The interesting part is where that derivation stops being valid. It assumes the optimum is interior, and in dry soil it often is not: the leaf ends up at the very edge of the range of collar potentials it can physically occupy. There dprofit/dψ is not zero at the answer, and the formula is not merely imprecise but wrong by orders of magnitude, without failing loudly. The solve records which kind of point it found, and at a pinned one leaf_gradient() differentiates the bound that holds the leaf instead:

sweep <- lapply(seq(1.0, 5.0, by = 0.125), function(p) {
  g <- leaf_gradient(psi_soil = p, PPFD = 900, atm_vpd = 3.0, pars = "stem_P50",
                     supply = collar$supply, route = collar$route)
  tibble(psi_soil = p, dA_dstem_P50 = g$gradient["stem_P50", "A"],
         status = g$status)
}) |> bind_rows()

ggplot(sweep, aes(psi_soil, dA_dstem_P50)) +
  geom_line(colour = "grey60") +
  geom_point(aes(colour = status), size = 2) +
  labs(x = expression(psi[soil]~(MPa)),
       y = expression(dA/d*italic(stem_P50)),
       colour = NULL,
       title = "sensitivity to stem vulnerability, and where the optimum pins")

On the leaf route this sweep stays interior throughout: the potential upstream of the stem is pinned at psi_soil, so there is no collar interval to pin against.

Read that left to right. In wet soil, a less vulnerable stem is worth about two units of assimilation per MPa of stem_P50. The sensitivity rises as the soil dries, because vulnerability matters more when water is scarce, peaks near 2.9 MPa, eases off, and then falls away sharply. Past about 3.9 MPa the optimum is a bound: the leaf is stuck at the wettest collar potential it can physically reach, and a stem trait it has no freedom to exploit stops mattering at all. The tail is exactly zero.

The transition itself is a genuine kink in trait space, and a real non-differentiable set. It belongs to the constrained problem rather than to this implementation. A gradient-based fit that wanders across it is optimising through a corner, and no amount of care in the derivative changes that. What leaf_gradient() can do is tell you when you are there, which a naive finite difference cannot: status marks these points, from the kind of point the solve recorded.

A faster method: the closed form

Everything above reaches the operating point the same way: a root-find of the first-order condition over the feasible interval. That is the method = "exact" path, and it is the default. set_model() takes a third argument, and "closed" selects a different route to the same optimum.

The idea is the one vignette("the-models") builds the whole comparison on. Every model here satisfies dA/dE = lambda at its optimum, and given lambda the answer collapses to the Medlyn relation ci/ca = xi/(xi + sqrt(D)) with xi = sqrt(Q/lambda). The difficulty is that lambda depends on the psi you are solving for — but where it is a power law in psi at the wet end, that circularity inverts explicitly. One Newton step on the supply-minus-demand residual finishes the job, and a search over the bracket is not needed at all.

cf <- leaf_model(supply = leaf_supply_singlelayer())
set_drivers(cf, psi_soil = 0.5, PPFD = 1500, atm_vpd = 1,
            root_network = series_resistance(1500))

cf$set_model("TF24", "leaf", "closed")
cf$optimise()
tibble(method = cf$model_method(), psi_leaf = cf$opt_psi_leaf_,
       A = cf$assim_colimited_, fell_back = cf$last_solve_fell_back_)
#> # A tibble: 1 × 4
#>   method psi_leaf     A fell_back
#>   <chr>     <dbl> <dbl> <lgl>    
#> 1 closed     2.43  19.7 FALSE

When it applies, and when it refuses

It is refused rather than approximated wherever it does not hold, and the refusals carry their reasons:

TF24, CF77 and TF24_floor only, of the nine cost curves. The three-way split is worth knowing, because only one part of it is a gap in the code. A closed form exists exactly when h'(A), the derivative of the benefit link, does not depend on the solution. For JS22, CMax and ProfitMax it does not — so the inversion exists for them and is simply not written yet. For SOX and JW26 it does: their link is log, so h'(A) = 1/A, and lambda therefore carries the assimilation the solve is for. Given lambda -> xi -> ci -> A the dependency closes into a fixed point, which can be iterated but not inverted. Those two cannot have a closed form.

⚠️ TF24_floor needs $TF24_floor_lambda_o > 0, and its tail is worse than the other two’s. Its wet-end lambda is lambda_o — a constant, so the same n = 0 case as CF77 — and at lambda_o = 0 that limit is zero, xi = sqrt(Q/lambda) diverges, and there is nothing to expand about; that case is TF24, which has its own start, and the refusal points there. Where it does apply, the leading order discards the whole Θ̃'(ψ)/K(ψ) term, which at the DEFAULT TF24_cost_scale dominates lambda_o at the operating point. Measured over a 150-row driver grid, the worst served error is 2.9e-01 at TF24_cost_scale = 0.5 (against CF77’s own 2.4e-01) and 5.2e-01 at the default 7.5. The ci/ca guard does not catch that, because it reads the closed form’s own lambda and stays plausible while psi* is wrong — so on this curve the fast arm is a measurement tool rather than a production path.
The leaf route only, and one soil layer only. The inversion is against a supply measured from psi_soil; it says nothing about the soil-to-collar path, and the multi-layer closed form does not exist.
Not with the energy balance on. Fick’s law divides by the leaf-to-air deficit, which depends on leaf temperature, which depends on transpiration — which is what the inversion solves for. So D is not known at the moment of inversion and the form is not closed. This is a statement about the form rather than about its accuracy: measured with the refusal bypassed, the error reaches 2.3 MPa in psi_leaf even in the wettest bin, against 0.04 MPa on the prescribed-temperature path.

cf$set_model("SOX", "leaf", "closed")
cf$optimise()
#> Error:
#> ! the closed form does not exist for the SOX cost curve, and it is not a matter of unwritten algebra. SOX, JW26 and LeastCost are the three curves with the LOG benefit link, where h'(A) = 1/A -- so lambda carries the assimilation the solve is for, and given lambda -> xi -> ci -> A the dependency closes into a fixed point rather than an inversion. Use method "exact".

What it costs, and what it buys

The guard on validity tests an output — ci/ca at the answer — so it can only be applied after solving, and a rejected row falls back to the exact solve. That makes the fallback fraction phi part of the price: the realised speedup is 1/[phi + (1 - phi)/s] rather than the ceiling s.

Measured in tests/cpp/bench_solve at 0.16.2, with both arms in one process on a 96-point single-layer driver grid:

curve	realised speedup	`phi`	ceiling at `phi = 0`
`TF24`	2.2×	0.27	4.0×
`CF77`	2.0×	0.39	5.3×

Two-fold, in other words — not the order of magnitude an earlier version of this code advertised. Those figures predated the current solver, which root-finds the first-order condition where it used to search the objective, and they also priced the inversion alone rather than a solve a caller can read outputs from.

phi is reported per leaf, so you can measure it on your own drivers rather than taking the grid’s:

sweep_leaf <- leaf_model(supply = leaf_supply_singlelayer())
sweep_leaf$set_model("TF24", "leaf", "closed")
for (p in seq(0.25, 5, length.out = 20)) {
  set_drivers(sweep_leaf, psi_soil = p, PPFD = 1500, atm_vpd = 1,
              root_network = series_resistance(1500))
  sweep_leaf$optimise()
}
tibble(calls = sweep_leaf$closed_form_calls_,
       fell_back = sweep_leaf$closed_form_fallbacks_,
       phi = sweep_leaf$closed_form_fallback_fraction())
#> # A tibble: 1 × 3
#>   calls fell_back   phi
#>   <int>     <int> <dbl>
#> 1    20         6   0.3

Accuracy, and one failure that is not a tolerance

Three approximations sit inside the closed form: the power law for lambda, the single Newton step, and the Medlyn collapse itself — whose 3*Gstar group is the electron-transport-limited dA/dci, so it is not exact under colimitation. The last of those is the binding one, and the error is organised by ci/ca:

grid <- expand.grid(psi_soil = c(0.25, 0.5, 1, 2, 3, 4, 5),
                    PPFD = c(100, 500, 1500),
                    atm_vpd = c(0.5, 1, 2, 3))

solve_both <- function(psi_soil, PPFD, atm_vpd) {
  one <- function(method) {
    l <- leaf_model(supply = leaf_supply_singlelayer())
    set_drivers(l, psi_soil = psi_soil, PPFD = PPFD, atm_vpd = atm_vpd,
                root_network = series_resistance(1500))
    l$set_model("TF24", "leaf", method)
    l$optimise()
    l
  }
  ex <- one("exact")
  cl <- one("closed")
  tibble(ci_ca = ex$ci_ / ex$ca_, fell_back = cl$last_solve_fell_back_,
         err_A = abs(cl$assim_colimited_ - ex$assim_colimited_) /
           max(1, abs(ex$assim_colimited_)),
         err_psi = abs(cl$opt_psi_leaf_ - ex$opt_psi_leaf_))
}

res <- Map(solve_both, grid$psi_soil, grid$PPFD, grid$atm_vpd) |>
  bind_rows() |>
  filter(is.finite(ci_ca))

res |>
  mutate(bin = cut(ci_ca, c(0, 0.25, 0.4, 0.5, 0.6, 0.8, 1))) |>
  group_by(bin) |>
  summarise(n = n(), fell_back = sum(fell_back),
            mean_err_A = mean(err_A), worst_err_A = max(err_A),
            worst_err_psi = max(err_psi), .groups = "drop")
#> # A tibble: 6 × 6
#>   bin            n fell_back mean_err_A worst_err_A worst_err_psi
#>   <fct>      <int>     <int>      <dbl>       <dbl>         <dbl>
#> 1 (0,0.25]      11        11    0           0              0     
#> 2 (0.25,0.4]    13        13    0           0              0     
#> 3 (0.4,0.5]      6         3    0.0106      0.0245         0.164 
#> 4 (0.5,0.6]     14         0    0.0219      0.0659         1.50  
#> 5 (0.6,0.8]     25         0    0.00740     0.0364         0.265 
#> 6 (0.8,1]       15         0    0.00299     0.00679        0.0862

Read it as three regimes. Above ci/ca about 0.6 the closed form agrees with the exact solve to well under a percent in assimilation on average, and a few percent at worst. Between 0.5 and 0.6 the mean is a couple of percent. And below about 0.4 the guard has fired on nearly every row, so the tabulated error collapses to the solver’s own floor — those rows are the exact solve.

⚠️ The failure below ci/ca about 0.3 is topological, not a matter of tolerance. There the exact optimum turns back toward shutdown while the inversion runs on to psi_crit: the fixed point the closed form is looking for has ceased to exist, so the answer is wrong rather than imprecise, and no amount of refinement recovers it. That is why there is a guard and a fallback at all rather than a tighter convergence criterion — and it is why $last_solve_fell_back_ is worth reading if you are going to trust an individual row.

The caveat that would stop you using this in a demographic model

plant differentiates the operating point with respect to plant state, and that needs the argmax to vary smoothly with its inputs. Sweeping kmax — which is how height reaches the leaf — over ±30% and taking the mean absolute second difference of opt_psi_leaf_:

the sweep	exact	closed
never crosses the guard	6.1e-06	7.2e-06
crosses the guard	1.0e-05	1.0e-03

The closed form’s own arithmetic is barely rougher than the solve, at about 1.2×. What costs two orders of magnitude is the fallback boundary: the guard tests an output, so the set of drivers on which the method switches is a surface in driver space, and the argmax jumps as a smooth input carries the leaf across it. A mixture of two methods is discontinuous even where both are individually smooth.

So this is usable for a large batch of independent operating points, and it is not usable underneath a derivative. leaf_gradient() reflects that: seating a model resets the method to "exact", so a gradient always differentiates the exact solve whatever was seated before.

What the model assumes

Every model simplifies. The assumptions most worth keeping in mind when reading output from this one are:

Stomatal behaviour is optimal and instantaneous. The leaf sits at the maximum of its objective at every point, with no cost of adjustment and no lag.
The leaf is intensive. Every input is per unit leaf area, or is a dimensionless driver. Nothing scales with plant size, so whole-plant allometry belongs to the caller.
The hydraulic path is two Weibull curves in steady state. There is no capacitance and no refilling.
The default cost curve is one hypothesis among several. TF24 is the default because it is what plant runs, not because it has been shown to be the right account of what a leaf is spending. vignette("the-models") has the other six and how they differ.
The default traits are Eucalyptus saligna. They are a starting point rather than a recommendation, and the same is true of both default supply networks.

Where to go next

I want to know what is being solved. vignette("the-models") is the reference: the objective, the nine cost curves, the first-order conditions, and the two optional gates.
I want to fit this to data. vignette("fitting") works through a calibration, and covers when the exact gradient earns its complexity.
I want to call it without R. The model is a set of self-contained C++ headers, so the same solve is available from a C++ program or a Python extension. vignette("cpp-interface") covers it.
I want to know how this differs from what I already use. COMPARISON.md sets it against plantecophys, tealeaves and bigleaf.

One loose end worth flagging if you are supplying your own hydraulics. Both supply paths take resistances through the same root_network argument, so a caller with measured or fitted resistances never has to go through a root carbon model. The single exception is gravity_head, which leaf_supply_singlelayer() still takes as configuration, because that path has no depth profile to derive it from.