The two supply paths

Water reaches the leaf through soil, then roots, then the root collar — the point where the root system joins the stem — and then up the stem. The model can be set up to optimise the potential at either of two points on that path, and the two set-ups are called routes:

Those names are used throughout this page. The literature calls the first “non-root-based”, which is also what the package’s refusal message means by the term, and find_root_collar_psi() is the second under its own name.

These are the same optimality problem in two coordinates, over different supply. Neither is a check on the other. Over 30 driver rows they disagree on 20, and the disagreement is the soil-to-collar drop, which only the root-based route charges for. Remove that drop and the two converge. Sweeping the series resistance down through three decades closes the gap by about a factor of ten per decade, reaching 4.0e-04 MPa.

One point in that sweep looks like a structural gap and is not. At a resistance small enough that the whole range of flux fits inside a collar interval a few times the resistance wide, a collar known to double precision no longer resolves the flux, and once the interval closes the root-based solve reports hydraulic shutdown (#176). Real resistances are around 1e3, far from this, so read the trend rather than the endpoint.

One problem, several coordinates

Different papers in this literature optimise with respect to different things. Cowan and Farquhar (1977) solve for stomatal conductance. Sperry et al. (2017) solves for water potential. Others work in terms of transpiration, or of internal CO₂. It is easy to read that as a substantive disagreement about what the plant controls.

For the class of models here, it is not. These variables are interchangeable, and any of the models could be solved in any of them. The choice is one of convenience rather than of biology.

The reason is that, given the drivers, the quantities are not independent. They are linked by a chain of monotone maps:

\[\psi \;\xrightarrow{\;\text{supply}\;}\; E \;\xrightarrow{\;\text{diffusion}\;}\; g_c \;\xrightarrow{\;\text{supply}=\text{demand}\;}\; c_i \;\xrightarrow{\;A(\cdot)\;}\; A\]

Fix any one of these and the rest follow. Because each map is monotone, it carries the maximum to the maximum, so an optimum located in one coordinate is the optimum in all of them. “Optimise over leaf potential”, “optimise over stomatal conductance” and “optimise over c_i” therefore name one problem written three ways, and a model’s operating point does not depend on which was chosen.

This package solves in water potential, because that is the coordinate in which the hydraulic supply constraint is simplest to state. That is an implementation choice and nothing more. What actually distinguishes the models is the cost function, which is the next section.

The maps are monotone but not linear, and that limits how far the equivalence goes. A derivative does depend on the coordinate. dA/dE and dA/dψ are different numbers, related by dE/dψ. So a marginal cost only means something once you say what it is per unit of, and comparing a λ from one paper with a λ from another means checking that they share a denominator. A leaf’s $lambda is per unit E, which is the convention the λ literature uses.

Solving for the optimum by its first-order condition

There are two ways to find ψ*: search J for its maximum, or solve

\[\frac{dJ}{d\psi}=0\]

for a stationary point. Every solver in the package takes the second route, and the three subsections here are why.

The comparison below is against golden section, the standard bracketing search, because it is the natural alternative and because the difference is instructive on all three axes that can be measured. It is also the reason the rest of this page can promise exact gradients at all, so the argument is worth following even if you never touch a solver.

Precision. A bracketing search terminates on the width of an interval, so it never locates the argmax better than its own tolerance. Over the reference grid, golden section run to a 1e-3 bracket leaves |dJ/dψ| at a median of 7.8e-04 at the point it returns. The root-find leaves 5.6e-15. That gap is eleven orders of magnitude, and it is not a matter of tightening a tolerance. The two methods answer different questions. One says “the maximum is somewhere in here”. The other says “the derivative is zero at this point”.

Smoothness. This is what a gradient actually needs. An argmax determined only to a bracket width is piecewise constant at fine scales, so differentiating it with respect to a trait yields noise dressed as a derivative. Over eleven steps in a trait, golden section returns 6 distinct argmax values out of 11, with second differences of 3.9e-04. That figure is the step size itself, which is to say nothing at all. The root-find returns 11 distinct values of 11, with second differences of 3.4e-07.

The noisy answer is worse than merely noisy, because it is not obviously broken. One trait derivative taken that way comes back smooth, plausible, and sign-inverted: −2.6e-03 against a true +2.6e-04.

Speed. Golden section shrinks its bracket by a constant factor per evaluation. That is linear convergence, about 0.62 per step, and it cannot be improved on, because the method uses only comparisons. A Brent-class root-finder interpolates through the points it already has and converges superlinearly. On the collar solve, with both binaries built from one tree and run interleaved, golden section costs 3.51 µs and the root-find 2.65: a saving of 24%. The root-find is therefore faster and eleven orders more accurate, from the same model.

The three claims are one picture. Both methods in the figure below are given the same sign-changing bracket, which flatters the search: it is handed the basin rather than having to find it, and finding the basin is what the scan exists for.

sc <- plot_solver_convergence()

e_search <- abs(sc$search$x - sc$star)
e_root <- abs(sc$rootfind$x - sc$star)
reach <- function(e, tol) if (any(e < tol)) which(e < tol)[1] else NA_integer_
for (tol in c(1e-4, 1e-8, 1e-12)) {
  cat(sprintf("to reach %5.0e MPa:  root-find %2s evaluations, search %s\n", tol,
              reach(e_root, tol),
              ifelse(is.na(reach(e_search, tol)), "never", reach(e_search, tol))))
}
#> to reach 1e-04 MPa:  root-find  4 evaluations, search 7
#> to reach 1e-08 MPa:  root-find  5 evaluations, search never
#> to reach 1e-12 MPa:  root-find  5 evaluations, search never
cat(sprintf("search stalls at %.2e MPa; sqrt(4.eps.|J|/k) predicts %.2e\n",
            min(e_search), sc$search_floor))
#> search stalls at 2.08e-08 MPa; sqrt(4.eps.|J|/k) predicts 2.61e-08

The search does not merely converge more slowly. It stops. After about thirty evaluations it stops improving at all, and no further spending moves it.

That floor is not a tolerance anyone chose. It is where comparing values of J runs out of information. Near its maximum J is flat, so ψ can move while J changes by less than a rounding error, and the comparison driving the search becomes noise. The limit is sqrt(4ε|J|/k) for curvature k, drawn as the dashed line in the figure. The observed stall sits within a factor of two of it, which is close enough to identify the mechanism; the remaining factor is the constant in that expression rather than anything about the model.

A root-find has no such floor, because it is looking for a zero of a function that is steep exactly where J is flat. The quantity it tests still carries information at 1e-15, which is why panel (b) keeps descending.

This is the same flat maximum that Appendix C is about, seen from the inside. There it decides what a difference of the model can mean. Here it decides what a search can find.

Having the derivative in hand is what makes the trade one-sided. dJ/dψ is not a second model to evaluate: it comes out of the same forward pass as J. So one iteration of the root-find costs about what one iteration of the search costs, and the superlinear convergence is pure gain.

Take the derivative away and the comparison closes up again. Differencing dJ/dψ would cost two objective evaluations per iteration, and would bring a step-size problem of its own: a difference through a nested root-find has a noise floor of tolerance divided by the step, so the obvious small step is the wrong one. The efficiency argument here is therefore not really about Brent versus golden section. It is about what a cheap exact derivative buys, which is also why the trait gradients in Appendix B are worth their complexity.

The scan, which does not run by default

A root-find finds a stationary point, not the best one, so where the objective has two humps something has to choose between them. That job falls to a scan: a sweep of the interval on a coarse grid, picking the cell that contains the taller hump, after which the root-find refines inside that cell.

The sweep costs 64 objective evaluations before refinement begins, which is expensive in a package built around a solve of a few microseconds. So it is switched on only where the objective is measurably two-humped, which over a 1728-row driver sweep is only ever with the energy balance on. With that gate off, which is the default, there is no scan and the solver is endpoints plus a root-find.

When the answer is a bound, not a stationary point

Where profit is still climbing as it reaches the end of the feasible interval, dJ/dψ never crosses zero inside that interval and the optimum is the endpoint itself. This is a constrained optimum, and it is not a corner case: it happens on 51 of the 288 feasible rows of the reference grid at 25 °C.

A solver therefore has to test the derivative’s sign at each end and say explicitly that the answer is a bound. Handed no sign change, a bracketing root-find cannot on its own distinguish “the answer is the endpoint” from “I failed”, so the operating-point kind is what carries the distinction.

The feasible interval runs from the wet bound, where the stomata are closed and the stem carries only the cuticle’s residual flux (ψ_soil itself, with no flow, at the default g_min = 0), to the dry end. For SOX, JW26, CF77 and LeastCost the dry end is ψ_crit, where their own published forms end. For every other curve ψ_crit is a labelled point on a curve that continues past it, and the dry end is the end of the stem curve, where it carries all it ever can: their marginal cost of water diverges as the stem dries, and that, rather than a bound, keeps the optimum inside. Constrained answers are named for the bound that binds:

The solve records eight kinds in all; the others are rarer or are refusals. vignette("derivative-surface") lists them with what each means for a derivative.

boundary-soil is the stressed case, despite sitting at the wet end of the interval. That is why the two are named for the bound rather than for the moisture: ψ_soil is the least-tension end, so “wet” reads as the comfortable case, whereas in fact it describes what a leaf does when carbon cannot pay for water. The reference grid makes the direction plain: warming the leaf from 25 °C to 40 °C takes boundary-soil from 45 rows to 116, and boundary-dry from 11 rows to none.

boundary-dry is rarer than it looks. Under TF24 it is reached only where the soil is dry and the air humid enough that the leaf takes all its stem can still carry, and any product objective is immune to it at ψ_crit by construction, because g(ψ_crit) = 0 makes the cost diverge before the bound and water becomes infinitely expensive first.

A bound that binds is also not the same as a bound that is hit exactly. The bracket steps in by around one part in a million of its width, so testing ψ* == ψ_soil to detect closure is a test that never fires. Compare against the bound with a tolerance instead, or read the sign of the gradient there.

hydraulic-shutdown is separate from the other three and easily confused with them. It accounts for 96 of the 672 reference points, 48 at each leaf temperature, and it is a statement about feasibility rather than about a choice the leaf made.

Finally, a stationary point need not be unique. Where the marginal cost is peaked, which the next section explains, the first-order condition can have two roots and a bracketing solve finds one of them. This is the one real argument for keeping a scan: searching the objective over a grid finds the taller hump but cannot feed a gradient, while solving the condition feeds a gradient but can land on the wrong hump.

The two properties are not in conflict. What reconciles them is multi-start on the first-order condition, using several sub-brackets, each yielding an exact stationary point, compared on profit. Every candidate is then a genuine stationary point, so each one is smooth in the traits, and the comparison between them is a comparison of maxima, which is the well-conditioned quantity. A grid argmax is neither.

The benefit, and the root-find inside it

B(ψ) = A(c_i) is colimited assimilation, and c_i is not free. It is fixed by supply meeting demand, which is a root-find at every candidate ψ. There is therefore a solver inside the solver, and both its cost and its precision multiply through everything above.

The same move works one level down: differentiate the condition rather than the search. dA/dψ cannot be taken directly, because c_i is implicit. The chain runs A'(c_i) · dc_i/dψ, and dc_i/dψ comes from applying the implicit function theorem to the supply-equals-demand residual. The derivative therefore never passes through the iteration that located c_i.

That buys two things. Differencing the inner solve would multiply its iteration count into every outer step, and it would inherit a noise floor of the inner tolerance divided by the step, which is what makes an apparently sensible step return an apparently sensible wrong derivative. Differentiating the residual is exact and costs one solve.

The inner tolerance sets what every reported output of this model means, so it was chosen by measurement rather than by taste. Perturbing the model identically at each tolerance and reading the reference grid gives:

The model runs at 1e-13. 1e-10 looks like the knee of that table, and it is cheaper by about 5% per solve and per gradient, but what the table measures is the worst case at fixed drivers. A caller that differences the solve sees something else: where the root-finder stops inside its tolerance band depends on the state, so at 1e-10 the outputs are noisy along any smooth path through the drivers. Along plant height at fixed light that noise reached 3.7e-11 relative, which a finite difference at a step of 1e-7 turns into an error of 1e-4. At 1e-13 the same paths are smooth to 2.3e-14, and the trait gradient’s own finite-difference partials stop inheriting the noise: on the reference gradient file they moved by up to 2.3e-3 relative, toward a Richardson-extrapolated difference of the solve, which they now match to 1e-8.

The plateaus in that table matter. 1e-10 and 1e-11 return the same number, because they are the same count of solver iterations, so tightening between plateaus costs something and buys nothing. 1e-12 still left 1.2e-12 of noise at some states.

A general point follows, and it is the reason this subsection sits directly after the last one: a nested solver’s precision is only as visible as the loosest tolerance above it. An outer search terminating on a 1e-3 bracket would hide this inner tolerance completely, whatever it was set to. Solving the outer first-order condition is what makes it the model’s dominant amplifier. If you tighten one tolerance in a stack, expect the next one down to become the answer.

The tolerance you can set from R is not this one. The inner c_i root-find is fixed at 1e-13 on both routes. ci_abs_tol in leaf_control() is read in exactly one place, solve_medlyn_ci_numerical, which is the empirical Medlyn-conductance route. It reaches neither family of optimality solver, so tightening it buys no precision anywhere in the optimality model. The name promises otherwise, which is why it is worth stating.

There is one regime where the implicit function theorem does not apply. When the leaf is dark enough, or hot enough, that assimilation is negative across the whole [Γ*, c_a] bracket, there is no supply-equals-demand root. c_i is then put at c_a and the leaf nets A(c_a), as a leaf moving no gas does, so the open branch meets the shut one continuously as E → 0. There is no residual to differentiate there, so the derivatives take A(c_a)’s own slope instead. Only a genuine absence of a root reaches that branch; a failure of the root-finder stops rather than land there, because c_a would net the largest A on offer.

The benefit link, and the proof it costs nothing

Three of those nine rows do not subtract a cost: two are products and one is a ratio. This section shows exactly how they belong in one table with the six that do, and that nothing is lost in getting them there.

Write the objective as

\[J(\psi)=h\big(A(\psi)\big)-C(\psi)\]

with h a benefit link, meaning a transform applied to assimilation before the cost is subtracted. Three links cover every model here.

For a product objective, take the cost to be the logarithm of the reduction factor,

\[C(\psi):=-\log g(\psi) \qquad\Longrightarrow\qquad J=\log A+\log g=\log\big(A\,g\big)=\log P\]

so the additive form is not an analogy. J is the logarithm of the product. The link therefore applies to the whole objective rather than to half of it: the cost’s share of the transform is absorbed into which cost curve you name, which is the −log g row of the table above.

The first-order condition is unchanged, which can be shown by differentiating both descriptions of the same maximum:

\[\frac{dP}{d\psi}=A'g+Ag', \qquad \frac{dJ}{d\psi}=\frac{A'}{A}+\frac{g'}{g}=\frac{A'g+Ag'}{A\,g}=\frac{1}{A\,g}\frac{dP}{d\psi}\]

The two differ by the factor A·g, which is strictly positive wherever the log form is defined. They therefore vanish at exactly the same potentials. The framework’s stationary points are the product’s stationary points, and putting the link on the benefit drops nothing. It is the product rule, divided through.

That claim is checkable rather than decorative. In the chunk below, the left-hand side comes from the model’s own dJ/dψ, the single expression it uses for all nine curves. The right-hand side is A'g + Ag' assembled from parts: A' as the derivative of a zero-cost curve, and g' differentiated by hand from the closed form for g quoted above. So the check does not simply ask the implementation to agree with itself.

psi_soil <- 1.0
tr <- leaf_traits(JS22_gamma = 0)          # a zero unit cost: dC/dpsi is exactly 0,
l  <- leaf_model(tr, supply = leaf_supply_singlelayer())   # so JS22's derivative IS dA/dpsi
set_drivers(l, psi_soil = psi_soil, PPFD = 900)

f      <- function(p) 2^-((p / tr$stem_P50)^tr$stem_c)
fprime <- function(p) -f(p) * log(2) * tr$stem_c * p^(tr$stem_c - 1) / tr$stem_P50^tr$stem_c
model <- list(
  SOX  = list(P  = function(p) l$profit_psi_leaf_SOX(p, psi_soil),
              g  = function(p) (f(p) - 0.05) / 0.95,
              gp = function(p) fprime(p) / 0.95),
  JW26 = list(P  = function(p) l$profit_psi_leaf_JW26(p, psi_soil),
              g  = function(p) 1 - p / l$psi_crit,
              gp = function(p) -1 / l$psi_crit))

cat(sprintf("%-5s %6s %14s %14s %9s\n", "curve", "psi", "(logP)'*P", "A'g + Ag'", "rel"))
#> curve    psi      (logP)'*P      A'g + Ag'       rel
for (nm in names(model)) {
  m <- model[[nm]]
  for (q in c(1.5, 2.5, 3.5, 4.5)) {
    P  <- m$P(q)                                   # the UNtransformed objective, A*g
    A  <- l$assim_colimited_
    # The derivative reads whichever model is seated, so the seat moves with the
    # curve being checked. No index anywhere: a name in, a name out.
    l$set_model(nm, "leaf")
    lhs <- l$dprofit_dpsi_leaf_checked(q)[1] * P
    l$set_model("JS22", "leaf")
    Ap  <- l$dprofit_dpsi_leaf_checked(q)[1]
    rhs <- Ap * m$g(q) + A * m$gp(q)
    cat(sprintf("%-5s %6.1f %14.6f %14.6f %9.1e\n", nm, q, lhs, rhs,
                abs(lhs - rhs) / max(abs(Ap * m$g(q)), abs(A * m$gp(q)))))
  }
}
#> SOX      1.5       9.072038       9.072038   0.0e+00
#> SOX      2.5      -0.754568      -0.754568   9.6e-17
#> SOX      3.5      -4.010238      -4.010238   0.0e+00
#> SOX      4.5      -3.486972      -3.486972   5.0e-16
#> JW26     1.5       6.897052       6.897052   0.0e+00
#> JW26     2.5      -0.180299      -0.180299   1.2e-17
#> JW26     3.5      -2.214832      -2.214832   1.6e-16
#> JW26     4.5      -2.765833      -2.765833   0.0e+00

Agreement is at the last bit. The relative error is measured against the scale of the two terms rather than against their sum, and that is deliberate: A'g + Ag' cancels to exactly zero at the optimum, so an error relative to the sum would diverge there even for a perfectly correct implementation.

One consequence matters when reading profit_ across models. It now carries three kinds of number: carbon for the four difference objectives, a dimensionless normalised profit for ProfitMax, and carbon times a dimensionless factor for the two products. Comparing it between curves is not meaningful. Comparing each curve’s argmax is.

The wet-end price, and the one curve that has one

Read down the marginal-cost column and every entry but two carries a factor of f or of |f'|, both of which vanish as the drop closes. So for those curves

\[\lambda(\psi)=\frac{dC/d\psi}{dE/d\psi}\;\longrightarrow\;0 \qquad\text{as }\psi\to\psi_s\]

— water is free precisely when it is abundant. That is not a modelling choice anyone made; it is what “the cost of water is the conductivity you lose to move it” implies, and every conductance-loss cost in the literature inherits it. CF77 escapes it by not being a conductivity cost at all: its λ is a prescribed constant everywhere.

TF24_floor is the repair, and it takes the form the decomposition already suggests. Any cost splits into a part depending on the potential alone and a part linear in the flux,

\[\Theta(E)=\tilde\Theta(\psi)+\lambda_o E, \qquad \lambda=\frac{\tilde\Theta'(\psi)}{K(\psi)}+\lambda_o\]

with K = k_max f the stem conductance. The second term is what every curve above sets to zero.

\(\tilde\Theta\) here is TF24’s own cost, read at TF24’s own traits — so λ_o is this curve’s only parameter, and TF24 is the same curve at λ_o = 0. That is deliberate: it makes “is TF24 missing a price of water?” a one-restriction question rather than a comparison between two differently-parameterised fits. The two reductions are λ_o = 0 → TF24 and TF24_cost_scale = 0 → CF77, both bit-for-bit, which is how the implementation is pinned.

l <- leaf_model(supply = leaf_supply_singlelayer())
l$TF24_floor_lambda_o <- 1.5e5
set_drivers(l, psi_soil = 1.5, PPFD = 1500,
            root_network = series_resistance(1e4))

# The price at one micro-MPa of tension: essentially no flow at all. The limit is
# at zero POTENTIAL rather than at zero drop, because these two costs read the
# absolute potential.
cat(sprintf("at psi = 1e-06 MPa:  TF24 %10.3e    TF24_floor %9.1f\n",
            l$lambda_TF24(1e-6),
            l$lambda_TF24_floor(1e-6)))
#> at psi = 1e-06 MPa:  TF24  0.000e+00    TF24_floor  150000.0

λ_o has CF77_lambda_’s units and meaning — carbon per unit transpiration. Its default is 0, as in plant, where this curve is TF24 bit for bit, so TF24_floor is the default model of [leaf_solve()] and the gradients without changing a number. [leaf_solve()] and [leaf_batch()] take it as an argument; on a bare Leaf, assign the field.

Pricing the drop is not pricing the potential

JS22 is the only cost here that reads Δψ. The cleanest way to see what that changes is to ask which γ_J would reproduce TF24’s cost at TF24’s own optimum. If the two curves were reparameterisations of one another, the answer would be a constant.

out <- do.call(rbind, lapply(c(0.5, 1, 1.5, 2, 3), function(ps) {
  l <- leaf_model()
  set_drivers(l, psi_soil = ps, PPFD = 1500)
  l$set_model("TF24", "leaf"); l$optimise()
  dpsi <- l$opt_psi_leaf_ - ps
  data.frame(psi_soil = ps, psi_leaf = l$opt_psi_leaf_, dpsi = dpsi,
             cost_TF24 = l$hydraulic_cost_,
             # JS22 prices the DROP: C = gamma_J * dpsi^2.
             gamma_J = l$hydraulic_cost_ / dpsi^2,
             # CMax prices the ABSOLUTE potential, so with CMax_b = 0 the
             # matching slope solves C = CMax_a*(psi^2 - psi_soil^2)/2.
             CMax_a = 2 * l$hydraulic_cost_ / (l$opt_psi_leaf_^2 - ps^2))
}))
print(out, row.names = FALSE, digits = 4)
#>  psi_soil psi_leaf  dpsi cost_TF24 gamma_J CMax_a
#>       0.5    2.917 2.417     1.679  0.2872 0.4064
#>       1.0    3.215 2.215     2.260  0.4604 0.4840
#>       1.5    3.554 2.054     2.990  0.7091 0.5763
#>       2.0    3.931 1.931     3.844  1.0306 0.6711
#>       3.0    4.782 1.782     5.596  1.7614 0.8068
spread <- function(x) max(x) / min(x)
cat(sprintf("\ngamma_J (prices the drop)      spans %.3f to %.3f -- %.1fx",
            min(out$gamma_J), max(out$gamma_J), spread(out$gamma_J)))
#> 
#> gamma_J (prices the drop)      spans 0.287 to 1.761 -- 6.1x
cat(sprintf("\nCMax_a  (prices the potential) spans %.3f to %.3f -- %.1fx\n",
            min(out$CMax_a), max(out$CMax_a), spread(out$CMax_a)))
#> 
#> CMax_a  (prices the potential) spans 0.406 to 0.807 -- 2.0x

Neither is a constant, so neither curve is TF24 reparameterised. But the two spreads are not the same size, and the gap between them is the finding.

As the soil dries, TF24’s cost rises while Δψ shrinks. The leaf cannot pull as far, so the drop narrows even as the water gets dearer, and a cost quadratic in the drop has to be scaled up several-fold to keep pace. CMax is under no such pressure, because it reads the same absolute potential TF24 does, and its required spread is correspondingly smaller: roughly a third.

The sharpest division in the cost table is therefore not any of the individual formulae. It is whether a cost prices the potential the leaf has reached, or the distance it had to pull. That is the one difference no parameter can absorb.

It is also why JS22_gamma’s default carries a warning rather than a citation. The default matches TF24 at one reference soil potential (see Comparable at their defaults), and no single value would be matched across this band.

λ is not water use efficiency

The two are easy to conflate. They share units, and both read as “carbon per water”. They are different quantities, and the difference is between a marginal rate and an average.

\[\text{WUE}=\frac{A}{E}\quad\text{(an average)}, \qquad \lambda=\frac{dA}{dE}\quad\text{(a marginal rate)}\]

WUE is the slope of the chord from the origin to the operating point. λ is the slope of the tangent there. On a saturating A(E) those differ, and the gap is not a fixed factor:

cat(" psi_soil   lambda      A/E (WUE)   lambda/WUE\n")
#>  psi_soil   lambda      A/E (WUE)   lambda/WUE
for (ps in c(0.5, 1, 2, 3)) {
  lw <- leaf_model()
  set_drivers(lw, psi_soil = ps, PPFD = 1500)
  lw$set_model("TF24", "leaf"); lw$optimise()
  cat(sprintf(" %7.1f  %10.4g  %10.4g  %9.3f\n", ps, lw$lambda,
              lw$assim_colimited_ / lw$transpiration_,
              lw$lambda / (lw$assim_colimited_ / lw$transpiration_)))
}
#>      0.5    9.19e+04   2.606e+05      0.353
#>      1.0   1.195e+05   2.883e+05      0.414
#>      2.0   1.999e+05   3.549e+05      0.563
#>      3.0    3.15e+05   4.322e+05      0.729

So λ runs between a third and three-quarters of WUE here, and the ratio rises as the soil dries. Neither can be inferred from the other, in either direction.

What is true is narrower, and holds only at the optimum. The first-order condition is dA/dE = λ, so at the operating point λ is the marginal water use efficiency. Away from the optimum it is not even that. And “WUE” unqualified almost always means either the average A/E or the intrinsic A/g_s, so a λ from this model is not comparable with a WUE from the literature without checking which was meant.

λ is a common currency only across difference objectives

The claim above needs one qualification, and it is structural rather than incidental. Every model in the cost table except the last two maximises a difference, A(ψ) − C(ψ). Its first-order condition is A' = C', so

\[\lambda=\frac{dA}{dE}=\frac{C'(\psi)}{E'(\psi)}\]

in which A does not appear at all. That is what makes λ comparable across those models. Each prices water in absolute carbon, at a price set by the cost curve and the supply curve alone.

A product objective, A(ψ)·g(ψ), does not behave that way. Its first-order condition is A'g + Ag' = 0, so A' = -A(g'/g) and

\[\lambda=-\,\frac{A\,g'(\psi)}{g(\psi)\,E'(\psi)}\]

which is proportional to A. A model of that kind prices water relative to current assimilation rather than in absolute carbon. It is one line of algebra, and it holds for any g.

SOX is of this kind, with g = (f − 0.05)/(1 − 0.05), and so is JW26 (Jones et al. 2026), with g = 1 − ψ/ψ_crit giving λ = A / [(ψ_crit − ψ)·k_max·f(ψ)]. The two differ only in g, and everything in this section follows from the product form and holds for both.

Two consequences are worth carrying away.

The first is that a product objective’s optimum is invariant to a uniform rescaling of assimilation, where a difference objective’s is not. Replacing A by κA leaves A'g + Ag' = 0 unchanged, while κA' = C' moves. That is the same statement as λ ∝ A.

That is a property of the objective rather than a test you can run here, which is worth stating because the reverse is the obvious assumption. No trait uniformly rescales A(ψ). R_d is subtracted rather than scaled, and more fundamentally c_i sits at a root-find at fixed g_c, so scaling the capacities moves c_i and A with it. What is observable is a reduced sensitivity, and the reduction is large.

The second is that a product objective decouples stomatal behaviour from assimilation for free. As A falls, under warming for instance, λ falls in proportion, water gets cheaper, and the stomata do not close. Under a difference objective C' is unchanged while A' falls, so A' = C' is met at lower tension and the stomata do close. This is a first-order argument rather than an exact one, since it holds to the extent that a thermal decline acts as a scale reduction in A.

SOX is the product objective in the package, so the sensitivity claim can be measured rather than asserted. Scaling vcmax_25 and jmax_25 together is the closest the trait interface gets to a κ on A:

shift <- function(method, kappa) {
  l <- leaf_model(leaf_traits(vcmax_25 = 96 * kappa, jmax_25 = 157.44 * kappa),
                  supply = leaf_supply_singlelayer())
  set_drivers(l, psi_soil = 1.5, PPFD = 1500)
  l$set_model(if (method == "TF") "TF24" else method, "leaf")
  l$optimise()
  l$opt_psi_leaf_
}
for (m in c("SOX", "TF", "CMax")) {
  lo <- shift(m, 1)
  hi <- shift(m, 2)
  cat(sprintf("%-5s psi* %.4f -> %.4f  (%+.0f%%)\n", m, lo, hi,
              100 * (hi - lo) / lo))
}
#> SOX   psi* 2.7424 -> 2.8782  (+5%)
#> TF    psi* 3.5536 -> 4.3617  (+23%)
#> CMax  psi* 3.5410 -> 4.1144  (+16%)

The product objective moves by a few per cent where the two difference objectives move an order of magnitude more. That is the algebra showing through an imperfect experiment. The movement is not zero, because doubling the capacities is not a uniform rescaling of A, for the two reasons given above.

Comparable at their defaults

A comparison across curves at each curve’s own published defaults mixes two things: how each curve responds to the environment, and where its parameters happen to sit. So the defaults here are aligned. Each curve with a level parameter takes the value at which its leaf optimum has TF24’s conductance at one reference environment: a single soil potential of 1 MPa, PPFD 900, a 2 kPa deficit, 25 °C with the energy balance off, and every trait at its default. TF24 is the anchor and does not move. Because the supply is shared, matching the conductance there also matches ψ_leaf, E, A and ci.

[align_curves()] does the alignment, and the defaults are its answer at its own defaults. Four curves have a level parameter: CF77 (its price), JS22 (JS22_gamma), CMax (CMax_a) and LeastCost (LeastCost_beta). CF77 aligned to TF24 is an identity rather than a fit, since at TF24’s emergent λ the two first-order conditions coincide.

align_curves()
#>    CF77_lambda     JS22_gamma         CMax_a LeastCost_beta 
#>   1.077509e+05   4.717888e-01   6.385597e-01   5.159667e+01 
#> attr(,"target")
#> [1] 0.09042876
sp <- leaf_supply_singlelayer()
curves <- c("TF24", "CF77", "JS22", "CMax", "LeastCost", "SOX", "JW26",
            "ProfitMax")
gc <- sapply(curves, function(m) {
  leaf_solve(psi_soil = c(0.5, 1, 2, 3), model = m, route = "leaf",
             supply = sp)$gc
})
rownames(gc) <- paste("psi_soil", c(0.5, 1, 2, 3))
print(signif(gc, 3))
#>                TF24   CF77   JS22   CMax LeastCost    SOX   JW26 ProfitMax
#> psi_soil 0.5 0.1060 0.0904 0.1000 0.1050    0.0904 0.0776 0.0764    0.0770
#> psi_soil 1   0.0904 0.0904 0.0904 0.0904    0.0904 0.0650 0.0678    0.0651
#> psi_soil 2   0.0613 0.0804 0.0643 0.0600    0.0804 0.0415 0.0470    0.0421
#> psi_soil 3   0.0354 0.0416 0.0358 0.0307    0.0416 0.0218 0.0257    0.0225

The first five columns agree at psi_soil 1 and part either side of it, which is the comparison an alignment is for. Four curves are left unaligned:

• SOX, JW26 and ProfitMax have no level parameter; the supply puts their optimum. Matching them through kmax would work, but they would then no longer share the one supply every other curve uses.
• TF24_floor adds a price to TF24’s own cost at TF24’s own traits, so the only λ_o that reaches TF24’s conductance is zero, where the curve is TF24. That zero is its default.

Call [align_curves()] with a different reference or target to align for a different study, and pass its answer back as traits (JS22_gamma, CMax_a) or as [leaf_solve()] arguments (CF77_lambda, LeastCost_beta).

1. Photosynthetic demand

Farquhar-type assimilation (Farquhar et al. 1980). A Rubisco-limited rate and an electron-transport-limited rate, each a saturating function of c_i:

\[A_c=\frac{V_{c\max}\,(c_i-\Gamma^{*})}{c_i+K_m}, \qquad A_j=\frac{J}{4}\cdot\frac{c_i-\Gamma^{*}}{c_i+2\Gamma^{*}}\]

Electron transport itself is a smooth minimum of light supply and capacity:

\[J=\frac{aQ+J_{\max}-\sqrt{(aQ+J_{\max})^{2}-4\,\theta_J\,aQ\,J_{\max}}}{2\theta_J}\]

and the two rates are colimited by the same smooth-minimum form before respiration is subtracted:

\[A=\frac{A_c+A_j-\sqrt{(A_c+A_j)^{2}-4\,\theta_c A_c A_j}}{2\theta_c}-R_d\]

θ_J and θ_c are curvature factors. At 1 the expression is an exact minimum; below 1 it rounds the corner. The Michaelis-Menten composite carries the competitive inhibition by oxygen:

\[K_m=K_c\Big(1+\frac{O_2}{K_o}\Big)\]

A here is net assimilation. R_d is inside it, which is why A can be negative. A leaf with its stomata closed exchanges CO₂ only through its cuticle, at ρ g_min; at the default g_min = 0 it nets zero in light, where the transpiring leaf tends to A = 0 as E → 0, and −R_d in the dark.

2. Temperature responses

V_cmax, J_max, Γ*, K_c, K_o and R_d are all defined at 25 °C and moved from there. The first two take a peaked Arrhenius form, rising and then falling as denaturation takes over:

\[X(T)=X_{25}\exp\!\Big[\frac{E_a\,(T-T_{25})}{T_{25}\,R\,T}\Big] \cdot\frac{1+\exp\!\big[(\Delta S\,T_{25}-H_d)/(R\,T_{25})\big]} {1+\exp\!\big[(\Delta S\,T-H_d)/(R\,T)\big]}\]

with T in kelvin. Γ*, K_c and K_o take the first factor alone. Respiration is different again, following a Q₁₀ whose own value declines with temperature:

\[R_d=R_{d,25}\;q^{(T-25)/10},\qquad q=q_0-q_1\frac{T+25}{2}\]

A grid at 25 °C is blind to every one of these by construction, because each response is defined to be 1 at its reference. Any test of a temperature response needs a second temperature.

The thermal-cost gate reaches into this block. When it is on, J_max is multiplied by (1 − TC), which is why it moves the benefit and not only the cost.

3. Diffusive supply

CO₂ in and water out share the stomata, so one conductance sets both. Converting the mass flux E to a molar CO₂ conductance:

\[g_c=\frac{P\,E\,M_{\rm H_2O}}{1.6\;D_{\rm leaf}}\]

where 1.6 is the ratio of diffusivities of water vapour and CO₂, P is atmospheric pressure, and D_leaf is the leaf-to-air vapour pressure deficit. Fixing E fixes g_c, which is the first link of the coordinate chain above.

$stom_cond_CO2_ is a total leaf-to-air conductance rather than a stomatal one: the relation above is the whole path from leaf interior to bulk air, paired with a leaf-to-bulk-air deficit, which is the right quantity for the c_i solve below. With the energy balance on, the leaf has a boundary layer ($gb_H2O, its conductance to vapour), and that part of the path diffuses CO₂ at a different ratio from the stoma’s, r_b = r_s^{2/3} for a laminar layer, as in Sabot et al. (2022)’s equation 10. So the CO₂ conductance is the series of the two, with the stomatal water conductance taken from the total g_w = $gs_H2O:

\[\frac{1}{g_c}=\frac{r_s}{g_w}-\frac{r_s-r_b}{g_{bw}}\]

The transpiration itself is unchanged, because the hydraulics set it; the boundary layer changes only how much CO₂ that water buys: a median +0.3% in A at 5 m/s of wind and +0.8% at 0.5 m/s, at most +1.4%, measured on the default traits (#122). A measured g_s has the boundary layer removed from it, so compare it with $gs_stomatal_H2O, which takes $gb_H2O out of $gs_H2O in series, not with $gs_H2O. With the energy balance off there is no boundary layer, so $gb_H2O is Inf, the two are identical, and the relation above is exact.

4. Closing the loop: supply equals demand

c_i is not free. Diffusive supply must equal biochemical demand, which is a scalar root-find at every candidate ψ:

\[g(c_i)\;=\;A(c_i)\;-\;\frac{g_c\,(c_a-c_i)}{P}\;=\;0\]

solved over [Γ*, c_a] by a bracketing method. This is the innermost loop of the whole model. It runs at every objective evaluation, and its tolerance sets what every reported output means.

# The residual really is zero at the solved ci -- the statement above, checked.
l2 <- leaf_model()
set_drivers(l2, psi_soil = 1.0, PPFD = 1500)
l2$set_model("TF24", "leaf"); l2$optimise()
ci <- l2$ci_
gc <- l2$stom_cond_CO2_
resid <- l2$assim_colimited(ci) * 1e-6 - gc * (l2$ca_ - ci) / (l2$atm_kpa_ * 1000)
cat(sprintf("ci = %.4f Pa, supply-demand residual = %.2e mol m-2 s-1\n",
            ci, resid))
#> ci = 22.6517 Pa, supply-demand residual = -2.71e-20 mol m-2 s-1

5. Hydraulic supply, and why the problem is non-trivial

Water arrives through a path whose conductivity falls as tension rises — the Weibull vulnerability curve of Appendix A. Transpiration at a given leaf potential is that curve integrated along the path:

\[E(\psi_{\rm stem};\psi_{\rm up})=k_{\max}\int_{\psi_{\rm up}}^{\psi_{\rm stem}}\! f(u)\,du\]

This integral is what makes the optimisation non-trivial. Because f falls, E saturates: pulling harder recruits progressively less conductance, so the carbon return on further tension diminishes. Without that saturation, more tension would always buy more carbon and there would be no interior optimum to find.

6. Optionally, energy balance and thermal damage

A transpiring leaf is cooler than the air. Linearised, the leaf temperature is

\[T_{\rm leaf}=T_{\rm air}+\frac{R_{n,\rm iso}-\lambda_v E}{\rho c_p/r_a+g_R},\qquad g_R=4\varepsilon_l\sigma T_{\rm air}^3\]

with no clamp: set_drivers() refuses drivers under which a leaf the solve can reach would be past what the model can evaluate. \(R_{n,\rm iso}\) is the isothermal net radiation, absorbed shortwave less \((1-\varepsilon_a)\sigma T_{\rm air}^4\) with a Brutsaert sky emissivity \(\varepsilon_a\) that rises with air temperature and humidity and is held below 1 by a soft minimum of width 0.002, so the longwave term goes to zero smoothly where the fit passes unity. The leaf’s own emission, which goes as \(T_{\rm leaf}^4\), is linearised about \(T_{\rm air}\) into the radiative conductance \(g_R\), which is what keeps the balance explicit: \(T_{\rm leaf}\) is affine in E, so dT/dE = −λ_v/(\rho c_p/r_a + g_R) everywhere. Warming the leaf also raises its own vapour pressure, so the deficit driving diffusion is not the atmospheric one:

\[D_{\rm leaf}=\max\!\big(D_{\rm air}+e_s(T_{\rm leaf})-e_s(T_{\rm air}),\;D_{\min}\big)\]

That closes a feedback loop: E sets T_leaf, which sets D_leaf, which sets g_c, which sets E. It is why the energy-balance path is solved rather than evaluated, and why every derivative on it carries a dT/dE term.

Finally, the thermal cost is a sigmoid in leaf temperature, between a critical and a half-lethal value (Sicangco et al. 2026):

\[TC(T)=\frac{1}{1+\exp\!\big[-r\,(T-T_{50})\big]},\qquad r=\frac{2}{T_{50}-T_{\rm crit}}\]

With this on, a leaf has a reason to transpire that is not photosynthetic at all: evaporative cooling becomes a benefit in its own right.

What the models disagree about

Blocks 1 to 5 are shared. Every published model in the table above accepts them, and differs only in the cost. The options are a carbon penalty rising with the fraction of conductivity lost (Towers et al. 2024), the water itself priced at a constant marginal rate (Cowan and Farquhar 1977), the normalised loss of remaining conductance (Sperry et al. 2017), the square of the water-potential drop (Joshi et al. 2022), a rate rising linearly with the potential (Wolf et al. 2016), or the log of the conductance remaining, following either the vulnerability curve (Eller et al. 2020) or a straight line (Jones et al. 2026). These are not refinements of one another. They disagree about what is being spent.

Two divisions cut across that list, and both are sharper than any of the individual formulae.

The first is what the marginal cost reads, which is the third column of the table and the quantity that actually locates the optimum. For seven of the nine it is a function of ψ alone. JS22 is the exception that matters, because 2γ_J Δψ reads the drop; LeastCost’s reads the flux, and so the drop, only through its denominator. That is the one difference no parameter can absorb, since as the soil dries the drop shrinks while every absolute-potential cost rises. The gamma-match chunk measures both spreads side by side.

The second is whether the cost is subtracted or the objective is a product. Five subtract; SOX and JW26 multiply. That decides whether λ carries a factor of A, and so whether water is priced in absolute carbon or as a share of what the leaf is currently earning.

Why not (P50, b)

b is a scale parameter. It is the 63% loss point, in the same units and the same role as P50. Specifying both means specifying two points on one curve and recovering the shape from their ratio:

\[c=\frac{\ln(\ln 2)}{\ln(P_{50}/b)}\]

This diverges as P50 → b and is undefined at P50 = b, so small errors in either input make the shape explode. The pair is also badly correlated for estimation, which is why fitplc (Duursma and Choat 2017) fits Px and the slope at Px rather than b and c. c is the only genuinely independent second axis.

Why c is not a constant

c carries real biology. At c ≈ 1 the curve is exponential, with conductivity bleeding away from ψ = 0. At c > 2 it is sigmoid, with a safe plateau followed by rapid failure. Conifers are much steeper than angiosperms, and leaves and roots are flatter than stems in the same plant.

Fitted directly in (P50, c) from raw PLC curves for eight species:

The four Choat 2003 curves are not trustworthy for c. They are 13 to 14 points each, and already sit at 25 to 30% PLC by 0.8 MPa, which is a native-embolism or open-vessel artefact rather than a trait.

The defensible comparison is the three Blackman angiosperms, measured the same way. Across those, c runs 1.8 to 2.6 while |P50| spans 2.9 to 4.3 MPa: shape roughly conserved, scale free. That is the evidence for treating c as a parameter with a stable default rather than as a constant.

The default is c = 2.0 for angiosperm stems. Conifers need a substantially higher value; the 13.0 fitted above is inflated by the cavitron technique, and 4 to 5 is a safer literature-style default. Guard c > 0, and treat anything below about 1 with suspicion, since it means no safe plateau at all and in practice usually flags a measurement artefact.

Converting from published numbers

Published curves arrive in almost any form except (P50, c). Everything inverts in closed form. For a remaining conductivity fraction f, so that PLC = 1 − f:

\[P_f = P_{50}\big(\log_2(1/f)\big)^{1/c} \qquad\Longleftrightarrow\qquad P_{50} = P_f\big(\log_2(1/f)\big)^{-1/c}\]

From any two quantiles, both parameters follow:

\[c=\frac{\ln(L_2/L_1)}{\ln(P_{f_2}/P_{f_1})},\qquad L_i=\log_2(1/f_i)\]

And from P50 plus the slope at P50 in % MPa⁻¹, which is the other commonly published form:

\[S_{50}=\frac{50\,c\ln 2}{|P_{50}|}=\frac{34.657\,c}{|P_{50}|}\]

S50 should not be used as the internal parameter, because it has P50 built into it by construction.

L <- function(f) log2(1 / f)                    # f = fraction REMAINING
p_at <- function(P50, c, f) P50 * L(f)^(1 / c)  # any quantile from (P50, c)
c_from_two <- function(P1, f1, P2, f2) log(L(f2) / L(f1)) / log(P2 / P1)
s50_from <- function(P50, c) 50 * c * log(2) / P50

# round-trip: the package defaults, recovered from a (P50, P88) pair
tr <- leaf_traits()
l  <- leaf_model()
P88 <- p_at(tr$stem_P50, tr$stem_c, 0.12)
c_hat <- c_from_two(tr$stem_P50, 0.50, P88, 0.12)
P50_hat <- P88 * L(0.12)^(-1 / c_hat)
cat(sprintf("defaults: P50 = %.6f MPa, c = %.6f, so P88 = %.6f, P95 = %.6f\n",
            tr$stem_P50, tr$stem_c, P88, p_at(tr$stem_P50, tr$stem_c, 0.05)))
#> defaults: P50 = 3.400000 MPa, c = 2.680147, so P88 = 5.160000, P95 = 5.870283
cat(sprintf("recovered from (P50, P88): c = %.6f, P50 = %.6f  (errors %.1e, %.1e)\n",
            c_hat, P50_hat, abs(c_hat - tr$stem_c), abs(P50_hat - tr$stem_P50)))
#> recovered from (P50, P88): c = 2.680147, P50 = 3.400000  (errors 4.4e-16, 4.4e-16)
cat(sprintf("S50 at the defaults: %.1f %%/MPa\n", s50_from(tr$stem_P50, tr$stem_c)))
#> S50 at the defaults: 27.3 %/MPa

# the derived pair, at FULL precision -- b and psi_crit are read off the model,
# the right-hand sides are this appendix's own formulae
cat(sprintf("stem_b   = %.17g  vs  P50/(ln2)^(1/c) = %.17g\n",
            l$stem_b, tr$stem_P50 / log(2)^(1 / tr$stem_c)))
#> stem_b   = 3.8982451221145307  vs  P50/(ln2)^(1/c) = 3.8982451221145307
cat(sprintf("psi_crit = %.17g  vs  P95            = %.17g\n",
            l$psi_crit, p_at(tr$stem_P50, tr$stem_c, 0.05)))
#> psi_crit = 5.8702827267723245  vs  P95            = 5.8702827267723245
cat(sprintf("f(psi_crit) = %.17g   (the critical fraction, exactly)\n",
            l$proportion_of_conductivity(l$psi_crit)))
#> f(psi_crit) = 0.05000000000000001   (the critical fraction, exactly)

Those are equalities rather than agreements. stem_b and psi_crit are computed from (P50, c) by the expressions printed beside them, so there is no gap to quote.

That matters more than it looks, because rounding a derived quantity is amplified by the solve. If you set stem_b and psi_crit by hand to a rounded P95 instead of deriving them, the error does not stay where you put it. Rounding psi_crit to seven figures displaces it by 4.7e-08 and stem_b by 3.1e-08, and that propagates to a worst case of 4.4e-05 on assimilation across 576 of the reference grid’s points. A thousandfold amplification of a rounding-level input, which is still below the 1e-4 at which this project calls a difference real, but well above where it would be if the amplification were 1.

So specify (P50, c) at whatever precision you have, and let everything else be derived.

psi_crit is a quantile, not a parameter

psi_crit is P_f at f = 0.05, which is to say P95. Sperry’s reference conductivity k_crit = f·kmax uses the same f. So psi_crit, P95 and k_crit are one number wearing three names. Derived from a single named constant they cannot disagree; held as a free trait plus a hard-coded 0.05, they needed a consistency check to police them. f = 0.05 is a convention rather than a knob, and it is 0.05 everywhere it appears across this family.

Two things follow, and both are easy to miss.

• The domain check disappears, because it becomes vacuous. A derived psi_crit describes the curve it was derived from, so there is no consistency left to validate.
• For the curves that can pin at psi_crit (CF77 and LeastCost; SOX and JW26 end there too, but shut a leaf with A <= 0 rather than pin it), a leaf pinned there has a hydraulic cost that does not depend on the curve’s scale: the remaining conductivity is exactly f whatever P50 is, so dprofit/dP50 at that pin is exactly zero rather than merely small.

Three consequences for gradients

All three are asserted in the suite rather than described here, and the numbers below are what those assertions compare against.

1. P50 has a homogeneity shortcut. At fixed c, b is a constant multiple of P50, so the cumulative-vulnerability integral is homogeneous of degree 1 and d/dP50 = (ln 2)^{-1/c} · d/db, a factor of 1.1465 at the defaults. It is also why the stem’s vulnerability table is built in ψ/b and depends on c alone: moving stem_P50 rebuilds nothing. psi_crit scales by the same factor.

2. The c gradient splits. Once P50 is the scale parameter, b depends on c as well:

\[\frac{d}{dc}\bigg|_{P_{50}}=\frac{d}{dc}\bigg|_{b}+\frac{db}{dc}\frac{d}{db}, \qquad \frac{db}{dc}=\frac{b\ln(\ln 2)}{c^{2}}\]

The second term is nearly free, because it reuses the cheap b-derivative. But omitting it silently converts a total derivative into a partial, which produces a plausible wrong gradient rather than an error.

3. A derived dry end adds a term. Moving P50 moves the dry end, so d/dP50 picks up a term through it: dpsi_crit/dP50 = (log₂(1/f))^{1/c} (1.7266 at the defaults) for the curves ψ_crit bounds, and dpsi_end/dP50 = (40/ln 2)^{1/c} (4.5409) for the rest, both constant at fixed c. The term is exactly zero at an interior optimum and of order 1 at one sitting on a bound, so the clean identity in (1) holds only where the optimum is interior.

The size of that term is the reason to care. At a ψ_crit pin the two halves were 1.1465 × 4.6991 from the curve and 1.7266 × 1.2605 through the bound, so dropping the second gave 5.39 where the answer was 7.56. That is a 29% error with nothing about it that looks wrong.

The root curve obeys the identical homogeneity identity and takes the same (root_P50, root_c) treatment. Two curves on different conventions is the shape of the mistake that once carried a root parameter into a stem cost.

The implicit function theorem

The reported outputs are evaluated at the optimum ψ, so a parameter moves them in two ways: directly, and by moving ψ itself. At an interior optimum dJ/dψ = 0, and differentiating that condition gives ψ*’s response without re-solving:

\[\frac{d\psi^*}{d\theta}=-\frac{M}{H},\qquad M=\frac{\partial^2 J}{\partial\psi\,\partial\theta},\qquad H=\frac{\partial^2 J}{\partial\psi^2}\]

\[\frac{dY}{d\theta}=\left.\frac{\partial Y}{\partial\theta}\right|_{\psi} +\frac{\partial Y}{\partial\psi}\,\frac{d\psi^*}{d\theta}\]

The second term is not a correction. For the cost parameters it is the whole answer.

Every term is exact. The outputs, dJ/dψ and the bounds are closed forms in the parameters, and one pass of forward-mode automatic differentiation carries their derivatives in every parameter at once, so a gradient costs the same whatever it is asked for. Only H is a difference, of the exact dJ/dψ at ψ*. vignette("derivative-surface") has the machinery.

The link makes one derivative serve every model

M and H are second derivatives of J = h(A) − C, so both need dJ/dψ. That is where the benefit link earns its place:

\[\frac{dJ}{d\psi}=h'(A)\,\frac{dA}{d\psi}-\frac{dC}{d\psi}\]

That is one expression, with h' the only per-model term: 1, 1/A, or 1/|A|max. A product objective is therefore not a special case needing its own derivative, and neither is ProfitMax.

The envelope theorem, and where it stops

At a stationary point dJ/dψ = 0 identically, so the objective’s own indirect term vanishes and dJ/dθ is the direct partial alone. The gradient reads J at the held ψ* there, so that column is the direct partial exactly.

The envelope theorem does not survive a boundary. Where the optimum is a bound, ψ* is a parameter-dependent bound, dJ/dψ is not zero there, and −M/H is not the bound’s derivative. So the premise is not tested numerically: the solve records which kind of point it found, and at a bound the gradient differentiates the bound itself, the wet end where uptake meets the residual flux or the dry end of the stem.

Every route root-finds dJ/dψ = 0, so a point classified interior is a root. Measured over 427 interior rows across six curves, the residual |dJ/dψ| at the returned ψ* has a median of 1.2e-15, against 4.4e-05 when those routes refined on bracket width instead.

ProfitMax is the one curve where the two routes answer different questions, and the next section is about that.

Prescribed versus optimised

Pass psi and the outputs are evaluated at that potential instead of at the argmax. This is for a caller whose model tracks the optimum rather than finding it. There ψ does not move with θ, so there is no indirect term for any output and no envelope identity to invoke: the answer is the direct partial.

The exceptions are at the dry end. A psi outside the feasible interval is refused rather than moved onto the end, whose derivative would be the bound’s; one inside it whose leaf potential lies in the stem curve’s tail has a value and no gradient, because every row there would be rounding. The returned status says which.

A partial is not always a total

ProfitMax normalises by |A|max, the largest assimilation available on the supply stream. That is itself an optimum, found by endpoints plus a root-find on dA/dψ = 0, and it sits at the dry bound on 1318 of 1320 driver rows.

Because it is an optimum, there are two things dY/dθ could mean for this curve, and they are different numbers. Holding |A|max fixed across a perturbation gives a partial derivative at a fixed normaliser. Letting it move with the parameters gives the total. At psi_soil 1.5 and PPFD 1500, for vcmax_25, the partial is 0.0766 and the total 0.0579: the difference is 32% of the answer, so this is not a rounding distinction.

The gradient returns the total, which is what a fit needs. The model re-solves |A|max at every parameter vector, so an optimiser fed the partial would stop where the partial vanishes rather than where the objective does. The gradient carries |A|max as a function of the parameters, through its own first-order condition, giving 0.057885, which a difference of the whole solve reproduces to seven figures at relative steps from 1e-06 to 1e-03. No other curve carries this distinction, because no other curve normalises by a quantity that is itself solved for.

The energy-balance term

With the energy balance on, ψ reaches profit by two routes beyond the two the gate-off derivative already carries, both through the leaf temperature: directly, ψ → ψ_leaf → E → T_leaf → the temperature responses → A; and indirectly, through the c_i root-find, whose residual reads the temperature-dependent members and so gains an explicit T argument. dprofit_energy_balance_term() adds both, and subtracts the temperature slope of the two costs that read the leaf temperature (ProfitMax’s thermal cost and LeastCost’s capacity term); dprofit_dpsi_leaf() and dprofit_at_collar_psi() call it, and it is exactly zero with the gate off. Below, \(A_T=\partial A/\partial T\) at fixed \(c_i\) and \(A'=\partial A/\partial c_i\) at fixed temperature (not the \(d A/d\psi\) written \(A'\) earlier on this page).

The c_i solve finds the root of

\[g(c_i;\psi_{\rm stem},\psi,T)=A(c_i,T)\,u-g_c(T)\,(c_a-c_i)\,a=0,\]

with \(u\) the µmol-to-mol factor and \(a\) the inverse of atmospheric pressure. So \(g_T=A_T u-g_c'(T)(c_a-c_i)a\), \(g_{c_i}=A' u+g_c a\), and \(dc_i/dT=-g_T/g_{c_i}\). The total temperature slope of assimilation is then

\[\frac{dA}{dT}=A_T+A'\frac{dc_i}{dT}=\frac{A_T\,g_c a+A'\,g_c'(T)\,(c_a-c_i)\,a}{g_{c_i}},\]

and the term added to \(dJ/d\psi\) is that slope times \(\tau=(dT/dE)(dE/d\psi)\), through the curve’s link \(h'(A)\).

Two checks come free. At \(g_c'(T)=0\) the slope is \(A_T\) times the damping factor \(g_c a/g_{c_i}\), which lies strictly in (0, 1): the direct thermal effect is partly undone as \(c_i\) re-equilibrates, so an implementation returning more than \(|A_T\tau|\) is wrong. And \(dE/d\psi>0\) while \(dT/dE<0\), so \(\tau<0\): above the thermal optimum (\(A_T<0\)) the term is positive and moves the optimum drier, toward more transpiration, because the extra flux buys evaporative cooling. That is the decoupling mechanism, and its sign is the sharpest test of the term.

\(g_c'(T)\) is not zero, and a form that puts it under \(A_T\) rather than under \(A'\) agrees with the one above only at \(g_c'(T)=0\), so it passes every test written while that slope was zero.

The maximum is flat, so its location is the delicate quantity

J is flat at its optimum, with curvature k ≈ 1 in J ≈ J* − k(ψ − ψ*)². An error dJ in the value therefore displaces its location by √(dJ/k). The two are not interchangeable as things to check, and that has two practical consequences.

The first is that you should check the residual |dJ/dψ| rather than the profit when judging a change to a solver. Profit is the wrong instrument precisely because it is flat: a change that locates the optimum ten orders of magnitude more precisely can still report a profit 6e-07 lower at a handful of points, simply because the two points differ by less than the objective can resolve. A change that improves the answer can make the objective look worse.

The second is that a relative disagreement in either quantity needs care. Comparing this model’s outputs between platforms, where glibc and Apple’s libm do not agree in the last bits of exp and pow, the reference grid disagrees by 9.05e-13 µmol C m⁻² s⁻¹ in profit and 3.78e-11 Pa in c_i. Those figures are set by the solver’s tolerances rather than by libm: with the inner c_i root-find at 1e-10 rather than 1e-13, both were about 1e-9.

The worst relative figures, 2.13e-11 in profit and 8.6e-08 among the argmax-evaluated fields, land on rows where the quantity concerned is small: a profit of 0.029 against a grid median of 2.47, and an assimilation of 5.4e-06 against a median of 1.44, a leaf whose assimilation barely clears respiration. A relative difference divides by the value, so at such rows the metric is cancelling rather than the model disagreeing.

A finite difference through a solve has a noise floor, so sweep the step

Differencing anything computed by a solver divides that solver’s tolerance by the step, so the relative error grows as 1/h until truncation takes over. The error-versus-h curve is a V, and only its floor is a derivative. A small step is not the safe choice; it is the noise side.

The V’s position is a property of the solver rather than of the model, and that is the part that catches people. How well a solver locates its argmax sets how small a step can still move it.

A root-find on dJ/dψ = 0 locates ψ* to about 1e-15, so a 1e-06 relative parameter step moves it far more than the solve’s own precision, and the floor sits near 1e-06. A search refining on bracket width locates ψ* only to about 1e-06 of the bracket, which is 4e-06 MPa on a 4.4 MPa one, and the same step then differences pure quantisation. Measured on SOX for vcmax_25, against an exact gradient of 0.0551:

At 1e-06 the width-refined difference also gets the sign wrong for JW26. The lesson is that a single step is not a verification: two bugs can be reported against a correct gradient on the strength of one, and were.

Every route in the package root-finds, so a step of 1e-06 sits on the floor for all of them, and the package’s own referee of the gradient sweeps it anyway. Re-sweep if a solver changes, because the number is not a constant of the model.

A scan’s argmax cannot feed a parameter gradient

The section above is about resolution. This one is about a harder obstruction.

A grid argmax is piecewise constant in the parameters. It sits at one grid point until the winner changes cells, so its derivative is not merely imprecise: it is zeros interrupted by jumps. No tolerance fixes that. It is why a bracketing search cannot support a trait gradient however tightly it is run, and why every solver here solves a first-order condition instead. Over eleven parameter steps, golden section gives 6 distinct argmax values where the root-find gives 11, with second differences of 3.9e-04 against 3.4e-07. The larger figure is the step size itself, which is to say nothing.

This is also why ProfitMax’s |A|max is found by a root-find rather than by a sweep of the supply stream: it is differentiated, so it has to be differentiable. Appendix B has what that buys.