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/ext/CreateParametersExt.jl b/ext/CreateParametersExt.jl
index 514495916b..77553f39ec 100644
--- a/ext/CreateParametersExt.jl
+++ b/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
 
 
diff --git a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Ha1/US-Ha1_parameters.jl b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Ha1/US-Ha1_parameters.jl
index a7fc8a915e..2a474c98e2 100644
--- a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Ha1/US-Ha1_parameters.jl
+++ b/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
 
diff --git a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-MOz/US-MOz_parameters.jl b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-MOz/US-MOz_parameters.jl
index db9110af70..5c308cf482 100644
--- a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-MOz/US-MOz_parameters.jl
+++ b/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
 
diff --git a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-NR1/US-NR1_parameters.jl b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-NR1/US-NR1_parameters.jl
index 311552a439..ed031790b2 100644
--- a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-NR1/US-NR1_parameters.jl
+++ b/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
 
diff --git a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Var/US-Var_parameters.jl b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Var/US-Var_parameters.jl
index 7228ccab2a..f572aaa25c 100644
--- a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Var/US-Var_parameters.jl
+++ b/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
 
diff --git a/lib/ClimaLandSimulations/src/utilities/climaland_output_dataframe.jl b/lib/ClimaLandSimulations/src/utilities/climaland_output_dataframe.jl
index 3db7826318..bc6d2fa8b4 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, :photosynthesis, :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..4f828b2086 100644
--- a/src/standalone/Vegetation/Canopy.jl
+++ b/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,
diff --git a/src/standalone/Vegetation/canopy_parameterizations.jl b/src/standalone/Vegetation/canopy_parameterizations.jl
index 6522c50d2f..1b44386b81 100644
--- a/src/standalone/Vegetation/canopy_parameterizations.jl
+++ b/src/standalone/Vegetation/canopy_parameterizations.jl
@@ -1041,3 +1041,47 @@ 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(
+        APAR::FT,
+        Tc::FT,
+        Vcmax25::FT,
+        R::FT,
+        photosynthesis_parameters,
+    ) 
+
+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.
+Lee et al, 2015. Global Change Biology 21, 3469-3477, doi:10.1111/gcb.12948
+"""
+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
diff --git a/src/standalone/Vegetation/optimality_farquhar.jl b/src/standalone/Vegetation/optimality_farquhar.jl
index a2212269b7..b4e6ca3728 100644
--- a/src/standalone/Vegetation/optimality_farquhar.jl
+++ b/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)
diff --git a/src/standalone/Vegetation/photosynthesis.jl b/src/standalone/Vegetation/photosynthesis.jl
index c9282c423d..b86f1f83f4 100644
--- a/src/standalone/Vegetation/photosynthesis.jl
+++ b/src/standalone/Vegetation/photosynthesis.jl
@@ -1,4 +1,34 @@
-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 
+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
 
 abstract type AbstractPhotosynthesisMechanism end
 """
@@ -17,7 +47,7 @@ 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)
@@ -25,6 +55,7 @@ $(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 +93,8 @@ Base.@kwdef struct FarquharParameters{
     pc::FT
     "Photosynthesis mechanism: C3 or C4"
     mechanism::MECH
+    "Parameters of the sif model, Lee et al 2015"
+    sif_parameters::SP
 end
 
 Base.eltype(::FarquharParameters{FT}) where {FT} = FT
@@ -78,12 +111,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 +159,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 +181,7 @@ function update_photosynthesis!(
     Rd,
     An,
     Vcmax25field,
+    SIF,
     model::FarquharModel,
     T,
     APAR,
@@ -170,8 +204,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")
diff --git a/test/standalone/Vegetation/test_bigleaf_parameterizations.jl b/test/standalone/Vegetation/test_bigleaf_parameterizations.jl
index 5859cab915..472516fcea 100644
--- a/test/standalone/Vegetation/test_bigleaf_parameterizations.jl
+++ b/test/standalone/Vegetation/test_bigleaf_parameterizations.jl
@@ -41,13 +41,15 @@ for FT in (Float32, Float64)
         params = OptimalityFarquharParameters(FT)
         @test params.mechanism == C3()
         model = OptimalityFarquharModel(params)
-        @test ClimaLand.auxiliary_vars(model) == (:An, :GPP, :Rd, :Vcmax25)
-        @test ClimaLand.auxiliary_types(model) == (FT, FT, FT, FT)
+        @test ClimaLand.auxiliary_vars(model) ==
+              (:An, :GPP, :Rd, :Vcmax25, :SIF)
+        @test ClimaLand.auxiliary_types(model) == (FT, FT, FT, FT, FT)
         @test ClimaLand.auxiliary_domain_names(model) ==
-              (:surface, :surface, :surface, :surface)
+              (:surface, :surface, :surface, :surface, :surface)
         Rd = zeros(FT, 1)
         An = similar(Rd)
         Vcmax25 = similar(An)
+        SIF = similar(An)
         T = FT(280)
         β = FT(1)
         medlyn_factor = FT(10.0)
@@ -57,6 +59,7 @@ for FT in (Float32, Float64)
             Rd,
             An,
             Vcmax25,
+            SIF,
             model,
             T,
             APAR,
@@ -68,6 +71,7 @@ for FT in (Float32, Float64)
         @test Rd[1] != 0.0
         @test Vcmax25[1] != 0.0
         @test An[1] != 0.0
+        @test SIF[1] != 0.0
     end
 
     @testset "Big Leaf Parameterizations, FT = $FT" begin