fix pmodel gpp for c4

Author: AlexisRenchon

src/standalone/Vegetation/Canopy.jl

@@ -222,7 +222,8 @@ end
         toml_dict::CP.ParamDict;
         is_c3 = clm_photosynthesis_parameters(domain.space.surface).is_c3,
         cstar = toml_dict["pmodel_cstar"],
-        β = toml_dict["pmodel_β"],
+        β_c3 = toml_dict["pmodel_β_c3"],
+        β_c4 = toml_dict["pmodel_β_c4"],
         temperature_dep_yield = true,
         ϕ0_c3 = toml_dict["pmodel_ϕ0_c3"],
         ϕ0_c4 = toml_dict["pmodel_ϕ0_c4"],
@@ -239,7 +240,8 @@ Constructs a P-model (an optimality model for photosynthesis) using default para
 
 The following default parameters (from the TOML file) are used:
 - cstar = 0.41 (unitless) - 4 * dA/dJmax, assumed to be a constant marginal cost (Wang 2017, Stocker 2020)
-- β = 146 (unitless) - Unit cost ratio of Vcmax to transpiration (Stocker 2020)
+- β_c3 = 146 (unitless) - Unit cost ratio of Vcmax to transpiration for C3 plants (Stocker 2020)
+- β_c4 = 146/9 ≈ 16.222 (unitless) - Unit cost ratio of Vcmax to transpiration for C4 plants (pyrealm)
 - ϕ0_c3 = 0.052 (unitless) - constant intrinsic quantum yield. Skillman (2008)
 - ϕ0_c4 = 0.057 (unitless) - constant intrinsic quantum yield. Skillman (2008)
 - ϕa0_c3 = 0.352*0.087 (unitless) - constant term in quadratic intrinsic quantum yield (Stocker 2020)
@@ -255,7 +257,8 @@ function PModel{FT}(
     toml_dict::CP.ParamDict;
     is_c3 = clm_photosynthesis_parameters(domain.space.surface).is_c3,
     cstar = toml_dict["pmodel_cstar"],
-    β = toml_dict["pmodel_β"],
+    β_c3 = toml_dict["pmodel_β_c3"],
+    β_c4 = toml_dict["pmodel_β_c4"],
     temperature_dep_yield = true,
     ϕ0_c3 = toml_dict["pmodel_ϕ0_c3"],
     ϕ0_c4 = toml_dict["pmodel_ϕ0_c4"],
@@ -269,7 +272,8 @@ function PModel{FT}(
 ) where {FT <: AbstractFloat}
     parameters = ClimaLand.Canopy.PModelParameters(
         cstar,
-        β,
+        β_c3,
+        β_c4,
         temperature_dep_yield,
         ϕ0_c3,
         ϕ0_c4,

src/standalone/Vegetation/optimal_lai.jl

@@ -557,6 +557,7 @@ function call_update_optimal_LAI(p, Y, t, current_date; canopy, dt, local_noon)
             P_air,
             p.canopy.biomass.vpd_gs,
             ca,
+            is_c3,
         ),
     )
 

src/standalone/Vegetation/pmodel.jl

@@ -17,8 +17,10 @@ $(DocStringExtensions.FIELDS)
 Base.@kwdef struct PModelParameters{FT <: AbstractFloat}
     "Constant describing cost of maintaining electron transport (unitless)"
     cstar::FT
-    "Ratio of unit costs of transpiration and carboxylation (unitless)"
-    β::FT
+    "Ratio of unit costs of transpiration and carboxylation for C3 plants (unitless)"
+    β_c3::FT
+    "Ratio of unit costs of transpiration and carboxylation for C4 plants (unitless)"
+    β_c4::FT
     "A boolean flag indicating if the quantum yield is a function of temperature or not"
     temperature_dep_yield::Bool
     "Temp-independent intrinsic quantum yield. (unitless); C3"
@@ -112,6 +114,13 @@ Base.eltype(::PModelConstants{FT}) where {FT} = FT
 Base.broadcastable(x::PModelParameters) = tuple(x)
 Base.broadcastable(x::PModelConstants) = tuple(x)
 
+"""
+    get_beta(is_c3, β_c3, β_c4)
+
+Returns the appropriate β (cost ratio) based on whether the plant is C3 or C4.
+"""
+get_beta(is_c3, β_c3, β_c4) = ifelse(is_c3 == one(is_c3), β_c3, β_c4)
+
 """
     PModelConstants(toml_dict::CP.ParamDict;
                     Kc25 = toml_dict["pmodel_Kc25"],
@@ -336,7 +345,8 @@ function compute_full_pmodel_outputs(
     is_c3 = FT(1),
 ) where {FT}
     # Unpack parameters
-    (; cstar, β) = parameters
+    (; cstar, β_c3, β_c4) = parameters
+    β = get_beta(is_c3, β_c3, β_c4)
 
     # Unpack constants
     (;
@@ -507,7 +517,8 @@ function update_optimal_EMA(
 ) where {FT}
     if local_noon_mask == FT(1.0)
         # Unpack parameters
-        (; cstar, β, α) = parameters
+        (; cstar, β_c3, β_c4, α) = parameters
+        β = get_beta(is_c3, β_c3, β_c4)
 
         # Unpack constants
         (;
@@ -670,7 +681,8 @@ function set_historical_cache!(p, Y0, model::PModel, canopy)
 
     parameters_init = PModelParameters(
         cstar = parameters.cstar,
-        β = parameters.β,
+        β_c3 = parameters.β_c3,
+        β_c4 = parameters.β_c4,
         temperature_dep_yield = parameters.temperature_dep_yield,
         ϕ0_c4 = parameters.ϕ0_c4,
         ϕ0_c3 = parameters.ϕ0_c3,
@@ -1400,7 +1412,8 @@ function compute_A0_daily(
     PPFD::FT,
     βm::FT,
 ) where {FT}
-    (; cstar, β) = parameters
+    (; cstar, β_c3, β_c4) = parameters
+    β = get_beta(is_c3, β_c3, β_c4)
     (; R, Kc25, Ko25, To, ΔHkc, ΔHko, Drel, ΔHΓstar, Γstar25, Mc, oi, ρ_water) =
         constants
 
@@ -1429,17 +1442,18 @@ function compute_A0_daily(
 end
 
 """
