add simple SIF model to photosynthesis

ext/CreateParametersExt.jl

@@ -6,7 +6,6 @@ import Insolation.Parameters.InsolationParameters
 import SurfaceFluxes.Parameters.SurfaceFluxesParameters
 import SurfaceFluxes.UniversalFunctions as UF
 import ClimaParams as CP
-
 import ClimaLand
 import ClimaLand.Soil
 # Parameter structs
@@ -15,6 +14,7 @@ import ClimaLand.Soil.EnergyHydrologyParameters
 import ClimaLand.Canopy.AutotrophicRespirationParameters
 import ClimaLand.Canopy.FarquharParameters
 import ClimaLand.Canopy.OptimalityFarquharParameters
+import ClimaLand.Canopy.SIFParameters
 import ClimaLand.Canopy.MedlynConductanceParameters
 import ClimaLand.Canopy.BeerLambertParameters
 import ClimaLand.Canopy.TwoStreamParameters
@@ -133,16 +133,24 @@ toml_dict = CP.create_toml_dict(Float32);
 ClimaLand.Canopy.FarquharParameters(toml_dict, ClimaLand.Canopy.C3(); Vcmax25 = 99999999, pc = 444444444)
 ```
 """
-FarquharParameters(
+function FarquharParameters(
     ::Type{FT},
     mechanism;
     kwargs...,
-) where {FT <: AbstractFloat} =
-    FarquharParameters(CP.create_toml_dict(FT), mechanism; kwargs...)
+) where {FT <: AbstractFloat}
+    sif_parameters = SIFParameters{FT}()
+    FarquharParameters(
+        CP.create_toml_dict(FT),
+        mechanism,
+        sif_parameters;
+        kwargs...,
+    )
+end
 
 function FarquharParameters(
     toml_dict::CP.AbstractTOMLDict,
-    mechanism;
+    mechanism,
+    sif_parameters;
     Vcmax25 = 5e-5,
     kwargs...,
 )
@@ -167,9 +175,11 @@ function FarquharParameters(
     parameters = CP.get_parameter_values(toml_dict, name_map, "Land")
     FT = CP.float_type(toml_dict)
     MECH = typeof(mechanism)
-    return FarquharParameters{FT, MECH}(;
+    SP = typeof(sif_parameters)
+    return FarquharParameters{FT, MECH, SP}(;
         mechanism,
         Vcmax25,
+        sif_parameters,
         parameters...,
         kwargs...,
     )
@@ -199,13 +209,19 @@ toml_dict = CP.create_toml_dict(Float32);
 ClimaLand.Canopy.OptimalityFarquharParameters(toml_dict; pc = 444444444)
 ```
 """
-OptimalityFarquharParameters(
+function OptimalityFarquharParameters(
     ::Type{FT};
     kwargs...,
-) where {FT <: AbstractFloat} =
-    OptimalityFarquharParameters(CP.create_toml_dict(FT); kwargs...)
+) where {FT <: AbstractFloat}
+    sif_parameters = SIFParameters{FT}()
+    OptimalityFarquharParameters(
+        CP.create_toml_dict(FT),
+        sif_parameters;
+        kwargs...,
+    )
+end
 
-function OptimalityFarquharParameters(toml_dict; kwargs...)
+function OptimalityFarquharParameters(toml_dict, sif_parameters; kwargs...)
     name_map = (;
         :Jmax_activation_energy => :ΔHJmax,
         :intercellular_O2_concentration => :oi,
@@ -229,7 +245,12 @@ function OptimalityFarquharParameters(toml_dict; kwargs...)
     params = CP.get_parameter_values(toml_dict, name_map, "Land")
     FT = CP.float_type(toml_dict)
     mechanism = ClimaLand.Canopy.C3()
-    return OptimalityFarquharParameters{FT}(; params..., kwargs..., mechanism)
+    return OptimalityFarquharParameters{FT, typeof(sif_parameters)}(;
+        sif_parameters,
+        params...,
+        kwargs...,
+        mechanism,
+    )
 end
 
 

lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Ha1/US-Ha1_parameters.jl

@@ -158,6 +158,7 @@ function photosynthesis_harvard(;
     ΔHΓstar = FT(37830),
     ΔHJmax = FT(43540),
     ΔHRd = FT(46390),
+    sif_parameters = SIFParameters{FT}(),
 )
     return FarquharParameters(
         Vcmax25,
@@ -178,6 +179,7 @@ function photosynthesis_harvard(;
         sc,
         pc,
         mechanism,
+        sif_parameters,
     )
 end
 

lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-MOz/US-MOz_parameters.jl

@@ -158,6 +158,7 @@ function photosynthesis_ozark(;
     ΔHΓstar = FT(37830),
     ΔHJmax = FT(43540),
     ΔHRd = FT(46390),
+    sif_parameters = SIFParameters{FT}(),
 )
     return FarquharParameters(
         Vcmax25,
@@ -178,6 +179,7 @@ function photosynthesis_ozark(;
         sc,
         pc,
         mechanism,
+        sif_parameters,
     )
 end
 

lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-NR1/US-NR1_parameters.jl

@@ -158,6 +158,7 @@ function photosynthesis_niwotridge(;
     ΔHΓstar = FT(37830),
     ΔHJmax = FT(43540),
     ΔHRd = FT(46390),
+    sif_parameters = SIFParameters{FT}(),
 )
     return FarquharParameters(
         Vcmax25,
@@ -178,6 +179,7 @@ function photosynthesis_niwotridge(;
         sc,
         pc,
         mechanism,
+        sif_parameters,
     )
 end
 

lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Var/US-Var_parameters.jl

@@ -158,6 +158,7 @@ function photosynthesis_vairaranch(;
     ΔHΓstar = FT(37830),
     ΔHJmax = FT(43540),
     ΔHRd = FT(46390),
+    sif_parameters = SIFParameters{FT}(),
 )
     return FarquharParameters(
         Vcmax25,
@@ -178,6 +179,7 @@ function photosynthesis_vairaranch(;
         sc,
         pc,
         mechanism,
+        sif_parameters,
     )
 end
 

src/standalone/Vegetation/Canopy.jl

@@ -543,10 +543,12 @@ function ClimaLand.make_update_aux(
 
         # Update Rd, An, Vcmax25 (if applicable to model) in place
         Vcmax25 = p.canopy.photosynthesis.Vcmax25
+        SIF = p.canopy.photosynthesis.SIF
         update_photosynthesis!(
             Rd,
             An,
             Vcmax25,
+            SIF,
             canopy.photosynthesis,
             T_canopy,
             p.canopy.radiative_transfer.par.abs,

src/standalone/Vegetation/canopy_parameterizations.jl

@@ -1041,3 +1041,40 @@ function plant_respiration_growth(f2::FT, An::FT, Rpm::FT) where {FT}
     Rg = f2 * (An - Rpm)
     return Rg
 end
+
+# 4 Solar Induced Fluorescence (SIF)
+
+# call function below inside photosynthesis.jl p
+
+"""
+    compute_SIF_at_a_point(Tc::FT,APAR::FT, Vcmax25::FT)
+
+Computes 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.
+"""
+function compute_SIF_at_a_point(
+    APAR::FT,
+    Tc::FT,
+    Vcmax25::FT,
+    R::FT,
+    photosynthesis_parameters,
+) where {FT}
+
+    (; ΔHJmax, To, θj, ϕ, sif_parameters) = 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
+    Tf = FT(273.15)
+    kd = max(kd_p1 * (Tc - Tf) + 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

src/standalone/Vegetation/optimality_farquhar.jl

@@ -8,7 +8,7 @@ The required parameters for the optimality Farquhar photosynthesis model.
 Currently, only C3 photosynthesis is supported.
 $(DocStringExtensions.FIELDS)
 """
