Kd/rb/ar/basic sif refactor

docs/src/APIs/canopy/Photosynthesis.md

@@ -7,6 +7,7 @@ CurrentModule = ClimaLand.Canopy
 ## Parameters
 
 ```@docs
+ClimaLand.Canopy.SIFParameters
 ClimaLand.Canopy.FarquharParameters
 ClimaLand.Canopy.OptimalityFarquharParameters
 ```
@@ -31,4 +32,4 @@ ClimaLand.Canopy.compute_GPP
 ClimaLand.Canopy.MM_Kc
 ClimaLand.Canopy.MM_Ko
 ClimaLand.Canopy.compute_Vcmax
-```
+```

lib/ClimaLandSimulations/src/utilities/climaland_output_dataframe.jl

@@ -57,6 +57,7 @@ function make_output_df(
         collect(map(i -> (i, :soil, :T), 1:20)), # 20 shouldn't be hard-coded, but an arg, equal to n layers
         collect(map(i -> (i, :soil, :θ_l), 1:20)),
         (1, :soil, :turbulent_fluxes, :vapor_flux),
+        (1, :canopy, :sif, :SIF),
     )
 
     output_vectors = [getoutput(sv, args...) for args in output_list]

lib/ClimaLandSimulations/src/utilities/makie_plots.jl

@@ -401,8 +401,9 @@ function diurnals_fig(inputs, climaland, earth_param_set; dashboard = false) # w
         ylabel = L"\text{CO}_{2} \, (\mu\text{mol m}^{-2} \, \text{s}^{-1})",
     ) # C fluxes
     ax_W = Axis(fig[2, 1], ylabel = L"\text{H}_{2}\text{O} \, \text{(mm)}") # h2o fluxes
+    ax_SIF = Axis(fig[3, 1], ylabel = L"\text{SIF}") # SIF
     ax_E = Axis(
-        fig[3, 1],
+        fig[4, 1],
         ylabel = L"\text{Radiation} \, (\text{W} \, \text{m}^{-2})",
         xlabel = L"\text{Hour of the day}",
         xgridvisible = false,
@@ -456,6 +457,11 @@ function diurnals_fig(inputs, climaland, earth_param_set; dashboard = false) # w
         linestyle = :dot,
     )
 
+    # SIF
+    p_SIF_m =
+        diurnal_plot!(fig, ax_SIF, climaland.DateTime, climaland.SIF, :black)
+
+
     # Energy fluxes
     # model
     # diurnal_plot!(fig, ax_E, climaland.DateTime, climaland.LW_out, :red)
@@ -472,7 +478,7 @@ function diurnals_fig(inputs, climaland, earth_param_set; dashboard = false) # w
         linestyle = :dot,
     )
 
-    [xlims!(axes, (0, 24)) for axes in [ax_C, ax_W, ax_E]]
+    [xlims!(axes, (0, 24)) for axes in [ax_C, ax_W, ax_SIF, ax_E]]
 
     axislegend(
         ax_C,
@@ -490,6 +496,14 @@ function diurnals_fig(inputs, climaland, earth_param_set; dashboard = false) # w
         position = :rt,
         orientation = :vertical,
     )
