diff --git a/NEWS.md b/NEWS.md
index 2e40365078..a2c29ef5e6 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -4,6 +4,8 @@ ClimaLand.jl Release Notes
 main
 --------
  - Add regional simulation example PR [#757](https://github.com/CliMA/ClimaLand.jl/pull/757)
+ - ![][badge-💥breaking] Extend photosynthesis mechanism parameter to support fields.
+  PR[#774](https://github.com/CliMA/ClimaLand.jl/pull/774)
 
 v0.14.3
 --------
diff --git a/docs/src/standalone/apps/leaf_an.jl b/docs/src/standalone/apps/leaf_an.jl
index 5960996a6c..618a24284a 100644
--- a/docs/src/standalone/apps/leaf_an.jl
+++ b/docs/src/standalone/apps/leaf_an.jl
@@ -38,10 +38,11 @@ function ParamViz.parameterisation(
     Γstar = co2_compensation(Γstar25, ΔHΓstar, T, To, R)
     medlynt = 1 + g1 / sqrt(VPD)
     ci = intercellular_co2(ca, Γstar, medlynt)
-    Aj = light_assimilation(Canopy.C3(), J, ci, Γstar)
+    is_c3 = 1.0
+    Aj = light_assimilation(is_c3, J, ci, Γstar)
     Kc = MM_Kc(Kc25, ΔHkc, T, To, R)
     Ko = MM_Ko(Ko25, ΔHko, T, To, R)
-    Ac = rubisco_assimilation(Canopy.C3(), Vcmax, ci, Γstar, Kc, Ko, oi)
+    Ac = rubisco_assimilation(is_c3, Vcmax, ci, Γstar, Kc, Ko, oi)
     Rd = dark_respiration(Vcmax25, β, f, ΔHRd, T, To, R)
     An = net_photosynthesis(Ac, Aj, Rd, β)
     return An
@@ -83,7 +84,7 @@ function An_app_f()
             "Ω", # Clumping index
             "Γstar25", # co2 compensation at 25c
             "ΔHΓstar", # a constant energy of activation
-            "To", # a standard temperature 
+            "To", # a standard temperature
             "R", # the universal gas constant
             "Vcmax25", # the maximum rate of carboxylation of Rubisco
             "ΔHJmax", # a constant
diff --git a/docs/src/standalone/apps/leaf_an_ci.jl b/docs/src/standalone/apps/leaf_an_ci.jl
index 351338a131..417d658e4c 100644
--- a/docs/src/standalone/apps/leaf_an_ci.jl
+++ b/docs/src/standalone/apps/leaf_an_ci.jl
@@ -36,8 +36,9 @@ function ParamViz.parameterisation(
     Γstar = co2_compensation(Γstar25, ΔHΓstar, T, To, R)
     Kc = MM_Kc(Kc25, ΔHkc, T, To, R)
     Ko = MM_Ko(Ko25, ΔHko, T, To, R)
-    Ac = rubisco_assimilation(Canopy.C3(), Vcmax, ci, Γstar, Kc, Ko, oi)
-    Aj = light_assimilation(Canopy.C3(), J, ci, Γstar)
+    is_c3 = 1.0
+    Ac = rubisco_assimilation(is_c3, Vcmax, ci, Γstar, Kc, Ko, oi)
+    Aj = light_assimilation(is_c3, J, ci, Γstar)
     Rd = dark_respiration(Vcmax25, β, f, ΔHRd, T, To, R)
     An = net_photosynthesis(Ac, Aj, Rd, β)
     return An
@@ -79,7 +80,7 @@ function An_ci_app_f()
             "Ω", # Clumping index
             "Γstar25", # co2 compensation at 25c
             "ΔHΓstar", # a constant energy of activation
-            "To", # a standard temperature 
+            "To", # a standard temperature
             "R", # the universal gas constant
             "ΔHJmax", # a constant
             "θj", # an empirical "curvature parameter"
diff --git a/docs/tutorials/integrated/soil_canopy_tutorial.jl b/docs/tutorials/integrated/soil_canopy_tutorial.jl
index ff4ece3849..ce6c4acc83 100644
--- a/docs/tutorials/integrated/soil_canopy_tutorial.jl
+++ b/docs/tutorials/integrated/soil_canopy_tutorial.jl
@@ -231,8 +231,9 @@ radiative_transfer_args = (;
 
 conductance_args = (; parameters = MedlynConductanceParameters(FT; g1 = 141))
 
+is_c3 = FT(1) # set the photosynthesis mechanism to C3
 photosynthesis_args =
-    (; parameters = FarquharParameters(FT, Canopy.C3(); Vcmax25 = FT(5e-5)));
+    (; parameters = FarquharParameters(FT, is_c3; Vcmax25 = FT(5e-5)));
 
 K_sat_plant = FT(1.8e-8)
 RAI = (SAI + maxLAI) * f_root_to_shoot;
diff --git a/docs/tutorials/standalone/Canopy/canopy_tutorial.jl b/docs/tutorials/standalone/Canopy/canopy_tutorial.jl
index f6efdcf314..ac29d233fe 100644
--- a/docs/tutorials/standalone/Canopy/canopy_tutorial.jl
+++ b/docs/tutorials/standalone/Canopy/canopy_tutorial.jl
@@ -98,7 +98,7 @@ land_domain = Point(; z_sfc = FT(0.0))
 #   the simulation with atmospheric and radiative flux models. We also
 # read in the observed LAI and let that vary in time in a prescribed manner.
 
-# Use the data tools for reading FLUXNET data sets 
+# Use the data tools for reading FLUXNET data sets
 include(
     joinpath(pkgdir(ClimaLand), "experiments/integrated/fluxnet/data_tools.jl"),
 );
@@ -174,8 +174,8 @@ cond_params = MedlynConductanceParameters(FT; g1 = FT(141.0))
 stomatal_model = MedlynConductanceModel{FT}(cond_params);
 
 # Arguments for photosynthesis model:
-
-photo_params = FarquharParameters(FT, Canopy.C3(); Vcmax25 = FT(5e-5))
+is_c3 = FT(1) # set the photosynthesis mechanism to C3
+photo_params = FarquharParameters(FT, is_c3; Vcmax25 = FT(5e-5))
 
 photosynthesis_model = FarquharModel{FT}(photo_params);
 
diff --git a/experiments/benchmarks/land.jl b/experiments/benchmarks/land.jl
index 6024e0af00..cc809ef0c4 100644
--- a/experiments/benchmarks/land.jl
+++ b/experiments/benchmarks/land.jl
@@ -380,6 +380,14 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
         regridder_kwargs = (; extrapolation_bc,),
         file_reader_kwargs = (; preprocess_func = (data) -> data / 1_000_000,),
     )
+    # photosynthesis mechanism is read as a float, where 1.0 indicates c3 and 0.0 c4
+    is_c3 = SpaceVaryingInput(
+        joinpath(clm_artifact_path, "mechanism_map.nc"),
+        "c3_dominant",
+        surface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+    )
 
     # Plant Hydraulics and general plant parameters
     SAI = FT(0.0) # m2/m2
@@ -465,13 +473,8 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
     conductance_args =
         (; parameters = Canopy.MedlynConductanceParameters(FT; g1))
     # Set up photosynthesis
-    photosynthesis_args = (;
-        parameters = Canopy.FarquharParameters(
-            FT,
-            Canopy.C3();
-            Vcmax25 = Vcmax25,
-        )
-    )
+    photosynthesis_args =
+        (; parameters = Canopy.FarquharParameters(FT, is_c3; Vcmax25 = Vcmax25))
     # Set up plant hydraulics
     # Note that we clip all values of LAI below 0.05 to zero.
     # This is because we currently run into issues when LAI is
diff --git a/experiments/integrated/fluxnet/ozark_pft.jl b/experiments/integrated/fluxnet/ozark_pft.jl
index 4a3cd83d89..4d7be37565 100644
--- a/experiments/integrated/fluxnet/ozark_pft.jl
+++ b/experiments/integrated/fluxnet/ozark_pft.jl
@@ -184,8 +184,9 @@ radiative_transfer_args = (;
 conductance_args = (; parameters = MedlynConductanceParameters(FT; g1))
 # Set up photosynthesis
 # Set up photosynthesis
+is_c3 = FT(1) # set the photosynthesis mechanism to C3
 photosynthesis_args =
-    (; parameters = FarquharParameters(FT, Canopy.C3(); Vcmax25 = Vcmax25))
+    (; parameters = FarquharParameters(FT, is_c3; Vcmax25 = Vcmax25))
 # Set up plant hydraulics
 ai_parameterization = PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI)
 
diff --git a/experiments/integrated/fluxnet/run_fluxnet.jl b/experiments/integrated/fluxnet/run_fluxnet.jl
index 4ceed3ca3d..e41095c316 100644
--- a/experiments/integrated/fluxnet/run_fluxnet.jl
+++ b/experiments/integrated/fluxnet/run_fluxnet.jl
@@ -141,8 +141,9 @@ radiative_transfer_args = (;
 # Set up conductance
 conductance_args = (; parameters = MedlynConductanceParameters(FT; g1))
 # Set up photosynthesis
+is_c3 = FT(1) # set the photosynthesis mechanism to C3
 photosynthesis_args =
-    (; parameters = FarquharParameters(FT, Canopy.C3(); Vcmax25 = Vcmax25))
+    (; parameters = FarquharParameters(FT, is_c3; Vcmax25 = Vcmax25))
 # Set up plant hydraulics
 ai_parameterization = PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI)
 
diff --git a/experiments/integrated/global/global_soil_canopy.jl b/experiments/integrated/global/global_soil_canopy.jl
index 010c7c663d..977656bbfd 100644
--- a/experiments/integrated/global/global_soil_canopy.jl
+++ b/experiments/integrated/global/global_soil_canopy.jl
@@ -222,9 +222,9 @@ radiative_transfer_args = (;
 # Set up conductance
 conductance_args = (; parameters = Canopy.MedlynConductanceParameters(FT; g1))
 # Set up photosynthesis
-photosynthesis_args = (;
-    parameters = Canopy.FarquharParameters(FT, Canopy.C3(); Vcmax25 = Vcmax25)
-)
+is_c3 = FT(1) # set the photosynthesis mechanism to C3
+photosynthesis_args =
+    (; parameters = Canopy.FarquharParameters(FT, is_c3; Vcmax25 = Vcmax25))
 # Set up plant hydraulics
 # Not ideal
 LAIfunction = TimeVaryingInput(
diff --git a/experiments/integrated/performance/conservation/ozark_conservation_setup.jl b/experiments/integrated/performance/conservation/ozark_conservation_setup.jl
index 39c7e87a32..c195862b7b 100644
--- a/experiments/integrated/performance/conservation/ozark_conservation_setup.jl
+++ b/experiments/integrated/performance/conservation/ozark_conservation_setup.jl
@@ -129,8 +129,9 @@ radiative_transfer_args = (;
 # Set up conductance
 conductance_args = (; parameters = MedlynConductanceParameters(FT; g1))
 # Set up photosynthesis
+is_c3 = FT(1) # set the photosynthesis mechanism to C3
 photosynthesis_args =
-    (; parameters = FarquharParameters(FT, Canopy.C3(); Vcmax25 = Vcmax25))
+    (; parameters = FarquharParameters(FT, is_c3; Vcmax25 = Vcmax25))
 # Set up plant hydraulics
 ai_parameterization = PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI)
 
diff --git a/experiments/integrated/performance/profile_allocations.jl b/experiments/integrated/performance/profile_allocations.jl
index 22c41dfea1..2ec19f9192 100644
--- a/experiments/integrated/performance/profile_allocations.jl
+++ b/experiments/integrated/performance/profile_allocations.jl
@@ -266,8 +266,9 @@ radiative_transfer_args = (; parameters = TwoStreamParameters(FT))
 # Set up conductance
 conductance_args = (; parameters = MedlynConductanceParameters(FT; g1))
 # Set up photosynthesis
+is_c3 = FT(1) # set the photosynthesis mechanism to C3
 photosynthesis_args =
-    (; parameters = FarquharParameters(FT, Canopy.C3(); Vcmax25 = Vcmax25))
+    (; parameters = FarquharParameters(FT, is_c3; Vcmax25 = Vcmax25))
 # Set up plant hydraulics
 ai_parameterization = PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI)
 function root_distribution(z::T) where {T}
diff --git a/experiments/long_runs/land.jl b/experiments/long_runs/land.jl
index d044ff25dc..676b13136b 100644
--- a/experiments/long_runs/land.jl
+++ b/experiments/long_runs/land.jl
@@ -389,7 +389,14 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         regridder_kwargs = (; extrapolation_bc,),
         file_reader_kwargs = (; preprocess_func = (data) -> data / 1_000_000,),
     )
-
+    # photosynthesis mechanism is read as a float, where 1.0 indicates c3 and 0.0 c4
+    is_c3 = SpaceVaryingInput(
+        joinpath(clm_artifact_path, "mechanism_map.nc"),
+        "c3_dominant",
+        surface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+    )
     # Plant Hydraulics and general plant parameters
     SAI = FT(0.0) # m2/m2
     f_root_to_shoot = FT(3.5)
@@ -474,13 +481,8 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
     conductance_args =
         (; parameters = Canopy.MedlynConductanceParameters(FT; g1))
     # Set up photosynthesis
-    photosynthesis_args = (;
-        parameters = Canopy.FarquharParameters(
-            FT,
-            Canopy.C3();
-            Vcmax25 = Vcmax25,
-        )
-    )
+    photosynthesis_args =
+        (; parameters = Canopy.FarquharParameters(FT, is_c3; Vcmax25 = Vcmax25))
     # Set up plant hydraulics
 
     # Note that we clip all values of LAI below 0.05 to zero.
diff --git a/experiments/long_runs/land_region.jl b/experiments/long_runs/land_region.jl
index 92fa8afa87..f9ac565f4a 100644
--- a/experiments/long_runs/land_region.jl
+++ b/experiments/long_runs/land_region.jl
@@ -368,8 +368,18 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (10, 10, 15))
 
     # Energy Balance model
     ac_canopy = FT(2.5e4) # this will likely be 10x smaller!
-
     clm_artifact_path = ClimaLand.Artifacts.clm_data_folder_path(; context)
+    # Conductance Model
+    # g1 is read in units of sqrt(kPa) and then converted to sqrt(Pa)
+    g1 = SpaceVaryingInput(
+        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
+        "medlynslope",
+        surface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+        file_reader_kwargs = (; preprocess_func = (data) -> data * 10^(3 / 2),),
+    )
+
     # vcmax is read in units of umol CO2/m^2/s and then converted to mol CO2/m^2/s
     Vcmax25 = SpaceVaryingInput(
         joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
@@ -380,14 +390,13 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (10, 10, 15))
         file_reader_kwargs = (; preprocess_func = (data) -> data / 1_000_000,),
     )
 
-    # g1 is read in units of sqrt(kPa) and then converted to sqrt(Pa)
-    g1 = SpaceVaryingInput(
-        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
-        "medlynslope",
+    # c3 is read in as a float
+    is_c3 = SpaceVaryingInput(
+        joinpath(clm_artifact_path, "mechanism_map.nc"),
+        "c3_dominant",
         surface_space;
         regridder_type,
         regridder_kwargs = (; extrapolation_bc,),
-        file_reader_kwargs = (; preprocess_func = (data) -> data * 10^(3 / 2),),
     )
 
     # Plant Hydraulics and general plant parameters
@@ -474,13 +483,8 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (10, 10, 15))
     conductance_args =
         (; parameters = Canopy.MedlynConductanceParameters(FT; g1))
     # Set up photosynthesis
-    photosynthesis_args = (;
-        parameters = Canopy.FarquharParameters(
-            FT,
-            Canopy.C3();
-            Vcmax25 = Vcmax25,
-        )
-    )
+    photosynthesis_args =
+        (; parameters = Canopy.FarquharParameters(FT, is_c3; Vcmax25 = Vcmax25))
     # Set up plant hydraulics
 
     # Note that we clip all values of LAI below 0.05 to zero.
diff --git a/experiments/standalone/Vegetation/no_vegetation.jl b/experiments/standalone/Vegetation/no_vegetation.jl
index 570e456825..b71240098e 100644
--- a/experiments/standalone/Vegetation/no_vegetation.jl
+++ b/experiments/standalone/Vegetation/no_vegetation.jl
@@ -82,8 +82,8 @@ cond_params = MedlynConductanceParameters(FT; g1 = FT(141.0))
 
 stomatal_model = MedlynConductanceModel{FT}(cond_params);
 
-
-photo_params = FarquharParameters(FT, Canopy.C3(); Vcmax25 = FT(5e-5))
+is_c3 = FT(1) # set the photosynthesis mechanism to C3
+photo_params = FarquharParameters(FT, is_c3; Vcmax25 = FT(5e-5))
 
 photosynthesis_model = FarquharModel{FT}(photo_params);
 
diff --git a/experiments/standalone/Vegetation/varying_lai.jl b/experiments/standalone/Vegetation/varying_lai.jl
index 883b424af3..ecf9dbec3a 100644
--- a/experiments/standalone/Vegetation/varying_lai.jl
+++ b/experiments/standalone/Vegetation/varying_lai.jl
@@ -82,8 +82,8 @@ cond_params = MedlynConductanceParameters(FT; g1 = FT(141.0))
 
 stomatal_model = MedlynConductanceModel{FT}(cond_params);
 
-
-photo_params = FarquharParameters(FT, Canopy.C3(); Vcmax25 = FT(5e-5))
+is_c3 = FT(1) # set the photosynthesis mechanism to C3
+photo_params = FarquharParameters(FT, is_c3; Vcmax25 = FT(5e-5))
 
 photosynthesis_model = FarquharModel{FT}(photo_params);
 
diff --git a/experiments/standalone/Vegetation/varying_lai_with_stem.jl b/experiments/standalone/Vegetation/varying_lai_with_stem.jl
index ac403f8fe1..d1068db1ec 100644
--- a/experiments/standalone/Vegetation/varying_lai_with_stem.jl
+++ b/experiments/standalone/Vegetation/varying_lai_with_stem.jl
@@ -82,8 +82,8 @@ cond_params = MedlynConductanceParameters(FT; g1 = FT(141.0))
 
 stomatal_model = MedlynConductanceModel{FT}(cond_params);
 
-
-photo_params = FarquharParameters(FT, Canopy.C3(); Vcmax25 = FT(5e-5))
+is_c3 = FT(1) # set the photosynthesis mechanism to C3
+photo_params = FarquharParameters(FT, is_c3; Vcmax25 = FT(5e-5))
 
 photosynthesis_model = FarquharModel{FT}(photo_params);
 
diff --git a/ext/CreateParametersExt.jl b/ext/CreateParametersExt.jl
index f9abc5359b..6fcda3cad1 100644
--- a/ext/CreateParametersExt.jl
+++ b/ext/CreateParametersExt.jl
@@ -107,14 +107,14 @@ end
 """
     function FarquharParameters(
         FT,
-        mechanism::AbstractPhotosynthesisMechanism;
+        is_c3::Union{FT, ClimaCore.Fields.Field};
         Vcmax25 = FT(5e-5),
         kwargs...  # For individual parameter overrides
     )
 
     function FarquharParameters(
         toml_dict::CP.AbstractTOMLDict,
-        mechanism::AbstractPhotosynthesisMechanism;
+        is_c3::Union{AbstractFloat, ClimaCore.Fields.Field};
         Vcmax25 = FT(5e-5),
         kwargs...  # For individual parameter overrides
     )
@@ -123,25 +123,25 @@ Constructors for the FarquharParameters struct. Two variants:
 1. Pass in the float-type and retrieve parameter values from the default TOML dict.
 2. Pass in a TOML dictionary to retrieve parameter values.Possible calls:
 ```julia
-ClimaLand.Canopy.FarquharParameters(Float64, ClimaLand.Canopy.C3())
+ClimaLand.Canopy.FarquharParameters(Float64, 1.0)
 # Kwarg overrides
-ClimaLand.Canopy.FarquharParameters(Float64, ClimaLand.Canopy.C3(); Vcmax25 = 99999999, pc = 444444444)
+ClimaLand.Canopy.FarquharParameters(Float64, 1.0; Vcmax25 = 99999999, pc = 444444444)
 # TOML Dictionary:
 import ClimaParams as CP
 toml_dict = CP.create_toml_dict(Float32);
-ClimaLand.Canopy.FarquharParameters(toml_dict, ClimaLand.Canopy.C3(); Vcmax25 = 99999999, pc = 444444444)
+ClimaLand.Canopy.FarquharParameters(toml_dict, 1.0f0; Vcmax25 = 99999999, pc = 444444444)
 ```
 """
 FarquharParameters(
     ::Type{FT},
-    mechanism;
+    is_c3::Union{FT, ClimaCore.Fields.Field};
     kwargs...,
 ) where {FT <: AbstractFloat} =
-    FarquharParameters(CP.create_toml_dict(FT), mechanism; kwargs...)
+    FarquharParameters(CP.create_toml_dict(FT), is_c3; kwargs...)
 
 function FarquharParameters(
     toml_dict::CP.AbstractTOMLDict,
-    mechanism;
+    is_c3::Union{AbstractFloat, ClimaCore.Fields.Field};
     Vcmax25 = 5e-5,
     kwargs...,
 )
@@ -165,11 +165,24 @@ function FarquharParameters(
     )
     parameters = CP.get_parameter_values(toml_dict, name_map, "Land")
     FT = CP.float_type(toml_dict)
-    MECH = typeof(mechanism)
+    if maximum(is_c3) > 1
+        error(
+            "is_c3 has maximum of $(maximum(is_c3)). is_c3 should be between 0 and 1",
+        )
+    end
+    if minimum(is_c3) < 0
+        error(
+            "is_c3 has minimum of $(minimum(is_c3)). is_c3 should be between 0 and 1",
+        )
+    end
+    # if is_c3 is a field, is_c3 may contain values between 0.0 and 1.0 after regridding
+    # this deals with that possibility by rounding to the closest int
+    is_c3 = round.(is_c3)
+    MECH = typeof(is_c3)
     Vcmax25 = FT.(Vcmax25)
     VC = typeof(Vcmax25)
     return FarquharParameters{FT, MECH, VC}(;
-        mechanism,
+        is_c3,
         Vcmax25,
         parameters...,
         kwargs...,
@@ -229,8 +242,8 @@ 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)
+    is_c3 = FT(1)
+    return OptimalityFarquharParameters{FT, FT}(; params..., kwargs..., is_c3)
 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 fde43e2b8d..cc6b100104 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
@@ -139,7 +139,7 @@ function conductance_harvard(;
 end
 
 function photosynthesis_harvard(;
-    mechanism = Canopy.C3(),
+    is_c3 = FT(1),
     oi = FT(0.209),
     ϕ = FT(0.6),
     θj = FT(0.9),
@@ -176,7 +176,7 @@ function photosynthesis_harvard(;
         f,
         sc,
         pc,
-        mechanism,
+        is_c3,
     )
 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 c51d2dbd80..6326334d43 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
@@ -139,7 +139,7 @@ function conductance_ozark(;
 end
 
 function photosynthesis_ozark(;
-    mechanism = Canopy.C3(),
+    is_c3 = FT(1),
     oi = FT(0.209),
     ϕ = FT(0.6),
     θj = FT(0.9),
@@ -176,7 +176,7 @@ function photosynthesis_ozark(;
         f,
         sc,
         pc,
-        mechanism,
+        is_c3,
     )
 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 2b30e3eb11..bafdb21c53 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
@@ -139,7 +139,7 @@ function conductance_niwotridge(;
 end
 
 function photosynthesis_niwotridge(;
-    mechanism = Canopy.C3(),
+    is_c3 = FT(1),
     oi = FT(0.209),
     ϕ = FT(0.6),
     θj = FT(0.9),
@@ -176,7 +176,7 @@ function photosynthesis_niwotridge(;
         f,
         sc,
         pc,
-        mechanism,
+        is_c3,
     )
 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 4e8947bf8b..6fac7afa48 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
@@ -139,7 +139,7 @@ function conductance_vairaranch(;
 end
 
 function photosynthesis_vairaranch(;
-    mechanism = Canopy.C3(),
+    is_c3 = FT(1),
     oi = FT(0.209),
     ϕ = FT(0.6),
     θj = FT(0.9),
@@ -176,7 +176,7 @@ function photosynthesis_vairaranch(;
         f,
         sc,
         pc,
-        mechanism,
+        is_c3,
     )
 end
 
diff --git a/lib/ClimaLandSimulations/src/Fluxnet/run_fluxnet.jl b/lib/ClimaLandSimulations/src/Fluxnet/run_fluxnet.jl
index 1efe21c2be..e584770469 100644
--- a/lib/ClimaLandSimulations/src/Fluxnet/run_fluxnet.jl
+++ b/lib/ClimaLandSimulations/src/Fluxnet/run_fluxnet.jl
@@ -127,7 +127,7 @@ function run_fluxnet(
     photosynthesis_args = (;
         parameters = FarquharParameters(
             FT,
-            Canopy.C3();
+            FT(1);
             oi = params.photosynthesis.oi,
             ϕ = params.photosynthesis.ϕ,
             θj = params.photosynthesis.θj,
@@ -253,7 +253,7 @@ function run_fluxnet(
             land.soil.parameters.earth_param_set,
         )
     Y.soilco2.C .= FT(0.000412) # set to atmospheric co2, mol co2 per mol air
-    ψ_stem_0 = FT(-1e5 / 9800) # pressure in the leaf divided by rho_liquid*gravitational acceleration [m] 
+    ψ_stem_0 = FT(-1e5 / 9800) # pressure in the leaf divided by rho_liquid*gravitational acceleration [m]
     ψ_leaf_0 = FT(-2e5 / 9800)
     ψ_comps = setup.n_stem > 0 ? [ψ_stem_0, ψ_leaf_0] : ψ_leaf_0
 
diff --git a/src/standalone/Vegetation/canopy_parameterizations.jl b/src/standalone/Vegetation/canopy_parameterizations.jl
index b90ceaabdd..9926496089 100644
--- a/src/standalone/Vegetation/canopy_parameterizations.jl
+++ b/src/standalone/Vegetation/canopy_parameterizations.jl
@@ -45,8 +45,8 @@ end
         θs::FT,
     )
 
-Returns the leaf angle distribution value for CLM G function as a function of the 
-solar zenith angle and the leaf orientation index. See section 3.1 of 
+Returns the leaf angle distribution value for CLM G function as a function of the
+solar zenith angle and the leaf orientation index. See section 3.1 of
 https://www2.cesm.ucar.edu/models/cesm2/land/CLM50_Tech_Note.pdf
 """
 function compute_G(G::CLMGFunction{FT}, θs::FT) where {FT}
@@ -72,8 +72,8 @@ end
     )
 
 Computes the PAR and NIR absorbances, reflectances, and tranmittances
-for a canopy in the case of the 
-Beer-Lambert model. The absorbances are a function of the radiative transfer 
+for a canopy in the case of the
+Beer-Lambert model. The absorbances are a function of the radiative transfer
 model, as well as the magnitude of incident PAR and NIR radiation in W/m^2,
  the leaf area index, the extinction coefficient, and the
 soil albedo in the PAR and NIR bands. Returns a
@@ -130,11 +130,11 @@ end
     )
 
 Computes the PAR and NIR absorbances, reflectances, and tranmittances
-for a canopy in the case of the 
-Beer-Lambert model. The absorbances are a function of the radiative transfer 
-model, as well as the magnitude of incident PAR and NIR radiation in W/m^2, 
+for a canopy in the case of the
+Beer-Lambert model. The absorbances are a function of the radiative transfer
+model, as well as the magnitude of incident PAR and NIR radiation in W/m^2,
 the leaf area index, the extinction coefficient, and the
-soil albedo in the PAR and NIR bands. 
+soil albedo in the PAR and NIR bands.
 
 This model also depends on the diffuse fraction and the zenith angle.
 Returns a
@@ -191,14 +191,14 @@ end
         α_soil::FT
     )
 
-Computes the absorbed, reflected, and transmitted photon flux density 
-in terms of mol photons per m^2 per 
+Computes the absorbed, reflected, and transmitted photon flux density
+in terms of mol photons per m^2 per
 second for a radiation band.
 
-This applies the Beer-Lambert law, which is a function of incident 
+This applies the Beer-Lambert law, which is a function of incident
 radiation (`SW_IN`; moles of photons/m^2/), leaf reflectance
 (`α_leaf`), the extinction coefficient (`K`), leaf area index (`LAI`),
-and the albedo of the soil (`α_soil`). 
+and the albedo of the soil (`α_soil`).
 
 Returns a tuple of reflected, absorbed, and transmitted radiation in
 mol photons/m^2/s.
@@ -230,13 +230,13 @@ end
         α_soil::FT,
     )
 
-Computes the absorbed, transmitted, and reflected  photon flux density 
-in terms of mol photons  per m^2 per second for a radiation band. 
+Computes the absorbed, transmitted, and reflected  photon flux density
+in terms of mol photons  per m^2 per second for a radiation band.
 
 This applies the two-stream radiative transfer solution which takes into account
-the impacts of scattering within the canopy. The function takes in all 
-parameters from the parameter struct of a TwoStreamModel, along with the 
-incident radiation, LAI, extinction coefficient K, soil albedo from the 
+the impacts of scattering within the canopy. The function takes in all
+parameters from the parameter struct of a TwoStreamModel, along with the
+incident radiation, LAI, extinction coefficient K, soil albedo from the
 canopy soil_driver, solar zenith angle, and τ.
 
 Returns a tuple of reflected, absorbed, and transmitted radiation in
@@ -276,7 +276,7 @@ function plant_absorbed_pfd(
     c²θ̄ = pi * G / 4
     β = 0.5 * (ω + diff * c²θ̄) / ω
 
-    # Compute coefficients for two-stream solution 
+    # Compute coefficients for two-stream solution
     b = 1 - ω + ω * β
     c = ω * β
     d = ω * β₀ * μ̄ * K
@@ -312,7 +312,7 @@ function plant_absorbed_pfd(
             p₁ * s₂ * (d - c - h₁ / σ * (u₁ + μ̄ * K))
         )
 
-    # h coefficients for direct downward flux 
+    # h coefficients for direct downward flux
     h₄ = -f * p₃ - c * d
     h₅ =
         -1 / d₂ *
@@ -351,7 +351,7 @@ function plant_absorbed_pfd(
         # Compute cumulative LAI at this layer
         L = i * Lₗ
 
-        # Compute the direct fluxes into/out of the layer 
+        # Compute the direct fluxes into/out of the layer
         I_dir_up =
             h₁ * exp(-K * L * Ω) / σ +
             h₂ * exp(-h * L * Ω) +
@@ -368,7 +368,7 @@ function plant_absorbed_pfd(
         I_dif_up = h₇ * exp(-h * L * Ω) + h₈ * exp(h * L * Ω)
         I_dif_dn = h₉ * exp(-h * L * Ω) + h₁₀ * exp(h * L * Ω)
 
-        # Energy balance giving radiation absorbed in the layer 
+        # Energy balance giving radiation absorbed in the layer
         if i == 0
             I_dir_abs = 0
             I_dif_abs = 0
@@ -385,7 +385,7 @@ function plant_absorbed_pfd(
         # Add radiation absorbed in the layer to total absorbed radiation
         F_abs += (1 - frac_diff) * I_dir_abs + (frac_diff) * I_dif_abs
 
-        # Save input/output values to compute energy balance of next layer 
+        # Save input/output values to compute energy balance of next layer
         I_dir_up_prev = I_dir_up
         I_dir_dn_prev = I_dir_dn
         I_dif_up_prev = I_dif_up
@@ -411,7 +411,7 @@ end
 Computes the fraction of diffuse radiation (`diff_frac`) as a function
 of the solar zenith angle (`θs`), the total surface incident shortwave radiation (`SW_IN`),
 the air temperature (`T`), air pressure (`P`), specific humidity (`q`), and the day of the year
-(`td`). 
+(`td`).
 
 See Appendix A of Braghiere, "Evaluation of turbulent fluxes of CO2, sensible heat,
 and latent heat as a function of aerosol optical depth over the course of deforestation
@@ -419,10 +419,10 @@ in the Brazilian Amazon" 2013.
 
 Note that cos(θs) is equal to zero when θs = π/2, and this is a coefficient
 of k₀, which we divide by in this expression. This can amplify small errors
-when θs is near π/2. 
+when θs is near π/2.
 
-This formula is empirical and can yied negative numbers depending on the 
-input, which, when dividing by something very near zero, 
+This formula is empirical and can yied negative numbers depending on the
+input, which, when dividing by something very near zero,
 can become large negative numbers.
 
 Because of that, we cap the returned value to lie within [0,1].
@@ -481,7 +481,7 @@ end
     intercellular_co2(ca::FT, Γstar::FT, medlyn_factor::FT) where{FT}
 
 Computes the intercellular CO2 concentration (mol/mol) given the
-atmospheric concentration (`ca`, mol/mol), the CO2 compensation (`Γstar`, 
+atmospheric concentration (`ca`, mol/mol), the CO2 compensation (`Γstar`,
  mol/mol), and the Medlyn factor (unitless).
 """
 function intercellular_co2(ca::FT, Γstar::FT, medlyn_term::FT) where {FT}
@@ -516,8 +516,21 @@ function co2_compensation(
 end
 
 """
-    rubisco_assimilation(::C3,
-                         Vcmax::FT,
+    rubisco_assimilation(is_c3::AbstractFloat, args...)
+
+Calls the correct rubisco assimilation function based on the `is_c3`.
+
+A `is_c3` value of 1.0 corresponds to C3 photosynthesis and calls
+`c3_rubisco_assimilation`, while 0.0 corresponds to C4 photsynthesis and calls
+`c4_rubisco_assimilation`.
+"""
+function rubisco_assimilation(is_c3::AbstractFloat, args...)
+    is_c3 > 0.5 ? c3_rubisco_assimilation(args...) :
+    c4_rubisco_assimilation(args...)
+end
+
+"""
+    c3_rubisco_assimilation(Vcmax::FT,
                          ci::FT,
                          Γstar::FT,
                          Kc::FT,
@@ -526,16 +539,15 @@ end
 
 Computes the Rubisco limiting rate of photosynthesis for C3 plants (`Ac`),
 in units of moles CO2/m^2/s,
-as a function of the maximum rate of carboxylation of Rubisco (`Vcmax`), 
-the leaf internal carbon dioxide partial pressure (`ci`), 
+as a function of the maximum rate of carboxylation of Rubisco (`Vcmax`),
+the leaf internal carbon dioxide partial pressure (`ci`),
 the CO2 compensation point (`Γstar`), and Michaelis-Menten parameters
 for CO2 and O2, respectively, (`Kc`) and (`Ko`).
 
 The empirical parameter oi is equal to 0.209 (mol/mol).
 See Table 11.5 of G. Bonan's textbook, Climate Change and Terrestrial Ecosystem Modeling (2019).
 """
-function rubisco_assimilation(
-    ::C3,
+function c3_rubisco_assimilation(
     Vcmax::FT,
     ci::FT,
     Γstar::FT,
@@ -548,19 +560,32 @@ function rubisco_assimilation(
 end
 
 """
-    rubisco_assimilation(::C4, Vcmax::FT,_...) where {FT}
+    c4_rubisco_assimilation(Vcmax::FT,_...) where {FT}
 
 Computes the Rubisco limiting rate of photosynthesis for C4 plants (`Ac`)
 in units of moles CO2/m^2/s,
 as equal to the maximum rate of carboxylation of Rubisco (`Vcmax`).
 """
-function rubisco_assimilation(::C4, Vcmax::FT, _...) where {FT}
+function c4_rubisco_assimilation(Vcmax::FT, _...) where {FT}
     Ac = Vcmax
     return Ac
 end
 
 """
-    light_assimilation(::C3,
+    light_assimilation(is_c3::AbstractFloat, args...)
+
+Calls the correct light assimilation function based on the `is_c3`.
+
+A `is_c3` value of 1.0 corresponds to C3 photosynthesis and calls
+`c3_light_assimilation`, while 0.0 corresponds to C4 photsynthesis and calls
+`c4_light_assimilation`.
+"""
+function light_assimilation(is_c3::AbstractFloat, args...)
+    is_c3 > 0.5 ? c3_light_assimilation(args...) :
+    c4_light_assimilation(args...)
+end
+"""
+    c3_light_assimilation(
                        J::FT,
                        ci::FT,
                        Γstar::FT) where {FT}
@@ -572,19 +597,19 @@ and the CO2 compensation point (`Γstar`).
 
 See Table 11.5 of G. Bonan's textbook, Climate Change and Terrestrial Ecosystem Modeling (2019).
 """
-function light_assimilation(::C3, J::FT, ci::FT, Γstar::FT) where {FT}
+function c3_light_assimilation(J::FT, ci::FT, Γstar::FT) where {FT}
     Aj = J * (ci - Γstar) / (4 * (ci + 2 * Γstar))
     return Aj
 end
 
 """
-    light_assimilation(::C4, J::FT, _...) where {FT}
+    light_assimilation(J::FT, _...) where {FT}
 
 Computes the electron transport limiting rate (`Aj`),
 in units of moles CO2/m^2/s, for C4 plants, as equal to
 the rate of electron transport (`J`).
 """
-function light_assimilation(::C4, J::FT, _...) where {FT}
+function c4_light_assimilation(J::FT, _...) where {FT}
     Aj = J
     return Aj
 end
@@ -593,12 +618,12 @@ end
     max_electron_transport(Vcmax::FT) where {FT}
 
 Computes the maximum potential rate of electron transport (`Jmax`),
-in units of mol/m^2/s, 
+in units of mol/m^2/s,
 as a function of Vcmax at 25 °C (`Vcmax25`),
 a constant (`ΔHJmax`), a standard temperature (`To`),
 the unversal gas constant (`R`), and the temperature (`T`).
 
-See Table 11.5 of G. Bonan's textbook, 
+See Table 11.5 of G. Bonan's textbook,
 Climate Change and Terrestrial Ecosystem Modeling (2019).
 """
 function max_electron_transport(
@@ -624,7 +649,7 @@ in units of mol/m^2/s, as a function of
 the maximum potential rate of electron transport (`Jmax`),
 absorbed photosynthetically active radiation (`APAR`),
 an empirical "curvature parameter" (`θj`; Bonan Eqn 11.21)
-and the quantum yield of photosystem II (`ϕ`). 
+and the quantum yield of photosystem II (`ϕ`).
 
 See Ch 11, G. Bonan's textbook, Climate Change and Terrestrial Ecosystem Modeling (2019).
 """
@@ -701,9 +726,9 @@ end
                        β::FT) where {FT}
 
 Computes the total net carbon assimilation (`An`),
-in units of mol CO2/m^2/s, as a function of 
+in units of mol CO2/m^2/s, as a function of
 the Rubisco limiting factor (`Ac`), the electron transport limiting rate (`Aj`),
-dark respiration (`Rd`), and the moisture stress factor (`β`). 
+dark respiration (`Rd`), and the moisture stress factor (`β`).
 
 See Table 11.5 of G. Bonan's textbook, Climate Change and Terrestrial Ecosystem Modeling (2019).
 """
@@ -719,10 +744,10 @@ end
 
 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) . 
+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, 
+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}
@@ -745,7 +770,7 @@ and the moisture stress factor (`β`), an empirical factor `f` is equal to 0.015
 a constant (`ΔHRd`), a standard temperature (`To`),
 the unversal gas constant (`R`), and the temperature (`T`).
 
-See Table 11.5 of G. Bonan's textbook, 
+See Table 11.5 of G. Bonan's textbook,
 Climate Change and Terrestrial Ecosystem Modeling (2019).
 """
 function dark_respiration(
@@ -767,7 +792,7 @@ end
              LAI::FT,
              Ω::FT) where {FT}
 
-Computes the total canopy photosynthesis (`GPP`) as a function of 
+Computes the total canopy photosynthesis (`GPP`) as a function of
 the total net carbon assimilation (`An`), the extinction coefficient (`K`),
 leaf area index (`LAI`) and the clumping index (`Ω`).
 """
@@ -785,7 +810,7 @@ and (1) converts it to m/s, (2) upscales to the entire canopy, by assuming
 the leaves in the canopy are in parallel and hence multiplying
 by LAI.
 
-TODO: Check what CLM does, and check if we can use the same function 
+TODO: Check what CLM does, and check if we can use the same function
 for GPP from An, and make more general.
 """
 function upscale_leaf_conductance(
@@ -803,10 +828,10 @@ end
     arrhenius_function(T::FT, To::FT, R::FT, ΔH::FT)
 
 Computes the Arrhenius function at temperature `T` given
-the reference temperature `To=298.15K`, the universal 
+the reference temperature `To=298.15K`, the universal
 gas constant `R`, and the energy activation `ΔH`.
 
-See Table 11.5 of G. Bonan's textbook, 
+See Table 11.5 of G. Bonan's textbook,
 Climate Change and Terrestrial Ecosystem Modeling (2019).
 """
 function arrhenius_function(T::FT, To::FT, R::FT, ΔH::FT) where {FT}
@@ -826,7 +851,7 @@ as a function of its value at 25 °C (`Kc25`),
 a constant (`ΔHkc`), a standard temperature (`To`),
 the unversal gas constant (`R`), and the temperature (`T`).
 
-See Table 11.5 of G. Bonan's textbook, 
+See Table 11.5 of G. Bonan's textbook,
 Climate Change and Terrestrial Ecosystem Modeling (2019).
 """
 function MM_Kc(Kc25::FT, ΔHkc::FT, T::FT, To::FT, R::FT) where {FT}
@@ -847,7 +872,7 @@ as a function of its value at 25 °C (`Ko25`),
 a constant (`ΔHko`), a standard temperature (`To`),
 the universal gas constant (`R`), and the temperature (`T`).
 
-See Table 11.5 of G. Bonan's textbook, 
+See Table 11.5 of G. Bonan's textbook,
 Climate Change and Terrestrial Ecosystem Modeling (2019).
 """
 function MM_Ko(Ko25::FT, ΔHko::FT, T::FT, To::FT, R::FT) where {FT}
@@ -863,11 +888,11 @@ end
            ep5::FT) where {FT}
 
 Computes the maximum rate of carboxylation of Rubisco (`Vcmax`),
-in units of mol/m^2/s, 
+in units of mol/m^2/s,
 as a function of temperature (`T`), Vcmax at the reference temperature 25 °C (`Vcmax25`),
 the universal gas constant (`R`), and the reference temperature (`To`).
 
-See Table 11.5 of G. Bonan's textbook, 
+See Table 11.5 of G. Bonan's textbook,
 Climate Change and Terrestrial Ecosystem Modeling (2019).
 """
 function compute_Vcmax(
@@ -911,13 +936,13 @@ end
                        An::FT,
                        ca::FT) where {FT}
 
-Computes the stomatal conductance according to Medlyn, as a function of 
-the minimum stomatal conductance (`g0`), 
+Computes the stomatal conductance according to Medlyn, as a function of
+the minimum stomatal conductance (`g0`),
 the relative diffusivity of water vapor with respect to CO2 (`Drel`),
 the Medlyn term (unitless), the biochemical demand for CO2 (`An`), and the
 atmospheric concentration of CO2 (`ca`).
 
-This returns the conductance in units of mol/m^2/s. It must be converted to 
+This returns the conductance in units of mol/m^2/s. It must be converted to
 m/s using the molar density of water prior to use in SurfaceFluxes.jl.
 """
 function medlyn_conductance(
@@ -933,11 +958,11 @@ end
 
 """
     penman_monteith(
-        Δ::FT, # Rate of change of saturation vapor pressure with air temperature. (Pa K−1)  
+        Δ::FT, # Rate of change of saturation vapor pressure with air temperature. (Pa K−1)
         Rn::FT, # Net irradiance (W m−2)
         G::FT, # Ground heat flux (W m−2)
         ρa::FT, # Dry air density (kg m−3)
-        cp::FT, # Specific heat capacity of air (J kg−1 K−1) 
+        cp::FT, # Specific heat capacity of air (J kg−1 K−1)
         VPD::FT, # vapor pressure deficit (Pa)
         ga::FT, # atmospheric conductance (m s−1)
         γ::FT, # Psychrometric constant (γ ≈ 66 Pa K−1)
@@ -945,7 +970,7 @@ end
         Lv::FT, # Volumetric latent heat of vaporization (J m-3)
         ) where {FT}
 
-Computes the evapotranspiration in m/s using the Penman-Monteith equation.       
+Computes the evapotranspiration in m/s using the Penman-Monteith equation.
 """
 function penman_monteith(
     Δ::FT,
@@ -970,11 +995,11 @@ end
                      LAI::FT, # Leaf area index
                      SAI::FT,
                      RAI::FT,
-                     ηsl::FT, # live stem  wood coefficient (kg C m-3) 
+                     ηsl::FT, # live stem  wood coefficient (kg C m-3)
                      h::FT, # canopy height (m)
                      σl::FT # Specific leaf density (kg C m-2 [leaf])
                      μr::FT, # Ratio root nitrogen to top leaf nitrogen (-), typical value 1.0
-                     μs::FT, # Ratio stem nitrogen to top leaf nitrogen (-), typical value 0.1 
+                     μs::FT, # Ratio stem nitrogen to top leaf nitrogen (-), typical value 0.1
                     ) where {FT}
 
 Computes the nitrogen content of leafs (Nl), roots (Nr) and stems (Ns).
@@ -985,11 +1010,11 @@ function nitrogen_content(
     LAI::FT, # Leaf area index
     SAI::FT,
     RAI::FT,
-    ηsl::FT, # live stem  wood coefficient (kg C m-3) 
+    ηsl::FT, # live stem  wood coefficient (kg C m-3)
     h::FT, # canopy height (m)
     σl::FT, # Specific leaf density (kg C m-2 [leaf])
     μr::FT, # Ratio root nitrogen to top leaf nitrogen (-), typical value 1.0
-    μs::FT, # Ratio stem nitrogen to top leaf nitrogen (-), typical value 0.1 
+    μs::FT, # Ratio stem nitrogen to top leaf nitrogen (-), typical value 0.1
 ) where {FT}
     Sc = ηsl * h * LAI * ClimaLand.heaviside(SAI)
     Rc = σl * RAI
@@ -1010,7 +1035,7 @@ end
         ) where {FT}
 
 Computes plant maintenance respiration as a function of dark respiration (Rd),
-the nitrogen content of leafs (Nl), roots (Nr) and stems (Ns), 
+the nitrogen content of leafs (Nl), roots (Nr) and stems (Ns),
 and the soil moisture factor (β).
 """
 function plant_respiration_maintenance(
diff --git a/src/standalone/Vegetation/optimality_farquhar.jl b/src/standalone/Vegetation/optimality_farquhar.jl
index 1422b20271..ee8228b8fa 100644
--- a/src/standalone/Vegetation/optimality_farquhar.jl
+++ b/src/standalone/Vegetation/optimality_farquhar.jl
@@ -8,9 +8,12 @@ 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,
+    MECH <: Union{FT, ClimaCore.Fields.Field},
+}
     "Photosynthesis mechanism: C3 only"
-    mechanism::C3
+    is_c3::MECH
     "Γstar at 25 °C (mol/mol)"
     Γstar25::FT
     "Michaelis-Menten parameter for CO2 at 25 °C (mol/mol)"
@@ -91,7 +94,7 @@ ClimaLand.auxiliary_domain_names(::OptimalityFarquharModel) =
 
 Computes the net photosynthesis rate `An` for the Optimality Farquhar model, along with the
 dark respiration `Rd`, and the value of `Vcmax25`, and updates them in place.
-        
+
  To do so, we require the canopy leaf temperature `T`, Medlyn factor, `APAR` in
 photons per m^2 per second, CO2 concentration in the atmosphere,
 moisture stress factor `β` (unitless), and the universal gas constant
@@ -118,7 +121,7 @@ function update_photosynthesis!(
         To,
         θj,
         ϕ,
-        mechanism,
+        is_c3,
         oi,
         Kc25,
         Ko25,
@@ -146,8 +149,8 @@ function update_photosynthesis!(
     Jmax = rates.:1
     Vcmax = rates.:2
     J = electron_transport.(APAR, Jmax, θj, ϕ)
-    Aj = light_assimilation.(mechanism, J, ci, Γstar)
-    Ac = rubisco_assimilation.(mechanism, Vcmax, ci, Γstar, Kc, Ko, oi)
+    Aj = light_assimilation.(is_c3, J, ci, Γstar)
+    Ac = rubisco_assimilation.(is_c3, Vcmax, ci, Γstar, Kc, Ko, oi)
 
     @. Vcmax25 = Vcmax / arrhenius_function(T, To, R, ΔHVcmax)
     @. Rd = dark_respiration(Vcmax25, β, f, ΔHRd, T, To, R)
diff --git a/src/standalone/Vegetation/photosynthesis.jl b/src/standalone/Vegetation/photosynthesis.jl
index 3257fd894f..8d220f3c12 100644
--- a/src/standalone/Vegetation/photosynthesis.jl
+++ b/src/standalone/Vegetation/photosynthesis.jl
@@ -1,25 +1,11 @@
-export FarquharParameters, FarquharModel, C3, C4
-
-abstract type AbstractPhotosynthesisMechanism end
-"""
-    C3 <: AbstractPhotosynthesisMechanism
-
-Helper struct for dispatching between C3 and C4 photosynthesis.
-"""
-struct C3 <: AbstractPhotosynthesisMechanism end
-
-"""
-    C4 <: AbstractPhotosynthesisMechanism
-
-Helper struct for dispatching between C3 and C4 photosynthesis.
-"""
-struct C4 <: AbstractPhotosynthesisMechanism end
+export FarquharParameters, FarquharModel
 
 abstract type AbstractPhotosynthesisModel{FT} <: AbstractCanopyComponent{FT} end
+
 """
     FarquharParameters{
         FT<:AbstractFloat,
-        MECH <: AbstractPhotosynthesisMechanism,
+        MECH <: Union{FT, ClimaCore.Fields.Field},
         VC <: Union{FT, ClimaCore.Fields.Field},
     }
 
@@ -28,7 +14,7 @@ $(DocStringExtensions.FIELDS)
 """
 Base.@kwdef struct FarquharParameters{
     FT <: AbstractFloat,
-    MECH <: AbstractPhotosynthesisMechanism,
+    MECH <: Union{FT, ClimaCore.Fields.Field},
     VC <: Union{FT, ClimaCore.Fields.Field},
 }
     "Vcmax at 25 °C (mol CO2/m^2/s)"
@@ -65,8 +51,8 @@ Base.@kwdef struct FarquharParameters{
     sc::FT
     "Reference water pressure for the moisture stress factor (Pa) [Tuzet et al. (2003)]"
     pc::FT
-    "Photosynthesis mechanism: C3 or C4"
-    mechanism::MECH
+    "Photosynthesis mechanism: 1.0 indicates C3, 0.0 indicates C4"
+    is_c3::MECH
 end
 
 Base.eltype(::FarquharParameters{FT}) where {FT} = FT
@@ -108,7 +94,7 @@ function photosynthesis_at_a_point_Farquhar(
     To,
     θj,
     ϕ,
-    mechanism,
+    is_c3,
     oi,
     Kc25,
     Ko25,
@@ -120,10 +106,10 @@ function photosynthesis_at_a_point_Farquhar(
     Vcmax = compute_Vcmax(Vcmax25, T, To, R, ΔHVcmax)
     Γstar = co2_compensation(Γstar25, ΔHΓstar, T, To, R)
     ci = intercellular_co2(c_co2, Γstar, medlyn_factor)
-    Aj = light_assimilation(mechanism, J, ci, Γstar)
+    Aj = light_assimilation(is_c3, J, ci, Γstar)
     Kc = MM_Kc(Kc25, ΔHkc, T, To, R)
     Ko = MM_Ko(Ko25, ΔHko, T, To, R)
-    Ac = rubisco_assimilation(mechanism, Vcmax, ci, Γstar, Kc, Ko, oi)
+    Ac = rubisco_assimilation(is_c3, Vcmax, ci, Γstar, Kc, Ko, oi)
     return net_photosynthesis(Ac, Aj, Rd, β)
 end
 
@@ -172,7 +158,7 @@ function update_photosynthesis!(
         To,
         θj,
         ϕ,
-        mechanism,
+        is_c3,
         oi,
         Kc25,
         Ko25,
@@ -198,7 +184,7 @@ function update_photosynthesis!(
         To,
         θj,
         ϕ,
-        mechanism,
+        is_c3,
         oi,
         Kc25,
         Ko25,
@@ -207,7 +193,6 @@ function update_photosynthesis!(
     )
     Vcmax25field .= Vcmax25
 end
-Base.broadcastable(m::AbstractPhotosynthesisMechanism) = tuple(m)
 Base.broadcastable(m::FarquharParameters) = tuple(m)
 
 include("./optimality_farquhar.jl")
diff --git a/test/standalone/Vegetation/canopy_model.jl b/test/standalone/Vegetation/canopy_model.jl
index 9f52df0fb2..0446693a06 100644
--- a/test/standalone/Vegetation/canopy_model.jl
+++ b/test/standalone/Vegetation/canopy_model.jl
@@ -21,13 +21,26 @@ import ClimaParams
 
 @testset "Canopy software pipes" begin
     for FT in (Float32, Float64)
-        domain = Point(; z_sfc = FT(0.0))
-        for g1 in (FT(790), fill(FT(790), domain.space.surface))
-            # create new field with constant value everywhere
-            Vcmax25 = fill(FT(9e-5), domain.space.surface)
+        domain = ClimaLand.Domains.SphericalSurface(;
+            radius = FT(100.0),
+            nelements = 10,
+            npolynomial = 1,
+        )
+        # create a field with both 1.0s and 0.0s
+        mechanism_field = ClimaCore.Fields.Field(FT, domain.space.surface)
+        ClimaCore.Fields.set!(
+            x -> x.coordinates.lat > 0 ? 0.0 : 1.0,
+            mechanism_field,
+        )
+        # create one case where parameters are spatially varying and one where not
+        g1_cases = (FT(790), fill(FT(790), domain.space.surface))
+        Vcmax25_cases = (FT(9e-5), fill(FT(9e-5), domain.space.surface))
+        mechanism_cases = (FT(1), mechanism_field)
+        zipped = zip(g1_cases, Vcmax25_cases, mechanism_cases)
+        for (g1, Vcmax25, is_c3) in zipped
             AR_params = AutotrophicRespirationParameters(FT)
             RTparams = BeerLambertParameters(FT)
-            photosynthesis_params = FarquharParameters(FT, C3(); Vcmax25)
+            photosynthesis_params = FarquharParameters(FT, is_c3; Vcmax25)
             stomatal_g_params = MedlynConductanceParameters(FT; g1)
 
             AR_model = AutotrophicRespirationModel{FT}(AR_params)
@@ -211,7 +224,7 @@ import ClimaParams
             # Check that structure of Y is value (will error if not)
             @test !isnothing(zero(Y))
             @test typeof(canopy.energy) == PrescribedCanopyTempModel{FT}
-            @test propertynames(p) == (:canopy, :drivers)
+            @test propertynames(p) == (:canopy, :dss_buffer_2d, :drivers)
             for component in ClimaLand.Canopy.canopy_components(canopy)
                 # Only hydraulics has a prognostic variable
                 if component == :hydraulics
@@ -348,22 +361,26 @@ import ClimaParams
             @test ClimaLand.surface_evaporative_scaling(canopy, Y, p) == FT(1.0)
             @test ClimaLand.surface_height(canopy, Y, p) == compartment_faces[1]
             T_sfc = FT.(T_atmos(t0))
-            @test Array(
-                parent(ClimaLand.surface_temperature(canopy, Y, p, t0)),
-            ) == [T_sfc]
+            @test all(
+                Array(
+                    parent(ClimaLand.surface_temperature(canopy, Y, p, t0)),
+                ) .== [T_sfc],
+            )
             @test ClimaLand.surface_temperature(canopy, Y, p, t0) isa
                   ClimaCore.Fields.Field
-            @test Array(
-                parent(
-                    ClimaLand.Canopy.canopy_temperature(
-                        canopy.energy,
-                        canopy,
-                        Y,
-                        p,
-                        t0,
+            @test all(
+                Array(
+                    parent(
+                        ClimaLand.Canopy.canopy_temperature(
+                            canopy.energy,
+                            canopy,
+                            Y,
+                            p,
+                            t0,
+                        ),
                     ),
-                ),
-            ) == [T_sfc]
+                ) .== [T_sfc],
+            )
             @test ClimaLand.Canopy.canopy_temperature(
                 canopy.energy,
                 canopy,
@@ -533,12 +550,25 @@ end
 
 @testset "Canopy software pipes with energy model" begin
     for FT in (Float32, Float64)
-        domain = Point(; z_sfc = FT(0.0))
-        for g1 in (FT(790), fill(FT(790), domain.space.surface))
-            # create new field with constant value everywhere
-            Vcmax25 = fill(FT(9e-5), domain.space.surface)
+        domain = ClimaLand.Domains.SphericalSurface(;
+            radius = FT(100.0),
+            nelements = 10,
+            npolynomial = 1,
+        )
+        # create a field with both 1.0s and 0.0s
+        mechanism_field = ClimaCore.Fields.Field(FT, domain.space.surface)
+        ClimaCore.Fields.set!(
+            x -> x.coordinates.lat > 0 ? 0.0 : 1.0,
+            mechanism_field,
+        )
+        # create one case where parameters are spatially varying and one where not
+        g1_cases = (FT(790), fill(FT(790), domain.space.surface))
+        Vcmax25_cases = (FT(9e-5), fill(FT(9e-5), domain.space.surface))
+        mechanism_cases = (FT(1), mechanism_field)
+        zipped = zip(g1_cases, Vcmax25_cases, mechanism_cases)
+        for (g1, Vcmax25, is_c3) in zipped
             RTparams = BeerLambertParameters(FT)
-            photosynthesis_params = FarquharParameters(FT, C3(); Vcmax25)
+            photosynthesis_params = FarquharParameters(FT, is_c3; Vcmax25)
             stomatal_g_params = MedlynConductanceParameters(FT; g1)
 
             stomatal_model = MedlynConductanceModel{FT}(stomatal_g_params)
@@ -723,7 +753,7 @@ end
 
             # Check that structure of Y is value (will error if not)
             @test !isnothing(zero(Y))
-            @test propertynames(p) == (:canopy, :drivers)
+            @test propertynames(p) == (:canopy, :dss_buffer_2d, :drivers)
             for component in ClimaLand.Canopy.canopy_components(canopy)
                 # Only hydraulics has a prognostic variable
                 if component == :hydraulics
@@ -786,10 +816,23 @@ end
 
 @testset "Zero LAI;" begin
     for FT in (Float32, Float64)
-        domain = Point(; z_sfc = FT(0.0))
-        for g1 in (FT(790), fill(FT(790), domain.space.surface))
-            # create new field with constant value everywhere
-            Vcmax25 = fill(FT(9e-5), domain.space.surface)
+        domain = ClimaLand.Domains.SphericalSurface(;
+            radius = FT(100.0),
+            nelements = 10,
+            npolynomial = 1,
+        )
+        # create a field with both 1.0s and 0.0s
+        mechanism_field = ClimaCore.Fields.Field(FT, domain.space.surface)
+        ClimaCore.Fields.set!(
+            x -> x.coordinates.lat > 0 ? 0.0 : 1.0,
+            mechanism_field,
+        )
+        # create one case where parameters are spatially varying and one where not
+        g1_cases = (FT(790), fill(FT(790), domain.space.surface))
+        Vcmax25_cases = (FT(9e-5), fill(FT(9e-5), domain.space.surface))
+        mechanism_cases = (FT(1), mechanism_field)
+        zipped = zip(g1_cases, Vcmax25_cases, mechanism_cases)
+        for (g1, Vcmax25, is_c3) in zipped
             BeerLambertparams = BeerLambertParameters(FT)
             # TwoStreamModel parameters
             Ω = FT(0.69)
@@ -811,7 +854,7 @@ end
                 τ_NIR_leaf,
                 G_Function,
             )
-            photosynthesis_params = FarquharParameters(FT, C3(); Vcmax25)
+            photosynthesis_params = FarquharParameters(FT, is_c3; Vcmax25)
             stomatal_g_params = MedlynConductanceParameters(FT; g1)
 
             stomatal_model = MedlynConductanceModel{FT}(stomatal_g_params)
diff --git a/test/standalone/Vegetation/plant_hydraulics_test.jl b/test/standalone/Vegetation/plant_hydraulics_test.jl
index afbc7457c1..481b78770f 100644
--- a/test/standalone/Vegetation/plant_hydraulics_test.jl
+++ b/test/standalone/Vegetation/plant_hydraulics_test.jl
@@ -115,7 +115,8 @@ for FT in (Float32, Float64)
 
         AR_params = AutotrophicRespirationParameters(FT)
         RTparams = BeerLambertParameters(FT)
-        photosynthesis_params = FarquharParameters(FT, C3())
+        is_c3 = FT(1) # set the photosynthesis mechanism to C3
+        photosynthesis_params = FarquharParameters(FT, is_c3)
         stomatal_g_params = MedlynConductanceParameters(FT)
 
         AR_model = AutotrophicRespirationModel{FT}(AR_params)
@@ -403,7 +404,8 @@ for FT in (Float32, Float64)
 
         AR_params = AutotrophicRespirationParameters(FT)
         RTparams = BeerLambertParameters(FT)
-        photosynthesis_params = FarquharParameters(FT, C3())
+        is_c3 = FT(1) # set the photosynthesis mechanism to C3
+        photosynthesis_params = FarquharParameters(FT, is_c3)
         stomatal_g_params = MedlynConductanceParameters(FT)
 
         AR_model = AutotrophicRespirationModel{FT}(AR_params)
diff --git a/test/standalone/Vegetation/test_bigleaf_parameterizations.jl b/test/standalone/Vegetation/test_bigleaf_parameterizations.jl
index a5e5e69a23..45893cf504 100644
--- a/test/standalone/Vegetation/test_bigleaf_parameterizations.jl
+++ b/test/standalone/Vegetation/test_bigleaf_parameterizations.jl
@@ -39,7 +39,7 @@ for FT in (Float32, Float64)
         @test typeof(Jmax) == FT
         @test typeof(Vcmax) == FT
         params = OptimalityFarquharParameters(FT)
-        @test params.mechanism == C3()
+        @test params.is_c3 ≈ 1.0
         model = OptimalityFarquharModel(params)
         @test ClimaLand.auxiliary_vars(model) == (:An, :GPP, :Rd, :Vcmax25)
         @test ClimaLand.auxiliary_types(model) == (FT, FT, FT, FT)
@@ -76,7 +76,8 @@ for FT in (Float32, Float64)
         ARparams = AutotrophicRespirationParameters(FT)
         RTparams = BeerLambertParameters(FT)
         RT = BeerLambertModel{FT}(RTparams)
-        photosynthesisparams = FarquharParameters(FT, C3())
+        is_c3 = FT(1) # set the photosynthesis mechanism to C3
+        photosynthesisparams = FarquharParameters(FT, is_c3)
         stomatal_g_params = MedlynConductanceParameters(FT)
 
         LAI = FT(5.0) # m2 (leaf) m-2 (ground)
@@ -165,7 +166,7 @@ for FT in (Float32, Float64)
         @test intercellular_co2(ca, FT(1), m_t) == FT(1)
 
         Ac = rubisco_assimilation(
-            photosynthesisparams.mechanism,
+            photosynthesisparams.is_c3,
             Vcmax,
             ci,
             Γstar,
@@ -202,13 +203,7 @@ for FT in (Float32, Float64)
             )
         )
 
-        Aj =
-            light_assimilation.(
-                Ref(photosynthesisparams.mechanism),
-                J,
-                ci,
-                Γstar,
-            )
+        Aj = light_assimilation.(Ref(photosynthesisparams.is_c3), J, ci, Γstar)
         @test all(@.(Aj == J * (ci - Γstar) / (4 * (ci + 2 * Γstar))))
         β = moisture_stress(
             p_l,
@@ -220,8 +215,9 @@ for FT in (Float32, Float64)
             1 + exp(photosynthesisparams.sc * (p_l - photosynthesisparams.pc))
         )
         #    C4 tests
+        is_c3 = 0.0
         @test rubisco_assimilation(
-            C4(),
+            is_c3,
             Vcmax,
             ci,
             Γstar,
@@ -229,7 +225,7 @@ for FT in (Float32, Float64)
             Ko,
             photosynthesisparams.oi,
         ) == Vcmax
-        @test light_assimilation(C4(), J, ci, Γstar) == J
+        @test light_assimilation(is_c3, J, ci, Γstar) == J
 
         Rd = dark_respiration(
             photosynthesisparams.Vcmax25,