from bonan_hydraulics.utils import satvap
from bonan_hydraulics.leaf_boundary_layer import *
from bonan_hydraulics.utils import calc_radiative_forcing_qa
from bonan_hydraulics.parameter_classes import Params, PhysCon, Atmos, LeafLeaf temperature
Start with initial guess: Tleaf = Tair
│
▼
┌──────────────────────────────┐
│ Calculate at current Tleaf: │
│ • esat(Tleaf), desat │ ← How much moisture the leaf "wants" to release
│ • lwrad = 2·ε·σ·T⁴ │ ← Infrared glow (both sides)
│ • shflx = 2·cp·(T-Ta)·gbh │ ← Heat to air (both sides)
│ • lhflx = λ·(esat-ea)·glw │ ← Evaporative cooling
└──────────────────────────────┘
│
▼
f0 = Qa - lwrad - shflx - lhflx ← Energy imbalance
│
▼
Is |f0| < 0.000001 W/m² ?
/ \
YES NO
│ │
▼ ▼
DONE dtleaf = -f0/df0
│ Tleaf += dtleaf
│ │
│ └──→ Go back to top
▼
Calculate final:
• rnet = Qa - lwrad
• etflx = lhflx / λ
• Verify: rnet - shflx - lhflx ≈ 0leaf_temperature
def leaf_temperature(
physcon:PhysCon, # Physical constants:
- tfrz : float
Freezing point of water (K).
- mmh2o : float
Molecular mass of water (kg/mol).
- sigma : float
Stefan-Boltzmann constant (W/m2/K4).
atmos:Atmos, # Atmospheric forcing variables:
- patm : float
Atmospheric pressure (Pa).
- cpair : float
Specific heat of air at constant pressure (J/mol/K).
- tair : float
Air temperature (K).
- eair : float
Vapor pressure of air (Pa).
leaf:Leaf, # Leaf parameters:
- emiss : float
Leaf emissivity (-).
- tleaf : float
Leaf temperature (K). Used as initial guess on input
flux:Flux, # Flux variables with the following inputs:
- gbh : float
Leaf boundary layer conductance for heat (mol/m2 leaf/s).
- gbv : float
Leaf boundary layer conductance for H2O (mol H2O/m2 leaf/s).
- gs : float
Leaf stomatal conductance (mol H2O/m2 leaf/s).
- qa : float
Leaf radiative forcing (W/m2 leaf).
)->Flux: # Updated flux object with the following attributes:
- tleaf : float
Leaf temperature (K). Updated to energy-balance solution on output.
- rnet : float
Leaf net radiation (W/m2 leaf).
- lwrad : float
Longwave radiation emitted from leaf (W/m2 leaf).
- shflx : float
Leaf sensible heat flux (W/m2 leaf).
- lhflx : float
Leaf latent heat flux (W/m2 leaf).
- etflx : float
Leaf transpiration flux (mol H2O/m2 leaf/s).
Calculate leaf temperature and energy fluxes.
Solve the leaf energy balance using Newton-Raphson iteration to find the temperature that balances radiative forcing against longwave emission, sensible heat flux, and latent heat flux.
Parameters:
PhysCon: Physical constants
- tfrz: Freezing point of water (K).
- mmh2o: Molecular mass of water (kg/mol).
- sigma: Stefan-Boltzmann constant (W/m2/K4).
Atmos: Atmospheric forcing variables:
- patm: Atmospheric pressure (Pa).
- cpair: Specific heat of air at constant pressure (J/mol/K).
- tair: Air temperature (K)
- eair: Vapor pressure of air (Pa).
Leaf: Leaf parameters
- emiss: Leaf emissivity (-).
- tleaf: Leaf temperature (K). Used as initial guess on input
Flux: variables with the following inputs:
- gbh: Leaf boundary layer conductance for heat (mol/m2 leaf/s).
- gbv: Leaf boundary layer conductance for H2O (mol H2O/m2 leaf/s).
- gs: Leaf stomatal conductance (mol H2O/m2 leaf/s).
- qa: Leaf radiative forcing (W/m2 leaf).
Returns:
Flux: Updated flux object with the following attributes:
- tleaf: Leaf temperature (K). Updated to energy-balance solution on output.
- rnet: Leaf net radiation (W/m2 leaf).
- lwrad: Longwave radiation emitted from leaf (W/m2 leaf).
- shflx: Leaf sensible heat flux (W/m2 leaf).
- lhflx: Leaf latent heat flux (W/m2 leaf).
- etflx: Leaf transpiration flux (mol H2O/m2 leaf/s).
Example leaf_temperature()
Reproducing the results from table 10.2 of the book
Load packages
Call objects
params = Params()
phys_constants = PhysCon
atmos = Atmos()
leaf = Leaf()Edit object parameters
# Atmospheric pressure
atmos.patm = 101325.0
# Ambient temp in Kelvin
atmos.tair = 35.0 + phys_constants.tfrz
relative_humidity = 50
# Saturation vapor pressure at 35°C. Ignoring dsat that is why the underscore (_)
esat_air, _ = satvap(35.0)
# Actual vapor pressure (Pa)
atmos.eair = (relative_humidity / 100.0) * esat_air
# Molar density of air: ρm = P / (R * T)
atmos.rhomol = atmos.patm / (phys_constants.rgas * atmos.tair)
# Specific heat of moist air (J/mol/K)
# cp = cpd * Md + cpw * Mw * (eair/patm), simplified to ~29.2 J/mol/K
atmos.cpair = phys_constants.cpd * phys_constants.mmdry
# Longwave radiation
atmos.irsky = 300.0# Leaf dimension (m)
leaf.dleaf = 0.05
# Emissivity
leaf.emiss = 0.98
# Reflectance [vis, nir]
leaf.rho = [0.1, 0.4]
# Transmittance [vis, nir]
leaf.tau = [0.1, 0.4]# Ground albedo [vis, nir]
ground_albedo = [0.1, 0.2]
# Ground longwave (W/m2)
ground_lw = 545.0# Calculate Qa for clear and cloudy sky
qa_clear = calc_radiative_forcing_qa(
1100.0, leaf, ground_albedo, ground_lw, atmos.irsky
)
qa_cloudy = calc_radiative_forcing_qa(
550.0, leaf, ground_albedo, ground_lw, atmos.irsky
)# Loop over conditions and compute leaf temperature
#
# For each combination of:
# - Sky condition (clear, cloudy)
# - Stomatal conductance (0, 0.4 mol/m2/s)
# - Wind speed (0.01, 0.1, 1, 5 m/s)
#
# (a) Compute boundary layer conductances (gbh, gbv, gbc)
# (b) Set stomatal conductance
# (c) Solve for leaf temperature via Newton-Raphson
wind_speeds = [0.01, 0.1, 1.0, 5.0]
gsw_values = [0.0, 0.4]
sky_conditions = [
("Clear sky", qa_clear, 1444),
("Cloudy sky", qa_cloudy, 1136),
]for sky_name, qa, qa_book in sky_conditions:
for each_gsw in gsw_values:
for each_wind in wind_speeds:
print(f"\ngs = {each_gsw}, wind_speed = {each_wind}")
# Set wind speed for boundary layer calculation
atmos.wind = each_wind
# Create fresh flux object for each case
flux = Flux()
# Initial guess = air temperature
flux.tleaf = atmos.tair
# Radiative forcing
flux.qa = qa
# (a) Set stomatal conductance
# Note: gs = 0 causes division by zero in gleaf = gs*gbv/(gs+gbv)
# Use a tiny value instead to represent closed stomata
flux.gs = max(each_gsw, 1e-06)
# (b) Iterate: boundary layer depends on Tleaf (via Grashof number)
# and Tleaf depends on boundary layer conductance.
# At low wind speeds, free convection matters and this
# coupling must be resolved iteratively.
for iteration in range(20):
tleaf_old = flux.tleaf
# Compute boundary layer conductances at current Tleaf
flux = leaf_boundary_layer(phys_constants, atmos, leaf, flux)
# Solve energy balance for leaf temperature
flux = leaf_temperature(phys_constants, atmos, leaf, flux)
# Check convergence
if abs(flux.tleaf - tleaf_old) < 0.01:
break
# Convert to Celsius and print
tleaf_C = flux.tleaf - phys_constants.tfrz
print(f"Leaf temperature = {tleaf_C:6.1f}\n", end="")
gs = 0.0, wind_speed = 0.01
Leaf temperature = 49.4
gs = 0.0, wind_speed = 0.1
Leaf temperature = 45.8
gs = 0.0, wind_speed = 1.0
Leaf temperature = 40.9
gs = 0.0, wind_speed = 5.0
Leaf temperature = 38.2
gs = 0.4, wind_speed = 0.01
Leaf temperature = 39.8
gs = 0.4, wind_speed = 0.1
Leaf temperature = 37.7
gs = 0.4, wind_speed = 1.0
Leaf temperature = 35.8
gs = 0.4, wind_speed = 5.0
Leaf temperature = 35.2
gs = 0.0, wind_speed = 0.01
Leaf temperature = 39.9
gs = 0.0, wind_speed = 0.1
Leaf temperature = 38.5
gs = 0.0, wind_speed = 1.0
Leaf temperature = 36.9
gs = 0.0, wind_speed = 5.0
Leaf temperature = 36.0
gs = 0.4, wind_speed = 0.01
Leaf temperature = 34.1
gs = 0.4, wind_speed = 0.1
Leaf temperature = 33.5
gs = 0.4, wind_speed = 1.0
Leaf temperature = 32.9
gs = 0.4, wind_speed = 5.0
Leaf temperature = 33.4