Optimal LAI model

Author: kmdeck

docs/list_standalone_models.jl

@@ -20,6 +20,7 @@ standalone_models = [
         ]
         "Canopy structure" => [
             "Prescribed structure" => "standalone/pages/vegetation/canopy_structure/prescribed_structure.md",
+            "Optimal LAI model" => "standalone/pages/vegetation/canopy_structure/optimal_lai.md",
         ]
         "Canopy energy" => [
             "Physics" => "standalone/pages/vegetation/canopy_energy/canopy_energy.md",

docs/src/standalone/pages/vegetation/canopy_structure/optimal_lai.md

@@ -0,0 +1,101 @@
+# Optimal LAI Model
+
+The Optimal LAI model predicts seasonal to decadal dynamics of leaf area index based on optimality principles, balancing energy and water constraints.
+
+This model is based on Zhou et al. (2025), which presents a general model for the seasonal to decadal dynamics of leaf area that combines predictions from both the light use efficiency (LUE) framework and optimization theory.
+
+## Model Overview
+
+The optimal LAI model computes LAI dynamically by:
+1. Calculating the seasonal maximum LAI (LAI$_{max}$) based on energy and water limitations
+2. Computing steady-state LAI from daily meteorological conditions
+3. Updating actual LAI using an exponential moving average to represent the lag in leaf development
+
+## Seasonal Maximum LAI
+
+The seasonal maximum LAI (LAI$_{max}$) is determined by the minimum of energy-limited and water-limited constraints:
+
+```math
+\begin{align}
+\text{fAPAR}_{max} &= \min\left(\text{fAPAR}_{energy}, \text{fAPAR}_{water}\right) \\
+\text{fAPAR}_{energy} &= 1 - \frac{z}{k \cdot A_{0,annual}} \\
+\text{fAPAR}_{water} &= \frac{c_a(1-\chi)}{1.6 \cdot D_{growing}} \cdot \frac{f_0 \cdot P_{annual}}{A_{0,annual}} \\
+\text{LAI}_{max} &= -\frac{1}{k} \ln(1 - \text{fAPAR}_{max})
+\end{align}
+```
+
+where:
+- $A_{0,annual}$ is the annual total potential GPP (mol m⁻² yr⁻¹)
+- $P_{annual}$ is the annual total precipitation (mol m⁻² yr⁻¹)
+- $D_{growing}$ is the mean vapor pressure deficit during the growing season (Pa)
+- $k$ is the light extinction coefficient (dimensionless)
+- $z$ is the unit cost of constructing and maintaining leaves (mol m⁻² yr⁻¹)
+- $c_a$ is the ambient CO₂ partial pressure (Pa)
+- $\chi$ is the ratio of leaf-internal to ambient CO₂ partial pressure (dimensionless)
+- $f_0$ is the fraction of annual precipitation used by plants (dimensionless)
+
+## Daily Steady-State LAI
+
+Given daily meteorological conditions, the steady-state LAI ($L_s$) represents the LAI that would be in equilibrium with GPP if conditions were held constant:
+
+```math
+\begin{align}
+\mu &= m \cdot A_{0,daily} \\
+L_s &= \min\left\{\mu + \frac{1}{k} W_0[-k\mu \exp(-k\mu)], \text{LAI}_{max}\right\}
+\end{align}
+```
+
+where:
+- $A_{0,daily}$ is the daily potential GPP (mol m⁻² day⁻¹)
+- $m$ is a parameter relating steady-state LAI to steady-state GPP, computed as:
+
+```math
+m = \frac{\sigma \cdot \text{GSL} \cdot \text{LAI}_{max}}{A_{0,annual} \cdot \text{fAPAR}_{max}}
+```
+
+where GSL is the growing season length (days) and $\sigma$ is a dimensionless parameter representing departure from square-wave LAI dynamics.
+
+- $W_0$ is the principal branch of the Lambert W function
+
+## LAI Update
+
+The actual LAI is updated using an exponential weighted moving average to represent the time lag for photosynthate allocation to leaves:
+
+```math
+\text{LAI}_{new} = \alpha \cdot L_s + (1-\alpha) \cdot \text{LAI}_{prev}
+```
+
+where $\alpha$ is a smoothing factor (dimensionless, 0-1). Setting $\alpha = 0.067$ corresponds to approximately 15 days of memory.
+
+## Parameters
+
+| Parameter | Symbol | Unit | Typical Value | Description |
+| :--- | :---: | :---: | :---: | :--- |
+| Light extinction coefficient | $k$ | - | 0.5 | Controls light attenuation through canopy |
+| Leaf construction cost | $z$ | mol m⁻² yr⁻¹ | 12.227 | Unit cost of building and maintaining leaves |
+| CO₂ concentration ratio | $\chi$ | - | 0.7 | Ratio of leaf-internal to ambient CO₂ |
+| Precipitation fraction | $f_0$ | - | 0.62 | Fraction of precipitation used by plants |
+| LAI dynamics parameter | $\sigma$ | - | 0.771 | Departure from square-wave dynamics |
+| Smoothing factor | $\alpha$ | - | 0.067 | Controls LAI response time (~15 days) |
+
+## Drivers
+
+| Driver | Symbol | Unit | Description |
+| :--- | :---: | :---: | :--- |
+| Daily potential GPP | $A_{0,daily}$ | mol m⁻² day⁻¹ | GPP with fAPAR = 1 |
+| Annual potential GPP | $A_{0,annual}$ | mol m⁻² yr⁻¹ | Yearly integral of $A_0$ |
+| Annual precipitation | $P_{annual}$ | mol m⁻² yr⁻¹ | Total yearly precipitation |
+| Growing season VPD | $D_{growing}$ | Pa | Mean VPD when T > 0°C |
+| Growing season length | GSL | days | Length of continuous period T > 0°C |
+| CO₂ partial pressure | $c_a$ | Pa | Ambient CO₂ concentration |
+
+## Output
+
+| Output | Symbol | Unit | Range |
+| :--- | :---: | :---: | :---: |
+| Leaf Area Index | LAI | m² m⁻² | 0-10 |
+
+## References
+
+Zhou et al. (2025) "A General Model for the Seasonal to Decadal Dynamics of Leaf Area" 
+Global Change Biology. https://onlinelibrary.wiley.com/doi/pdf/10.1111/gcb.70125

experiments/integrated/fluxnet/ozark_pmodel.jl

@@ -163,11 +163,11 @@ radiation_parameters = (;
     Ω,
     G_Function = CLMGFunction(χl),
     α_PAR_leaf,
-    τ_PAR_leaf,
+    #τ_PAR_leaf,
     α_NIR_leaf,
-    τ_NIR_leaf,
+    #τ_NIR_leaf,
 )
-radiative_transfer = Canopy.TwoStreamModel{FT}(
+radiative_transfer = Canopy.BeerLambertModel{FT}(
     surface_domain,
     toml_dict;
     radiation_parameters,
@@ -181,6 +181,11 @@ conductance = PModelConductance{FT}(toml_dict)
 is_c3 = FT(1)
 photosynthesis = PModel{FT}(surface_domain, toml_dict; is_c3)
 
+# Set up optimal LAI model
+lai_model = Canopy.OptimalLAIModel{FT}(
+    Canopy.OptimalLAIParameters{FT}(toml_dict)
+)
+
 # Set up soil moisture stress using soil retention parameters
 soil_moisture_stress = PiecewiseMoistureStressModel{FT}(
     land_domain,
@@ -230,6 +235,7 @@ canopy = Canopy.CanopyModel{FT}(
     radiative_transfer,
     photosynthesis,
     conductance,
+    lai_model,
     soil_moisture_stress,
     hydraulics,
     energy,
@@ -255,7 +261,7 @@ set_ic! = FluxnetSimulations.make_set_fluxnet_initial_conditions(
     land,
 )
 # Callbacks
-output_vars = ["gpp", "shf", "lhf", "swu", "lwu", "swc", "swe", "tsoil"]
+output_vars = ["gpp", "shf", "lhf", "swu", "lwu", "swc", "swe", "tsoil", "lai", "lai_pred"]
 diags = ClimaLand.default_diagnostics(
     land,
     start_date;
@@ -277,6 +283,7 @@ simulation = LandSimulation(
 
 @time solve!(simulation)
 
+#=
 comparison_data = FluxnetSimulations.get_comparison_data(site_ID, time_offset)
 savedir = joinpath(
     pkgdir(ClimaLand),
@@ -301,3 +308,57 @@ LandSimVis.make_timeseries(
     spinup_date = start_date + Day(20),
     comparison_data,
 )
+=#
+
+# Need to solve!
+# can also look at simulation._integrator.p
+
+short_names = [d.variable.short_name for d in simulation.diagnostics] # short_name_X_average e.g.
+diag_names = [d.output_short_name for d in simulation.diagnostics] # short_name_X_average e.g.
+diag_units = [d.variable.units for d in simulation.diagnostics]
+i = 10
+    dn = diag_names[i]
+    unit = diag_units[i]
+    sn = short_names[i]
+    model_time, LAI_opt = ClimaLand.Diagnostics.diagnostic_as_vectors(
+                                                                           simulation.diagnostics[1].output_writer,
+                                                                           dn,
+                                                                          )
+LAI_opt
+
+
+i = 9
+    dn = diag_names[i]
+    unit = diag_units[i]
+    sn = short_names[i]
+    model_time, LAI_obs = ClimaLand.Diagnostics.diagnostic_as_vectors(
+                                                                           simulation.diagnostics[1].output_writer,
+                                                                           dn,
+                                                                          )
+LAI_obs
+
+        save_Δt = model_time[2] - model_time[1] # in seconds since the start_date. if model_time is an Itime, the epoch should be start_date
+        import ClimaUtilities.TimeManager: ITime, date
+        function time_to_date(t::ITime, start_date)
+    start_date != t.epoch &&
+        @warn("$(start_date) is different from the simulation time epoch.")
+    return isnothing(t.epoch) ? start_date + t.counter * t.period : date(t)
+end
+        model_dates = time_to_date.(model_time, start_date)
+
+
+        fig = Figure()
+        ax = Axis(fig[1,1])
+        p = lines!(ax, model_dates, LAI_obs, label = "LAI obs")
+        p2 = lines!(ax, model_dates, LAI_opt, label = "LAI opt")
+        ax.ylabel= "LAI m2 m-2"
+        axislegend()
+        save("test.png", fig)
+
+
+# Thoughts:
+# it makes no sense to have LAI < 0
+# maybe params (α, z, m) are not adequate with our units of A (GPP)
+# Solutions:
+# we could convert A internally
+# and force L to be 0 min

src/diagnostics/default_diagnostics.jl

@@ -397,6 +397,7 @@ function get_possible_diagnostics(model::CanopyModel)
         # "fa", # return a Tuple
         "far",
         "lai",
+        "lai_pred",
         "msf",
         "rai",
         "sai",
@@ -528,7 +529,7 @@ function get_short_diagnostics(model::SoilCO2Model)
     return ["sco2"]
 end
 function get_short_diagnostics(model::CanopyModel)
-    return ["gpp", "ct", "lai", "trans", "er", "sif"]
+    return ["gpp", "ct", "lai", "lai_pred","trans", "er", "sif"]
 end
 function get_short_diagnostics(model::SnowModel)
     return get_possible_diagnostics(model)

src/diagnostics/define_diagnostics.jl

@@ -309,6 +309,16 @@ function define_diagnostics!(land_model)
             compute_leaf_area_index!(out, Y, p, t, land_model),
     )
 
+        add_diagnostic_variable!(
+        short_name = "lai_pred",
+        long_name = "Leaf area Index",
+        standard_name = "leaf_area_index",
+        units = "m^2 m^-2",
+        comments = "The area index of leaves, expressed in surface area of leaves per surface area of ground.",
+        compute! = (out, Y, p, t) ->
+            compute_leaf_area_index_pred!(out, Y, p, t, land_model),
+    )
+
     # Moisture stress factor
     add_diagnostic_variable!(
         short_name = "msf",

src/diagnostics/land_compute_methods.jl

@@ -257,6 +257,11 @@ end
     LandModel,
     CanopyModel,
 } p.canopy.biomass.area_index.leaf
+@diagnostic_compute "leaf_area_index_pred" Union{
+    SoilCanopyModel,
+    LandModel,
+    CanopyModel,
+} p.canopy.lai_model.LAI
 
 # Canopy - Soil moisture stress
 @diagnostic_compute "moisture_stress_factor" Union{

src/standalone/Vegetation/Canopy.jl

@@ -48,6 +48,7 @@ include("./stomatalconductance.jl")
 include("./photosynthesis.jl")
 include("./photosynthesis_farquhar.jl")
 include("./pmodel.jl")
+include("./optimal_lai.jl")
 include("./radiation.jl")
 include("./solar_induced_fluorescence.jl")
 include("./pfts.jl")
@@ -98,8 +99,8 @@ and wilting point (θ_low) from the soil parameters.
 The low and high thresholds for the piecewise soil moisture stress function are given by
 the residual soil water content and the soil porosity, respectively.
 
-The soil parameters should be a named tuple with keys of `ν` and `θ_r` 
-(additional keys may be present but they will not be used). These may 
+The soil parameters should be a named tuple with keys of `ν` and `θ_r`
+(additional keys may be present but they will not be used). These may
 be ClimaCore fields or floats, but must be consistently one or the other.
 
 Note that unlike other Canopy components, this model is defined on the subsurface domain.
@@ -606,7 +607,7 @@ treated differently.
 
 $(DocStringExtensions.FIELDS)
 """
-struct CanopyModel{FT, AR, RM, PM, SM, SMSM, PHM, EM, SIFM, BM, B, PS, D} <:
+struct CanopyModel{FT, AR, RM, PM, SM, SMSM, PHM, EM, SIFM, LAI, BM, B, PS, D} <:
        ClimaLand.AbstractImExModel{FT}
     "Autotrophic respiration model, a canopy component model"
     autotrophic_respiration::AR
@@ -624,6 +625,8 @@ struct CanopyModel{FT, AR, RM, PM, SM, SMSM, PHM, EM, SIFM, BM, B, PS, D} <:
     energy::EM
     "SIF model, a canopy component model"
     sif::SIFM
+    "LAI model, a canopy component model"
+    lai_model::LAI
     "Biomass parameterization, a canopy component model"
     biomass::BM
     "Boundary Conditions"
@@ -644,6 +647,7 @@ end
         hydraulics::AbstractPlantHydraulicsModel{FT},
         energy::AbstractCanopyEnergyModel{FT},
         sif::AbstractSIFModel{FT},
+        lai_model::AbstractLAIModel{FT},
         biomass::AbstractBiomassModel{FT},
         boundary_conditions::B,
         parameters::SharedCanopyParameters{FT, PSE},
@@ -669,6 +673,7 @@ function CanopyModel{FT}(;
     hydraulics::AbstractPlantHydraulicsModel{FT},
     soil_moisture_stress::AbstractSoilMoistureStressModel{FT},
     sif::AbstractSIFModel{FT},
+    lai_model::AbstractLAIModel{FT},
     energy = PrescribedCanopyTempModel{FT}(),
     biomass::PrescribedBiomassModel{FT},
     boundary_conditions::B,
@@ -708,6 +713,7 @@ function CanopyModel{FT}(;
         hydraulics,
         energy,
         sif,
+        lai_model,
         biomass,
         boundary_conditions,
         parameters,
@@ -780,6 +786,7 @@ function CanopyModel{FT}(
     energy = BigLeafEnergyModel{FT}(toml_dict),
     biomass = PrescribedBiomassModel{FT}(domain, LAI, toml_dict),
     sif = Lee2015SIFModel{FT}(toml_dict),
+    lai_model = OptimalLAIModel{FT}(toml_dict),
 ) where {FT}
     (; atmos, radiation, ground) = forcing
 
@@ -799,6 +806,7 @@ function CanopyModel{FT}(
         energy,
         biomass,
         sif,
+        lai_model,
     ]
         # For component models without parameters, skip the check
         !hasproperty(component, :parameters) && continue
@@ -831,6 +839,7 @@ function CanopyModel{FT}(
         hydraulics,
         energy,
         sif,
+        lai_model,
         biomass,
         boundary_conditions,
         parameters,
@@ -859,6 +868,7 @@ canopy_components(::CanopyModel) = (
     :autotrophic_respiration,
     :energy,
     :sif,
+    :lai_model,
     :soil_moisture_stress,
     :biomass,
 )