-Base.@kwdef struct OptimalityFarquharParameters{FT <: AbstractFloat}
+Base.@kwdef struct OptimalityFarquharParameters{FT <: AbstractFloat, SP}
     "Photosynthesis mechanism: C3 only"
     mechanism::C3
     "Γstar at 25 °C (mol/mol)"
@@ -45,6 +45,8 @@ Base.@kwdef struct OptimalityFarquharParameters{FT <: AbstractFloat}
     pc::FT
     "Constant describing cost of maintaining electron transport (unitless)"
     c::FT
+    "SIF Parameters"
+    sif_parameters::SP
 end
 
 Base.eltype(::OptimalityFarquharParameters{FT}) where {FT} = FT
@@ -72,14 +74,14 @@ function OptimalityFarquharModel{FT}(
 end
 
 ClimaLand.auxiliary_vars(model::OptimalityFarquharModel) =
-    (:An, :GPP, :Rd, :Vcmax25)
+    (:An, :GPP, :Rd, :Vcmax25, :SIF)
 ClimaLand.auxiliary_types(model::OptimalityFarquharModel{FT}) where {FT} =
-    (FT, FT, FT, FT)
+    (FT, FT, FT, FT, FT)
 ClimaLand.auxiliary_domain_names(::OptimalityFarquharModel) =
-    (:surface, :surface, :surface, :surface)
+    (:surface, :surface, :surface, :surface, :surface)
 
 """
-    update_photosynthesis!(Rd, An, Vcmax25,
+    update_photosynthesis!(Rd, An, Vcmax25,SIF,
         model::OptimalityFarquharModel,
         T,
         APAR,
@@ -101,6 +103,7 @@ function update_photosynthesis!(
     Rd,
     An,
     Vcmax25,
+    SIF,
     model::OptimalityFarquharModel,
     T,
     APAR,
@@ -152,4 +155,6 @@ function update_photosynthesis!(
     @. Vcmax25 = Vcmax / arrhenius_function(T, To, R, ΔHVcmax)
     @. Rd = dark_respiration(Vcmax25, β, f, ΔHRd, T, To, R)
     @. An = net_photosynthesis(Ac, Aj, Rd, β)
+    @. SIF = compute_SIF_at_a_point(APAR, T, Vcmax25, R, model.parameters)
 end
+Base.broadcastable(m::OptimalityFarquharParameters) = tuple(m)

src/standalone/Vegetation/photosynthesis.jl

@@ -1,4 +1,33 @@
-export FarquharParameters, FarquharModel, C3, C4
+export SIFParameters, FarquharParameters, FarquharModel, C3, C4
+
+abstract type AbstractPhotosynthesisModel{FT} <: AbstractCanopyComponent{FT} end
+
+"""
+    SIFParameters{FT<:AbstractFloat}
+
+The required parameters for the SIF parameterisation [give citation].
+$(DocStringExtensions.FIELDS)
+"""
+@kwdef struct SIFParameters{FT <: AbstractFloat}
+    ""
+    kf::FT = FT(0.05)
+    ""
+    kd_p1::FT = FT(0.03)
+    ""
+    kd_p2::FT = FT(0.0273)
+    ""
+    min_kd::FT = FT(0.087)
+    ""
+    kn_p1::FT = FT(6.2473)
+    ""
+    kn_p2::FT = FT(0.5944)
+    ""
+    kp::FT = FT(4.0)
+    ""
+    kappa_p1::FT = FT(0.045)
+    ""
+    kappa_p2::FT = FT(7.85)
+end
 
 abstract type AbstractPhotosynthesisMechanism end
 """
@@ -17,14 +46,15 @@ struct C4 <: AbstractPhotosynthesisMechanism end
 
 abstract type AbstractPhotosynthesisModel{FT} <: AbstractCanopyComponent{FT} end
 """
-    FarquharParameters{FT<:AbstractFloat, MECH <: AbstractPhotosynthesisMechanism}
+    FarquharParameters{FT<:AbstractFloat, MECH <: AbstractPhotosynthesisMechanism, SP <: SIFParameters}
 
 The required parameters for the Farquhar photosynthesis model.
 $(DocStringExtensions.FIELDS)
 """
 Base.@kwdef struct FarquharParameters{
     FT <: AbstractFloat,
     MECH <: AbstractPhotosynthesisMechanism,
+    SP <: SIFParameters,
 }
     "Vcmax at 25 °C (mol CO2/m^2/s)"
     Vcmax25::FT
@@ -62,6 +92,8 @@ Base.@kwdef struct FarquharParameters{
     pc::FT
     "Photosynthesis mechanism: C3 or C4"
     mechanism::MECH
+    ""
+    sif_parameters::SP
 end
 
 Base.eltype(::FarquharParameters{FT}) where {FT} = FT
@@ -78,12 +110,12 @@ function FarquharModel{FT}(
 end
 
 ClimaLand.name(model::AbstractPhotosynthesisModel) = :photosynthesis
-ClimaLand.auxiliary_vars(model::FarquharModel) = (:An, :GPP, :Rd, :Vcmax25)
+ClimaLand.auxiliary_vars(model::FarquharModel) =
+    (:An, :GPP, :Rd, :Vcmax25, :SIF)
 ClimaLand.auxiliary_types(model::FarquharModel{FT}) where {FT} =
-    (FT, FT, FT, FT)
+    (FT, FT, FT, FT, FT)
 ClimaLand.auxiliary_domain_names(::FarquharModel) =
-    (:surface, :surface, :surface, :surface)
-
+    (:surface, :surface, :surface, :surface, :surface)
 
 function photosynthesis_at_a_point_Farquhar(
     T,
@@ -126,7 +158,7 @@ function photosynthesis_at_a_point_Farquhar(
 end
 
 """
-    update_photosynthesis!(Rd, An, Vcmax25,
+    update_photosynthesis!(Rd, An, Vcmax25, SIF,
         model::FarquharModel,
         T,
         APAR,
@@ -148,6 +180,7 @@ function update_photosynthesis!(
     Rd,
     An,
     Vcmax25field,
+    SIF,
     model::FarquharModel,
     T,
     APAR,
@@ -170,8 +203,11 @@ function update_photosynthesis!(
         model.parameters,
     )
     Vcmax25field .= Vcmax25
+    @. SIF = compute_SIF_at_a_point(APAR, T, Vcmax25, R, model.parameters)
+
 end
 Base.broadcastable(m::AbstractPhotosynthesisMechanism) = tuple(m)
 Base.broadcastable(m::FarquharParameters) = tuple(m)
+Base.broadcastable(m::SIFParameters) = tuple(m)
 
 include("./optimality_farquhar.jl")