Skip to content

[TF24 hydraulics] Penman-Monteith leaf energy balance (opt-in, default off) - #567

Open
dfalster wants to merge 6 commits into
developfrom
feature/penman-monteith-leaf-energy-balance
Open

[TF24 hydraulics] Penman-Monteith leaf energy balance (opt-in, default off)#567
dfalster wants to merge 6 commits into
developfrom
feature/penman-monteith-leaf-energy-balance

Conversation

@dfalster

Copy link
Copy Markdown
Member

Summary

Adds a Penman-Monteith leaf energy balance to the shared TF24/TF24f Leaf
submodel (issue #523), opt-in and default off so it is fully
backward-compatible. Leaf temperature is solved from the hydraulically-pinned
transpiration and fed into the Farquhar temperature scaling, correcting the
Tleaf = Tair assumption that biases both carbon gain A and water cost E
under hot, high-radiation conditions.

Key design point (see notes/penman-monteith/implementation-plan.md): E stays
the Sperry hydraulic supply — PM is not used to re-predict E. Because
A(Tleaf) doesn't feed back into E, Tleaf is an explicit algebraic function
of E: one non-iterative forward pass per operating point, no new optimisation
dimension.

What's here (Steps 1–6 of the plan)

  • Core (1–3): Arrhenius block extracted to update_temperature_dependent_params; leaf_temp_from_E() + es/Δ helpers; use_energy_balance_ gate wired into set_leaf_states_rates_from_psi_stem. TF24f uses its finite-difference gradient when PM is on (the exact AD gradient assumes temperature-independent kinetics — deferred).
  • Interface (4): TF24_Pars gains use_energy_balance + d; new wind_speed driver; ra = C_ra·√(d/U₀); the flag is R-settable. RcppR6/Rcpp bindings regenerated; NEWS breaking-changes entry; model-version snapshot re-accepted (no scientific_version bump — PM off is unchanged).
  • Robustness (PM-only, gated): psi_stem_to_ci falls back to the compensation point when extreme heating leaves no feasible ci; Tleaf clamped to a physical range so a non-equilibrium probe can't make the Arrhenius block non-finite.
  • Evaluation (5): factorial Fick-vs-PM (notes/penman-monteith/). Verdict: keep the PM core — ΔTleaf up to 22 °C (median ~5 °C), A collapses 90–100 % in the hot/high-radiation/high-VPD corner, negligible (<2 °C, <2 % A) in cool/low-light cells.
  • Demo (6): staging leaf-level demo in overstorey_staging/ (.Rbuildignore'd) with a fast smoke test so it can't rot.

Backward compatibility / testing

  • With the gate off, output is bit-identical to before (verified against a leaf-level grid and a full SCM run — offspring to 17 digits).
  • Full devtools::test() suite green (2377 pass, 4 skip). New/updated: test-pm-leaf-demo.R, the two TF24 Defaults tests, and the model-version snapshot.

Not yet done (before PM outputs are used in publications)

Doc §7.5.2–7.5.4: fixed-ra vs wind-profile sensitivity, Rn ±30 % sensitivity, and the annual SCM integration. The wind profile and longwave term remain candidates to drop if they fail their own gates.

Related

A companion PR re-blesses the drifted scenario_baseline.rds and adds a regenerable reference-run generator on chore/refresh-pre-pm-baselines (→ develop); that drift is from earlier merged work (NSC #517, soil #558), not this PR.

🤖 Generated with Claude Code

dfalster and others added 3 commits July 22, 2026 10:51
…n plan

Design note (penman_monteith_implementation.md) plus an implementation plan
scoping how PM / computed leaf temperature would wire into the TF24 and TF24f
strategies via the shared Leaf submodel. Companion to epic #523; the plan also
records the ATLS thermal-damage follow-on (epic #566) and its dependency on PM.

Docs only, no code changes.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…#523)

Adds a Penman-Monteith leaf energy balance to the shared TF24/TF24f Leaf
submodel: leaf temperature is solved from the hydraulically-pinned transpiration
(Tleaf = Tair + (Rn - lambda*E)*ra/(rho*cp)) and fed to the Farquhar temperature
scaling, correcting the Tleaf=Tair assumption that biases A and E under hot,
high-radiation conditions.

Design: E stays the hydraulic supply (no PM inversion), and A(Tleaf) does not
feed back into E, so Tleaf is an explicit algebraic function of E -- a single
non-iterative forward pass per operating point, no new optimisation dimension.

Implementation (minimal core, Steps 1-3):
- extract the Arrhenius block into update_temperature_dependent_params(Tleaf);
- leaf_temp_from_E() + es/Delta helpers + physical constants;
- gate use_energy_balance_ (default off) wired into
  set_leaf_states_rates_from_psi_stem, defeating the photo_temp cache on the PM
  path only. TF24f runs on its finite-difference gradient with PM on (the exact
  AD gradient assumes temperature-independent kinetics; deferred).

Interface (full cut, Step 4): TF24_Pars gains use_energy_balance + d; a
wind_speed extrinsic driver; ra = C_ra*sqrt(d/U0) (fixed-ra fallback); the flag
is R-settable. Regenerated RcppR6/Rcpp bindings.

Robustness (PM path only, gated): psi_stem_to_ci falls back to the compensation
point when extreme heating leaves no feasible ci root; Tleaf is clamped to a
physical range so a non-equilibrium probe cannot drive the Arrhenius block
non-finite.

Backward compatibility: with the gate off, runtime behaviour is bit-identical
(verified vs leaf-level and SCM references) and scientific_version is unchanged;
new pars/driver fields are additive (NEWS breaking-changes entry + model-version
snapshot re-accepted). Full test suite green.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Step 5 (doc 7.5.1): factorial Fick-vs-PM at the leaf level across Tair x VPD x
PAR, reporting dTleaf/dA/dE/dgs, with a regenerable generator + data. Verdict:
KEEP the PM core -- Tleaf departs from Tair by up to 22C (median ~5C) and A
collapses 90-100% in the hot/high-radiation/high-VPD corner, while cool/low-light
cells barely move (<2C, <2% A). Still to run before publication: fixed-ra vs wind
profile, Rn sensitivity, and the annual SCM integration (doc 7.5.2-7.5.4).

Step 6: a staging leaf-level demo (overstorey_staging/, .Rbuildignore'd) --
profit anatomy, optimal outcomes across environments, and a where-PM-matters
heatmap, non-PM vs PM. Leaf-driving code is shared with a fast smoke test
(test-pm-leaf-demo.R) so the demo cannot silently rot.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
dfalster added a commit that referenced this pull request Jul 22, 2026
… + regenerable reference run (#568)

## Summary

Housekeeping companion to #567 (Penman-Monteith), kept separate so the
baseline
refresh is **not** attributed to PM (PM is off / bit-identical here).

- **Re-bless `scenario_baseline.rds`.** The recorded gateway baseline
was stamped
at `5879d73` (pkg 2.0.0.9001) and had drifted: 5 scenarios it recorded
as
`failure` now **persist**, from work merged since (NSC storage
#517/#554,
layered soil #558). That left the opt-in gateway test
(`PLANT_RUN_SCENARIOS=1`)
red for reasons unrelated to any one change. Regenerated on `develop`
via
  `make bless-scenarios` (mpl=100) so the diff is meaningful again.
**Note for reviewers:** the scenarios' a-priori `expected` column (from
`inst/scenarios/`) still encodes the old expectations, so `match_rate`
is now
  0.375 — worth revisiting those expectations separately.

- **Add `scripts/generate_reference_run.R`.** Replaces the former
hand-saved
`results_high.RDS`, which predated the `pars` refactor (#410), could no
longer
be loaded/replayed, and was in fact gitignored (never a committed
artifact).
The reproducible deliverable is the generator: a canonical
single-species TF24
  SCM run whose output round-trips against the current model.

🤖 Generated with [Claude Code](https://claude.com/claude-code)

---------

Co-authored-by: Claude Fable 5 <noreply@anthropic.com>

@itowers1 itowers1 left a comment

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Couple of comments, but ready to go

// shortwave absorbed ~= 2 * absorbed PAR (PAR ~= 50% of shortwave; doc 3.3)
static const double sw_abs_per_par = 2.0;
// fixed net longwave (cooling) offset, W m^-2 (clear-sky approximation; doc 3.3)
static const double longwave_net_offset = -40.0;

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Need to investigate this further

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Outgoing longwave radiation is actually a function of temperature, but this would cause another root-finding operation (I believe)

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Agreed — the fixed −40 W/m² offset is a Tleaf-independent approximation, and a faithful outgoing-longwave term makes Rn implicit in Tleaf (root-find), which would break this PR's single non-iterative forward pass. Raised as #581 with a closed-form Stefan-Boltzmann linearisation option (adds a radiative-conductance term to the balance, no iteration), to be gated by the Step 5 longwave sensitivity check per the plan rather than expanding this PR. Noted at the constant in 76df716.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Agreed — the fixed −40 W/m² offset is a Tleaf-independent approximation, and a faithful outgoing-longwave term makes Rn implicit in Tleaf (root-find), which would break this PR's single non-iterative forward pass. Raised as #581 with a closed-form Stefan-Boltzmann linearisation option (adds a radiative-conductance term to the balance, no iteration), to be gated by the Step 5 longwave sensitivity check per the plan rather than expanding this PR. Noted at the constant in 76df716.

Comment thread src/leaf_model.cpp
// (doc 4.1), with the per-strategy leaf dimension d_ and the above-canopy wind
// wind_speed_. Fall back to the fixed value when the wind model is unusable
// (non-finite / non-positive inputs), e.g. a bare Leaf that set neither.
ra_ = (std::isfinite(d_) && std::isfinite(wind_speed_) &&

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This will allow the model to work when wind speed is 0 OR NA, I understand the 0 case, but possibly just want it to fail if wind speed is NA.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good catch — fixed in 76df716. set_physiology now errors on a non-finite wind_speed_ (or d_) on the energy-balance path rather than silently using the fixed ra; zero wind/d (physically ra → ∞) stays a legitimate fallback. Added a test covering all three cases (NA → error, PM-off NA → ok, zero-wind → fallback).


In practice:
```
LW_net ≈ ε·σ·(εa·Ta⁴ - Ta⁴) = ε·σ·Ta⁴·(εa - 1)

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This could potentially be implemented in this PR or raised as an issue.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Agreed — the fixed −40 W/m² offset is a Tleaf-independent approximation, and a faithful outgoing-longwave term makes Rn implicit in Tleaf (root-find), which would break this PR's single non-iterative forward pass. Raised as #581 with a closed-form Stefan-Boltzmann linearisation option (adds a radiative-conductance term to the balance, no iteration), to be gated by the Step 5 longwave sensitivity check per the plan rather than expanding this PR. Noted at the constant in 76df716.

…ongwave (#523)

Addresses review feedback from @itowers1 on #567:

- set_physiology now errors when wind_speed_ or d_ is non-finite on the
  energy-balance path (a broken driver / unset trait), instead of silently
  falling back to the fixed ra. Zero wind or zero d (physically ra -> infinity)
  stays a legitimate case and keeps the fixed-ra fallback. Gated on
  use_energy_balance_, so the non-PM path is unchanged. New test covers all three
  cases (NA -> error, PM-off NA -> ok, zero-wind -> fallback).
- The fixed longwave offset is a Tleaf-independent approximation; a faithful
  temperature-dependent treatment would make Rn implicit in Tleaf. Documented at
  the constant and tracked as #581 (with a non-iterative linearisation option),
  to be gated by the Step 5 sensitivity protocol rather than expanding this PR.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants