SurEau plant hydraulics

Pressure-Volume curve functions

---

source

compute_Rs

def compute_Rs(
    pi_FT, # Osmotic potential at full turgor π₀ (MPa). Negative value,
    epsilon, # Bulk modulus of elasticity ε of the symplasm (MPa). 
    psi, # Current water potential ψ of the symplasm (MPa). 
): # Relative symplasmic water deficit (dimensionless, 0 to 1).

Compute the relative symplasmic water deficit Rs = 1 − RWC.

Inverts the pressure–volume curve (Ruffault et al. 2022, Eq. 43) to obtain the fractional water loss of the symplasm at a given water potential ψ. The function handles both regimes of the p–v curve (above and below the turgor loss point)

Above the turgor loss point (Rs1), Rs is found by substituting RWC = 1 - Rs into Equation 46 (the product form of Eq. 43), yielding a quadratic in Rs:

ε·Rs² + (ψ + π₀ − ε)·Rs − ψ = 0

which is solved with the quadratic formula (minus root).

Derivation of the quadratic solved by compute_Rs()

STEP 1 — Equation 43 (first case): the pressure-volume curve
─────────────────────────────────────────────────────────────
    ψ = −π₀ - ε(1 - RWC) + π₀ / RWC

            │
            │  Multiply both sides by RWC to clear the fraction
            ▼

STEP 2 — Equation 46: product form
───────────────────────────────────
    ψ·RWC + π₀·RWC + ε(1 - RWC)·RWC - π₀ = 0

            │
            │  Substitute  RWC = 1 - Rs
            ▼

STEP 3 — After substitution
────────────────────────────
    ψ(1 - Rs) + π₀(1 - Rs) + ε·Rs·(1 - Rs) - π₀ = 0

            │
            │  Expand every product
            ▼

STEP 4 — Expanded form
───────────────────────
    (ψ - ψ·Rs) + (π₀ - π₀·Rs) + (ε·Rs - ε·Rs² - π₀) = 0

            │
            │  Cancel  +π₀  and  −π₀ , then group by powers of Rs
            ▼

STEP 5 — Collected by powers of Rs
───────────────────────────────────
    -ε·Rs²  +  (−ψ - π₀ + ε)·Rs  +  ψ  =  0

            │
            │  Multiply everything by -1
            ▼

STEP 6 — Standard quadratic form  (a·x² + b·x + c = 0)
────────────────────────────────────────────────────────
    ε·Rs²  +  (ψ + π₀ - ε)·Rs  -  ψ  =  0

    where   a = ε
            b = (ψ + π₀ - ε)
            c = −ψ

            │
            │  Apply quadratic formula, take the MINUS root
            │  (the plus root gives Rs > 1, physically impossible)
            ▼

STEP 7 — Quadratic formula (minus root)
────────────────────────────────────────

            -b  -  √(b² - 4ac)
    Rs  =  ─────────────────────
                2a

            -(ψ + π₀ - ε)  -  √( (ψ + π₀ - ε)²  +  4·ε·ψ )
        =   ──────────────────────────────────────────────────────
                                    2·ε

    Note: the  +4·ε·ψ  under the square root comes from
        -4ac = -4·ε·(−ψ) = +4·ε·ψ

            │
            │  Direct translation to Python
            ▼

STEP 8 — Python implementation
──────────────────────────────

Rs1 = (
-1 * (psi + pi_FT - epsilon)    # ← −b
- np.sqrt((psi + pi_FT - epsilon) ** 2 + 4 * psi * epsilon) # ← −√(b²−4ac)
) / (2 * epsilon)  # ← ÷ 2a

Below the turgor loss point (Rs2), turgor is zero and only the osmotic term remains (Eq. 43, second case: ψ = π₀/RWC), giving:

Rs = 1 − π₀/ψ

Physical analogy: Think of a balloon filled with sugar water. Above TLP the balloon is inflated — both the elastic wall (turgor) and the sugar concentration (osmotic pressure) resist water loss. Below TLP the balloon has deflated — only the sugar concentration matters.

Parameters:

• pi_FT: Osmotic potential at full turgor π₀ (MPa). Negative value, e.g. −2.1.

• epsilon: Bulk modulus of elasticity ε of the symplasm (MPa). Positive value, e.g. 10.

• psi: Current water potential ψ of the symplasm (MPa). Negative value, e.g. −1.5. Corresponds to psi_LSym or psi_SSym

Returns:

• Rs: Relative symplasmic water deficit (dimensionless, 0 to 1). Rs = 0 means fully hydrated (RWC = 1). Rs = 1 means completely desiccated (RWC = 0).

---

source

compute_turgor

def compute_turgor(
    pi_FT, # Osmotic potential at full turgor π₀ (MPa).  
    epsilon, # Bulk modulus of elasticity ε of the symplasm (MPa).  
    Rs, # Relative symplasmic water deficit (dimensionless, 0 to 1). 
): # Turgor pressure (MPa). 

Compute turgor pressure from the elastic component of the pressure–volume curve.

Turgor is the positive internal pressure a plant cell exerts against its cell wall when hydrated. It is extracted from Equation 43 (first case) of Ruffault et al. (2022) by decomposing the water potential into its turgor and osmotic components:

ψ  =  P  +  π

where P is turgor pressure and π is osmotic potential. Equation 43 can be split as:

ψ  =  [ −π₀ − ε(1 − RWC) ]  +  [ π₀ / RWC ]
       ___________________      ___________
            turgor  P            osmotic  π

Substituting Rs = 1 − RWC into the turgor part gives:

P = −π₀ − ε · Rs

which is what this function computes.

Physical analogy: Think of the cell as a water balloon filled with sugar water. The term −π₀ is the “natural inflation” from osmosis — even with no external pressure, the sugar solution draws water in and inflates the balloon. The term −ε·Rs is the “deflation penalty” — each unit of water lost costs ε MPa of turgor. When Rs reaches −π₀/ε, the balloon has fully deflated (turgor = 0 = the turgor loss point). Beyond that, the formula returns negative values, which are physically impossible, the calling code ( compute_turgor_from_psi ) clamps the result to zero with max(0.0, turgor).

Parameters:

• pi_FT: Osmotic potential at full turgor π₀ (MPa). Negative value, e.g. −2.1. Equivalent to pi_full_turgor_leaf or pi_full_turgor_stem in SurEauVegetationParams.

• epsilon: Bulk modulus of elasticity ε of the symplasm (MPa). Positive value, e.g. 10. Controls how fast turgor drops per unit of water lost — a stiff cell wall (high ε) loses turgor rapidly; a soft wall (low ε) loses turgor gradually. Equivalent to epsilon_sym_leaf or epsilon_sym_stem in SurEauVegetationParams.

• Rs: Relative symplasmic water deficit (dimensionless, 0 to 1). Rs = 0 means fully hydrated (RWC = 1); Rs = 1 means completely desiccated (RWC = 0). Obtained from compute_Rs().

Returns:

• P: Turgor pressure (MPa). Positive when the cell is inflated (above turgor loss point), zero or negative when the cell wall is flaccid (at or below turgor loss point). The caller should clamp to zero if a physically meaningful turgor is needed.

---

source

compute_TLP

def compute_TLP(
    pi_FT, epsilon
):

Turgor loss point [MPa].

---

source

compute_turgor_from_psi

def compute_turgor_from_psi(
    pi_FT, epsilon, psi
):

Turgor directly from water potential.

Vulnerability curve

---

source

compute_PLC

def compute_PLC(
    psi, # Water potential of the apoplasm ψ (MPa). 
    slope, # Slope of the linear rate of embolism spread at the inflection point P50. 
    P50, # Water potential causing 50 % loss of hydraulic conductance (MPa). 
):

Compute the percent loss of hydraulic conductivity from a sigmoidal vulnerability curve.

Implements Equation 15 of Ruffault et al. (2022):

                  100
PLC  =  ─────────────────────────────────
         1 + exp( slope/25 · (ψ − P50) )

The logistic (S-shaped) curve models the progressive failure of xylem conduits as water potential becomes more negative. A population of vessels has a distribution of breaking points, wide conduits cavitate first, narrow ones last, producing a cumulative sigmoid failure curve?

Physical analogy: Imagine a bundle of 100 drinking straws of varying thickness. As suction increases, the thickest straws collapse first. P50 is the suction at which half the straws are blocked. slope controls how uniform the straws are: high slope means all similar thickness (sharp transition); low slope means a wide range (gradual transition).

The PLC value feeds directly into the conductance reduction of Equations 14 and 16:

K_actual = K_max × (100 − PLC) / 100

This creates the positive feedback loop driving hydraulic failure: more cavitation → lower conductance → less water supply → more negative ψ → more cavitation.

Parameters:

• psi: Water potential of the apoplasm ψ (MPa). Negative value, e.g. −3.0. Corresponds to psi_LApo (leaf, Eq. 15) or psi_SApo (stem, Eq. 16) in SurEauPlantState.

• slope: Slope of the linear rate of embolism spread at the inflection point P50 (% MPa⁻¹). Positive value, e.g. 60. Higher values produce a steeper (more abrupt) transition. Equivalent to slope_VC_leaf or slope_VC_stem in SurEauVegetationParams.

