Add a AbstractSoilMoistureStressModel, 3 methods

This commits adds a portion of Yuchen Li PR #1330,
a new Canopy component: AbstractSoilMoistureStressModel,
with 3 methods: the existing one (Tuzet) which was previously
in PlantHydraulics, a Piecewise method using dry to wet soil moisture,
and no moisture stress.

Co-authored-by: AlexisRenchon <[email protected]>
Co-authored-by: ychnli <[email protected]>

[skip ci]

Author: AlexisRenchon

src/integrated/land.jl

@@ -184,6 +184,17 @@ function LandModel{FT}(;
         energy_model = PrescribedCanopyTempModel{FT}()
     end
 
+    if :soil_moisture_stress in propertynames(canopy_component_types)
+        soil_moisture_stress_model =
+        canopy_component_types.soil_moisture_stress(
+                                                    canopy_component_args.soil_moisture_stress...,
+                                                   )
+    else
+        @info "No soil moisture stress model provided, using NoMoistureStressModel (βm = 1.0)"
+        soil_moisture_stress_model =
+        ClimaLand.Canopy.NoMoistureStressModel{FT}()
+    end
+
     canopy = Canopy.CanopyModel{FT}(;
         autotrophic_respiration = canopy_component_types.autotrophic_respiration(
             canopy_component_args.autotrophic_respiration...,
@@ -197,6 +208,7 @@ function LandModel{FT}(;
         conductance = canopy_component_types.conductance(
             canopy_component_args.conductance...,
         ),
+        soil_moisture_stress = soil_moisture_stress_model,
         hydraulics = canopy_component_types.hydraulics(;
             transpiration = transpiration,
             canopy_component_args.hydraulics...,

src/standalone/Vegetation/Canopy.jl

@@ -40,6 +40,7 @@ export SharedCanopyParameters, CanopyModel, set_canopy_prescribed_field!
 include("./component_models.jl")
 include("./PlantHydraulics.jl")
 using .PlantHydraulics
+include("./soil_moisture_stress.jl")
 include("./stomatalconductance.jl")
 include("./photosynthesis.jl")
 include("./photosynthesis_farquhar.jl")
@@ -137,8 +138,6 @@ end
         ϕa1 = FT(0.022),
         ϕa2 = FT(-0.00034),
         α = FT(0.933),
-        sc = LP.get_default_parameter(FT, :low_water_pressure_sensitivity),
-        pc = LP.get_default_parameter(FT, :moisture_stress_ref_water_pressure),
     ) where {FT <: AbstractFloat}
 
 Constructs a P-model (an optimality model for photosynthesis) using default parameters.
@@ -153,10 +152,7 @@ The following default parameters are used:
 - ϕa1 = 0.022 (K^-1) - first order term in quadratic intrinsic quantum yield (Stocker 2020)
 - ϕa2 = -0.00034 (K^-2) - second order term in quadratic intrinsic quantum yield (Stocker 2020)
 - α = 0.933 (unitless) - 1 - 1/T where T is the timescale of Vcmax, Jmax acclimation. Here T = 15 days. (Mengoli 2022)
-- sc = 5e-6 (Pa^{-1}) - sensitivity to low water pressure in the moisture stress factor [Tuzet et al. (2003)]
-- pc = -2e6 (Pa) - reference water pressure for the moisture stress factor [Tuzet et al. (2003)]
 """
-
 function PModel{FT}(;
     cstar = FT(0.41),
     β = FT(146),
@@ -166,8 +162,6 @@ function PModel{FT}(;
     ϕa1 = FT(0.022),
     ϕa2 = FT(-0.00034),
     α = FT(0.933),
-    sc = LP.get_default_parameter(FT, :low_water_pressure_sensitivity),
-    pc = LP.get_default_parameter(FT, :moisture_stress_ref_water_pressure),
 ) where {FT <: AbstractFloat}
     parameters = ClimaLand.Canopy.PModelParameters(
         cstar = cstar,
@@ -178,8 +172,6 @@ function PModel{FT}(;
         ϕa1 = ϕa1,
         ϕa2 = ϕa2,
         α = α,
-        sc = sc,
-        pc = pc,
     )
 
     return PModel{FT}(parameters)
@@ -402,6 +394,22 @@ function PModelConductance{FT}(; Drel = FT(1.6)) where {FT <: AbstractFloat}
     return PModelConductance{FT}(cond_params)
 end
 
+"""
+    TuzetMoistureStressModel{FT}() where {FT <: AbstractFloat}
+Creates a TuzetMoistureStressModel using default parameters of type FT.
+The parameters are set to:
+- sc = 5e-6 (Pa^{-1}) - sensitivity to low water pressure
+- pc = -2e6 (Pa) - reference water pressure for the moisture stress factor
+"""
+function TuzetMoistureStressModel{FT}(;
+    sc = LP.get_default_parameter(FT, :low_water_pressure_sensitivity),
+    pc = LP.get_default_parameter(FT, :moisture_stress_ref_water_pressure),
+) where {FT <: AbstractFloat}
+    parameters = TuzetMoistureStressParameters{FT}(sc, pc)
+    return TuzetMoistureStressModel{eltype(parameters), typeof(parameters)}(
+        parameters,
+    )
+end
 
 ########################################################
 # End component model convenience constructors
@@ -423,19 +431,23 @@ struct SharedCanopyParameters{FT <: AbstractFloat, PSE}
 end
 
 """
-     CanopyModel{FT, AR, RM, PM, SM, PHM, EM, SM, A, R, S, PS, D} <: ClimaLand.AbstractImExModel{FT}
+     CanopyModel{FT, AR, RM, PM, SM, PHM, EM, SM, SMSM, A, R, S, PS, D} <: ClimaLand.AbstractImExModel{FT}
 
 The model struct for the canopy, which contains
 - the canopy model domain (a point for site-level simulations, or
 an extended surface (plane/spherical surface) for regional or global simulations.
 - subcomponent model type for radiative transfer. This is of type
 `AbstractRadiationModel`.
 - subcomponent model type for photosynthesis. This is of type
-`AbstractPhotosynthesisModel`, and currently only the `FarquharModel`
-is supported.
+`AbstractPhotosynthesisModel` and supports `FarquharModel`, and `PModel`.
 - subcomponent model type for stomatal conductance. This is of type
- `AbstractStomatalConductanceModel` and currently only the `MedlynModel`
-is supported
+ `AbstractStomatalConductanceModel` and supports `MedlynConductanceModel` and
+ `PModelConductance`. Note if `PModel` is used for photosynthesis, then you
+ must also use `PModelConductance` for stomatal conductance, since these two models
+ are derived from the same set of conditions.
+- subcomponent model type for soil moisture stress. This is of type
+ `AbstractSoilMoistureStressModel`. Currently we support `TuzetMoistureStressModel` (default)
+ `PiecewiseMoistureStressModel`, and `NoMoistureStressModel` (stress factor = 1).
 - subcomponent model type for plant hydraulics. This is of type
  `AbstractPlantHydraulicsModel` and currently only a version which
 prognostically solves Richards equation in the plant is available.
@@ -465,7 +477,7 @@ treated differently.
 
 $(DocStringExtensions.FIELDS)
 """
-struct CanopyModel{FT, AR, RM, PM, SM, PHM, EM, SIFM, B, PS, D} <:
+struct CanopyModel{FT, AR, RM, PM, SM, SMSM, PHM, EM, SIFM, B, PS, D} <:
        ClimaLand.AbstractImExModel{FT}
     "Autotrophic respiration model, a canopy component model"
     autotrophic_respiration::AR
@@ -475,6 +487,8 @@ struct CanopyModel{FT, AR, RM, PM, SM, PHM, EM, SIFM, B, PS, D} <:
     photosynthesis::PM
     "Stomatal conductance model, a canopy component model"
     conductance::SM
+    "Soil moisture stress parameterization, a canopy component model"
+    soil_moisture_stress::SMSM
     "Plant hydraulics model, a canopy component model"
     hydraulics::PHM
     "Energy balance model, a canopy component model"
@@ -495,6 +509,7 @@ end
         radiative_transfer::AbstractRadiationModel{FT},
         photosynthesis::AbstractPhotosynthesisModel{FT},
         conductance::AbstractStomatalConductanceModel{FT},
+        soil_moisture_stress::AbstractSoilMoistureStressModel{FT} = TuzetMoistureStressModel{FT}(),
         hydraulics::AbstractPlantHydraulicsModel{FT},
         energy::AbstractCanopyEnergyModel{FT},
         sif::AbstractSIFModel{FT},
@@ -519,6 +534,9 @@ function CanopyModel{FT}(;
     radiative_transfer::AbstractRadiationModel{FT},
     photosynthesis::AbstractPhotosynthesisModel{FT},
     conductance::AbstractStomatalConductanceModel{FT},
+    soil_moisture_stress::AbstractSoilMoistureStressModel{FT} = TuzetMoistureStressModel{
+                                                                                         FT,
+                                                                                        }(),
     hydraulics::AbstractPlantHydraulicsModel{FT},
     energy = PrescribedCanopyTempModel{FT}(),
     sif = Lee2015SIFModel{FT}(),
@@ -548,6 +566,7 @@ function CanopyModel{FT}(;
         radiative_transfer,
         photosynthesis,
         conductance,
+        soil_moisture_stress,
         hydraulics,
         energy,
         sif,
@@ -575,6 +594,7 @@ end
         radiative_transfer = TwoStreamModel{FT}(domain),
         photosynthesis = FarquharModel{FT}(domain),
         conductance = MedlynConductanceModel{FT}(domain),
+            soil_moisture_stress = TuzetMoistureStressModel{FT}(),
         hydraulics = PlantHydraulicsModel{FT}(domain, LAI, toml_dict),
         energy = BigLeafEnergyModel{FT}(),
         sif = Lee2015SIFModel{FT}(),
@@ -614,6 +634,7 @@ function CanopyModel{FT}(
     radiative_transfer = TwoStreamModel{FT}(domain),
     photosynthesis = FarquharModel{FT}(domain),
     conductance = MedlynConductanceModel{FT}(domain),
+        soil_moisture_stress = TuzetMoistureStressModel{FT}(),
     hydraulics = PlantHydraulicsModel{FT}(domain, LAI, toml_dict),
     energy = BigLeafEnergyModel{FT}(),
     sif = Lee2015SIFModel{FT}(),
@@ -631,6 +652,7 @@ function CanopyModel{FT}(
         radiative_transfer,
         photosynthesis,
         conductance,
+        soil_moisture_stress,
         hydraulics,
         energy,
         sif,
@@ -691,6 +713,7 @@ canopy_components(::CanopyModel) = (
     :autotrophic_respiration,
     :energy,
     :sif,
+    :soil_moisture_stress,
 )
 
 """
@@ -878,13 +901,18 @@ function initialize_boundary_vars(model::CanopyModel{FT}, coords) where {FT}
 end
 
 """
-     ClimaLand.make_update_aux(canopy::CanopyModel{FT,
-                                                  <:AutotrophicRespirationModel,
-                                                  <:AbstractRadiationModel,
-                                                  <:AbstractPhotosynthesisModel,
-                                                  <:AbstractStomatalConductanceModel,
-                                                  <:PlantHydraulicsModel,},
-                              ) where {FT}
+     ClimaLand.make_update_aux(
+        canopy::CanopyModel{
+            FT,
+        <:AutotrophicRespirationModel,
+        <:AbstractRadiationModel,
+        <:AbstractPhotosynthesisModel,
+        <:AbstractStomatalConductanceModel,
+        <:PlantHydraulicsModel,
+        <:AbstractCanopyEnergyModel,
+        <:AbstractSoilMoistureStressModel,
+        },
+    ) where {FT}
 
 Creates the `update_aux!` function for the `CanopyModel`; a specific
 method for `update_aux!` for the case where the canopy model components
@@ -909,6 +937,7 @@ function ClimaLand.make_update_aux(
         <:AbstractStomatalConductanceModel,
         <:PlantHydraulicsModel,
         <:AbstractCanopyEnergyModel,
+        <:AbstractSoilMoistureStressModel,
     },
 ) where {FT}
     function update_aux!(p, Y, t)
@@ -1015,6 +1044,9 @@ function ClimaLand.make_update_aux(
         end
         # We update the fa[n_stem+n_leaf] element once we have computed transpiration
 
+        # Update soil moisture stress, used in photosynthesis and conductance
+        update_soil_moisture_stress!(p, Y, canopy.soil_moisture_stress, canopy)
+
         # Update Rd, An, Vcmax25 (if applicable to model) in place, GPP
         update_photosynthesis!(p, Y, canopy.photosynthesis, canopy)
 
@@ -1049,6 +1081,7 @@ function make_compute_exp_tendency(
         <:AbstractStomatalConductanceModel,
         <:PlantHydraulicsModel,
         <:AbstractCanopyEnergyModel,
+        <:AbstractSoilMoistureStressModel,
     },
 ) where {FT}
     components = canopy_components(canopy)
@@ -1079,6 +1112,7 @@ function make_compute_imp_tendency(
         <:AbstractStomatalConductanceModel,
         <:PlantHydraulicsModel,
         <:AbstractCanopyEnergyModel,
+        <:AbstractSoilMoistureStressModel,
     },
 ) where {FT}
     components = canopy_components(canopy)
@@ -1216,4 +1250,41 @@ function set_historical_cache!(p, Y0, m::AbstractPhotosynthesisModel, canopy)
     return nothing
 end
 
+"""
+    required_model_callbacks(start_date, t0, dt, model::CanopyModel)
+For some canopy components (namely the P-model), we need a callback to ensure regular
+updates at a preset frequency. This method dispatches off of the type of the photosynthesis
+model.
+"""
+function required_model_callbacks(start_date, t0, dt, model::CanopyModel)
+    return required_photosynthesis_model_callbacks(
+        start_date,
+        t0,
+        dt,
+        model,
+        model.photosynthesis,
+    )
+end
+
+
+"""
+    required_photosynthesis_model_callbacks(
+        start_date,
+        t0,
+        dt,
+        canopy,
+        photo_model::AbstractPhotosynthesisModel,
+    )
+For an AbstractPhotosynthesisModel in general, add no additional model callbacks.
+"""
+function required_photosynthesis_model_callbacks(
+    start_date,
+    t0,
+    dt,
+    canopy,
+    photo_model::AbstractPhotosynthesisModel,
+)
+    return ()
+end
+
 end

src/standalone/Vegetation/autotrophic_respiration.jl

@@ -6,7 +6,7 @@ abstract type AbstractAutotrophicRespirationModel{FT} <:
 """
     AutotrophicRespirationParameters{FT<:AbstractFloat}
 
-The required parameters for the autrophic respiration model, which is based 
+The required parameters for the autrophic respiration model, which is based
 off of the JULES model.
 Clark, D. B., et al. "The Joint UK Land Environment Simulator (JULES), model description–Part 2: carbon fluxes and vegetation dynamics." Geoscientific Model Development 4.3 (2011): 701-722.
 $(DocStringExtensions.FIELDS)
@@ -29,7 +29,7 @@ end
 Base.eltype(::AutotrophicRespirationParameters{FT}) where {FT} = FT
 
 """
-    AutotrophicRespirationModel{FT, ARP <: AutotrophicRespirationParameters{FT},} <: AbstractAutotrophicRespirationModel{FT} 
+    AutotrophicRespirationModel{FT, ARP <: AutotrophicRespirationParameters{FT},} <: AbstractAutotrophicRespirationModel{FT}
 
 The JULES autotrophic respiration model.
 
@@ -65,7 +65,7 @@ ClimaLand.auxiliary_domain_names(::AutotrophicRespirationModel) = (:surface,)
 Computes the autotrophic respiration rate  (mol co2 m^-2 s^-1) as the sum of the plant maintenance
 and growth respirations, according to the JULES model.
 
-Clark, D. B., et al. "The Joint UK Land Environment Simulator (JULES), model 
+Clark, D. B., et al. "The Joint UK Land Environment Simulator (JULES), model
 description–Part 2: carbon fluxes and vegetation dynamics." Geoscientific Model Development 4.3 (2011): 701-722.
 """
 function update_autotrophic_respiration!(
@@ -87,10 +87,9 @@ function update_autotrophic_respiration!(
     earth_param_set = canopy.parameters.earth_param_set
     grav = LP.grav(earth_param_set)
     ρ_l = LP.ρ_cloud_liq(earth_param_set)
-    (; sc, pc) = canopy.photosynthesis.parameters
     (; G_Function, Ω) = canopy.radiative_transfer.parameters
     cosθs = p.drivers.cosθs
-    β = @. lazy(moisture_stress(ψ.:($$i_end) * ρ_l * grav, sc, pc))
+    β = p.canopy.soil_moisture_stress.βm
     Vcmax25_leaf = get_Vcmax25_leaf(p, canopy.photosynthesis)
     Rd_leaf = get_Rd_leaf(p, canopy.photosynthesis)
     An_leaf = get_An_leaf(p, canopy.photosynthesis)

src/standalone/Vegetation/canopy_boundary_fluxes.jl

@@ -161,7 +161,8 @@ end
                                 <:AbstractPhotosynthesisModel,
                                 <:AbstractStomatalConductanceModel,
                                 <:PlantHydraulicsModel,
-                                <:AbstractCanopyEnergyModel}
+                                <:AbstractCanopyEnergyModel,
+                                <:AbstractSoilMoistureStressModel,
                             },
                             Y::ClimaCore.Fields.FieldVector,
                             t,
@@ -188,6 +189,7 @@ function canopy_boundary_fluxes!(
         <:AbstractStomatalConductanceModel,
         <:PlantHydraulicsModel,
         <:AbstractCanopyEnergyModel,
+        <:AbstractSoilMoistureStressModel,
     },
     Y::ClimaCore.Fields.FieldVector,
     t,

src/standalone/Vegetation/canopy_parameterizations.jl

@@ -3,7 +3,6 @@ export canopy_sw_rt_beer_lambert, # Radiative transfer
     canopy_sw_rt_two_stream,
     extinction_coeff,
     compute_G,
-    moisture_stress,
     # Conductance
     medlyn_term,
     medlyn_conductance,
@@ -392,24 +391,6 @@ function extinction_coeff(G_Function, cosθs::FT) where {FT}
     return K
 end
 
-"""
-    moisture_stress(pl::FT,
-                    sc::FT,
-                    pc::FT) where {FT}
-
-Computes the moisture stress factor (`β`), which is unitless,
- as a function of
-a constant (`sc`, 1/Pa), a reference pressure (`pc`, Pa), and
-the leaf water pressure (`pl`, Pa) .
-
-See Eqn 12.57 of G. Bonan's textbook,
-Climate Change and Terrestrial Ecosystem Modeling (2019).
-"""
-function moisture_stress(pl::FT, sc::FT, pc::FT) where {FT}
-    β = min(FT(1), (1 + exp(sc * pc)) / (1 + exp(sc * (pc - pl))))
-    return β
-end
-
 """
     conductance_molar_flux_to_m_per_s(gs::FT, T::FT, R::FT, P::FT) where {FT}