+    axislegend(
+        ax_SIF,
+        [p_SIF_m],
+        ["SIF model"],
+        "",
+        position = :rt,
+        orientation = :vertical,
+    )
     axislegend(
         ax_E,
         [p_SWout_d, p_SWout_m],
@@ -501,6 +515,7 @@ function diurnals_fig(inputs, climaland, earth_param_set; dashboard = false) # w
 
     hidexdecorations!(ax_C)
     hidexdecorations!(ax_W)
+    hidexdecorations!(ax_SIF)
 
     fig
     return fig

src/standalone/Vegetation/Canopy.jl

@@ -31,6 +31,7 @@ using .PlantHydraulics
 include("./stomatalconductance.jl")
 include("./photosynthesis.jl")
 include("./radiation.jl")
+include("./solar_induced_fluorescence.jl")
 include("./pfts.jl")
 include("./canopy_energy.jl")
 include("./canopy_parameterizations.jl")
@@ -53,7 +54,7 @@ struct SharedCanopyParameters{FT <: AbstractFloat, PSE}
 end
 
 """
-     CanopyModel{FT, AR, RM, PM, SM, PHM, EM, A, R, S, PS, D} <: AbstractExpModel{FT}
+     CanopyModel{FT, AR, RM, PM, SM, PHM, EM, SM, A, R, S, PS, D} <: AbstractExpModel{FT}
 
 The model struct for the canopy, which contains
 - the canopy model domain (a point for site-level simulations, or
@@ -73,6 +74,8 @@ prognostically solves Richards equation in the plant is available.
  `AbstractCanopyEnergyModel` and currently we support a version where
   the canopy temperature is prescribed, and one where it is solved for
   prognostically.
+- subcomponent model type for canopy SIF.
+  prognostically.
 - canopy model parameters, which include parameters that are shared
 between canopy model components or those needed to compute boundary
 fluxes.
@@ -94,7 +97,7 @@ treated differently.
 
 $(DocStringExtensions.FIELDS)
 """
-struct CanopyModel{FT, AR, RM, PM, SM, PHM, EM, A, R, S, PS, D} <:
+struct CanopyModel{FT, AR, RM, PM, SM, PHM, EM, SIFM, A, R, S, PS, D} <:
        AbstractExpModel{FT}
     "Autotrophic respiration model, a canopy component model"
     autotrophic_respiration::AR
@@ -108,6 +111,8 @@ struct CanopyModel{FT, AR, RM, PM, SM, PHM, EM, A, R, S, PS, D} <:
     hydraulics::PHM
     "Energy balance model, a canopy component model"
     energy::EM
+    "SIF model, a canopy component model"
+    sif::SIFM
     "Atmospheric forcing: prescribed or coupled"
     atmos::A
     "Radiative forcing: prescribed or coupled"
@@ -128,6 +133,7 @@ end
         conductance::AbstractStomatalConductanceModel{FT},
         hydraulics::AbstractPlantHydraulicsModel{FT},
         energy::AbstractCanopyEnergyModel{FT},
+        sif::AbstractSIFModel{FT},
         atmos::AbstractAtmosphericDrivers{FT},
         radiation::AbstractRadiativeDrivers{FT},
         soil::AbstractSoilDriver,
@@ -153,6 +159,7 @@ function CanopyModel{FT}(;
     conductance::AbstractStomatalConductanceModel{FT},
     hydraulics::AbstractPlantHydraulicsModel{FT},
     energy = PrescribedCanopyTempModel{FT}(),
+    sif = Lee2015SIFModel{FT}(),
     atmos::AbstractAtmosphericDrivers{FT},
     radiation::AbstractRadiativeDrivers{FT},
     soil_driver::AbstractSoilDriver,
@@ -175,6 +182,7 @@ function CanopyModel{FT}(;
         conductance,
         hydraulics,
         energy,
+        sif,
         atmos,
         radiation,
         soil_driver,
@@ -203,6 +211,7 @@ canopy_components(::CanopyModel) = (
     :radiative_transfer,
     :autotrophic_respiration,
     :energy,
+    :sif,
 )
 
 """
@@ -429,6 +438,7 @@ function ClimaLand.make_update_aux(
         grav = FT(LP.grav(earth_param_set))
         ρ_l = FT(LP.ρ_cloud_liq(earth_param_set))
         R = FT(LP.gas_constant(earth_param_set))
+        T_freeze = FT(LP.T_freeze(earth_param_set))
         thermo_params = earth_param_set.thermo_params
         (; G_Function, Ω, λ_γ_PAR, λ_γ_NIR) =
             canopy.radiative_transfer.parameters
@@ -555,6 +565,18 @@ function ClimaLand.make_update_aux(
             c_co2_air,
             R,
         )
+        # update SIF
+        SIF = p.canopy.sif.SIF
+        update_SIF!(
+            SIF,
+            canopy.sif,
+            p.canopy.radiative_transfer.par.abs,
+            T_canopy,
+            Vcmax25,
+            R,
+            T_freeze,
+            canopy.photosynthesis.parameters,
+        )
         @. GPP = compute_GPP(An, K, LAI, Ω)
         @. gs = medlyn_conductance(g0, Drel, medlyn_factor, An, c_co2_air)
         # update autotrophic respiration

src/standalone/Vegetation/solar_induced_fluorescence.jl

@@ -0,0 +1,117 @@
+export SIFParameters, Lee2015SIFModel
+
+abstract type AbstractSIFModel{FT} <: AbstractCanopyComponent{FT} end
+
+"""
+    SIFParameters{FT<:AbstractFloat}
+
+The required parameters for the SIF parameterisation 
+Lee et al, 2015. Global Change Biology 21, 3469-3477, doi:10.1111/gcb.12948.
+$(DocStringExtensions.FIELDS)
+"""
+@kwdef struct SIFParameters{FT <: AbstractFloat}
+    "The rate coefficient for florescence, unitless"
+    kf::FT = FT(0.05)
+    "Parameter used to compute the rate coefficient for heat loss in dark-adapted conditions, Tol et al. 2014, unitless"
+    kd_p1::FT = FT(0.03)
+    "Parameter used to compute the rate coefficient for heat loss in dark-adapted conditions, Tol et al. 2014, unitless"
+    kd_p2::FT = FT(0.0273)
+    "Parameter used to compute the rate coefficient for heat loss in dark-adapted conditions, Tol et al. 2014, unitless"
+    min_kd::FT = FT(0.087)
+    "Parameter used to compute the rate coefficient for heat loss in light-adapted conditions, Lee et al 2013 (unitless)"
+    kn_p1::FT = FT(6.2473)
+    "Parameter used to compute the rate coefficient for heat loss in light-adapted conditions, Lee et al 2013 (unitless)"
+    kn_p2::FT = FT(0.5944)
+    "Rate coefficient for photochemical quenching"
+    kp::FT = FT(4.0)
+    "Slope of line relating leaf-level fluorescence to spectrometer-observed fluorescence as a function of Vcmax 25. Lee et al 2015."
+    kappa_p1::FT = FT(0.045)
+    "Intercept of line relating leaf-level fluorescence to spectrometer-observed fluorescence as a function of Vcmax 25.  Lee et al 2015."
+    kappa_p2::FT = FT(7.85)
+end
+
+Base.eltype(::SIFParameters{FT}) where {FT} = FT
+
+struct Lee2015SIFModel{FT, SP <: SIFParameters{FT}} <: AbstractSIFModel{FT}
+    parameters::SP
+    function Lee2015SIFModel{FT}() where {FT}
+        parameters = SIFParameters{FT}()
+        new{FT, typeof(parameters)}(parameters)
+    end
+end
+
+ClimaLand.name(model::AbstractSIFModel) = :sif
+ClimaLand.auxiliary_vars(model::Lee2015SIFModel) = (:SIF,)
+ClimaLand.auxiliary_types(model::Lee2015SIFModel{FT}) where {FT} = (FT,)
+ClimaLand.auxiliary_domain_names(::Lee2015SIFModel) = (:surface,)
+
+
+# 4 Solar Induced Fluorescence (SIF)
+
+# call function below inside photosynthesis.jl p
+
+"""
+    update_SIF!(
+        SIF::FT,
+        sif_model::Lee2015SIFModel,
+        APAR::FT,
+        Tc::FT,
+        Vcmax25::FT,
+        R::FT,
+        T_freeze::FT,
+        photosynthesis_parameters,
+    ) where {FT}
+
+Updates observed SIF at 755 nm in W/m^2. Note that Tc is in Kelvin, and photo
+synthetic rates are in mol/m^2/s, and APAR is in PPFD.
+Lee et al, 2015. Global Change Biology 21, 3469-3477, doi:10.1111/gcb.12948
+"""
+function update_SIF!(
+    SIF,
+    sif_model::Lee2015SIFModel,
+    APAR,
+    Tc,
+    Vcmax25,
+    R,
+    T_freeze,
+    photosynthesis_parameters,
+)
+    sif_parameters = sif_model.parameters
+    @. SIF = compute_SIF_at_a_point(
+        APAR,
+        Tc,
+        Vcmax25,
+        R,
+        T_freeze,
+        photosynthesis_parameters,
+        sif_parameters,
+    )
+end
+
+Base.broadcastable(m::SIFParameters) = tuple(m)
+function compute_SIF_at_a_point(
+    APAR::FT,
+    Tc::FT,
+    Vcmax25::FT,
+    R::FT,
+    T_freeze::FT,
+    photosynthesis_parameters,
+    sif_parameters,
+) where {FT}
+    (; ΔHJmax, To, θj, ϕ) = photosynthesis_parameters
+    Jmax = max_electron_transport(Vcmax25, ΔHJmax, Tc, To, R)
+    J = electron_transport(APAR, Jmax, θj, ϕ)
+    (; kf, kd_p1, kd_p2, min_kd, kn_p1, kn_p2, kp, kappa_p1, kappa_p2) =
+        sif_parameters
+    kd = max(kd_p1 * (Tc - T_freeze) + kd_p2, min_kd)
+    x = 1 - J / Jmax
+    kn = (kn_p1 * x - kn_p2) * x
+    ϕp0 = kp / (kf + kp + kn)
+    ϕp = J / Jmax * ϕp0
+    ϕf = kf / (kf + kd + kn) * (1 - ϕp)
+    κ = kappa_p1 * Vcmax25 * FT(1e6) + kappa_p2 # formula expects Vcmax25 in μmol/m^2/s
+    F = APAR * ϕf
+    SIF_755 = F / κ
+
+    return SIF_755
+end