• P50: Water potential causing 50 % loss of hydraulic conductance (MPa). Negative value, e.g. −3.4. The single most important trait for drought resistance — more negative P50 means more resistant. Equivalent to P50_VC_leaf or P50_VC_stem in SurEauVegetationParams.

Returns:

• PLC: Percent loss of conductivity (%, 0 to 100). PLC ≈ 0 when ψ >> P50 (safe). PLC = 50 when ψ = P50 (by definition). PLC ≈ 100 when ψ << P50 (hydraulic failure).

References:

• Ruffault et al. (2022), Eq. 15 (leaf) and Eq. 16 (stem, same functional form with stem-specific parameters).

Numerical example (Quercus petraea: P50 = −3.4, slope = 60):

ψ = −1.0 MPa  →  PLC =   0.3%   (essentially intact)
ψ = −2.5 MPa  →  PLC =  10.3%   (early damage)
ψ = −3.4 MPa  →  PLC =  50.0%   (P50 by definition)
ψ = −4.3 MPa  →  PLC =  89.7%   (severe damage)
ψ = −6.0 MPa  →  PLC =  99.8%   (near-total failure)

---

source

compute_PLC_prime

def compute_PLC_prime(
    PLC, # Current percent loss of conductivity (%, 0 to 100). 
    slope, # Slope of the linear rate of embolism spread at the inflection point P50  
): # Derivative dPLC/dψ (% MPa⁻¹). Always negative or zero.

Compute the derivative of the percent loss of conductivity with respect to water potential: dPLC/dψ.

Implements Equation 26 of Ruffault et al. (2022):

             slope     PLC        PLC
PLC'  =  − ─────── × ───── × (1 − ─────)
              25       100        100

This is the derivative of the sigmoidal vulnerability curve (Eq. 15) with respect to ψ. It can be verified by differentiating the logistic function directly: the derivative of 100 / (1 + exp(−x)) is (100 · exp(−x)) / (1 + exp(−x))²,

which simplifies to (PLC/100) · (1 − PLC/100) times the chain rule factor slope/25.

The denominator 250,000 in the code is the product of the three scaling constants:

25 × 100 × 100 = 250,000

where:  25   from the slope/25 convention (Pammenter & Vander Willigen 1998)
       100   from PLC/100  (first factor)
       100   from (100 − PLC)/100  (second factor, written as (100 − PLC)/100)

So:  −slope × PLC × (100 − PLC) / 250,000
   = −(slope/25) × (PLC/100) × (1 − PLC/100)

Physical meaning: PLC’ tells you how fast pipes are breaking per unit drop in water potential. It is steepest at P50 (where PLC = 50, so both factors PLC/100 and (1 − PLC/100) equal 0.5) and approaches zero at both extremes (PLC ≈ 0 or PLC ≈ 100).

PLC’ is used in the solver (sureau_solver) to compute the cavitation water-release conductance K_cav (Equation 25):

K_cav = −Q_sat_Apo × PLC' / dt

This equivalent conductance lets the water released from bursting pipes be treated as a Darcy-type flux, fitting into the same numerical framework as all other fluxes in the model.

References:

Pammenter, N. W. and Vander Willigen, C.: A mathematical and statistical analysis of the curves illustrating vulnerability of xylem to cavitation, Tree Physiology, 18, 589–593, 1998.

Parameters:

• PLC: Current percent loss of conductivity (%, 0 to 100). Obtained from compute_PLC(). Corresponds to PLC_leaf or PLC_stem in SurEauPlantState.

• slope: Slope of the linear rate of embolism spread at the inflection point P50 (% MPa⁻¹). Positive value, e.g. 60. Equivalent to slope_VC_leaf or slope_VC_stem in SurEauVegetationParams.

Returns:

• PLC_prime: Derivative dPLC/dψ (% MPa⁻¹). Always negative or zero, because PLC increases (more damage) as ψ decreases (more negative water potential).

Derivation:

Eq. 15   PLC = 100 / (1 + exp(slope/25 · (ψ − P50)))
      │  differentiate with respect to ψ using the chain rule
      │  d/dψ [100/(1+e^u)] = −100·e^u/(1+e^u)² · du/dψ
      │  where u = slope/25·(ψ−P50), so du/dψ = slope/25
      ▼

Simplify             note that e^u/(1+e^u)² = [1/(1+e^u)]·[e^u/(1+e^u)]
      │                                      = (PLC/100)·(1 − PLC/100)
      ▼

Eq. 26   PLC' = −(slope/25) · (PLC/100) · (1 − PLC/100)
      │  multiply out the denominators: 25 × 100 × 100 = 250,000
      ▼

    Code    return −slope × PLC × (100 − PLC) / 250,000

Numerical example (slope = 60):

PLC =   0%  →  PLC' = −60 ×   0 × 100 / 250000 =  0.0000  (no pipes breaking)
PLC =  10%  →  PLC' = −60 ×  10 ×  90 / 250000 = −0.0216  (slow onset)
PLC =  50%  →  PLC' = −60 ×  50 ×  50 / 250000 = −0.0600  (maximum rate)
PLC =  90%  →  PLC' = −60 ×  90 ×  10 / 250000 = −0.0216  (slowing down)
PLC = 100%  →  PLC' = −60 × 100 ×   0 / 250000 =  0.0000  (all pipes gone)

Conductance utilities

---

source

kseries

def kseries(
    k1, k2
):

Two conductances in series: 1/(1/k1 + 1/k2) = k1k2/(k1+k2).*

---

source

distribute_conductances

def distribute_conductances(
    k_plant_init:float, ri:ndarray, sym_frac_leaf:float=0.4
)->tuple:

Distribute whole-plant K among root, stem-leaf, and leaf symplasm.

---

source

compute_g_crown

def compute_g_crown(
    g_crown0:float, wind_speed:float
)->float:

Crown aerodynamic conductance.

Root distribution

---

source

compute_root_profile_BRP

def compute_root_profile_BRP(
    z:ndarray, # Cumulative depth of the bottom of each soil layer (m).  
    beta:float, # Species-specific root distribution parameter (dimensionless, 0 < β < 1). 
)->ndarray: # Root fraction in each soil layer (dimensionless). Same length as z. 

Compute per-layer root fractions using the Beta Root Profile model.

Implements Equation 19 of Ruffault et al. (2022), based on the exponential root distribution of Jackson et al. (1996):

r_j = β^{z_{j-1} · 100} − β^{z_j · 100}

which is equivalent to differencing the cumulative root fraction:

Y(d) = 1 − β^d              (Jackson et al. 1996)

r_j  = Y(z_j) − Y(z_{j-1})

where d is depth in centimeters, β is a species-specific root distribution parameter (0 < β < 1), and z_j is the cumulative depth of the bottom of layer j in meters (converted to cm internally by multiplying by 100).

The parameter β controls how deep roots extend. Think of it as a “root shallowness dial” :

• β close to 0 (e.g. 0.90): roots are very shallow, most mass concentrated near the surface.

• β close to 1 (e.g. 0.99): roots are spread deeply into the soil, like a deep-rooted tree reaching for groundwater.

The raw fractions from this function may not sum to exactly 1.0 because the exponential tail extends beyond the deepest layer.

