diff --git a/docs/src/APIs/canopy/Photosynthesis.md b/docs/src/APIs/canopy/Photosynthesis.md index 6ac46fee08..6361344835 100644 --- a/docs/src/APIs/canopy/Photosynthesis.md +++ b/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 -``` \ No newline at end of file +``` diff --git a/lib/ClimaLandSimulations/src/utilities/climaland_output_dataframe.jl b/lib/ClimaLandSimulations/src/utilities/climaland_output_dataframe.jl index 3db7826318..26323a00a8 100644 --- a/lib/ClimaLandSimulations/src/utilities/climaland_output_dataframe.jl +++ b/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] diff --git a/lib/ClimaLandSimulations/src/utilities/makie_plots.jl b/lib/ClimaLandSimulations/src/utilities/makie_plots.jl index 5036b64d47..d9cf9910ef 100644 --- a/lib/ClimaLandSimulations/src/utilities/makie_plots.jl +++ b/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 diff --git a/src/standalone/Vegetation/Canopy.jl b/src/standalone/Vegetation/Canopy.jl index b9d5c3a5db..41dd15ea71 100644 --- a/src/standalone/Vegetation/Canopy.jl +++ b/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 diff --git a/src/standalone/Vegetation/solar_induced_fluorescence.jl b/src/standalone/Vegetation/solar_induced_fluorescence.jl new file mode 100644 index 0000000000..f2ca8e6df7 --- /dev/null +++ b/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