SurEau vegetation processes
Phenology
compute_pheno
def compute_pheno(
state:SurEauPlantState, # Plant state object. The following fields are read/written:
# - LAI_pheno : float
# Current phenological LAI (m²/m²). Written.
# - LAI : float
# Current actual LAI (m²/m²). Read during senescence.
# - sum_temperature : float
# Cumulative forcing temperature (°C). Read/written.
# - budburst_date : float
# Day of year of budburst (int or NaN). Read/written.
params:SurEauVegetationParams, # Vegetation parameters. Relevant attributes:
# - foliage : str
# ``"Evergreen"``, ``"Deciduous"``, or ``"Forced"``.
# - LAI_max : float
# Maximum leaf area index (m²/m²).
# - T_base : float
# Minimum temperature for forcing accumulation (°C).
# - F_crit : float
# Cumulative forcing threshold for budburst (°C).
# - day_start : int
# First DOY for temperature accumulation.
# - nb_day_LAI : int
# Days from budburst to full canopy.
# - day_start_forced : int
# Fixed budburst DOY for Forced model.
# - day_end_forced : int
# Fixed senescence DOY for Forced model.
temperature:float, # Daily mean air temperature (°C).
DOY:int, # Day of year (1–365).
)->SurEauPlantState: # Updated plant state with modified ``LAI_pheno``,``sum_temperature``, ``budburst_date``.
Update the phenological leaf area fraction (LAI_pheno) based on the seasonal cycle of leaf growth and loss.
Implements the leaf phenology module described in Appendix A of Ruffault et al. (2022). The function controls ∅ (the phenological fraction, 0 to 1) which feeds into Equation A1:
LAI = ∅ × LAI_maxEq. A1 itself is computed downstream in update_LAI_and_stocks. This function determines ∅ by tracking the progression through three seasonal phases: dormancy → leaf-out → full canopy → senescence.
Three foliage types are supported:
Evergreen: LAI_pheno is set to LAI_max on DOY 1 and never changes.The canopy is always at full capacity.
Deciduous: Implements Equation A2 for budburst timing:
Σ_{t0}^{tf} R_f(T_d) ≥ F*
Each day after day_start (t0), if the temperature exceeds T_base (T_D), the temperature is added to a running sum.When the sum exceeds F_crit (F*), budburst is triggered. After budburst, LAI_pheno increases linearly over nb_day_LAI days at rate LAI_max / nb_day_LAI (the paper’s R_LAI). In autumn (DOY ≥ 280), LAI_pheno decreases at the same rate until zero.
Forced: Same as Deciduous but budburst is triggered by a fixed calendar date (day_start_forced) instead of temperature accumulation. Senescence is triggered by day_end_forced instead of DOY 280. Useful for sensitivity experiments where phenology should be decoupled from climate.
Physical analogy: Think of the tree as a solar panel factory with a seasonal schedule. Evergreen factories run 24/7. Deciduous factories shut down in winter and wait for enough accumulated warmth to restart — like a bank account for warmth that must reach a threshold before the factory reopens. Forced factories follow a fixed calendar regardless of weather.
Parameters:
state: SurEauPlantState object. The following fields are read and/or written:
LAI_pheno: Current phenological LAI (m²/m²). Written each call. Set to LAI_max for Evergreen on DOY 1, incremented during leaf-out, decremented during senescence.
LAI: Current actual LAI (m²/m²). Read during senescence to compute the daily decrement. May differ from LAI_pheno if cavitation-induced defoliation has occurred.
sum_temperature: Cumulative forcing temperature (°C). Read and written during the Deciduous accumulation phase. Reset to 0 on DOY 1.
budburst_date: Day of year when budburst occurred (int or NaN). Written when the temperature sum exceeds F_crit. Reset to NaN on DOY 1 for non-Evergreen species.
params: SurEauVegetationParams object with the following attributes:
foliage: Foliage type.
"Evergreen","Deciduous", or"Forced".LAI_max: Maximum leaf area index of the stand (m²/m²). E.g. 4.5 for Quercus ilex, 5.0 for Fagus sylvatica.
T_base: Minimum temperature to start cumulating forcing temperatures (°C). Equivalent to T_D in Eq. A2. E.g. 3.0.
F_crit: Amount of cumulative forcing temperature required to trigger budburst (°C). Equivalent to F* in Eq. A2. E.g. 450.
day_start: First day of year to begin temperature accumulation (DOY). Equivalent to t0 in Eq. A2. E.g. 55 (late February).
nb_day_LAI: Number of days from budburst to full canopy (days). Controls both the leaf-out growth rate and the senescence decline rate: rate = LAI_max / nb_day_LAI. Equivalent to R_LAI in text near Eq. A1. E.g. 21.
day_start_forced: Fixed budburst date for the Forced model (DOY). Only used when
foliage="Forced". E.g. 40.day_end_forced: Fixed senescence start date for the Forced model (DOY). Only used when
foliage="Forced". E.g. 220.
temperature: Daily mean air temperature (°C). Used for temperature accumulation in the Deciduous model (Eq. A2). Corresponds to
clim_day.T_air_mean.DOY: Day of year (1–365). Controls the yearly reset (DOY 1), the start of accumulation (
day_start), and the onset of senescence (DOY 280 for Deciduous,day_end_forcedfor Forced).
Returns:
- state: Updated SurEauPlantState with modified
LAI_pheno,sum_temperature, andbudburst_date.
Seasonal phases for Deciduous:
DOY 1 Reset: LAI_pheno=0, sum_T=0, budburst=NaN
│
▼
Phase 1: Dormancy (DOY < day_start or T < T_base)
│ Nothing happens, waiting for warmth
│
▼
Phase 1b: Accumulation (DOY ≥ day_start and T > T_base)
│ sum_T += temperature each day
│ (Eq. A2: Σ R_f(T_d))
│
▼ when sum_T > F_crit (Eq. A2: ≥ F*)
│
Phase 2: Leaf-out (budburst_date ≤ DOY < budburst + nb_day_LAI)
│ LAI_pheno += LAI_max / nb_day_LAI each day
│ (rate = R_LAI)
│
▼ when DOY = budburst + nb_day_LAI
│
Summer plateau (budburst + nb_day_LAI ≤ DOY < 280)
│ LAI_pheno stays constant (≈ LAI_max)
│ No conditions match — function does nothing
│
▼ when DOY ≥ 280 (≈ October 7)
│
Phase 3: Senescence LAI_pheno decreases by LAI_max / nb_day_LAI
│ each day until LAI_pheno = 0
▼
Back to DOY 1 next yearNumerical example (Deciduous, LAI_max = 5.0, T_base = 3.0, F_crit = 450, day_start = 55, nb_day_LAI = 21):
Daily growth rate = 5.0 / 21 = 0.238 m²/m²/day
DOY 55: accumulation begins (day_start)
DOY 56: T = 2°C → below T_base, no deposit. sum_T = 0
DOY 57: T = 7°C → deposit 7. sum_T = 7
DOY 58: T = 12°C → deposit 12. sum_T = 19
... (roughly 8°C/day average contribution)
DOY 111: sum_T = 452 > 450 → budburst! budburst_date = 111
DOY 112: LAI_pheno = 0.238 (first day of leaf-out)
DOY 113: LAI_pheno = 0.476
...
DOY 132: LAI_pheno = 5.0 (full canopy after 21 days)
DOY 133–279: summer plateau, LAI_pheno = 5.0 (no change)
DOY 280: senescence begins
DOY 280: LAI_pheno = 5.0 − 0.238 = 4.762
DOY 281: LAI_pheno = 4.762 − 0.238 = 4.524
...
DOY 301: LAI_pheno ≈ 0 (fully bare after 21 days)Capacitances
update_capacitances
def update_capacitances(
state:SurEauPlantState, # Plant state object.
params:SurEauVegetationParams, # Vegetation parameters. Relevant attributes:
# - pi_full_turgor_leaf, pi_full_turgor_stem : float
# Osmotic potential at full turgor π₀ (MPa).
# - epsilon_sym_leaf, epsilon_sym_stem : float
# Bulk elastic modulus ε (MPa).
# - psi_TLP_leaf, psi_TLP_stem : float
# Turgor loss point (MPa).
# - C_LApo_init, C_SApo_init : float
# Constant apoplasm capacitances (mmol m⁻² MPa⁻¹).
)->SurEauPlantState: # Updated plant state with modified capacitances.
Update symplasmic and apoplasmic capacitances for all four plant compartments.
Capacitance C answers the question: “If the water potential drops by a small amount dψ, how much water is released (dW)?” It is the core of Equation 3 (C = dQ/dψ) and enters the water balance ODEs Equations 6–9 as the coefficient on dψ/dt:
Eq. 6: C_LApo × dψ_LApo/dt = fluxes (leaf apoplasm)
Eq. 7: C_SApo × dψ_SApo/dt = fluxes (stem apoplasm)
Eq. 8: C_LSym × dψ_LSym/dt = fluxes (leaf symplasm)
Eq. 9: C_SSym × dψ_SSym/dt = fluxes (stem symplasm)A large C means the tissue is a good buffer, it releases water readily without crashing its water potential (like a soft sponge). A small C means the tissue is stiff, ψ drops quickly when water is lost (like a rigid pipe).
The symplasmic capacitances (C_LSym, C_SSym) are variable, they depend on the current water potential via the pressure-volume curve. This function implements Equation 42 to compute them:
C = Q_sat × dRWC/dψ (Eq. 42)where Q_sat is the total water at saturation (mmol per unit leaf area), and dRWC/dψ is the derivative of the pressure-volume curve (Equation 44), which has two regimes separated by the turgor loss point (TLP):
Above TLP (ψ > ψ_tlp):
RWC
RWC' = ────────────────────── (Eq. 44, case 1)
−π₀ − ψ − ε + 2ε·RWCBelow TLP (ψ < ψ_tlp):
−π₀
RWC' = ───── (Eq. 44, case 2)
ψ²The apoplasmic capacitances (C_LApo, C_SApo) are set to constants. In the full model, the apoplasm releases water primarily through cavitation (Eqs. 24–27), which is handled separately via compute_PLC and compute_PLC_prime. The constant here represents only the small residual elastic compliance of the cell walls.
Physical analogy: The symplasm is a flexible water balloon, squeeze it (lower ψ) and how much water comes out depends on how full the balloon is and how stiff the walls are (ε). The apoplasm is a rigid pipe network, its primary water release mechanism is cavitation (pipes breaking), handled elsewhere. The tiny elastic compliance of the pipe walls is the constant here.
Important: The derivative dRWC/dψ is discontinuous at the turgor loss point. Above TLP, dψ/dRWC = −π₀/RWC² + ε; below TLP, dψ/dRWC = −π₀/RWC². The jump is exactly ε. This means the capacitance jumps upward at TLP: below TLP, only osmotic pressure resists dehydration, so more water is released per unit ψ drop.
Parameters:
state: SurEauPlantState object. The following fields are read and/or written:
Read:
- psi_LSym: Water potential of the leaf symplasm (MPa).
- psi_SSym: Water potential of the stem symplasm (MPa).
- Q_LSym_sat_mmol_per_LA: Saturated water content of the leaf symplasm (mmol per unit leaf area). Q_sat in Eq. 42.
- Q_SSym_sat_mmol_per_LA: Saturated water content of the stem symplasm (mmol per unit leaf area). Q_sat in Eq. 42.
- LAI: Current leaf area index (m²/m²). Used to zero out leaf capacitance when no leaves are present.
Written:
- C_LSym: Leaf symplasm capacitance (mmol m⁻² MPa⁻¹).
- C_SSym: Stem symplasm capacitance (mmol m⁻² MPa⁻¹).
- C_LApo: Leaf apoplasm capacitance (mmol m⁻² MPa⁻¹). Set to constant.
- C_SApo: Stem apoplasm capacitance (mmol m⁻² MPa⁻¹). Set to constant.
params: SurEauVegetationParams object with the following attributes:
- pi_full_turgor_leaf: Osmotic potential at full turgor π₀ for leaves (MPa). Negative, e.g. −2.1.
- epsilon_sym_leaf: Bulk elastic modulus ε for leaf symplasm (MPa). Positive, e.g. 10.
- psi_TLP_leaf: Turgor loss point for leaves (MPa). Negative, e.g. −2.66. Computed from π₀ and ε.
- pi_full_turgor_stem: Same as above, for stem.
- epsilon_sym_stem: Same as above, for stem.
- psi_TLP_stem: Same as above, for stem.
- C_LApo_init: Constant leaf apoplasm capacitance (mmol m⁻² MPa⁻¹).
- C_SApo_init: Constant stem apoplasm capacitance (mmol m⁻² MPa⁻¹).
Returns:
- state: Updated SurEauPlantState with modified
C_LSym,C_SSym,C_LApo,C_SApo.
References:
- Ruffault et al. (2022), Eq. 3 (definition of capacitance), Eq. 42 (symplasmic capacitance = Q_sat × RWC’), Eqs. 43–47 (pressure-volume curve and its derivatives), Eqs. 6–9 (water balance ODEs where capacitances appear).
Derivation of the above-TLP derivative (Eq. 44, case 1):
Eq. 46 (product form) ψ·RWC + π₀·RWC + ε(1−RWC)·RWC − π₀ = 0
│ differentiate both sides w.r.t. ψ (implicit differentiation)
▼
Eq. 47 ψ·RWC' + RWC + π₀·RWC' + ε(1−RWC)·RWC'
− ε·RWC·RWC' = 0
│ factor out RWC'
▼
Rearrange RWC' × [ψ + π₀ + ε − 2ε·RWC] + RWC = 0
│ solve for RWC'
▼
Eq. 44 (case 1) RWC' = RWC / (−π₀ − ψ − ε + 2ε·RWC)
│ ▲ ▲
│ └── this is exactly the └── this is exactly the
│ code's RWC/denom code's denom
▼
Code denom = −π₀ − ψ − ε + 2ε·RWC
RWC_prime = RWC / denomEquivalence to the alternative form RWC²/(−π₀ + ε·RWC²):
Substitute ψ = −π₀ − ε(1−RWC) + π₀/RWC (Eq. 43) into denom:
denom = −π₀ − [−π₀ − ε(1−RWC) + π₀/RWC] − ε + 2ε·RWC
= −π₀ + π₀ + ε(1−RWC) − π₀/RWC − ε + 2ε·RWC
= ε − ε·RWC − π₀/RWC − ε + 2ε·RWC
= ε·RWC − π₀/RWC
= (ε·RWC² − π₀) / RWC
= (−π₀ + ε·RWC²) / RWC
So: RWC / denom = RWC / [(−π₀ + ε·RWC²)/RWC]
= RWC² / (−π₀ + ε·RWC²) Derivation of the below-TLP derivative (Eq. 44, case 2):
Eq. 43 (below TLP) ψ = π₀ / RWC (turgor = 0, osmotic only)
│ solve for RWC
▼
Invert RWC = π₀ / ψ
│ differentiate w.r.t. ψ
▼
Eq. 44 (case 2) RWC' = −π₀ / ψ²
│ since π₀ < 0, −π₀ > 0, so RWC' > 0
▼
Code RWC_prime = −π₀ / ψ²Computation flow:
For each organ (leaf, stem):
┌──────────────────────────────────────────────────────────┐
│ Step 1: Get RWC from PV curve │
│ RWC = 1 − compute_Rs(π₀, ε, ψ−dbxmin) │
│ (inverse of Eq. 43 via Eqs. 45/46) │
│ │
│ Step 2: Compute dRWC/dψ (Eq. 44) │
│ if ψ > ψ_TLP: │
│ denom = −π₀ − ψ − ε + 2ε·RWC │
│ RWC' = RWC / denom (Eq. 44, case 1) │
│ else: │
│ RWC' = −π₀ / ψ² (Eq. 44, case 2) │
│ │
│ Step 3: Scale to capacitance (Eq. 42) │
│ C = Q_sat × RWC' │
└──────────────────────────────────────────────────────────┘Numerical example (leaf symplasm, π₀ = −2.1, ε = 10, ψ_TLP = −2.66, Q_sat = 50 mmol/m²):
Above TLP (ψ = −1.5 MPa):
Rs = 0.1211, RWC = 0.8789
denom = −(−2.1) − (−1.5) − 10 + 2×10×0.8789
= 2.1 + 1.5 − 10 + 17.578 = 11.179
RWC' = 0.8789 / 11.179 = 0.0786
C = 50 × 0.0786 = 3.93 mmol/m²/MPa
Below TLP (ψ = −4.0 MPa):
RWC' = −(−2.1) / (−4.0)² = 2.1 / 16 = 0.1313
C = 50 × 0.1313 = 6.56 mmol/m²/MPaSpecial cases:
LAI = 0 (deciduous tree in winter):
C_LSym = 0 (no leaves → no leaf symplasm water storage)
C_SSym computed normally (stem always exists)
|denom| ≤ dbxmin (numerical guard):
RWC' = 0 (treat compartment as infinitely stiff for
one timestep — safe numerical default)
|ψ| ≤ dbxmin (near zero, below-TLP guard):
RWC' = 0 (avoid division by zero)LAI and water stocks
update_LAI_and_stocks
def update_LAI_and_stocks(
state:SurEauPlantState, # Plant state object.
params:SurEauVegetationParams, # Vegetation parameters object.
)->SurEauPlantState: # Updated plant state with all LAI-dependent fields recalculated.
Update actual LAI from phenology and cavitation-induced defoliation, then recompute all LAI-dependent water storage volumes and capacitances.
This function is the daily bookkeeper of the canopy. It sits between compute_pheno (which sets the phenological LAI target) and the hourly solver (which needs water storage capacities to solve the ODEs). Called once per day in run_sureau:
The function performs three tasks in sequence:
Compute actual LAI by subtracting drought-killed leaves from the phenological target (modified Eq. A1), plus derived canopy properties (fractional cover, interception capacity, dry matter).
Compute saturated water content Q_sat for all four plant compartments — leaf symplasm (Eq. 36), leaf apoplasm (Eq. 37), stem symplasm (Eq. 39), stem apoplasm (Eq. 40) — using leaf water content from the succulence relation (Eq. 38).
Update capacitances by delegating to
update_capacitances(Eqs. 42, 44), which needs the Q_sat values computed in step 2.
update_LAI_and_stocks is the daily inventory check, counting how many leaves are actually working (some may have been destroyed by drought), and recalculating the water storage capacity based on the current LAI.
Parameters:
state: SurEauPlantState object. The following fields are read and/or written:
Read:
- LAI_pheno: Phenological LAI target (m²/m²). Set by
compute_pheno. - PLC_leaf: Percent loss of leaf hydraulic conductivity (%). Used for cavitation-induced defoliation.
Written (LAI and canopy properties):
- LAI_dead: LAI of drought-killed leaves (m²/m²).
- LAI: Actual functional LAI (m²/m²).
LAI_pheno − LAI_dead. - FCC: Fractional canopy cover (0–1). Beer-Lambert law.
- canopy_storage_capacity: Rainfall interception capacity (mm).
- DM_live_canopy: Dry mass of living leaves (g/m² ground).
- DM_dead_canopy: Dry mass of dead leaves (g/m² ground).
Written (saturated water contents, Eqs. 36–40):
- Q_LSym_sat_L: Leaf symplasm water at saturation (L/m²).
- Q_LSym_sat_mmol: Same, in mmol/m² ground.
- Q_LSym_sat_mmol_per_LA: Same, in mmol/m² leaf area. Feeds Eq. 42 in
update_capacitances. - Q_SSym_sat_L, Q_SSym_sat_mmol, Q_SSym_sat_mmol_per_LA: Same three units for stem symplasm.
- Q_LApo_sat_L, Q_LApo_sat_mmol, Q_LApo_sat_mmol_per_LA: Same three units for leaf apoplasm.
- Q_SApo_sat_L, Q_SApo_sat_mmol, Q_SApo_sat_mmol_per_LA: Same three units for stem apoplasm.
Written (capacitances, via
update_capacitances):- C_LSym, C_SSym, C_LApo, C_SApo: Capacitances for all four compartments (mmol m⁻² MPa⁻¹).
- LAI_pheno: Phenological LAI target (m²/m²). Set by
params: SurEauVegetationParams object with the following attributes:
Defoliation:
- defoliation: Whether to enable cavitation-induced leaf shedding (bool). Default False.
Canopy optics and interception:
- K: Light extinction coefficient (dimensionless, e.g. 0.5). Used in the Beer-Lambert law for fractional canopy cover.
- canopy_storage_param: Rainfall interception capacity per unit LAI (mm per m²/m²).
Leaf traits:
- LMA: Leaf Mass per Area (g/m²). Converts LAI to dry mass.
- LDMC: Leaf Dry Matter Content (mg/g). Ratio of dry mass to fresh mass, Used to compute succulence (Eq. 38).
- apo_frac_leaf: Apoplasmic fraction of leaf water (dimensionless, 0–1). α_LApo in Eqs. 36–37.
Stem traits:
- vol_stem: Total stem water volume (L/m² ground). Combines V_S / M_H2O × α_Water from Eqs. 39–40.
- sym_frac_stem: Symplasmic fraction of stem water (dimensionless). α_SSym in Eq. 39.
- apo_frac_stem: Apoplasmic fraction of stem water (dimensionless). α_SApo in Eq. 40.
Passed through to
update_capacitances:- pi_full_turgor_leaf, pi_full_turgor_stem: Osmotic potential at full turgor π₀ (MPa).
- epsilon_sym_leaf, epsilon_sym_stem: Bulk elastic modulus ε (MPa).
- psi_TLP_leaf, psi_TLP_stem: Turgor loss point (MPa).
- C_LApo_init, C_SApo_init: Constant apoplasm capacitances.
Returns:
- state: Updated SurEauPlantState with all LAI-dependent fields recalculated.
Computation flow:
┌─────────────────────────────────────────────────────────┐
│ BLOCK 1: Defoliation (text, p. 5614) │
│ │
│ if defoliation enabled AND PLC_leaf > 10%: │
│ LAI_dead = LAI_pheno × PLC_leaf / 100 │
│ else: │
│ LAI_dead = 0 │
└───────────────────────┬─────────────────────────────────┘
│
▼
┌─────────────────────────────────────────────────────────┐
│ BLOCK 2: Actual LAI and canopy properties │
│ │
│ LAI = LAI_pheno − LAI_dead (modified Eq. A1) │
│ FCC = 1 − exp(−K × LAI) (Beer-Lambert) │
│ canopy_storage = param × LAI │
│ DM_live = LAI × LMA │
│ DM_dead = LAI_dead × LMA │
└───────────────────────┬─────────────────────────────────┘
│
▼
┌─────────────────────────────────────────────────────────┐
│ BLOCK 3: Saturated water contents (Eqs. 36–40) │
│ │
│ succulence = 1/LDMC − 1 (Eq. 38) │
│ Q_L_total = succulence × DM_live (Eq. 38) │
│ │
│ Q_LSym = Q_L_total × (1 − α_LApo) (Eq. 36) │
│ Q_LApo = Q_L_total × α_LApo (Eq. 37) │
│ Q_SSym = vol_stem × α_SSym (Eq. 39) │
│ Q_SApo = vol_stem × α_SApo (Eq. 40) │
│ │
│ All four converted: L → mmol → mmol/m²_leaf │
└───────────────────────┬─────────────────────────────────┘
│
▼
┌─────────────────────────────────────────────────────────┐
│ BLOCK 4: Capacitances (Eqs. 42, 44) │
│ │
│ update_capacitances(state, params) │
│ → C_LSym, C_SSym = Q_sat × dRWC/dψ (variable) │
│ → C_LApo, C_SApo = constants │
└─────────────────────────────────────────────────────────┘Unit conversion chain:
succulence × DM [g/m²] → Q_sat [g/m² ground]
│ ÷ 1000
▼
Q_sat [L/m² ground] (= mm water depth)
│ × 1e6 / 18
▼
Q_sat [mmol/m² ground]
│ ÷ LAI
▼
Q_sat [mmol/m² leaf] → feeds Eq. 42 (C = Q_sat × RWC')
The factor 1e6/18 converts liters to millimoles:
1 L = 1000 g H₂O
1 mol H₂O = 18 g
1 L = 1000/18 mol = 1e6/18 mmol ≈ 55,556 mmolNumerical example (Quercus ilex: LDMC = 515 mg/g, LMA = 110 g/m², apo_frac_leaf = 0.5, vol_stem = 20 L/m², sym_frac_stem = 0.5, apo_frac_stem = 0.5, LAI = 4.5):
Step 1 — Defoliation (defoliation=False):
LAI_dead = 0
LAI = 4.5
Step 2 — Canopy properties:
DM_live = 4.5 × 110 = 495 g/m²
Step 3 — Leaf symplasm (Eqs. 36, 38):
LDMC = 515/1000 = 0.515
succulence = 1/0.515 − 1 = 0.942 g_water/g_dry
Q_L_total = 0.942 × 495 = 466.3 g/m²
Q_LSym = 466.3 × (1 − 0.5) / 1000 = 0.2332 L/m²
Q_LSym_mmol = 0.2332 × 1e6/18 = 12,953 mmol/m²
Q_LSym_mmol_per_LA = 12,953 / 4.5 = 2,878 mmol/m²_leaf
Step 4 — Stem symplasm (Eq. 39):
Q_SSym = 20 × 0.5 = 10.0 L/m²
Q_SSym_mmol = 10.0 × 1e6/18 = 555,556 mmol/m²
Q_SSym_mmol_per_LA = 555,556 / 4.5 = 123,457 mmol/m²_leaf
Note: Stem water per leaf area (123,457) is ~43× larger than
leaf water per leaf area (2,878). The stem is the plant's main
water "savings account."Defoliation numerical example (PLC_leaf = 40%, LAI_pheno = 4.5):
LAI_dead = 4.5 × 40/100 = 1.8
LAI = 4.5 − 1.8 = 2.7 (40% of canopy shed)
All Q_sat values for leaves scale proportionally:
DM_live = 2.7 × 110 = 297 g/m² (was 495)
Q_LSym_sat_L = 0.1399 L/m² (was 0.2332)
But stem water volumes are UNCHANGED:
Q_SSym_sat_L = 10.0 L/m² (same)
And per-leaf-area stem water INCREASES:
Q_SSym_mmol_per_LA = 555,556 / 2.7 = 205,761 (was 123,457)
This is the survival benefit of defoliation: fewer leaves
means each remaining leaf has access to more of the stem's
water reserve, buying time during drought.k_plant update
update_kplant
def update_kplant(
state:SurEauPlantState, # Updated plant state with modified conductances.
soil:SurEauSoil, params:SurEauVegetationParams, # Vegetation parameters.
)->SurEauPlantState:
Update plant hydraulic conductances accounting for cavitation and soil-to-root series resistance.
Implements Equations 13, 14, 16, and 20 of Ruffault et al. (2022). These equations define the “pipe sizes” that control how easily water flows from soil to roots to stem to leaves. Cavitation progressively narrows the pipes as water potential drops, so this function is called every sub-daily timestep.
The plant is modeled as three pipe segments in series:
Soil layers ──┐
(parallel, │ k_soil_to_stem k_SLApo k_LSym
Eq. 20 each) ├──►[soil→root→stem]──►[stem→leaf]──►[apo→sym]
│ Eq. 16 Eq. 14
Soil layers ──┘
Total: k_plant = 1/(1/Σk_RSApo + 1/k_SLApo + 1/k_LSym) (Eq. 13)The root segments from different soil layers act in parallel (water enters from all layers simultaneously, conductances add), while the three pathway segments act in series (water must pass through all three sequentially, resistances add).
Cavitation reduces the xylem conductances via the vulnerability curve (Eq. 15, implemented in compute_PLC):
K_actual = K_max × (1 − PLC/100)This creates the positive feedback loop driving hydraulic failure: more negative ψ → higher PLC → lower conductances → less water supply → even more negative ψ → even higher PLC → mortality.
Parameters:
state: SurEauPlantState object. The following fields are read and/or written:
Read:
- PLC_stem: Percent loss of stem hydraulic conductivity (%). Used in Eq. 16 to reduce root-to-stem conductance.
- PLC_leaf: Percent loss of leaf hydraulic conductivity (%). Used in Eq. 14 to reduce stem-to-leaf conductance.
- k_LSym: Leaf symplasm conductance (mmol m⁻² s⁻¹ MPa⁻¹). Set during initialization, does not degrade with cavitation (membrane resistance, not xylem).
Written:
- k_RSApo: Root-to-stem apoplasmic conductance per soil layer (mmol m⁻² s⁻¹ MPa⁻¹). Array, one value per soil layer. Eq. 16.
- k_SLApo: Stem-to-leaf apoplasmic conductance (mmol m⁻² s⁻¹ MPa⁻¹). Scalar. Eq. 14.
- k_soil_to_stem: Combined soil-to-stem conductance per soil layer (mmol m⁻² s⁻¹ MPa⁻¹). Array, one value per soil layer. Eq. 20.
- k_plant: Total plant conductance from soil to leaf symplasm (mmol m⁻² s⁻¹ MPa⁻¹). Scalar. Eq. 13.
soil: SurEauSoil object with the following attribute:
- k_soil: Soil-to-root-surface conductance per soil layer (mmol m⁻² s⁻¹ MPa⁻¹). Array, one value per soil layer. Computed from Gardner-Cowan geometry and van Genuchten hydraulic conductivity (Eq. 21). Depends on current soil moisture.
params: SurEauVegetationParams object with the following attributes:
- k_RSApo_init: Maximum root-to-stem conductance per soil layer (mmol m⁻² s⁻¹ MPa⁻¹). Array. From
distribute_conductances(Eq. 17): K_max_j = RAI_j × K_R→SApo. Layers with more roots have higher values. - k_SLApo_init: Maximum stem-to-leaf apoplasmic conductance (mmol m⁻² s⁻¹ MPa⁻¹). Scalar. From
distribute_conductances.
- k_RSApo_init: Maximum root-to-stem conductance per soil layer (mmol m⁻² s⁻¹ MPa⁻¹). Array. From
Returns:
- state: Updated SurEauPlantState with modified
k_RSApo,k_SLApo,k_soil_to_stem,k_plant.
References:
- Ruffault et al. (2022), Eq. 13 (total plant conductance), Eq. 14 (stem-to-leaf with embolism), Eq. 16 (root-to-stem with embolism), Eq. 20 (soil-to-stem series combination).
- Eqs. 15, 17 are prerequisites computed elsewhere (
compute_PLC,distribute_conductances).
Equation-to-code mapping:
Eq. 16 K_{Rj→SApo} = K_{Rj→SApo,max} × (100 − PLC_S) / 100
│ k_RSApo = k_RSApo_init × (1 − PLC_stem/100)
│ Array: one value per soil layer
│ Uses STEM PLC (root xylem connects to stem xylem)
▼
Eq. 14 K_{SApo→LApo} = K_{SApo→LApo,max} × (100 − PLC_L) / 100
│ k_SLApo = k_SLApo_init × (1 − PLC_leaf/100)
│ Scalar: one trunk-to-canopy pathway
│ Uses LEAF PLC (leaf xylem vulnerability)
▼
Eq. 20 K_{soil_j→SApo} = 1 / (1/K_{Rj→SApo} + 1/K_{soil_j→Rj})
│ k_soil_to_stem = kseries(k_soil, k_RSApo)
│ Array: series combination per layer
│ Two bottlenecks: soil drying + xylem cavitation
▼
Eq. 13 K_Plant = 1 / (1/Σ_j K_{Rj→SApo} + 1/K_{SL} + 1/K_{LSym})
│ k_plant = 1/(1/Σk_RSApo + 1/k_SLApo + 1/k_LSym)
│ Scalar: total soil-to-leaf-symplasm throughput
▼
Feeds into Eqs. 6–9 (water balance ODEs)
and calculate_Ebound_Granier (transpiration bound)Numerical example (3 soil layers, PLC_stem = 20%, PLC_leaf = 30%):
Suppose k_RSApo_init = [0.6, 0.3, 0.1] (more roots shallow)
k_SLApo_init = 2.0
k_LSym = 3.0
k_soil = [1.5, 1.0, 0.5] (drier soil at depth)
Step 1 — Eq. 16 (root-to-stem, reduce by stem PLC):
k_RSApo = [0.6, 0.3, 0.1] × (1 − 20/100)
= [0.6, 0.3, 0.1] × 0.80
= [0.48, 0.24, 0.08]
Step 2 — Eq. 14 (stem-to-leaf, reduce by leaf PLC):
k_SLApo = 2.0 × (1 − 30/100)
= 2.0 × 0.70
= 1.40
Step 3 — Eq. 20 (soil-to-stem, series per layer):
k_soil_to_stem[0] = 1/(1/1.5 + 1/0.48) = 1/(0.667+2.083)
= 1/2.750 = 0.364
k_soil_to_stem[1] = 1/(1/1.0 + 1/0.24) = 1/(1.000+4.167)
= 1/5.167 = 0.194
k_soil_to_stem[2] = 1/(1/0.5 + 1/0.08) = 1/(2.000+12.50)
= 1/14.50 = 0.069
Note: layer 2 is severely bottlenecked by both dry soil
(k_soil=0.5) AND low root density (k_RSApo=0.08).
Step 4 — Eq. 13 (total plant conductance):
Σk_RSApo = 0.48 + 0.24 + 0.08 = 0.80
k_plant = 1/(1/0.80 + 1/1.40 + 1/3.0)
= 1/(1.250 + 0.714 + 0.333)
= 1/2.298
= 0.435
The root segment (1/0.80 = 1.250) dominates the total
resistance — it's the narrowest bottleneck.Interception
compute_interception
def compute_interception(
state:SurEauPlantState, # Plant state object.
fluxes:SurEauPlantFluxes, # Plant fluxes object.
ppt:float, # Hourly precipitation (mm).
)->SurEauPlantFluxes: # Updated fluxes with ``intercepted_water`` and ``ppt_soil`` set for this timestep.
Partition hourly rainfall into canopy interception and throughfall to soil.
Implements a simple bucket model for canopy rainfall interception, described on page 5601 of Ruffault et al. (2022). The canopy is treated as a reservoir with finite capacity (set by LAI in update_LAI_and_stocks). Each hour, a fraction of rainfall proportional to the fractional canopy cover (FCC) lands on leaves; the rest falls through gaps. Once the reservoir is full, all additional rainfall passes through to the soil.
The output ppt_soil feeds directly into the precipitation term of Equation 10 (top soil layer water balance). The accumulated intercepted_water is drawn down by compute_evapo_intercepted, which evaporates it back to the atmosphere — this evaporation consumes energy that would otherwise drive transpiration, so interception indirectly reduces plant water loss.
Two cases are handled:
Case 1 — Canopy has room (ppt × FCC ≤ available): Rain is partitioned by FCC. The fraction (1 − FCC) passes through gaps to the soil; the fraction FCC is caught by leaves and added to the interception reservoir.
Case 2 — Canopy overflows (ppt × FCC > available): The canopy absorbs only what it can (available), filling the reservoir to capacity. Everything else reaches the soil. FCC partitioning no longer applies because even rain landing on saturated leaves runs off immediately.
Parameters:
state: SurEauPlantState object. The following fields are read:
canopy_storage_capacity: Maximum water the canopy can hold (mm). Computed in
update_LAI_and_stocksascanopy_storage_param × LAI. E.g. 1.8 mm for LAI = 4.5 with param = 0.4 mm per LAI unit.FCC: Fractional canopy cover (dimensionless, 0–1). From the Beer-Lambert law:
1 − exp(−K × LAI). Determines what fraction of rainfall hits the canopy vs. falls through gaps.
fluxes: SurEauPlantFluxes object. The following fields are read and written:
intercepted_water: Running total of water currently stored on leaf surfaces (mm). Read to compute remaining capacity; written to update the total after this timestep. Accumulated across hourly timesteps within a day; reset as it is evaporated by
compute_evapo_intercepted.ppt_soil: Rainfall reaching the soil surface this timestep (mm). Written. Feeds into Eq. 10 (top soil layer water balance).
ppt: Hourly precipitation (mm). The raw rainfall before canopy interception.
Returns:
- fluxes: Updated SurEauPlantFluxes with modified
intercepted_waterandppt_soil.
Logic flow:
ppt (hourly rainfall)
│
▼
available = capacity − intercepted_water
│
├─── ppt × FCC ≤ available? ──── YES ──► Case 1
│ ppt_soil = ppt × (1−FCC)
│ intercepted += ppt × FCC
│
└─── ppt × FCC > available? ──── YES ──► Case 2
ppt_soil = ppt − available
intercepted = capacity (full)Water balance check:
In both cases, total water is conserved:
Case 1: ppt_soil + Δintercepted
= ppt×(1−FCC) + ppt×FCC
= ppt
Case 2: ppt_soil + Δintercepted
= (ppt − available) + available
= ppt Numerical example (canopy_storage_capacity = 1.8 mm, FCC = 0.9):
Hour 1: ppt = 0.5 mm, intercepted = 0.0 mm
available = 1.8 − 0.0 = 1.8
ppt×FCC = 0.45 ≤ 1.8 → Case 1
ppt_soil = 0.5 × 0.1 = 0.05 mm
intercepted = 0.0 + 0.45 = 0.45 mm
Hour 2: ppt = 1.0 mm, intercepted = 0.45 mm
available = 1.8 − 0.45 = 1.35
ppt×FCC = 0.90 ≤ 1.35 → Case 1
ppt_soil = 1.0 × 0.1 = 0.10 mm
intercepted = 0.45 + 0.90 = 1.35 mm
Hour 3: ppt = 2.0 mm, intercepted = 1.35 mm
available = 1.8 − 1.35 = 0.45
ppt×FCC = 1.80 > 0.45 → Case 2 (overflow!)
ppt_soil = 2.0 − 0.45 = 1.55 mm
intercepted = 1.8 mm (full)
Hour 4: ppt = 3.0 mm, intercepted = 1.8 mm
available = 1.8 − 1.8 = 0.0
ppt×FCC = 2.70 > 0.0 → Case 2
ppt_soil = 3.0 − 0.0 = 3.0 mm (all rain passes through)
intercepted = 1.8 mm (stays full)
Water balance check across all 4 hours:
Total ppt = 0.5 + 1.0 + 2.0 + 3.0 = 6.5 mm
Total ppt_soil = 0.05 + 0.10 + 1.55 + 3.0 = 4.7 mm
Δintercepted = 1.8 − 0.0 = 1.8 mm
Sum = 4.7 + 1.8 = 6.5 mmcompute_evapo_intercepted
def compute_evapo_intercepted(
fluxes:SurEauPlantFluxes, # Plant fluxes object.
)->SurEauPlantFluxes: # Updated fluxes with ``evaporation_intercepted``, ``ETP_r``, and ``intercepted_water`` set for this timestep.
Evaporate intercepted water from the canopy using available atmospheric demand, and compute the residual ETP available for transpiration.
Intercepted water (free water sitting on leaf surfaces from rainfall) gets first priority on the atmospheric evaporative demand (ETP). Only the leftover demand (ETP_r) is available to drive transpiration through the plant (Eqs. 28–34). Two active cases are handled:
Case 1 — Demand exceeds intercepted water (ETP ≥ intercepted): All intercepted water evaporates. The canopy dries completely. The residual ETP (ETP_r = ETP − intercepted) is passed to the transpiration calculations.
Case 2 — Intercepted water exceeds demand (ETP < intercepted): All ETP is consumed evaporating the wet canopy. No energy remains for transpiration (ETP_r = 0). Some water remains on the leaves for the next timestep.
Parameters:
fluxes: SurEauPlantFluxes object. The following fields are read and/or written:
Read:
ETP: Potential evapotranspiration this timestep (mm). The total atmospheric evaporative demand. Comes from climate forcing (Penman-Monteith or prescribed).
intercepted_water: Water currently stored on leaf surfaces (mm). Accumulated by
compute_interceptionduring rain.
Written:
evaporation_intercepted: Water evaporated from leaf surfaces this timestep (mm). Always ≤ min(ETP, intercepted_water).
ETP_r: Residual potential evapotranspiration after interception evaporation (mm). This is what drives transpiration in
compute_transpirationandcalculate_Ebound_Granier. ETP_r = ETP − evaporation_intercepted.intercepted_water: Updated after evaporation. Reduced by evaporation_intercepted, or set to 0 if fully dried.
Returns:
- fluxes: Updated SurEauPlantFluxes with modified
evaporation_intercepted,ETP_r, andintercepted_water.
Logic flow:
ETP (atmospheric demand)
│
├─── ETP < 0? ──────────────── YES ──► ETP_r = 0 (safety)
│
├─── ETP ≥ intercepted? ────── YES ──► Case 1
│ evap = intercepted (all of it)
│ ETP_r = ETP − evap
│ intercepted = 0 (dry canopy)
│
└─── ETP < intercepted? ────── YES ──► Case 2
evap = ETP (all demand consumed)
ETP_r = 0 (nothing left)
intercepted −= evapConservation checks:
Energy balance (both cases):
ETP = evaporation_intercepted + ETP_r
Case 1: ETP = intercepted + (ETP − intercepted)
Case 2: ETP = ETP + 0
Water balance (both cases):
intercepted_old = evaporation_intercepted + intercepted_new
Case 1: old = old + 0
Case 2: old = ETP + (old − ETP) Numerical examples:
Case 1 — Dry canopy after evaporation:
ETP = 0.5 mm, intercepted = 0.2 mm
evaporation_intercepted = 0.2 mm (all intercepted water)
ETP_r = 0.5 − 0.2 = 0.3 mm (leftover for transpiration)
intercepted = 0.0 mm (canopy now dry)
Check: 0.2 + 0.3 = 0.5 = ETP ✓
Case 2 — Canopy still wet, transpiration suppressed:
ETP = 0.3 mm, intercepted = 1.5 mm
evaporation_intercepted = 0.3 mm (all ETP consumed)
ETP_r = 0.0 mm (no transpiration this hour!)
intercepted = 1.5 − 0.3 = 1.2 mm (still wet)
Check: 0.3 + 0.0 = 0.3 = ETP ✓
Case 1 — Dry canopy, no interception:
ETP = 0.4 mm, intercepted = 0.0 mm
evaporation_intercepted = 0.0 mm (nothing to evaporate)
ETP_r = 0.4 mm (full ETP for transpiration)
intercepted = 0.0 mm
This is the typical clear-day scenario.Transpiration
compute_transpiration
def compute_transpiration(
state:SurEauPlantState, # Plant state object.
fluxes:SurEauPlantFluxes, # Plant fluxes object.
params:SurEauVegetationParams, # Vegetation parameters.
clim:dict, # Climate snapshot.
N_hours:float, # Hours in this sub-daily timestep. Used by Granier model only.
)->SurEauPlantFluxes: # Updated fluxes with all transpiration-related fields set.
Compute all transpiration-related fluxes for the current sub-daily timestep.
This is the central transpiration calculator that ties together the leaf energy balance, cuticular conductance (Eqs. 31– 32), stomatal regulation (Eq. 34), and the implicit solver’s transpiration linearization (Eq. 60). It implements Equations 28–34 from Ruffault et al. (2022).
Two transpiration model formulations are available:
Granier A simplified approach where the maximum transpiration bound comes from scaling potential evapotranspiration (ETP) by LAI. Stomatal regulation then reduces it. The leaf energy balance takes the total transpiration as given and solves for leaf temperature.
Jarvis A mechanistic approach that builds canopy conductance explicitly from stomatal conductance (light-dependent), boundary layer conductance, and crown aerodynamic conductance, all in series (Eq. 29). The leaf energy balance is solved twice per timestep (once with previous-timestep gs to bootstrap, once with updated gs to close the loop).
Parameters:
state: SurEauPlantState object. The following fields are read:
- psi_LSym: Leaf symplasm water potential (MPa). Drives stomatal regulation (Eq. 34) and the Kelvin equation correction in the leaf energy balance.
- LAI: Current functional LAI (m²/m²). Used by the Granier model to scale ETP to the canopy.
- LAI_pheno: Phenological LAI target (m²/m²). Used by the Jarvis model to detect leafless deciduous trees in winter.
fluxes: SurEauPlantFluxes object. Read and written. The following fields are read from the previous timestep (lagged feedback):
- E_lim: Previous stomatal transpiration (mmol/m²/s). Used with E_min to compute E_inst for the Granier energy balance.
- E_min: Previous cuticular transpiration (mmol/m²/s).
- gs_lim: Previous water-limited stomatal conductance (mmol/m²/s). Used by the Jarvis energy balance.
- gmin: Previous cuticular conductance (mmol/m²/s). Used by the Jarvis energy balance.
The following fields are written:
- leaf_temperature: Leaf temperature (°C). From CPRM21 energy balance.
- g_BL: Boundary layer conductance (mmol/m²/s).
- leaf_VPD: Leaf-to-air VPD (kPa). Drives all transpiration.
- gmin: Cuticular conductance (mmol/m²/s). Eqs. 31–32.
- E_min: Leaf cuticular transpiration (mmol/m²/s).
- E_min_S: Stem cuticular transpiration (mmol/m²/s). Eq. 30.
- E_bound: Unstressed maximum transpiration (mmol/m²/s).
- E_lim: Water-limited stomatal transpiration (mmol/m²/s). Eq. 33 applied to flux.
- E_prime: dE_lim/dψ (mmol/m²/s/MPa). For Eq. 60.
- regul_fact: Stomatal regulation factor γ (0–1). Eq. 34.
- g_crown: Crown aerodynamic conductance (Jarvis only).
- gs_lim: Water-limited stomatal conductance (Jarvis only).
- gs_bound: Light-limited stomatal conductance (Jarvis only). Written by
calculate_gs_jarvis(mutates fluxes). - g_canopy_bound: Unstressed canopy conductance (Jarvis only).
- g_canopy_lim: Water-limited canopy conductance (Jarvis only).
params: SurEauVegetationParams object with the following attributes:
- transpiration_model:
"Granier"or"Jarvis". - gmin20: Cuticular conductance at 20°C (mmol/m²/s).
- TPhase_gmin: Phase transition temperature (°C). Eqs. 31–32.
- Q10_1_gmin: Q10 below T_Phase. Eq. 31.
- Q10_2_gmin: Q10 above T_Phase. Eq. 32.
- gmin_S: Stem cuticular conductance (mmol/m²/s). Eq. 30.
- f_TRB_to_leaf: Stem-to-leaf area conversion factor.
- g_crown0: Reference crown conductance (Jarvis only).
- gs_max, jarvis_PAR, etc.: Jarvis stomatal params.
- P50_gs, slope_gs: Sigmoid stomatal regulation (Eq. 34).
- transpiration_model:
clim: Climate snapshot dictionary with keys:
T_air_mean,PAR,potential_PAR,WS,RH_air_mean,VPD,ETP,ETP_veg(optional).N_hours: Number of hours in this sub-daily timestep. Used by the Granier model to convert daily ETP to sub-daily rate.
Returns:
- fluxes: Updated SurEauPlantFluxes with all transpiration- related fields set for this timestep.
Granier branch computation flow:
Step 1 E_inst = E_lim + E_min (from previous timestep)
│ Total transpiration for energy balance bootstrap
▼
Step 2 compute_T_leaf(E_inst, ψ_LSym, ...)
│ → T_leaf, g_BL, VPD_leaf
▼
Step 3 E_bound = f(ETP_veg, LAI, N_hours)
│ Maximum unstressed transpiration from ETP scaling
▼
Step 4 gmin = compute_gmin(T_leaf, ...) (Eqs. 31–32)
│ E_min = gmin × VPD / P_atm (simplified Eq. 29)
│ E_min_S = f_TRB × gmin_S × VPD / P_atm (Eq. 30)
▼
Step 5 γ, dγ/dψ = compute_regul_fact(ψ_LSym) (Eq. 34)
│ E_lim = E_bound × γ (Eq. 33)
│ E_prime = E_bound × dγ/dψ + 1e-100 (for Eq. 60)
▼
Done → feeds into solver (Eqs. 6–9, 60)Jarvis branch computation flow:
Step 1 g_crown = f(g_crown0, WS)
│ E_min_S = f_TRB × compute_E_min(gmin_S, 2000, g_crown, VPD)
│ ↑ stem transpiration computed OUTSIDE LAI guard
▼
Step 2 if LAI_pheno = 0 → zero all leaf fluxes, return
│
▼
Step 3 compute_T_leaf(gs_lim_prev, gmin_prev, ψ_LSym, ...)
│ → T_leaf, g_BL, VPD_leaf (bootstrap with old gs)
▼
Step 4 gmin = compute_gmin(T_leaf, ...) (Eqs. 31–32)
│ E_min = compute_E_min(gmin, g_BL, g_crown, VPD_leaf)
▼
Step 5 γ, dγ/dψ = compute_regul_fact(ψ_LSym) (Eq. 34)
▼
Step 6 gs_bound = calculate_gs_jarvis(PAR, T)
│ Light-limited stomatal conductance
▼
Step 7 g_canopy_bound = 1/(1/g_crown + 1/gs_bound + 1/g_BL)
│ E_bound = g_canopy_bound × VPD / P_atm (Eq. 29)
▼
Step 8 gs_lim = gs_bound × γ (Eq. 33)
│ g_canopy_lim = 1/(1/g_crown + 1/gs_lim + 1/g_BL)
│ E_lim = g_canopy_lim × VPD / P_atm
▼
Step 9 E_prime = E_lim × gs_lim_prime / denom (for Eq. 60)
│ (derivative of series conductance formula)
▼
Step 10 compute_T_leaf(gs_lim_new, gmin_new, ψ_LSym, ...)
│ → T_leaf, g_BL, VPD_leaf (close the loop with updated gs)
▼
Done → feeds into solver (Eqs. 6–9, 60)Key differences between Granier and Jarvis:
Feature Granier Jarvis
──────────────────────────────────────────────────────────
E_bound source ETP × LAI scaling Explicit g_canopy × VPD
Energy balance Takes E as input Takes gs as input
compute_T_leaf Called once Called twice (bootstrap)
E_min Simplified g×VPD/P Full series resistance
E_prime E_bound × dγ/dψ Derivative of series formula
gs_lim tracked No Yes (used for 2nd T_leaf call)
Light response No Yes (calculate_gs_jarvis)
LAI guard None LAI_pheno > 0Note on the Jarvis double-call to compute_T_leaf:
The function calls compute_T_leaf twice in the Jarvis branch. The first call (step 3) uses previous-timestep gs_lim and gmin to calculate leaf temperature. Steps 4–8 then update gmin and gs_lim using the freshly computed T_leaf and γ. The second call (step 10) recomputes T_leaf with these updated values, closing the T_leaf-gs feedback loop. One iteration (not iterating to convergence) is used, which is a reasonable approximation since the coupling is typically weak over a single sub-hourly timestep.
Water storage
compute_water_storage
def compute_water_storage(
state:SurEauPlantState, # Plant state object.
diag:SurEauPlantDiagnostics, # Diagnostics object. Written: ``LFMC``, ``LFMC_symp``, ``LFMC_apo``, ``DFMC``, ``FMC_canopy``.
params:SurEauVegetationParams, # Vegetation parameters.
VPD:float, # Vapor pressure deficit (kPa).
)->tuple: # Updated with current water volumes in all compartments.
Compute actual water content in all four plant compartments and derive fuel moisture diagnostics for wildfire risk assessment.
Converts the water potentials (ψ) tracked by the solver into actual water volumes (L/m² ground) using two mechanisms:
Symplasm (leaf and stem living cells): Water content is determined by the relative water content (RWC) from the pressure-volume curve (Eqs. 43–46 via
compute_Rs). The sponge-like tissue releases water as ψ drops.Apoplasm (leaf and stem xylem pipes): Water content is determined by cavitation (Eq. 15 via PLC). Intact pipes are full; embolized pipes have released their water. Fraction remaining = (1 − PLC/100).
The function then computes three fuel moisture diagnostics:
- LFMC (Live Fuel Moisture Content): Total water in live leaves / leaf dry mass × 100. The key wildfire metric.
- DFMC (Dead Fuel Moisture Content): Equilibrium moisture of dead leaves as a function of VPD.
- FMC_canopy: Combined moisture of live + dead leaves. The most operationally relevant metric for fire prediction.
This is the diagnostic counterpart to update_capacitances — while update_capacitances computes the derivative dRWC/dψ (how much water is released per unit ψ drop), this function computes the integral (how much water is currently present).
Parameters:
state: SurEauPlantState object. The following fields are read and/or written:
Read:
- psi_LSym: Leaf symplasm water potential (MPa).
- psi_SSym: Stem symplasm water potential (MPa).
- PLC_leaf: Percent loss of leaf hydraulic conductivity (%).
- PLC_stem: Percent loss of stem hydraulic conductivity (%).
- Q_LSym_sat_L: Leaf symplasm saturated water (L/m²). From
update_LAI_and_stocks(Eqs. 36, 38). - Q_LApo_sat_L: Leaf apoplasm saturated water (L/m²). From
update_LAI_and_stocks(Eqs. 37, 38). - Q_SSym_sat_L: Stem symplasm saturated water (L/m²). From
update_LAI_and_stocks(Eq. 39). - Q_SApo_sat_L: Stem apoplasm saturated water (L/m²). From
update_LAI_and_stocks(Eq. 40). - DM_live_canopy: Dry mass of living leaves (g/m²).
- DM_dead_canopy: Dry mass of dead leaves (g/m²).
Written:
- Q_LSym_L: Current leaf symplasm water (L/m²).
- Q_LApo_L: Current leaf apoplasm water (L/m²).
- Q_SSym_L: Current stem symplasm water (L/m²).
- Q_SApo_L: Current stem apoplasm water (L/m²).
diag: SurEauPlantDiagnostics object. The following fields are written:
- LFMC_symp: Symplasmic leaf fuel moisture content (%).
- LFMC_apo: Apoplasmic leaf fuel moisture content (%).
- LFMC: Total live fuel moisture content (%). The key wildfire metric. Typically 40–150%.
- DFMC: Dead fuel moisture content (%). From VPD equilibrium via
compute_DFMC. - FMC_canopy: Canopy fuel moisture content (%, live + dead).
params: SurEauVegetationParams object with the following attributes:
- pi_full_turgor_leaf: Osmotic potential at full turgor π₀ for leaves (MPa). For PV curve inversion.
- epsilon_sym_leaf: Bulk elastic modulus ε for leaf symplasm (MPa). For PV curve inversion.
- pi_full_turgor_stem: Osmotic potential at full turgor π₀ for stem.
- epsilon_sym_stem: Bulk elastic modulus ε for stem symplasm (MPa).
- apo_frac_leaf: Apoplasmic fraction of leaf water (dimensionless, 0–1). Used to split DM for LFMC compartment diagnostics.
VPD: Vapor pressure deficit (kPa). Drives dead fuel moisture equilibrium via
compute_DFMC.
Returns:
state: Updated SurEauPlantState with modified
Q_LSym_L,Q_LApo_L,Q_SSym_L,Q_SApo_L.diag: Updated SurEauPlantDiagnostics with modified
LFMC,LFMC_symp,LFMC_apo,DFMC,FMC_canopy.
Computation flow:
┌─────────────────────────────────────────────────────────┐
│ LEAF SYMPLASM (Eqs. 43/46 via compute_Rs) │
│ │
│ RWC = 1 − compute_Rs(π₀_leaf, ε_leaf, ψ_LSym) │
│ Q_LSym = max(0, RWC) × Q_LSym_sat │
│ LFMC_symp = 100 × Q_LSym / (DM × (1−α_apo) / 1000) │
└───────────────────────┬─────────────────────────────────┘
│
▼
┌─────────────────────────────────────────────────────────┐
│ LEAF APOPLASM (Eq. 15 via PLC) │
│ │
│ Q_LApo = (1 − PLC_leaf/100) × Q_LApo_sat │
│ LFMC_apo = 100 × Q_LApo / (DM × α_apo / 1000) │
└───────────────────────┬─────────────────────────────────┘
│
▼
┌─────────────────────────────────────────────────────────┐
│ TOTAL LFMC │
│ │
│ LFMC = 100 × (Q_LApo + Q_LSym) / (DM / 1000) │
└───────────────────────┬─────────────────────────────────┘
│
▼
┌─────────────────────────────────────────────────────────┐
│ STEM (same logic, no fuel moisture diagnostics) │
│ │
│ Q_SSym = max(0, RWC_stem) × Q_SSym_sat (Eq. 43) │
│ Q_SApo = (1 − PLC_stem/100) × Q_SApo_sat (Eq. 15) │
└───────────────────────┬─────────────────────────────────┘
│
▼
┌─────────────────────────────────────────────────────────┐
│ CANOPY FUEL MOISTURE (live + dead) │
│ │
│ DFMC = compute_DFMC(VPD) (empirical) │
│ Q_LDead = DFMC/100 × DM_dead / 1000 │
│ FMC_canopy = 100 × (Q_LApo+Q_LSym+Q_LDead) │
│ / (DM_live+DM_dead) / 1000 │
└─────────────────────────────────────────────────────────┘Two mechanisms for water loss:
SYMPLASM (elastic sponge):
RWC = f(ψ) via PV curve (Eq. 43)
Q = RWC × Q_sat
Mechanism: osmotic + turgor pressure release water
as ψ drops
Reversible: if ψ recovers, RWC recovers
APOPLASM (rigid pipes):
Q = (1 − PLC/100) × Q_sat
Mechanism: cavitation bursts pipes, releasing water
Irreversible: PLC only increases (in this model)
PLC = 0% → all pipes full → Q = Q_sat
PLC = 50% → half pipes burst → Q = 0.5 × Q_sat
PLC = 100% → all pipes burst → Q = 0