-    compute_chi(pmodel_parameters, pmodel_constants, T, P_air, VPD, ca)
+    compute_chi(pmodel_parameters, pmodel_constants, T, P_air, VPD, ca, is_c3=FT(1))
 
 Compute the optimal ratio of intercellular to ambient CO2 (chi = ci/ca) using P-model.
 
 # Arguments
-- `pmodel_parameters`: P-model parameters (including beta)
+- `pmodel_parameters`: P-model parameters (including β_c3 and β_c4)
 - `pmodel_constants`: P-model constants (including Gamma_star25, Kc25, Ko25, etc.)
 - `T::FT`: Temperature (K)
 - `P_air::FT`: Atmospheric pressure (Pa)
 - `VPD::FT`: Vapor pressure deficit (Pa)
 - `ca::FT`: Ambient CO2 mixing ratio (mol/mol)
+- `is_c3::FT`: C3 (1) or C4 (0) flag (default: C3)
 
 # Returns
 - `chi::FT`: Optimal ci/ca ratio (dimensionless), typically 0.7-0.85
@@ -1451,8 +1465,10 @@ function compute_chi(
     P_air::FT,
     VPD::FT,
     ca::FT,
+    is_c3::FT = FT(1),
 ) where {FT}
-    (; β) = pmodel_parameters
+    (; β_c3, β_c4) = pmodel_parameters
+    β = get_beta(is_c3, β_c3, β_c4)
     (; R, Kc25, Ko25, To, ΔHkc, ΔHko, Drel, ΔHΓstar, Γstar25, oi, ρ_water) =
         pmodel_constants
 

test/standalone/Vegetation/test_pmodel.jl

@@ -62,7 +62,8 @@ Constructs the parameter structure for the P-Model
 """
 function PModelParameters(inputs::Dict{String, Any}, FT)
     # these are default values used in Stocker 2020, Scott and Smith 2022
-    β = FT(inputs["beta"])
+    β_c3 = FT(inputs["beta"])
+    β_c4 = FT(inputs["beta"]) / FT(9)
     cstar = FT(0.41)
 
     # handle temperature-dependent quantum yield
@@ -82,7 +83,8 @@ function PModelParameters(inputs::Dict{String, Any}, FT)
 
     return ClimaLand.Canopy.PModelParameters(;
         cstar,
-        β,
+        β_c3,
+        β_c4,
         temperature_dep_yield,
         ϕ0_c3,
         ϕ0_c4,
@@ -271,7 +273,8 @@ end
         toml_dict = LP.create_toml_dict(FT)
         parameters = ClimaLand.Canopy.PModelParameters(
             cstar = FT(0.41),
-            β = FT(146),
+            β_c3 = FT(146),
+            β_c4 = FT(146 / 9),
             temperature_dep_yield = true,
             ϕ0_c3 = FT(0.052),
             ϕ0_c4 = FT(0.057),

toml/default_parameters.toml

@@ -223,15 +223,21 @@ tag = "ConstantResetAlbedo"
 
 # P model parameters
 ["pmodel_cstar"]
-value = 0.43
+value = 0.41
 type = "float"
 description = "4 * dA/dJmax, assumed to be a constant marginal cost (Wang 2017, Stocker 2020)"
 tag = "PModelParameters"
 
-["pmodel_β"]
-value = 51
+["pmodel_β_c3"]
+value = 146
 type = "float"
-description = "Unit cost ratio of Vcmax to transpiration (Stocker 2020)"
+description = "Unit cost ratio of Vcmax to transpiration for C3 plants (Stocker 2020)"
+tag = "PModelParameters"
+
+["pmodel_β_c4"]
+value = 16.222222222222222
+type = "float"
+description = "Unit cost ratio of Vcmax to transpiration for C4 plants (146/9, pyrealm)"
 tag = "PModelParameters"
 
 ["pmodel_ϕ0_c3"]

toml/uncalibrated_parameters.toml

@@ -202,10 +202,16 @@ type = "float"
 description = "4 * dA/dJmax, assumed to be a constant marginal cost (Wang 2017, Stocker 2020)"
 tag = "PModelParameters"
 
-["pmodel_β"]
+["pmodel_β_c3"]
 value = 146
 type = "float"
-description = "Unit cost ratio of Vcmax to transpiration (Stocker 2020)"
+description = "Unit cost ratio of Vcmax to transpiration for C3 plants (Stocker 2020)"
+tag = "PModelParameters"
+
+["pmodel_β_c4"]
+value = 16.222222222222222
+type = "float"
+description = "Unit cost ratio of Vcmax to transpiration for C4 plants (146/9, pyrealm)"
 tag = "PModelParameters"
 
 ["pmodel_ϕ0_c3"]