The calling function [compute_root_profile](https://ecamo19.github.io/plant_hydraulics/sureau_plant_hydraulics.html#compute_root_profile) handles normalization.

Parameters:

• z: Cumulative depth of the bottom of each soil layer (m). NumPy array, e.g. np.array([0.2, 0.8, 2.0]). Corresponds to soil_params.depth in SurEauSoilParams.

• beta: Species-specific root distribution parameter (dimensionless, 0 < β < 1). Higher values produce deeper root profiles. Equivalent to beta_root_profile in SurEauVegetationParams. Reference values from the paper (Table 4): Quercus ilex β = 0.97, Fagus sylvatica β = 0.98.

Returns:

• rd: Root fraction in each soil layer (dimensionless array). Same length as z. Values are non-negative but may not sum to exactly 1.0 (see normalization note above).

References:

• Ruffault et al. (2022), Eq. 19.

• Jackson, R. B. et al. (1996). A global analysis of root distributions for terrestrial biomes. Oecologia, 108, 389–411.

Derivation:

Jackson et al. 1996    Y(d) = 1 − β^d     (cumulative, d in cm)
      │                evaluate at each layer bottom: 
      │                z in m → z×100 in cm
      ▼

Cumulative array    rd = 1 − β^{z·100}  = [Y(z₁), Y(z₂), ..., Y(zₙ)]
      │             difference consecutive entries (prepend 0 for 
      │             first layer)
      ▼
    Eq. 19          r_j = Y(z_j) − Y(z_{j-1})
      │             direct translation to Python
      ▼
     Code          rd - np.concatenate([[0.0], rd[:-1]])

Numerical example (β = 0.97, depths = [0.2, 0.8, 2.0] m):

• Cumulative:

            Y(20cm)  = 1 − 0.97^20  = 0.4562
  
            Y(80cm)  = 1 − 0.97^80  = 0.9124
  
            Y(200cm) = 1 − 0.97^200 = 0.9976

• Per-layer:

        r₁ = 0.4562 − 0.0    = 0.4562  (45.6% in top 20 cm)
  
        r₂ = 0.9124 − 0.4562 = 0.4562  (45.6% in 20–80 cm)
  
        r₃ = 0.9976 − 0.9124 = 0.0852  ( 8.5% in 80–200 cm)
  
        Total = 0.9976  (not exactly 1.0)

---

source

compute_root_profile_LDR

def compute_root_profile_LDR(
    z:ndarray, # Cumulative depth of the bottom of each soil layer (m).
    z50:float, # Depth above which 50 % of root density is found (m). 
    z95:float, # Depth above which 95 % of root density is found (m). 
)->ndarray: # Root fraction in each soil layer (dimensionless). Same length as z. 

Compute per-layer root fractions using the Linear Dose-Response model.

Implements the logistic root distribution of Schenk & Jackson (2002).

The cumulative root fraction at depth d is:

                1
Y(d) = ────────────────────
        1 + (d / Z50)^{pc}

where:

pc = 2.94 / ln(Z50 / Z95)

The constant 2.94 comes from the logistic curve properties: at d = Z50, Y = 0.5 (50% of roots above), and at d = Z95, Y ≈ 0.05 (95% of roots above). The shape parameter pc is derived so that these two constraints are satisfied simultaneously.

Since Z95 > Z50, the ratio Z50/Z95 < 1, so ln(Z50/Z95) < 0, and pc is negative. This makes (d/Z50)^pc a decreasing function of d, so Y(d) increases with depth — Y represents the cumulative fraction of roots that are ABOVE depth d (or equivalently, 1 − Y is the fraction below d). The differencing yz - concat([0], yz[:-1]) then yields positive per-layer fractions.

Per-layer fractions are obtained by differencing consecutive cumulative values, identical to the BRP approach:

r_j = Y(z_j) − Y(z_{j-1})

Physical meaning: Z50 and Z95 are field measurements, “half the roots are above this depth” and “almost all roots are above this depth.” This makes the LDR model convenient when root excavation data are reported as depth quantiles rather than as a fitted β parameter.

The raw fractions from this function may not sum to exactly 1.0. The calling function [compute_root_profile](https://ecamo19.github.io/plant_hydraulics/sureau_plant_hydraulics.html#compute_root_profile) handles normalization.

Parameters:

• z: Cumulative depth of the bottom of each soil layer (m). NumPy array, e.g. np.array([0.2, 0.8, 2.0]). Corresponds to soil_params.depth in SurEauSoilParams.

• z50: Depth above which 50 % of root density is found (m). Positive value, e.g. 0.3. Equivalent to root_Z50 in SurEauVegetationParams.

• z95: Depth above which 95 % of root density is found (m). Positive value, must be greater than z50, e.g. 1.2. Equivalent to root_Z95 in SurEauVegetationParams.

Returns:

• rd: Root fraction in each soil layer (dimensionless array). Same length as z. Values are non-negative but may not sum to exactly 1.0 (see normalization note above).

References:

• Schenk, H. J. and Jackson, R. B. (2002). The global biogeography of roots. Ecological Monographs, 72(3), 311–328.

Derivation:

Schenk & Jackson 2002   Y(d) = 1 / (1 + (d/Z50)^pc)
      │                 compute shape parameter from Z50, Z95
      │                 pc = 2.94 / ln(Z50 / Z95)
      │                 NOTE: pc is NEGATIVE because Z50 < Z95
      ▼
Cumulative array        yz = [Y(z₁), Y(z₂), ..., Y(zₙ)]
      │                 yz INCREASES with depth (because pc < 0)
      │                 difference consecutive entries
      │                 (prepend 0 for first layer)
      ▼
Per-layer               r_j = Y(z_j) − Y(z_{j-1})

Numerical example (Z50 = 0.3 m, Z95 = 1.2 m, depths = [0.2, 0.8, 2.0] m):

Step 1 — Shape parameter:
    pc = 2.94 / ln(0.3/1.2) = 2.94 / (−1.386) = −2.121

Step 2 — Cumulative fractions (yz INCREASES with depth):

    z = 0.2:  (0.2/0.3)^{−2.121} = (0.667)^{−2.121} = 2.404
              yz = 1/(1 + 2.404) = 0.294

    z = 0.8:  (0.8/0.3)^{−2.121} = (2.667)^{−2.121} = 0.128
              yz = 1/(1 + 0.128) = 0.887

    z = 2.0:  (2.0/0.3)^{−2.121} = (6.667)^{−2.121} = 0.018
              yz = 1/(1 + 0.018) = 0.982

Step 3 — Per-layer fractions (all positive):

    r₁ = 0.294 − 0.0   = 0.294  (29.4% in top 20 cm)
    r₂ = 0.887 − 0.294 = 0.593  (59.3% in 20–80 cm)
    r₃ = 0.982 − 0.887 = 0.095  ( 9.5% in 80–200 cm)
                           ─────
                   Total = 0.982  (not exactly 1.0)

Note: most roots (59.3%) are in the 20–80 cm layer, which
makes sense for Z50 = 0.3 m — half the roots are above 30 cm,
so the 0–20 cm layer captures less than half (29.4%), and the
20–80 cm layer captures the bulk of the remaining roots.

---

source

compute_root_profile

def compute_root_profile(
    soil_params:SurEauSoilParams, # Soil parameters object with the following attributes:
       # - depth : np.ndarray
       #   Cumulative depth of the bottom of each soil layer (m).
    veg_params:SurEauVegetationParams, # Vegetation parameters object with the following attributes:
        # - root_distribution_model : str
        #   "BRP" for Beta Root Profile (Eq. 19) or "LDR"
        # - beta_root_profile : float
        #   β parameter for BRP (dimensionless, 0 < β < 1). Only used when model is "BRP".
        # - root_Z50 : float or None
        #   Depth above which 50 % of roots are found (m). Only used when model is "LDR".
        # - root_Z95 : float or None
        #   Depth above which 95 % of roots are found (m). Only used when model is "LDR".
        # - root_depth_max : float or None
        #   Maximum rooting depth (m). If None, defaults to the depth of the deepest soil layer.
)->ndarray: # Normalized root fraction in each soil layer (dimensionless).

Compute and normalize the vertical root distribution across soil layers.

Dispatches to either the Beta Root Profile model (compute_root_profile_BRP, implementing Equation 19 of Ruffault et al. 2022) or the Linear Dose-Response model (compute_root_profile_LDR, Schenk & Jackson 2002), depending on veg_params.root_distribution_model. Then normalizes the resulting fractions so they sum to exactly 1.0 within the maximum rooting depth.

The normalization is necessary because:

1. The exponential (BRP) or logistic (LDR) tail extends mathematically beyond the deepest soil layer, leaving a small fraction unaccounted for (e.g. 0.24% for β = 0.97 at 2.0 m depth).
2. The user may specify a root_depth_max shallower than the deepest soil layer, requiring roots below that depth to be zeroed and redistributed.

The normalization strategy is: (a) find all layers within root_depth_max, (b) assign any “missing” root fraction (1.0 − sum) to the deepest valid layer, and (c) zero out all layers below root_depth_max.

Physical analogy: Think of the root system as a bundle of straws spread vertically through the soil. Each soil layer gets a certain number of straws. β (or Z50/Z95) controls how the bundle tapers with depth. The normalization step says: “All straws must be accounted for within the rooting depth. If some mathematically fell below the deepest layer, assign them to that deepest layer instead.”

The output feeds directly into:

• compute_root_length → Eqs. 17–18: per-layer root area RAI_j = RAI × r_j and root length density Lv.
• distribute_conductances → Eq. 17: per-layer root-to-stem conductance K_{Rj} = RAI_j × K_{R-SApo}.
• The Gardner-Cowan soil-to-root geometry factor B_GC (compute_soil_root_geometry), which depends on La and Lv derived from the root distribution.

Parameters:

• soil_params: SurEauSoilParams object with the following attributes:

  • depth: Cumulative depth of the bottom of each soil layer (m). NumPy array, e.g. np.array([0.2, 0.8, 2.0]).

• veg_params: SurEauVegetationParams object with the following attributes:

  • root_distribution_model: Which model to use. "BRP" (default) for Beta Root Profile (Eq. 19), or "LDR" for Linear Dose-Response (Schenk & Jackson 2002).
  • beta_root_profile: β parameter for BRP (dimensionless, 0 < β < 1). Only used when model is "BRP".
  • root_Z50: Depth above which 50 % of roots are found (m). Only used when model is "LDR".
  • root_Z95: Depth above which 95 % of roots are found (m). Only used when model is "LDR".
  • root_depth_max: Maximum rooting depth (m). If None, defaults to the depth of the deepest soil layer. Roots below this depth are zeroed and their fraction is reassigned to the deepest valid layer.

Returns:

• rd: Normalized root fraction in each soil layer (dimensionless NumPy array). Same length as soil_params.depth. Guaranteed to sum to exactly 1.0. Layers below root_depth_max have rd = 0.

References:

• Ruffault et al. (2022), Eq. 19 (BRP model).
• Jackson, R. B. et al. (1996). A global analysis of root distributions for terrestrial biomes. Oecologia, 108, 389–411.
• Schenk, H. J. and Jackson, R. B. (2002). The global biogeography of roots. Ecological Monographs, 72(3), 311–328.

Normalization logic:

1. Compute raw fractions via BRP (Eq. 19) or LDR
raw fractions may sum to < 1.0 (exponential tail)
      │         
      ▼
2. Find layers within root_depth_max
idx = layers where z ≤ zmax
      │  
      ▼
3. Compute missing fraction 
deficit = 1.0 − sum(rd[idx])
      │  
      ▼
4. Assign deficit to deepest valid layer
rd[last_valid] += deficit
      │  
      ▼
5. Zero out layers below root_depth_max 
rd[below_zmax] = 0.0
      │  
      ▼
6. Return rd (sums to exactly 1.0)

Numerical example (BRP, β = 0.97, depths = [0.2, 0.8, 2.0] m, root_depth_max = 2.0 m):

Raw BRP output:    [0.4562, 0.4562, 0.0852]   sum = 0.9976
Deficit:           1.0 − 0.9976 = 0.0024
Add to last layer: [0.4562, 0.4562, 0.0876]   sum = 1.0000 

Numerical example (BRP, β = 0.97, depths = [0.2, 0.8, 2.0] m, root_depth_max = 0.8 m):

Raw BRP output:    [0.4562, 0.4562, 0.0852]
Valid layers:      idx = [0, 1] (z ≤ 0.8)
Deficit:           1.0 − (0.4562 + 0.4562) = 0.0876
Add to last valid: [0.4562, 0.5438, 0.0000]  sum = 1.0000 ✓
Layer 3 zeroed:    roots below 0.8 m reassigned to layer 2

---

source

compute_root_length

def compute_root_length(
    soil_params:SurEauSoilParams, # Soil parameters object with the following attributes:
# - layer_thickness : np.ndarray
#    Thickness of each soil layer (m).    
# - RFC : np.ndarray
#    Rock fragment content of each soil layer (%). 
    veg_params:SurEauVegetationParams, # Vegetation parameters object with the following attributes:
# - LAI_max : float
#    Maximum leaf area index of the stand (m²_leaf / m²_soil).
# - f_root_to_leaf : float
#    Root-to-leaf area ratio RaLa (dimensionless). 
# - root_distribution : np.ndarray
#    Normalized root fraction in each soil layer (dimensionless, sums to 1.0). 
# - root_radius : float
#    Fine root radius r (m). E.g. 0.0002.
)->tuple: # Root length per unit soil area for each layer

Compute root length per soil area (La) and root length density per soil volume (Lv) for each soil layer.

Implements Equation 18 of Ruffault et al. (2022) plus two geometric conversions described in the text surrounding Eq. 21:

1. Distribute total root area index (RAI) across layers using the root fractions r_j from compute_root_profile (Eq. 18):

   RAI_j = RAI × r_j

2. Convert root surface area to root length using cylinder geometry (A = 2πrL, so L = A / 2πr):

   La_j = RAI_j / (2π r)

3. Convert root length per soil area to root length per soil volume by dividing by the effective layer thickness (corrected for rock fragments, same correction as Eq. 23):

   Lv_j = La_j / [ thickness_j × (1 − RFC_j/100) ]

Both outputs feed into Equation 21, the soil-to-root conductance via compute_soil_root_geometry, which computes the Gardner-Cowan geometry factor B_GC:

b    = 1 / √(π Lv_j)           (mean half-distance between roots)
B_GC = 2π La_j / ln(b / r)     (geometry factor in Eq. 21)

La_j represents the total length of roots, obtained by dividing surface area by the circumference of one straw (2πr).

Lv_j represents how densely packed the roots are per unit room volume, accounting for rocks taking upspace. From Lv, you can calculate how far apart the straws are on average (1/√(πLv)), which determines how far water must travel through the soil to reach the nearest straw.

Parameters:

• soil_params: SurEauSoilParams object with the following attributes:

  • layer_thickness: Thickness of each soil layer (m). NumPy array, e.g. np.array([0.2, 0.6, 1.2]). Computed by sureau_soil_params() from soil_params.depth.

  • RFC: Rock fragment content of each soil layer (%). NumPy array, e.g. np.array([75.0, 75.0, 75.0]). Reduces the effective soil volume available for roots.

• veg_params: SurEauVegetationParams object with the following attributes:

  • LAI_max: Maximum leaf area index of the stand (m²_leaf / m²_soil). E.g. 4.5.

  • f_root_to_leaf: Root-to-leaf area ratio RaLa (dimensionless). E.g. 1.0 means root surface area equals leaf area. Equivalent to RaLa in Eq. 18 text.

  • root_distribution: Normalized root fraction in each soil layer (dimensionless array summing to 1.0). Output of compute_root_profile(). Equivalent to r_j in Eq. 18.

  • root_radius: Fine root radius r (m). E.g. 0.0002 (= 0.2 mm). Equivalent to dR in Table 1.

Returns:

• La: Root length per unit soil area for each layer (m_root / m²_soil). NumPy array, same length as soil_params.depth. Stored as veg_params.La.

• Lv: Root length density per unit soil volume for each layer (m_root / m³_soil). NumPy array, same length as soil_params.depth. Stored as veg_params.Lv.

References:

• Ruffault et al. (2022), Eq. 18 (per-layer RAI), Eq. 21 (where La and Lv are used), Eq. 23 (same rock fragment correction).
• Jackson, R. B. et al. (1996) for the root distribution r_j.

Derivation:

 Eq. 18 RAI = LAI × RaLa        (total root area index)
      │                     distribute across layers using r_j 
      │                                 (Eq. 18)
      ▼
    Eq. 18               RAI_j = RAI × r_j    (per-layer root area)
      │                   cylinder geometry: A = 2πrL → L = A/(2πr)
      ▼
      La        La_j = RAI_j / (2π r)   (length per soil area)
      │         divide by effective layer volume
      │         V_eff = thickness × (1 − RFC/100)   (same as Eq. 23)
      ▼
      Lv         Lv_j = La_j / V_eff     (density per soil volume)
      │          feeds into Eq. 21 via compute_soil_root_geometry
      ▼
   Eq. 21               b = 1/√(πLv)       (half-distance between roots)
                     B_GC = 2πLa / ln(b/r)   (Gardner-Cowan factor)

Numerical example (LAI_max = 4.5, f_root_to_leaf = 1.0, root_radius = 0.0002 m, root_distribution = [0.4562, 0.4562, 0.0876], layer_thickness = [0.2, 0.6, 1.2] m, RFC = [75, 75, 75] %):

Step 1 — Total RAI:
    RAI = 4.5 × 1.0 = 4.5 m²_root / m²_soil

Step 2 — Per-layer RAI (Eq. 18):
    RAI₁ = 4.5 × 0.4562 = 2.053 m²/m²
    RAI₂ = 4.5 × 0.4562 = 2.053 m²/m²
    RAI₃ = 4.5 × 0.0876 = 0.394 m²/m²

Step 3 — Root length per soil area (La):
    La₁ = 2.053 / (2π × 0.0002) = 2053 / 0.001257 = 1634 m/m²
    La₂ = 2.053 / 0.001257                         = 1634 m/m²
    La₃ = 0.394 / 0.001257                         =  313 m/m²

Step 4 — Root length density per soil volume (Lv):
    V_eff₁ = 0.2 × (1 − 0.75) = 0.05 m
    Lv₁ = 1634 / 0.05 = 32,680 m/m³   (32.7 km of root per m³)

    V_eff₂ = 0.6 × 0.25 = 0.15 m
    Lv₂ = 1634 / 0.15 = 10,893 m/m³   (10.9 km per m³)

    V_eff₃ = 1.2 × 0.25 = 0.30 m
    Lv₃ = 313 / 0.30  =  1,043 m/m³   ( 1.0 km per m³)

These feed into compute_soil_root_geometry:
    b₁ = 1/√(π × 32680) = 0.0031 m  (3.1 mm between roots)
    b₂ = 1/√(π × 10893) = 0.0054 m  (5.4 mm between roots)
    b₃ = 1/√(π ×  1043) = 0.0175 m  (17.5 mm between roots)

Stomatal regulation

---

source

compute_regul_fact

def compute_regul_fact(
    psi, # Water potential of the leaf symplasm ψ_LSym (MPa). 
    params:SurEauVegetationParams, # Vegetation parameters object. Which attributes are read depends on the formulation:
# - stomatal_reg_formulation : str
#    ``"Sigmoid"`` (default), ``"PiecewiseLinear"``, or ``"Turgor"``.
# Sigmoid:
# - P50_gs : float
#    Water potential at 50 % stomatal closure (MPa). Derived
#    from ``(P12_gs + P88_gs) / 2``.
# - slope_gs : float
#    Slope of the stomatal sigmoid at P50_gs (% MPa⁻¹).
#    Derived from ``100 / (P12_gs − P88_gs)``.
# PiecewiseLinear:
# - psi_start_closing : float
#    ψ above which stomata are fully open (MPa).
# - psi_close : float
#    ψ below which stomata are fully closed (MPa).
# Turgor:
# - pi_full_turgor_leaf : float
#    Osmotic potential at full turgor π₀ (MPa).
# - epsilon_sym_leaf : float
#    Bulk elastic modulus ε (MPa).
# - turgor_pressure_at_gs_max : float
#    Turgor at which stomata are fully open (MPa).
): # Stomatal regulation factor γ (dimensionless, 0 to 1).

Compute the stomatal regulation factor γ and its derivative dγ/dψ.

The regulation factor γ is a dimensionless number between 0 and 1 that controls how open the stomata are, linking water potential to stomatal conductance via Equation 33 of Ruffault et al. (2022):

g_stom = γ × g_stom,max

where g_stom,max is the maximum stomatal conductance (set by light, temperature, CO₂). When γ = 1, stomata are fully open; when γ = 0, stomata are fully closed.

The derivative dγ/dψ is needed by the implicit solver for the transpiration linearization in Equation 60, which accounts for the quick adjustment of stomatal regulation to water potential changes within a single timestep.

Three formulations are available, selected by params.stomatal_reg_formulation:

Sigmoid (default, used in the paper) Implements Equation 34: γ follows a sigmoid (S-curve) that transitions smoothly from 1 (open) to 0 (closed) as ψ_LSym decreases. The shape is controlled by P50_gs (ψ at 50% closure) and slope_gs (steepness of the transition). This is the same logistic functional form as the xylem vulnerability curve (Eq. 15), but applied to stomata using leaf symplasm ψ rather than apoplasm ψ, and scaled to 0–1 instead of 0–100.

PiecewiseLinear (not in the paper) A simple three-zone model: fully open above ψ_start_closing, linearly decreasing between ψ_start_closing and ψ_close, fully closed below ψ_close.

Turgor (not in the paper) A mechanistic model that ties stomatal opening directly to guard cell turgor pressure, computed from the pressure–volume curve (Eqs. 43/46 via compute_Rs() and compute_turgor() ). Stomata open proportionally to turgor / turgor_at_gs_max.

Physical analogy: Think of γ as a valve controlled by water pressure. When the pipe pressure (ψ) is high (less negative), the valve is wide open. As pressure drops, the valve progressively shuts. The Sigmoid formulation describes a valve that shuts rapidly over a narrow pressure range (like a snap-action valve). The PiecewiseLinear formulation describes a valve that shuts at a constant rate. The Turgor formulation describes a valve directly driven by the inflation of the guard cells.

Parameters:

• psi: Water potential of the leaf symplasm ψ_LSym (MPa). Negative value, e.g. −2.0. Corresponds to psi_LSym in SurEauPlantState.

• params: SurEauVegetationParams object with the following attributes (which are read depends on the formulation):

  • stomatal_reg_formulation: Which model to use. "Sigmoid" (default, Eq. 34), "PiecewiseLinear", or "Turgor".

  • Sigmoid parameters:

    • P50_gs: Water potential at 50 % stomatal closure (MPa). Derived from (P12_gs + P88_gs) / 2 in sureau_vegetation_params .

    • slope_gs: Slope of the stomatal sigmoid at P50_gs (% MPa⁻¹). Derived from 100 / (P12_gs − P88_gs) in sureau_vegetation_params .

  • PiecewiseLinear parameters:

    • psi_start_closing: ψ above which stomata are fully open (MPa). E.g. −0.5.
    • psi_close: ψ below which stomata are fully closed (MPa). E.g. −2.0.

  • Turgor parameters:

    • pi_full_turgor_leaf: Osmotic potential at full turgor π₀ (MPa). E.g. −2.1.
    • epsilon_sym_leaf: Bulk elastic modulus ε (MPa). E.g. 10.
    • turgor_pressure_at_gs_max: Turgor at which stomata are fully open (MPa). E.g. 2.0.

Returns:

• rf: Stomatal regulation factor γ (dimensionless, 0 to 1). γ = 1 means stomata fully open (no water stress). γ = 0 means stomata fully closed (maximum water stress).

• rfp: Derivative dγ/dψ (MPa⁻¹). Always non-negative (γ increases as ψ increases). Used in the implicit solver for the transpiration linearization (Eq. 60). Returns 0.0 for the Turgor formulation (no analytical derivative provided).

References:

• Ruffault et al. (2022), Eq. 33 (stomatal conductance), Eq. 34 (sigmoid regulation factor), Eq. 60 (transpiration linearization using dγ/dψ).

• Pammenter, N. W. and Vander Willigen, C.: A mathematical and statistical analysis of the curves illustrating vulnerability of xylem to cavitation, Tree Physiology, 18, 589–593, 1998.

Derivation of the Sigmoid formulation (Eq. 34):

    Eq. 34      γ = 1 − 1 / (1 + exp(α(ψ − P50_gs)))
      │         define PL_gs = 1/(1+exp(α(ψ−P50)))  where α = slope_gs/25
      ▼
 Code line 1   PL_gs = 1/(1 + exp(α(ψ − P50_gs)))
      │             flip: γ = 1 − PL_gs
      ▼
Code line 2           rf = 1 − PL_gs

Derivation of the Sigmoid derivative:

  Start from            γ = 1 − PL_gs
      │             dγ/dψ = −d(PL_gs)/dψ
      │             d(PL_gs)/dψ = −α × PL_gs × (1−PL_gs)
      │             and (1−PL_gs) = γ = rf
      ▼
    Result                dγ/dψ = α × PL_gs × rf
      │                direct translation to Python
      ▼
    Code               al = slope_gs / 25.0
                      rfp = al * PL_gs * rf

Numerical example, Sigmoid (Quercus petraea: P12_gs = −2.07, P88_gs = −2.62, giving P50_gs = −2.345, slope_gs = 182):

ψ = −1.5 MPa (mild stress):
    α = 182/25 = 7.28
    PL_gs = 1/(1+exp(7.28×0.845)) = 0.002
    γ = 0.998          ← stomata 99.8% open
    dγ/dψ = 0.015

ψ = −2.345 MPa (at P50_gs):
    PL_gs = 0.5
    γ = 0.500          ← stomata half closed
    dγ/dψ = 1.82      ← maximum sensitivity

ψ = −2.8 MPa (severe stress):
    PL_gs = 0.965
    γ = 0.035          ← stomata 96.5% closed
    dγ/dψ = 0.246

PiecewiseLinear formulation:

Zone 1 (ψ > ψ_start):    γ = 1.0,  dγ/dψ = 0
      │     ψ decreasing
      ▼
Zone 2 (ψ_close < ψ < ψ_start):
                          γ = (ψ − ψ_close) / (ψ_start − ψ_close)
                          dγ/dψ = 1 / (ψ_start − ψ_close)
      │  ψ decreasing
      ▼
Zone 3 (ψ < ψ_close):    γ = 0.0,  dγ/dψ = 0

Numerical example, PiecewiseLinear (ψ_start = −0.5, ψ_close = −2.0):

ψ = −0.2 MPa  →  γ = 1.000,  dγ/dψ = 0.000  (Zone 1: fully open)
ψ = −0.5 MPa  →  γ = 1.000,  dγ/dψ = 0.000  (boundary)
ψ = −1.0 MPa  →  γ = 0.667,  dγ/dψ = 0.667  (Zone 2: linear ramp)
ψ = −1.5 MPa  →  γ = 0.333,  dγ/dψ = 0.667  (Zone 2: linear ramp)
ψ = −2.0 MPa  →  γ = 0.000,  dγ/dψ = 0.000  (boundary)
ψ = −2.5 MPa  →  γ = 0.000,  dγ/dψ = 0.000  (Zone 3: fully closed)

Turgor formulation:

Eq. 43 (via compute_Rs)      Rs = f(ψ, π₀, ε)
      │
      ▼
Turgor (via compute_turgor)   P = −π₀ − ε·Rs
      │
      ▼
Regulation factor             γ = clamp(P / P_at_gs_max, 0, 1)
                              dγ/dψ = 0   (no analytical derivative)

Numerical example, Turgor (π₀ = −2.1, ε = 10, turgor_at_gs_max = 2.0):

ψ = −0.5 MPa  →  Rs = 0.026, P = 1.84  →  γ = 0.920  (mostly open)
ψ = −1.5 MPa  →  Rs = 0.121, P = 0.89  →  γ = 0.445  (half closed)
ψ = −2.3 MPa  →  Rs = 0.193, P = 0.17  →  γ = 0.085  (nearly closed)
ψ = −2.7 MPa  →  Rs = 0.210, P = 0.00  →  γ = 0.000  (turgor loss)

---

source

calculate_gs_jarvis

def calculate_gs_jarvis(
    fluxes:SurEauPlantFluxes, params:SurEauVegetationParams, PAR:float
)->SurEauPlantFluxes:

Jarvis-type stomatal conductance with temperature correction.

Granier Ebound

---

source

calculate_Ebound_mm_Granier

def calculate_Ebound_mm_Granier(
    ETP, LAI, a:float=-0.006, b:float=0.134, c:float=0.0
):

Granier-type energy-limited transpiration [mm].

---

source

calculate_Ebound_Granier

def calculate_Ebound_Granier(
    ETP:float, LAI:float, time_step:float, params:SurEauVegetationParams
)->float:

Ebound in mmol/m²leaf/s from Granier formulation.

Cuticular conductance

---

source

compute_gmin

def compute_gmin(
    fluxes:SurEauPlantFluxes, veg_params:SurEauVegetationParams
)->float:

gmin as function of leaf temperature (Cochard et al. 2019).

---

source

compute_E_min

def compute_E_min(
    gmin:float, g_BL:float, g_crown:float, VPD:float, P_atm:float=101.3
)->float:

Minimal transpiration from cuticular/boundary layer/crown conductances.

Leaf energy balance

---

source

compute_T_leaf

def compute_T_leaf(
    state:SurEauPlantState, fluxes:SurEauPlantFluxes, veg_params:SurEauVegetationParams, clim:dict
)->tuple: # Leaf temperature (°C).

Compute leaf temperature, boundary-layer conductance, and leaf VPD by solving a linearized leaf energy balance.

Implements the leaf energy balance described in SurEau C version (Cochard et al. 2021):

The approach is a linearized single-step Penman-Monteith solution: the leaf temperature is found analytically as a departure from air temperature, rather than iteratively (as in the Newton-Raphson approach of leaf_temperature.py in the Bonan framework). Both solve the same physical problem, finding the temperature at which energy in (absorbed radiation) equals energy out (longwave emission + sensible heat + latent heat).

The core equation for the leaf-to-air temperature difference is:

         γ* × Rni × r_blr / (ρ·cp)  −  VPD
δT  =  ──────────────────────────────────────
                    Δ  +  γ*

where:

• Rni is the isothermal net radiation (W/m²), the radiation the leaf would absorb if it were at air temperature.
• r_blr is the combined boundary-layer + radiative resistance (s/m), representing the non-evaporative heat loss pathway.
• γ* = γ × (r_st / r_blr) is the modified psychrometric constant, weighting the stomatal resistance against the convective + radiative resistance. When stomata are open (r_st small), γ* is small and transpirational cooling dominates (δT < 0). When stomata are closed (r_st huge), γ* is large and the leaf heats up (δT > 0).
• Δ is the slope of the saturation vapor pressure curve (kPa/°C).
• VPD is the air vapor pressure deficit (kPa).

The function also computes the leaf VPD, the actual driving force for transpiration, which accounts for both the elevated leaf temperature (raising e_sat) and the reduced water activity at low water potentials (the Kelvin equation correction).This correction is critical for accurately computing residual transpiration under extreme drought.

Physical analogy: The leaf is a wet surface sitting in a breeze. Radiation heats it up; transpiration cools it down. The energy balance finds the temperature at which these effects balance. A leaf with open stomata (low r_st) is like a wet towel in the wind — much cooler than the air. A leaf with closed stomata (high r_st) is like a dry rock — warmer than the air because it can only cool by convection and longwave emission.

Parameters:

• T_air: Air temperature (°C).

• PAR: Photosynthetically active radiation (µmol/m²/s). E.g. 500. Used to estimate shortwave radiation via SWR = PAR × 0.5495 and cloud cover via PAR / potential_PAR.

• potential_PAR_val: Clear-sky potential PAR (µmol/m²/s). Used to estimate cloud cover for atmospheric longwave radiation. If 0 (nighttime), cloud cover is set to 0.

• WS: Wind speed (m/s). Floored at 0.1 m/s internally. Faster wind thins the boundary layer, increasing convective cooling.

• RH: Relative humidity (%). E.g. 65. Used to compute actual vapor pressure and VPD.

• gs: Stomatal conductance (mmol/m²/s). Only used when transpiration_model="Jarvis". Default None (treated as 0). Combined with g_cuti in parallel to give the total leaf conductance to water vapor, as in Eq. 29.

• g_cuti: Cuticular conductance (mmol/m²/s). Only used when transpiration_model="Jarvis". Default None (treated as 0). The residual conductance when stomata are fully closed.

• E_inst: Instantaneous transpiration rate (mmol/m²/s). Only used when transpiration_model="Granier". Default None. The effective conductance is back-calculated as g = E / (VPD/P_atm).

• psi_leaf: Leaf water potential (MPa). Default 0.0. Used in the Kelvin equation correction for vapor pressure inside the leaf. Only significant under severe drought (ψ < −5 MPa).

• leaf_size: Characteristic leaf dimension (mm). Default 50. Leaves > 3 mm use the flat-plate boundary layer formula; smaller leaves (needles) use the cylinder formula. Larger leaves have thicker boundary layers.

• leaf_angle: Leaf angle from horizontal (degrees). Default 45. Affects how much direct shortwave radiation the leaf intercepts: SWRabs = a_SWR × cos(angle) × SWR.

• turn_off_EB: If True, skip the energy balance and return T_air as the leaf temperature. Default False. Diagnostic option for testing without energy balance feedback.

• transpiration_model: "Jarvis" (default) computes stomatal resistance from gs + g_cuti conductances. "Granier" back-calculates effective conductance from E_inst.

Returns:

• T_leaf: Leaf temperature (°C).

• g_bl: Boundary layer conductance (mmol/m²/s). Used downstream in the conductance chain for computing E_min (Eqs. 29–30) and canopy-level conductances. Stored as fluxes.g_BL.

• VPD_leaf: Vapor pressure deficit at the leaf surface (kPa). Accounts for leaf temperature (via e_sat at T_leaf) and water activity correction (via Kelvin equation at ψ_leaf). Drives all transpiration calculations downstream.

References:

• Cochard, H. et al. (2021) — CPRM21 — SurEau: a mechanistic model of plant water relations under extreme drought. Ann. For. Sci., 78, 1–23. (Energy balance formulation.)

• Jones, H. G. (2013). Plants and microclimate. Cambridge University Press. (Boundary layer conductance, Eqs. 29–30 text.)

• Brutsaert, W. (1975). On a derivable formula for longwave radiation from clear skies. Water Resources Research, 11(5), 742–744. (Atmospheric emissivity.)

• Buck, A. L. (1981). New equations for computing vapor pressure and enhancement factor. J. Appl. Meteorol., 20, 1527–1532. (Saturation vapor pressure.)

Computation stages:

Stage 1 — Atmospheric properties
─────────────────────────────────
Tetens equation       e_sat = a·exp(b·T/(T+z))        (kPa)
Actual VP             ea = e_sat × RH/100              (kPa)
Slope                 Δ = e_sat·b·z / (T+z)²          (kPa/°C)
Cloud cover           cc = PAR / potential_PAR          (0–1)
                      NOTE:cc=1 means CLEAR sky, cc=0 means overcast,
                      opposite to meteorological convention.
Atm. emissivity       ε_air = f(ea, T, cc)             (Brutsaert 1975)
VPD                   VPDx = e_sat − ea                (kPa)

Stage 2 — Radiation balance
───────────────────────────
Absorbed SW           SWRabs = 0.5 × cos(angle) × SWR  (W/m²)
Incoming LW           LWRin = ε_air × σ × T⁴           (W/m²)
Outgoing LW           LWRout = ε_leaf × σ × T⁴         (W/m², at T_air)
Isothermal net rad    Rni = SWRabs + LWRin − LWRout     (W/m²)
Radiative resistance  r_r = ρ·cp / (4·ε·σ·T³)          (s/m)

Stage 3 — Boundary layer
────────────────────────
BL resistance         r_bl = f(WS, leaf_size)           (Jones 2013)
BL conductance        g_bl = (1/r_bl) × 1000 × 40      (mmol/m²/s)
Combined non-evap     r_blr = r_bl ∥ r_r                (parallel)

Stage 4 — Stomatal resistance
─────────────────────────────
Jarvis:  r_st = 1/(gs + g_cuti) × 1000 × 40            (s/m)
Granier: r_st = 1/(E/VPD × P_atm) × 1000 × 40         (s/m)

Stage 5 — Penman-Monteith temperature solution
───────────────────────────────────────────────
Modified γ            γ* = γ × (r_st / r_blr)
Leaf−air ΔT           δT = (γ*·Rni·r_blr/(ρ·cp) − VPD) / (Δ + γ*)
Leaf temperature      T_leaf = T_air + δT               (°C)

Stage 6 — Leaf VPD with Kelvin correction
──────────────────────────────────────────
e_sat at T_air        Buck (1981) equation               (Pa)
e_air                 e_sat_air × RH/100                 (Pa)
e_sat at T_leaf       Buck (1981) at T_leaf              (Pa)
Kelvin correction     e_leaf = e_sat_leaf × exp(ψ·Vw/(R·T))  (Pa)
Leaf VPD              VPD_leaf = (e_leaf − e_air) / 1000 (kPa)

Kelvin equation detail:

The vapor pressure inside the leaf is not simply e_sat(T_leaf) because dissolved solutes and the negative water potential reduce the water activity. The correction is:

e_leaf = e_sat(T_leaf) × exp(ψ × Vw / (R × T))

where Vw = 18.02 × 10⁻⁶ m³/mol (molar volume of water), R = 8.314 J/mol/K, T in Kelvin, and ψ in MPa. The pre-computed constant 2.16947115 ≈ Vw × 10⁶ / R handles the MPa-to-Pa conversion.

At ψ = −2 MPa, T = 25°C: correction = exp(−2 × 2.169/298) = 0.986 (1.4% reduction — small).

At ψ = −10 MPa, T = 25°C: correction = exp(−10 × 2.169/298) = 0.930 (7% reduction — significant for residual transpiration under extreme drought).

Numerical example (T_air = 25°C, PAR = 500, WS = 2 m/s, RH = 50%, gs = 100 mmol/m²/s, g_cuti = 4 mmol/m²/s):

Stage 1: e_sat = 3.17 kPa, ea = 1.58 kPa, VPD = 1.59 kPa
         Δ = 0.189 kPa/°C

Stage 2: SWR ≈ 275 W/m², SWRabs ≈ 97 W/m²
         LWRin ≈ 370 W/m², LWRout ≈ 382 W/m²
         Rni ≈ 85 W/m²

Stage 3: r_bl ≈ 90 s/m, r_r ≈ 220 s/m
         r_blr ≈ 64 s/m (parallel)

Stage 4: r_st = 40000 / 104 ≈ 385 s/m
         (gs+g_cuti = 104 mmol/m²/s → 0.0026 m/s → r = 385 s/m)

Stage 5: γ* = 0.066 × (385/64) = 0.397
         δT = (0.397×85×64/1306 − 1.59) / (0.189+0.397)
            = (1.66 − 1.59) / 0.586
            = 0.12°C
         T_leaf = 25.12°C (barely above air — transpiration
         cools effectively when stomata are open)

With stomata closed (gs = 0, g_cuti = 4 only):
         r_st ≈ 10,000 s/m
         γ* = 0.066 × (10000/64) = 10.3
         δT = (10.3×85×64/1306 − 1.59) / (0.189+10.3)
            = (42.9 − 1.59) / 10.49
            = 3.94°C
         T_leaf = 28.94°C (nearly 4°C above air —
         much less evaporative cooling)

 INPUTS
 ──────
   T_air              [°C]           — air temperature
   PAR                [µmol/m²/s]    — photosynthetically active radiation
   potential_PAR_val  [µmol/m²/s]    — clear-sky potential PAR
   WS                 [m/s]          — wind speed (floored at 0.1)
   RH                 [%]            — relative humidity
   gs                 [mmol/m²/s]    — stomatal conductance (Jarvis only)
   g_cuti             [mmol/m²/s]    — cuticular conductance (Jarvis only)
   E_inst             [mmol/m²/s]    — instantaneous transpiration (Granier only)
   psi_leaf           [MPa]          — leaf water potential (for Kelvin correction)
   leaf_size          [mm]           — characteristic leaf dimension
   leaf_angle         [degrees]      — leaf angle from horizontal
 
                         │
                         ▼
┌─────────────────────────────────────────────────────────┐
│                                                         │
│    STAGE 1: ATMOSPHERIC PROPERTIES                      │
│    "What are the thermodynamic conditions of the air?"  │
│                                                         │
│  Tetens equation (kPa):                                 │
│  e_sat = 0.61121 × exp(17.502 × T / (T + 240.97))       │
│                                                         │
│  Actual vapor pressure:                                 │
│  ea = e_sat × RH / 100                                  │
│                                                         │
│  Slope of saturation curve (kPa/°C):                    │
│  Δ = e_sat × 17.502 × 240.97 / (T + 240.97)²            │
│                                                         │
│  Clearness index (NOT cloud cover — see note):          │
│  cc = PAR / potential_PAR       clipped to [0, 1]       │
│  NOTE: cc=1 means CLEAR sky, cc=0 means overcast.       │
│        Opposite to meteorological convention!           │
│                                                         │
│  Atmospheric emissivity (Brutsaert 1975):               │
│  ε_air = (1−0.84·cc) × 1.31 × (10·ea/(T+273))^(1/7)     │
│        + 0.84 × cc                                      │
│                                                         │
│  Vapor pressure deficit:                                │
│  VPDx = e_sat − ea                          [kPa]       │
│                                                         │
│  T=25°C, RH=50%:  e_sat=3.17, ea=1.58, VPD=1.59 kPa     │
│                    Δ = 0.189 kPa/°C                     │
│                                                         │
└────────────────────────┬────────────────────────────────┘
                         │
                         ▼
┌─────────────────────────────────────────────────────────┐
│                                                         │
│    STAGE 2: RADIATION BALANCE                           │
│    "How much energy is the leaf absorbing?"             │
│                                                         │
│  Convert PAR → shortwave:                               │
│  SWR = PAR × 0.5495                          [W/m²]     │
│                                                         │
│  Absorbed shortwave (leaf absorptance × projection):    │
│  SWRabs = 0.5 × cos(leaf_angle) × SWR        [W/m²]     │
│                                                         │
│    0.5  = leaf absorptance (50% absorbed,               │
│           rest reflected + transmitted)                 │
│    cos() = projection onto leaf surface                 │
│    angle=45° → cos=0.707 → SWRabs = 0.354 × SWR         │
│                                                         │
│  Incoming longwave (Stefan-Boltzmann):                  │
│  LWRin = ε_air × σ × (T_air + 273.15)⁴      [W/m²]      │
│                                                         │
│  Outgoing longwave (linearized at T_air):               │
│  LWRout = ε_leaf × σ × (T_air + 273.15)⁴    [W/m²]      │
│  NOTE: uses T_air, not T_leaf — this is the             │
│  linearization. The δT correction comes in Stage 5.     │
│                                                         │
│  Isothermal net radiation:                              │
│  Rni = SWRabs + LWRin − LWRout               [W/m²]     │
│                                                         │
│  Radiative resistance (linearized LW derivative):       │
│  r_r = ρ·cp / (4 × ε_leaf × σ × T³)         [s/m]       │
│                                                         │
│  Constants: σ = 5.6704e−8 W/m²/K⁴                       │
│             ε_leaf = 0.97                               │
│             ρ = 1.292 kg/m³, cp = 1010 J/kg/K           │
│                                                         │
│  Example: SWR≈275, SWRabs≈97, LWRin≈370, LWRout≈382     │
│           Rni ≈ 85 W/m²,  r_r ≈ 220 s/m                 │
│                                                         │
└────────────────────────┬────────────────────────────────┘
                         │
                         ▼
┌─────────────────────────────────────────────────────────┐
│                                                         │
│    STAGE 3: BOUNDARY LAYER RESISTANCE                   │
│    "How thick is the still-air layer around the leaf?"  │
│                                                         │
│  Two regimes (Jones 2013):                              │
│                                                         │
│  Broad leaves (leaf_size > 3 mm):                       │
│  ┌────────────────────────────────────────────────┐     │
│  │  gflat = 0.00662,  jflat = 0.5                 │     │
│  │  r_bl = 1 / (1.5 × gflat × WS^0.5              │     │
│  │              / (leaf_size/1000)^0.5 )          │     │
│  │                                                │     │
│  │  Bigger leaf → thicker BL → higher resistance  │     │
│  │  Faster wind → thinner BL → lower resistance   │     │
│  └────────────────────────────────────────────────┘     │
│                                                         │
│  Needles (leaf_size ≤ 3 mm):                            │
│  ┌────────────────────────────────────────────────┐     │
│  │  gcyl = 0.00403,  jcyl = 0.6                   │     │
│  │  r_bl = 1 / (1.5 × gcyl × WS^0.6               │     │
│  │              / (leaf_size/1000)^0.4 )          │     │
│  │                                                │     │
│  │  Cylinder geometry for thin needles            │     │
│  └────────────────────────────────────────────────┘     │
│                                                         │
│  Convert to SurEau conductance units:                   │
│  g_bl = (1 / r_bl) × 1000 × 40       [mmol/m²/s]        │
│                                                         │
│  Combined non-evaporative resistance (parallel):        │
│  r_blr = 1 / (1/r_bl + 1/r_r)                [s/m]      │
│                                                         │
│  Physical meaning: two pathways for non-evaporative     │
│  heat loss operate simultaneously:                      │
│    • Convection through the boundary layer (r_bl)       │
│    • Longwave radiation to surroundings (r_r)           │
│  Parallel combination = total non-evap resistance       │
│                                                         │
│  Example: r_bl ≈ 90 s/m, r_r ≈ 220 s/m                  │
│           r_blr ≈ 64 s/m (parallel always < smallest)   │
│                                                         │
└────────────────────────┬────────────────────────────────┘
                         │
                         ▼
┌─────────────────────────────────────────────────────────┐
│                                                         │
│    STAGE 4: STOMATAL RESISTANCE                         │
│    "How hard is it for water vapor to escape?"          │
│                                                         │
│  Depends on transpiration_model:                        │
│                                                         │
│  ┌──── Jarvis ──────────────────────────────────────┐   │
│  │                                                  │   │
│  │  Total leaf conductance (stomata + cuticle       │   │
│  │  in parallel, as in Eq. 29):                     │   │
│  │  g_total = gs + g_cuti       [mmol/m²/s]         │   │
│  │                                                  │   │
│  │  Convert to resistance:                          │   │
│  │  r_st = (1 / g_total) × 1000 × 40      [s/m]     │   │
│  │                                                  │   │
│  │  If g_total = 0:  r_st = 9999.99 (≈ infinite)    │   │
│  │                                                  │   │
│  │  Open:   gs=100 + g_cuti=4 → r_st ≈ 385 s/m      │   │
│  │  Closed: gs=0   + g_cuti=4 → r_st ≈ 10000 s/m    │   │
│  │                                                  │   │
│  └──────────────────────────────────────────────────┘   │
│                                                         │
│  ┌──── Granier ─────────────────────────────────────┐   │
│  │                                                  │   │
│  │  Back-calculate conductance from flux:           │   │
│  │  g = E_inst / VPDx × P_atm                       │   │
│  │                                                  │   │
│  │  Convert to resistance:                          │   │
│  │  r_st = (1 / g) × 1000 × 40            [s/m]     │   │
│  │                                                  │   │
│  └──────────────────────────────────────────────────┘   │
│                                                         │
└────────────────────────┬────────────────────────────────┘
                         │
                         ▼
┌─────────────────────────────────────────────────────────┐
│                                                         │
│    STAGE 5: PENMAN-MONTEITH TEMPERATURE SOLUTION        │
│    "What temperature balances energy in vs energy out?" │
│                                                         │
│  Modified psychrometric constant:                       │
│  γ* = γ × (r_st / r_blr)                                │
│                                                         │
│    γ = 0.066 kPa/°C (psychrometric constant)            │
│    r_st / r_blr = ratio of stomatal to non-evaporative  │
│                   resistance                            │
│                                                         │
│    Stomata open:  r_st ≈ 385,  γ* = 0.066×385/64        │
│                                   = 0.40  (small)       │
│    Stomata closed: r_st ≈ 10000, γ* = 0.066×10000/64    │
│                                   = 10.3  (large)       │
│                                                         │
│  Leaf-to-air temperature difference:                    │
│  ┌─────────────────────────────────────────────────┐    │
│  │                                                 │    │
│  │         γ* × Rni × r_blr / (ρ·cp)  −  VPDx      │    │
│  │  δT  = ─────────────────────────────────────    │    │
│  │                    Δ  +  γ*                     │    │
│  │                                                 │    │
│  │  Numerator:                                     │    │
│  │    1st term: heating from absorbed radiation    │    │
│  │              (larger when stomata closed → γ*↑) │    │
│  │    2nd term: cooling from transpiration         │    │
│  │              (larger when air is dry → VPD↑)    │    │
│  │                                                 │    │
│  │  Denominator:                                   │    │
│  │    Δ + γ* = total system sensitivity            │    │
│  │                                                 │    │
│  └─────────────────────────────────────────────────┘    │
│                                                         │
│  T_leaf = T_air + δT                          [°C]      │
│                                                         │
│  Example (open stomata, γ*=0.40):                       │
│    δT = (0.40×85×64/1306 − 1.59) / (0.189+0.40)         │
│       = (1.66 − 1.59) / 0.59 = +0.12°C                  │
│    T_leaf = 25.12°C  (barely warmer — transpiration     │
│                       cools effectively)                │
│                                                         │
│  Example (closed stomata, γ*=10.3):                     │
│    δT = (10.3×85×64/1306 − 1.59) / (0.189+10.3)         │
│       = (42.9 − 1.59) / 10.49 = +3.94°C                 │
│    T_leaf = 28.94°C  (nearly 4°C warmer — much less     │
│                       evaporative cooling)              │
│                                                         │
└────────────────────────┬────────────────────────────────┘
                         │
                         ▼
┌─────────────────────────────────────────────────────────┐
│                                                         │
│    STAGE 6: LEAF VPD WITH KELVIN CORRECTION             │
│   "What is the actual driving force for transpiration?" │
│                                                         │
│  Saturation VP at T_air (Buck 1981, in Pa):             │
│  e_sat_air = 611.21 × exp((18.678 − T/234.5)            │
│                            × T / (257.14 + T))          │
│                                                         │
│  Actual VP of ambient air:                              │
│  e_air = e_sat_air × RH / 100                  [Pa]     │
│                                                         │
│  Saturation VP at T_leaf (Buck 1981, in Pa):            │
│  e_sat_leaf = 611.21 × exp((18.678 − Tl/234.5)          │
│                              × Tl / (257.14 + Tl))      │
│                                                         │
│  NOTE: e_sat_leaf > e_sat_air when T_leaf > T_air       │
│  (warmer leaf → higher vapor pressure inside)           │
│                                                         │
│  Kelvin equation correction:                            │
│  ┌─────────────────────────────────────────────────┐    │
│  │                                                 │    │
│  │  e_leaf = e_sat_leaf × exp(ψ × Vw / (R × T))    │    │
│  │                                                 │    │
│  │  In code:                                       │    │
│  │  e_leaf = e_sat_leaf × exp(psi_leaf             │    │
│  │                            × 2.16947115         │    │
│  │                            / (T_leaf + 273.15)) │    │
│  │                                                 │    │
│  │  where 2.16947115 ≈ Vw × 10⁶ / R                │    │
│  │    Vw = 18.02e−6 m³/mol  (molar volume H₂O)     │    │
│  │    R  = 8.314 J/mol/K    (gas constant)         │    │
│  │    10⁶ converts ψ from MPa to Pa                │    │
│  │                                                 │    │
│  │  Physical meaning: negative ψ reduces the       │    │
│  │  water activity inside the leaf, lowering the   │    │
│  │  effective vapor pressure below e_sat.          │    │
│  │                                                 │    │
│  │  Impact by drought severity:                    │    │
│  │  ψ =  −2 MPa → correction = 0.986 (−1.4%)       │    │
│  │  ψ =  −5 MPa → correction = 0.964 (−3.6%)       │    │
│  │  ψ = −10 MPa → correction = 0.930 (−7.0%)       │    │
│  │                                                 │    │
│  │  Only significant under extreme drought,        │    │
│  │  the regime SurEau-Ecos was designed for.       │    │
│  │                                                 │    │
│  └─────────────────────────────────────────────────┘    │
│                                                         │
│  Leaf VPD:                                              │
│  VPD_leaf = max(0, (e_leaf − e_air) / 1000)    [kPa]    │
│                                                         │
│  NOTE: VPD_leaf uses e_leaf (corrected for both         │
│  T_leaf and ψ_leaf) versus e_air (at T_air).            │
│  This is the ACTUAL gradient driving water out          │
│  of the leaf, not the ambient VPD.                      │
│                                                         │
└────────────────────────┬────────────────────────────────┘
                         │
                         ▼
┌─────────────────────────────────────────────────────────┐
│                                                         │
│    RETURN LOGIC                                         │
│                                                         │
│  if turn_off_EB:                                        │
│      return T_air, g_bl, VPD_leaf                       │
│      (skip energy balance — diagnostic mode)            │
│                                                         │
│  else:                                                  │
│      return T_leaf, g_bl, VPD_leaf                      │
│      (normal operation)                                 │
│                                                         │
└─────────────────────────────────────────────────────────┘
 
 OUTPUTS
 ───────
   T_leaf    [°C]          → stored as fluxes.leaf_temperature
                             used by compute_gmin() (Eqs. 31–32)
                             to compute temperature-dependent
                             cuticular conductance
 
   g_bl      [mmol/m²/s]   → stored as fluxes.g_BL
                             used by compute_E_min() (Eqs. 29–30)
                             in the conductance chain:
                             1/g_total = 1/g_stom + 1/g_bl + 1/g_crown
 
   VPD_leaf  [kPa]         → stored as fluxes.leaf_VPD
                             drives ALL transpiration calculations:
                             • E_min  = g_min × VPD_leaf / P_atm
                             • E_lim  = g_canopy_lim × VPD_leaf / P_atm
                             • E_bound = g_canopy_bound × VPD_leaf / P_atm
 
 CALLED BY
 ─────────
   compute_transpiration()  (sureau_vegetation.py)
     — called every sub-daily timestep
     — Jarvis model calls it TWICE per timestep:
         1st: with current gs → get initial T_leaf
         2nd: after updating gs from regulation → get final T_leaf
 
 REFERENCES
 ──────────
   Cochard et al. (2021)  — CPRM21 — energy balance formulation
   Jones (2013)           — boundary layer conductance
   Brutsaert (1975)       — atmospheric emissivity
   Buck (1981)            — saturation vapor pressure (Pa version)
   Tetens (via Buck 1981) — saturation vapor pressure (kPa version)
   Kelvin equation        — water activity correction at low ψ

Dead fuel moisture

---

source

compute_DFMC

def compute_DFMC(
    VPD, FM0:float=5.43, FM1:float=52.91, m:float=0.64
):

Dead fuel moisture content [%dry mass] (De Dios et al. 2015).