diff --git a/NEWS.md b/NEWS.md
index 47615846ce..3d3f7003bd 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -3,6 +3,7 @@ ClimaLand.jl Release Notes
 
 main
 --------
+- PR[#690] Use the soil parameters in creating the biogeochemistry SoilMet driver for consistency.
 
 v0.13.0
 --------
diff --git a/docs/tutorials/integrated/soil_canopy_tutorial.jl b/docs/tutorials/integrated/soil_canopy_tutorial.jl
index aed1c5bb6e..ed52d28940 100644
--- a/docs/tutorials/integrated/soil_canopy_tutorial.jl
+++ b/docs/tutorials/integrated/soil_canopy_tutorial.jl
@@ -174,7 +174,7 @@ soil_model_type = Soil.EnergyHydrology{FT}
 # The domain is defined similarly to the soil domain described above.
 soilco2_type = Soil.Biogeochemistry.SoilCO2Model{FT}
 
-soilco2_ps = SoilCO2ModelParameters(FT; ν = soil_ν);
+soilco2_ps = SoilCO2ModelParameters(FT);
 
 # soil microbes args
 Csom = (z, t) -> eltype(z)(5); # kg C m⁻³, this is a guess, not measured at the site
@@ -185,18 +185,11 @@ soilco2_sources = (MicrobeProduction{FT}(),);
 
 soilco2_boundary_conditions = (; top = soilco2_top_bc, bottom = soilco2_bot_bc);
 
-soilco2_drivers = Soil.Biogeochemistry.SoilDrivers(
-    Soil.Biogeochemistry.PrognosticMet{FT}(),
-    Soil.Biogeochemistry.PrescribedSOC{FT}(Csom),
-    atmos,
-);
-
 soilco2_args = (;
     boundary_conditions = soilco2_boundary_conditions,
     sources = soilco2_sources,
     domain = soil_domain,
     parameters = soilco2_ps,
-    drivers = soilco2_drivers,
 );
 
 # Next we need to set up the [`CanopyModel`](https://clima.github.io/ClimaLand.jl/dev/APIs/canopy/Canopy/#Canopy-Model-Structs).
@@ -306,7 +299,7 @@ canopy_model_args = (; parameters = shared_params, domain = canopy_domain);
 # atmospheric and radiative flux conditions from the observations at the Ozark
 # site as was done in the previous tutorial.
 
-land_input = (atmos = atmos, radiation = radiation)
+land_input = (atmos = atmos, radiation = radiation, soil_organic_carbon = Csom)
 
 land = SoilCanopyModel{FT}(;
     soilco2_type = soilco2_type,
diff --git a/experiments/benchmarks/land.jl b/experiments/benchmarks/land.jl
index 85cd4d727f..013ae3abaf 100644
--- a/experiments/benchmarks/land.jl
+++ b/experiments/benchmarks/land.jl
@@ -400,10 +400,7 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
     z0_b = FT(0.1) * z0_m
 
 
-    soilco2_ps = Soil.Biogeochemistry.SoilCO2ModelParameters(
-        FT;
-        ν = 1.0,# INCORRECT!
-    )
+    soilco2_ps = Soil.Biogeochemistry.SoilCO2ModelParameters(FT)
 
     soil_args = (domain = domain, parameters = soil_params)
     soil_model_type = Soil.EnergyHydrology{FT}
@@ -422,18 +419,11 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
     soilco2_boundary_conditions =
         (; top = soilco2_top_bc, bottom = soilco2_bot_bc)
 
-    soilco2_drivers = Soil.Biogeochemistry.SoilDrivers(
-        Soil.Biogeochemistry.PrognosticMet{FT}(),
-        Soil.Biogeochemistry.PrescribedSOC{FT}(Csom),
-        atmos,
-    )
-
     soilco2_args = (;
         boundary_conditions = soilco2_boundary_conditions,
         sources = soilco2_sources,
         domain = domain,
         parameters = soilco2_ps,
-        drivers = soilco2_drivers,
     )
 
     # Now we set up the canopy model, which we set up by component:
@@ -538,7 +528,12 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
     )
 
     # Integrated plant hydraulics and soil model
-    land_input = (atmos = atmos, radiation = radiation, runoff = runoff_model)
+    land_input = (
+        atmos = atmos,
+        radiation = radiation,
+        runoff = runoff_model,
+        soil_organic_carbon = Csom,
+    )
     land = SoilCanopyModel{FT}(;
         soilco2_type = soilco2_type,
         soilco2_args = soilco2_args,
diff --git a/experiments/integrated/fluxnet/ozark_pft.jl b/experiments/integrated/fluxnet/ozark_pft.jl
index 22d412e0c7..3ccde05c20 100644
--- a/experiments/integrated/fluxnet/ozark_pft.jl
+++ b/experiments/integrated/fluxnet/ozark_pft.jl
@@ -135,10 +135,7 @@ soil_model_type = Soil.EnergyHydrology{FT}
 # Soil microbes model
 soilco2_type = Soil.Biogeochemistry.SoilCO2Model{FT}
 
-soilco2_ps = SoilCO2ModelParameters(
-    FT;
-    ν = soil_ν, # same as soil
-)
+soilco2_ps = SoilCO2ModelParameters(FT)
 
 # soil microbes args
 Csom = (z, t) -> eltype(z)(5.0)
@@ -150,18 +147,11 @@ soilco2_sources = (MicrobeProduction{FT}(),)
 
 soilco2_boundary_conditions = (; top = soilco2_top_bc, bottom = soilco2_bot_bc)
 
-soilco2_drivers = Soil.Biogeochemistry.SoilDrivers(
-    Soil.Biogeochemistry.PrognosticMet{FT}(),
-    Soil.Biogeochemistry.PrescribedSOC{FT}(Csom),
-    atmos,
-)
-
 soilco2_args = (;
     boundary_conditions = soilco2_boundary_conditions,
     sources = soilco2_sources,
     domain = soil_domain,
     parameters = soilco2_ps,
-    drivers = soilco2_drivers,
 )
 
 # Now we set up the canopy model, which we set up by component:
@@ -242,7 +232,7 @@ shared_params = SharedCanopyParameters{FT, typeof(earth_param_set)}(
 canopy_model_args = (; parameters = shared_params, domain = canopy_domain)
 
 # Integrated plant hydraulics and soil model
-land_input = (atmos = atmos, radiation = radiation)
+land_input = (atmos = atmos, radiation = radiation, soil_organic_carbon = Csom)
 land = SoilCanopyModel{FT}(;
     soilco2_type = soilco2_type,
     soilco2_args = soilco2_args,
diff --git a/experiments/integrated/fluxnet/run_fluxnet.jl b/experiments/integrated/fluxnet/run_fluxnet.jl
index 6b89cd17e0..d802db5b5f 100644
--- a/experiments/integrated/fluxnet/run_fluxnet.jl
+++ b/experiments/integrated/fluxnet/run_fluxnet.jl
@@ -94,10 +94,7 @@ soil_model_type = Soil.EnergyHydrology{FT}
 # Soil microbes model
 soilco2_type = Soil.Biogeochemistry.SoilCO2Model{FT}
 
-soilco2_ps = SoilCO2ModelParameters(
-    FT;
-    ν = soil_ν, # same as soil
-)
+soilco2_ps = SoilCO2ModelParameters(FT)
 
 # soil microbes args
 Csom = (z, t) -> eltype(z)(5.0)
@@ -109,18 +106,11 @@ soilco2_sources = (MicrobeProduction{FT}(),)
 
 soilco2_boundary_conditions = (; top = soilco2_top_bc, bottom = soilco2_bot_bc)
 
-soilco2_drivers = Soil.Biogeochemistry.SoilDrivers(
-    Soil.Biogeochemistry.PrognosticMet{FT}(),
-    Soil.Biogeochemistry.PrescribedSOC{FT}(Csom),
-    atmos,
-)
-
 soilco2_args = (;
     boundary_conditions = soilco2_boundary_conditions,
     sources = soilco2_sources,
     domain = soil_domain,
     parameters = soilco2_ps,
-    drivers = soilco2_drivers,
 )
 
 # Now we set up the canopy model, which we set up by component:
@@ -199,7 +189,7 @@ shared_params = SharedCanopyParameters{FT, typeof(earth_param_set)}(
 canopy_model_args = (; parameters = shared_params, domain = canopy_domain)
 
 # Integrated plant hydraulics and soil model
-land_input = (atmos = atmos, radiation = radiation)
+land_input = (atmos = atmos, radiation = radiation, soil_organic_carbon = Csom)
 land = SoilCanopyModel{FT}(;
     soilco2_type = soilco2_type,
     soilco2_args = soilco2_args,
diff --git a/experiments/integrated/global/global_parameters.jl b/experiments/integrated/global/global_parameters.jl
index 79e3c5c485..56fe6fb824 100644
--- a/experiments/integrated/global/global_parameters.jl
+++ b/experiments/integrated/global/global_parameters.jl
@@ -188,8 +188,4 @@ z0_m = FT(0.13) * h_canopy
 z0_b = FT(0.1) * z0_m
 
 
-soilco2_ps = Soil.Biogeochemistry.SoilCO2ModelParameters(
-    FT;
-    ν = 1.0,# INCORRECT! This should be the same as the soil porosity, but
-    # currently, SoilCO2 does not support spatially varying parameters
-)
+soilco2_ps = Soil.Biogeochemistry.SoilCO2ModelParameters(FT)
diff --git a/experiments/integrated/global/global_soil_canopy.jl b/experiments/integrated/global/global_soil_canopy.jl
index adc3424b09..ef2bb9c0d9 100644
--- a/experiments/integrated/global/global_soil_canopy.jl
+++ b/experiments/integrated/global/global_soil_canopy.jl
@@ -190,18 +190,11 @@ soilco2_sources = (Soil.Biogeochemistry.MicrobeProduction{FT}(),)
 
 soilco2_boundary_conditions = (; top = soilco2_top_bc, bottom = soilco2_bot_bc)
 
-soilco2_drivers = Soil.Biogeochemistry.SoilDrivers(
-    Soil.Biogeochemistry.PrognosticMet{FT}(),
-    Soil.Biogeochemistry.PrescribedSOC{FT}(Csom),
-    atmos,
-)
-
 soilco2_args = (;
     boundary_conditions = soilco2_boundary_conditions,
     sources = soilco2_sources,
     domain = domain,
     parameters = soilco2_ps,
-    drivers = soilco2_drivers,
 )
 
 # Now we set up the canopy model, which we set up by component:
@@ -295,7 +288,12 @@ canopy_model_args = (;
 )
 
 # Integrated plant hydraulics and soil model
-land_input = (atmos = atmos, radiation = radiation, runoff = runoff_model)
+land_input = (
+    atmos = atmos,
+    radiation = radiation,
+    runoff = runoff_model,
+    soil_organic_carbon = Csom,
+)
 land = SoilCanopyModel{FT}(;
     soilco2_type = soilco2_type,
     soilco2_args = soilco2_args,
diff --git a/experiments/integrated/performance/conservation/ozark_conservation_setup.jl b/experiments/integrated/performance/conservation/ozark_conservation_setup.jl
index 74983ef2f8..b40f34c80c 100644
--- a/experiments/integrated/performance/conservation/ozark_conservation_setup.jl
+++ b/experiments/integrated/performance/conservation/ozark_conservation_setup.jl
@@ -79,7 +79,7 @@ soil_model_type = Soil.EnergyHydrology{FT}
 # Soil microbes model
 soilco2_type = Soil.Biogeochemistry.SoilCO2Model{FT}
 
-soilco2_ps = SoilCO2ModelParameters(FT; ν = soil_ν)
+soilco2_ps = SoilCO2ModelParameters(FT)
 
 # soil microbes args
 Csom = (z, t) -> eltype(z)(5.0)
@@ -91,18 +91,11 @@ soilco2_sources = (MicrobeProduction{FT}(),)
 
 soilco2_boundary_conditions = (; top = soilco2_top_bc, bottom = soilco2_bot_bc)
 
-soilco2_drivers = Soil.Biogeochemistry.SoilDrivers(
-    Soil.Biogeochemistry.PrognosticMet{FT}(),
-    Soil.Biogeochemistry.PrescribedSOC{FT}(Csom),
-    atmos,
-)
-
 soilco2_args = (;
     boundary_conditions = soilco2_boundary_conditions,
     sources = soilco2_sources,
     domain = soil_domain,
     parameters = soilco2_ps,
-    drivers = soilco2_drivers,
 )
 
 # Now we set up the canopy model, which we set up by component:
@@ -180,7 +173,7 @@ shared_params = SharedCanopyParameters{FT, typeof(earth_param_set)}(
 canopy_model_args = (; parameters = shared_params, domain = canopy_domain)
 
 # Integrated plant hydraulics and soil model
-land_input = (atmos = atmos, radiation = radiation)
+land_input = (atmos = atmos, radiation = radiation, soil_organic_carbon = Csom)
 land = SoilCanopyModel{FT}(;
     soilco2_type = soilco2_type,
     soilco2_args = soilco2_args,
diff --git a/experiments/integrated/performance/profile_allocations.jl b/experiments/integrated/performance/profile_allocations.jl
index 95b114da4d..e9b898383c 100644
--- a/experiments/integrated/performance/profile_allocations.jl
+++ b/experiments/integrated/performance/profile_allocations.jl
@@ -226,7 +226,7 @@ soil_model_type = Soil.EnergyHydrology{FT}
 # Soil microbes model
 soilco2_type = Soil.Biogeochemistry.SoilCO2Model{FT}
 
-soilco2_ps = SoilCO2ModelParameters(FT; ν = soil_ν)
+soilco2_ps = SoilCO2ModelParameters(FT)
 
 # soil microbes args
 Csom = (z, t) -> eltype(z)(5.0)
@@ -238,18 +238,11 @@ soilco2_sources = (MicrobeProduction{FT}(),)
 
 soilco2_boundary_conditions = (; top = soilco2_top_bc, bottom = soilco2_bot_bc)
 
-soilco2_drivers = Soil.Biogeochemistry.SoilDrivers(
-    Soil.Biogeochemistry.PrognosticMet{FT}(),
-    Soil.Biogeochemistry.PrescribedSOC{FT}(Csom),
-    atmos,
-)
-
 soilco2_args = (;
     boundary_conditions = soilco2_boundary_conditions,
     sources = soilco2_sources,
     domain = soil_domain,
     parameters = soilco2_ps,
-    drivers = soilco2_drivers,
 )
 
 # Now we set up the canopy model, which we set up by component:
@@ -315,7 +308,7 @@ shared_params = SharedCanopyParameters{FT, typeof(earth_param_set)}(
 canopy_model_args = (; parameters = shared_params, domain = canopy_domain)
 
 # Integrated plant hydraulics and soil model
-land_input = (atmos = atmos, radiation = radiation)
+land_input = (atmos = atmos, radiation = radiation, soil_organic_carbon = Csom)
 land = SoilCanopyModel{FT}(;
     soilco2_type = soilco2_type,
     soilco2_args = soilco2_args,
diff --git a/experiments/standalone/Biogeochemistry/experiment.jl b/experiments/standalone/Biogeochemistry/experiment.jl
index 44903485a7..a45fc4559f 100644
--- a/experiments/standalone/Biogeochemistry/experiment.jl
+++ b/experiments/standalone/Biogeochemistry/experiment.jl
@@ -71,7 +71,7 @@ for (FT, tf) in ((Float32, 2 * dt), (Float64, tf))
     # Make biogeochemistry model args
     Csom = (z, t) -> eltype(z)(5.0)
 
-    co2_parameters = Soil.Biogeochemistry.SoilCO2ModelParameters(FT; ν = 0.556)
+    co2_parameters = Soil.Biogeochemistry.SoilCO2ModelParameters(FT)
     C = FT(100)
 
     co2_top_bc = Soil.Biogeochemistry.SoilCO2StateBC((p, t) -> 0.0)
@@ -102,23 +102,17 @@ for (FT, tf) in ((Float32, 2 * dt), (Float64, tf))
         earth_param_set;
         c_co2 = TimeVaryingInput(atmos_co2),
     )
-
-    soil_drivers = Soil.Biogeochemistry.SoilDrivers(
-        Soil.Biogeochemistry.PrognosticMet{FT}(),
-        Soil.Biogeochemistry.PrescribedSOC{FT}(Csom),
-        atmos,
-    )
-
     soilco2_args = (;
         boundary_conditions = co2_boundary_conditions,
         sources = co2_sources,
         domain = lsm_domain,
         parameters = co2_parameters,
-        drivers = soil_drivers,
     )
 
     # Create integrated model instance
+    land_args = (atmos = atmos, soil_organic_carbon = Csom)
     model = LandSoilBiogeochemistry{FT}(;
+        land_args = land_args,
         soil_args = soil_args,
         soilco2_args = soilco2_args,
     )
diff --git a/ext/CreateParametersExt.jl b/ext/CreateParametersExt.jl
index 4a2770ae6b..514495916b 100644
--- a/ext/CreateParametersExt.jl
+++ b/ext/CreateParametersExt.jl
@@ -484,10 +484,6 @@ SoilCO2ModelParameters(::Type{FT}; kwargs...) where {FT <: AbstractFloat} =
     SoilCO2ModelParameters(CP.create_toml_dict(FT); kwargs...)
 
 function SoilCO2ModelParameters(toml_dict::CP.AbstractTOMLDict; kwargs...)
-    # These parameters have defaults that should not go in ClimaParams
-    θ_a100 = 0.1816
-    b = 4.547
-
     name_map = (;
         :CO2_diffusion_coefficient => :D_ref,
         :soil_C_substrate_diffusivity => :D_liq,
@@ -504,8 +500,6 @@ function SoilCO2ModelParameters(toml_dict::CP.AbstractTOMLDict; kwargs...)
     earth_param_set = LP.LandParameters(toml_dict)
     return SoilCO2ModelParameters{FT, typeof(earth_param_set)}(;
         earth_param_set,
-        θ_a100,
-        b,
         parameters...,
         kwargs...,
     )
diff --git a/lib/ClimaLandSimulations/experiments/ozark.jl b/lib/ClimaLandSimulations/experiments/ozark.jl
index 115ee5cca0..a2e6a2679c 100644
--- a/lib/ClimaLandSimulations/experiments/ozark.jl
+++ b/lib/ClimaLandSimulations/experiments/ozark.jl
@@ -4,7 +4,7 @@ using ClimaLandSimulations.Fluxnet
 # default parameters, except for custom parameter hetero_resp.b
 sv_test, sol_test, Y_test, p_test = run_fluxnet(
     "US-MOz";
-    params = ozark_default_params(; hetero_resp = hetero_resp_ozark(; b = 2)),
+    params = ozark_default_params(; hetero_resp = hetero_resp_ozark()),
 )
 
 # defaults, except start time
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 c1328a0c9a..a7fc8a915e 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
@@ -31,9 +31,7 @@ function harvard_default_params(;
 end
 
 function hetero_resp_harvard(;
-    θ_a100 = FT(0.1816),
     D_ref = FT(1.39e-5),
-    b = FT(4.547),
     D_liq = FT(3.17),
     α_sx = FT(194e3),
     Ea_sx = FT(61e3),
@@ -44,9 +42,7 @@ function hetero_resp_harvard(;
     p_sx = FT(0.024),
 )
     return HeteroRespP(
-        θ_a100,
         D_ref,
-        b,
         D_liq,
         α_sx,
         Ea_sx,
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 79068adaf4..db9110af70 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
@@ -31,9 +31,7 @@ function ozark_default_params(;
 end
 
 function hetero_resp_ozark(;
-    θ_a100 = FT(0.1816),
     D_ref = FT(1.39e-5),
-    b = FT(4.547),
     D_liq = FT(3.17),
     α_sx = FT(194e3),
     Ea_sx = FT(61e3),
@@ -44,9 +42,7 @@ function hetero_resp_ozark(;
     p_sx = FT(0.024),
 )
     return HeteroRespP(
-        θ_a100,
         D_ref,
-        b,
         D_liq,
         α_sx,
         Ea_sx,
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 bb64a87392..311552a439 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
@@ -31,9 +31,7 @@ function niwotridge_default_params(;
 end
 
 function hetero_resp_niwotridge(;
-    θ_a100 = FT(0.1816),
     D_ref = FT(1.39e-5),
-    b = FT(4.547),
     D_liq = FT(3.17),
     α_sx = FT(194e3),
     Ea_sx = FT(61e3),
@@ -44,9 +42,7 @@ function hetero_resp_niwotridge(;
     p_sx = FT(0.024),
 )
     return HeteroRespP(
-        θ_a100,
         D_ref,
-        b,
         D_liq,
         α_sx,
         Ea_sx,
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 7cbc37c9d1..7228ccab2a 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
@@ -31,9 +31,7 @@ function vairaranch_default_params(;
 end
 
 function hetero_resp_vairaranch(;
-    θ_a100 = FT(0.1816),
     D_ref = FT(1.39e-5),
-    b = FT(4.547),
     D_liq = FT(3.17),
     α_sx = FT(194e3),
     Ea_sx = FT(61e3),
@@ -44,9 +42,7 @@ function hetero_resp_vairaranch(;
     p_sx = FT(0.024),
 )
     return HeteroRespP(
-        θ_a100,
         D_ref,
-        b,
         D_liq,
         α_sx,
         Ea_sx,
diff --git a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_utilities/make_parameters.jl b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_utilities/make_parameters.jl
index 965495a5be..9f8913180f 100644
--- a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_utilities/make_parameters.jl
+++ b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_utilities/make_parameters.jl
@@ -27,10 +27,8 @@ end
 # or, another example, parameters calculated from other parameters (e.g., κ_solid)
 # the structs below only contain unique params, that are given as number (not computed FTwhere)
 
-struct HeteroRespP # differs from SoilCO2ModelParameters because of porosity
-    θ_a100::FT
+struct HeteroRespP
     D_ref::FT
-    b::FT
     D_liq::FT
     α_sx::FT
     Ea_sx::FT
diff --git a/lib/ClimaLandSimulations/src/Fluxnet/run_fluxnet.jl b/lib/ClimaLandSimulations/src/Fluxnet/run_fluxnet.jl
index 72ad5333ce..e5920a746c 100644
--- a/lib/ClimaLandSimulations/src/Fluxnet/run_fluxnet.jl
+++ b/lib/ClimaLandSimulations/src/Fluxnet/run_fluxnet.jl
@@ -46,10 +46,7 @@ function run_fluxnet(
 
     soilco2_ps = SoilCO2ModelParameters(
         FT;
-        ν = params.soil.ν,
-        θ_a100 = params.hetero_resp.θ_a100,
         D_ref = params.hetero_resp.D_ref,
-        b = params.hetero_resp.b,
         D_liq = params.hetero_resp.D_liq,
         # DAMM
         α_sx = params.hetero_resp.α_sx,
@@ -72,18 +69,11 @@ function run_fluxnet(
     soilco2_boundary_conditions =
         (; top = soilco2_top_bc, bottom = CO2 = soilco2_bot_bc)
 
-    soilco2_drivers = Soil.Biogeochemistry.SoilDrivers(
-        Soil.Biogeochemistry.PrognosticMet{FT}(),
-        Soil.Biogeochemistry.PrescribedSOC{FT}(Csom),
-        drivers.atmos,
-    )
-
     soilco2_args = (;
         boundary_conditions = soilco2_boundary_conditions,
         sources = soilco2_sources,
         domain = soil_domain,
         parameters = soilco2_ps,
-        drivers = soilco2_drivers,
     )
 
     # Now we set up the canopy model, which we set up by component:
@@ -218,7 +208,11 @@ function run_fluxnet(
         (; parameters = shared_params, domain = domain.canopy_domain)
 
     # Integrated plant hydraulics and soil model
-    land_input = (atmos = drivers.atmos, radiation = drivers.radiation)
+    land_input = (
+        atmos = drivers.atmos,
+        radiation = drivers.radiation,
+        soil_organic_carbon = Csom,
+    )
     land = SoilCanopyModel{FT}(;
         soilco2_type = soilco2_type,
         soilco2_args = soilco2_args,
diff --git a/src/integrated/soil_canopy_model.jl b/src/integrated/soil_canopy_model.jl
index 5c68aa3948..04577975d6 100644
--- a/src/integrated/soil_canopy_model.jl
+++ b/src/integrated/soil_canopy_model.jl
@@ -88,7 +88,7 @@ function SoilCanopyModel{FT}(;
     MM <: Soil.Biogeochemistry.SoilCO2Model{FT},
 }
 
-    (; atmos, radiation) = land_args
+    (; atmos, radiation, soil_organic_carbon) = land_args
     # These should always be set by the constructor.
     sources = (RootExtraction{FT}(), Soil.PhaseChange{FT}())
     if :runoff ∈ propertynames(land_args)
@@ -116,7 +116,7 @@ function SoilCanopyModel{FT}(;
     )
 
     transpiration = Canopy.PlantHydraulics.DiagnosticTranspiration{FT}()
-    soil_driver = PrognosticSoil{typeof(soil.parameters.PAR_albedo)}(
+    canopy_soil_driver = PrognosticSoil{typeof(soil.parameters.PAR_albedo)}(
         soil.parameters.PAR_albedo,
         soil.parameters.NIR_albedo,
     )
@@ -142,7 +142,7 @@ function SoilCanopyModel{FT}(;
             energy = canopy_component_types.energy(
                 canopy_component_args.energy.parameters,
             ),
-            soil_driver = soil_driver,
+            soil_driver = canopy_soil_driver,
             atmos = atmos,
             radiation = radiation,
             canopy_model_args...,
@@ -165,18 +165,20 @@ function SoilCanopyModel{FT}(;
                 transpiration = transpiration,
                 canopy_component_args.hydraulics...,
             ),
-            soil_driver = soil_driver,
+            soil_driver = canopy_soil_driver,
             atmos = atmos,
             radiation = radiation,
             canopy_model_args...,
         )
     end
 
-    soilco2 = Soil.Biogeochemistry.SoilCO2Model{FT}(; soilco2_args...)
-
-    if !(soilco2_args.drivers.met isa PrognosticMet)
-        throw(AssertionError("Must be of type PrognosticMet."))
-    end
+    co2_prognostic_soil = Soil.Biogeochemistry.PrognosticMet(soil.parameters)
+    soilco2_drivers = Soil.Biogeochemistry.SoilDrivers(
+        co2_prognostic_soil,
+        Soil.Biogeochemistry.PrescribedSOC{FT}(soil_organic_carbon),
+        atmos,
+    )
+    soilco2 = soilco2_type(; soilco2_args..., drivers = soilco2_drivers)
 
     return SoilCanopyModel{FT, typeof(soilco2), typeof(soil), typeof(canopy)}(
         soilco2,
diff --git a/src/integrated/soil_energy_hydrology_biogeochemistry.jl b/src/integrated/soil_energy_hydrology_biogeochemistry.jl
index 3172bd2666..e25d351f6c 100644
--- a/src/integrated/soil_energy_hydrology_biogeochemistry.jl
+++ b/src/integrated/soil_energy_hydrology_biogeochemistry.jl
@@ -38,35 +38,56 @@ Additional arguments, like parameters and driving atmospheric data, can be passe
 in as needed.
 """
 function LandSoilBiogeochemistry{FT}(;
+    land_args::NamedTuple,
     soil_args::NamedTuple = (;),
     soilco2_args::NamedTuple = (;),
 ) where {FT}
 
+    (; atmos, soil_organic_carbon) = land_args
     soil = Soil.EnergyHydrology{FT}(;
         soil_args..., # soil_args must have sources, boundary_conditions, domain, parameters
     )
-
-    soilco2 = Soil.Biogeochemistry.SoilCO2Model{FT}(; soilco2_args...)
-    if !(soilco2_args.drivers.met isa PrognosticMet)
-        throw(
-            AssertionError(
-                "To run with a prognostic soil energy and hydrology model, the met driver must be of type PrognosticMet.",
-            ),
-        )
-    end
-    if soil_args.parameters.ν != soilco2_args.parameters.ν
-        throw(
-            AssertionError(
-                "The provided porosity must be the same in both model parameter structures.",
-            ),
-        )
-    end
-
+    prognostic_soil = Soil.Biogeochemistry.PrognosticMet(soil.parameters)
+    soil_co2_drivers = Soil.Biogeochemistry.SoilDrivers(
+        prognostic_soil,
+        Soil.Biogeochemistry.PrescribedSOC{FT}(soil_organic_carbon),
+        atmos,
+    )
+    soilco2 = Soil.Biogeochemistry.SoilCO2Model{FT}(;
+        soilco2_args...,
+        drivers = soil_co2_drivers,
+    )
     args = (soil, soilco2)
     return LandSoilBiogeochemistry{FT, typeof.(args)...}(args...)
 end
 
-struct PrognosticMet{FT} <: Soil.Biogeochemistry.AbstractSoilDriver end
+struct PrognosticMet{FT, F <: Union{AbstractFloat, ClimaCore.Fields.Field}} <:
+       Soil.Biogeochemistry.AbstractSoilDriver
+    "Soil porosity (m³ m⁻³)"
+    ν::F
+    "Air-filled porosity at soil water potential of -100 cm H₂O (~ 10 Pa)"
+    θ_a100::F
+    "Absolute value of the slope of the line relating log(ψ) versus log(S) (unitless)"
+    b::F
+    function PrognosticMet(
+        soil_params::Soil.EnergyHydrologyParameters{FT},
+    ) where {FT}
+        ν = soil_params.ν
+        θ_r = soil_params.θ_r
+        hcm = soil_params.hydrology_cm
+        F = typeof(ν)
+        if F <: AbstractFloat
+            θ_a100 =
+                Soil.inverse_matric_potential(hcm, -FT(1)) * (ν - θ_r) + θ_r
+            b = Soil.approximate_ψ_S_slope(hcm)
+        else
+            θ_a100 =
+                @. Soil.inverse_matric_potential(hcm, -FT(1)) * (ν - θ_r) + θ_r
+            b = @. Soil.approximate_ψ_S_slope(hcm)
+        end
+        return new{FT, F}(ν, θ_a100, b)
+    end
+end
 
 """
     soil_temperature(driver::PrognosticSoil, p, Y, t, z)
@@ -96,6 +117,6 @@ function ClimaLand.get_drivers(model::LandSoilBiogeochemistry)
     }
         return (bc.atmos, bc.radiation)
     else
-        return (model.soilco2.driver.atmos, nothing)
+        return (model.soilco2.drivers.atmos, nothing)
     end
 end
diff --git a/src/standalone/Soil/Biogeochemistry/Biogeochemistry.jl b/src/standalone/Soil/Biogeochemistry/Biogeochemistry.jl
index e149532e6f..0e3b2c00e6 100644
--- a/src/standalone/Soil/Biogeochemistry/Biogeochemistry.jl
+++ b/src/standalone/Soil/Biogeochemistry/Biogeochemistry.jl
@@ -39,17 +39,13 @@ export SoilCO2ModelParameters,
     SoilCO2ModelParameters{FT <: AbstractFloat, PSE}
 
 A struct for storing parameters of the `SoilCO2Model`.
+
+All of these parameters are currently treated as global constants.
 $(DocStringExtensions.FIELDS)
 """
 Base.@kwdef struct SoilCO2ModelParameters{FT <: AbstractFloat, PSE}
-    "Soil porosity (m³ m⁻³)"
-    ν::FT
-    "Air-filled porosity at soil water potential of -100 cm H₂O (~ 10 Pa)"
-    θ_a100::FT
     "Diffusion coefficient for CO₂ in air at standard temperature and pressure (m² s⁻¹)"
     D_ref::FT
-    "Absolute value of the slope of the line relating log(ψ) versus log(θ) (unitless)"
-    b::FT
     "Diffusivity of soil C substrate in liquid (unitless)"
     D_liq::FT
     # DAMM
@@ -97,7 +93,7 @@ struct SoilCO2Model{FT, PS, D, BC, S, DT} <:
     "A tuple of sources, each of type AbstractSource"
     sources::S
     "Drivers"
-    driver::DT
+    drivers::DT
 end
 
 
@@ -312,28 +308,56 @@ without a prognostic soil model.
 
 $(DocStringExtensions.FIELDS)
 """
-struct PrescribedMet{FT, F1 <: Function, F2 <: Function} <: AbstractSoilDriver
+struct PrescribedMet{
+    FT,
+    F1 <: Function,
+    F2 <: Function,
+    F <: Union{AbstractFloat, ClimaCore.Fields.Field},
+} <: AbstractSoilDriver
     "The temperature of the soil, of the form f(z::FT,t) where FT <: AbstractFloat"
     temperature::F1
     "Soil moisture, of the form f(z::FT,t) FT <: AbstractFloat"
     volumetric_liquid_fraction::F2
+    "Soil porosity (m³ m⁻³)"
+    ν::F
+    "Air-filled porosity at soil water potential of -100 cm H₂O (~ 10 Pa)"
+    θ_a100::F
+    "Absolute value of the slope of the line relating log(ψ) versus log(S) (unitless)"
+    b::F
 end
 
 function PrescribedMet{FT}(
     temperature::Function,
     volumetric_liquid_fraction::Function,
-) where {FT <: AbstractFloat}
+    ν::F,
+    θ_r::F,
+    hcm::CF,
+) where {FT <: AbstractFloat, F, CF}
+    if F <: AbstractFloat
+        θ_a100 =
+            ClimaLand.Soil.inverse_matric_potential(hcm, -FT(1)) * (ν - θ_r) +
+            θ_r
+        b = ClimaLand.Soil.approximate_ψ_S_slope(hcm)
+    else
+        θ_a100 = @. ClimaLand.Soil.inverse_matric_potential(hcm, -FT(1)) *
+           (ν - θ_r) + θ_r
+        b = @. ClimaLand.Soil.approximate_ψ_S_slope(hcm)
+    end
+
     return PrescribedMet{
         FT,
         typeof(temperature),
         typeof(volumetric_liquid_fraction),
+        F,
     }(
         temperature,
         volumetric_liquid_fraction,
+        ν,
+        θ_a100,
+        b,
     )
 end
 
-
 """
     PrescribedSOC <: AbstractSoilDriver
 
@@ -383,14 +407,6 @@ function soil_SOM_C(driver::PrescribedSOC, p, Y, t, z)
     return driver.soil_organic_carbon.(z, t)
 end
 
-"""
-    air_pressure(driver::PrescribedAtmosphere, t)
-
-Returns the prescribed air pressure at the top boundary condition at time (t).
-"""
-function air_pressure(driver::PrescribedAtmosphere, p, Y, t) # not sure if/why p and Y are needed?
-    return p.drivers.P
-end
 
 """
     make_update_aux(model::SoilCO2Model)
@@ -404,14 +420,18 @@ function ClimaLand.make_update_aux(model::SoilCO2Model)
     function update_aux!(p, Y, t)
         params = model.parameters
         z = model.domain.fields.z
-        T_soil = soil_temperature(model.driver.met, p, Y, t, z)
-        θ_l = soil_moisture(model.driver.met, p, Y, t, z)
-        Csom = soil_SOM_C(model.driver.soc, p, Y, t, z)
-        P_sfc = air_pressure(model.driver.atmos, p, Y, t)
+        T_soil = soil_temperature(model.drivers.met, p, Y, t, z)
+        θ_l = soil_moisture(model.drivers.met, p, Y, t, z)
+        Csom = soil_SOM_C(model.drivers.soc, p, Y, t, z)
+        P_sfc = p.drivers.P
         θ_w = θ_l
+        ν = model.drivers.met.ν
+        θ_a100 = model.drivers.met.θ_a100
+        b = model.drivers.met.b
 
-        p.soilco2.D .= co2_diffusivity.(T_soil, θ_w, P_sfc, Ref(params))
-        p.soilco2.Sm .= microbe_source.(T_soil, θ_l, Csom, Ref(params))
+        p.soilco2.D .=
+            co2_diffusivity.(T_soil, θ_w, P_sfc, θ_a100, b, ν, params)
+        p.soilco2.Sm .= microbe_source.(T_soil, θ_l, Csom, ν, params)
     end
     return update_aux!
 end
@@ -581,9 +601,11 @@ function ClimaLand.boundary_flux(
 end
 
 function ClimaLand.get_drivers(model::SoilCO2Model)
-    return (model.driver.atmos, nothing)
+    return (model.drivers.atmos, nothing)
 end
 
+Base.broadcastable(ps::SoilCO2ModelParameters) = tuple(ps)
+
 include("./co2_parameterizations.jl")
 
 end # module
diff --git a/src/standalone/Soil/Biogeochemistry/co2_parameterizations.jl b/src/standalone/Soil/Biogeochemistry/co2_parameterizations.jl
index 719c54983e..cf4db435d3 100644
--- a/src/standalone/Soil/Biogeochemistry/co2_parameterizations.jl
+++ b/src/standalone/Soil/Biogeochemistry/co2_parameterizations.jl
@@ -1,4 +1,3 @@
-
 export volumetric_air_content, co2_diffusivity, microbe_source
 
 
@@ -6,6 +5,7 @@ export volumetric_air_content, co2_diffusivity, microbe_source
     microbe_source(T_soil::FT,
                    θ_l::FT,
                    Csom::FT,
+                   ν::FT,
                    params::SoilCO2ModelParameters{FT}
                    ) where {FT}
 
@@ -16,9 +16,10 @@ function microbe_source(
     T_soil::FT,
     θ_l::FT,
     Csom::FT,
+    ν::FT,
     params::SoilCO2ModelParameters{FT},
 ) where {FT}
-    (; α_sx, Ea_sx, kM_sx, kM_o2, ν, D_liq, p_sx, D_oa, O2_a, earth_param_set) =
+    (; α_sx, Ea_sx, kM_sx, kM_o2, D_liq, p_sx, D_oa, O2_a, earth_param_set) =
         params
     R = FT(LP.gas_constant(earth_param_set))
     Vmax = α_sx * exp(-Ea_sx / (R * T_soil)) # Maximum potential rate of respiration
@@ -33,18 +34,15 @@ end
 
 """
     volumetric_air_content(θ_w::FT,
-                           params::SoilCO2ModelParameters{FT}
+                           ν::FT,
                            ) where {FT}
 
 Computes the volumetric air content (`θ_a`) in the soil,
 which is related to the total soil porosity (`ν`) and
 volumetric soil water content (`θ_w = θ_l+θ_i`).
 """
-function volumetric_air_content(
-    θ_w::FT,
-    params::SoilCO2ModelParameters{FT},
-) where {FT}
-    θ_a = params.ν - θ_w
+function volumetric_air_content(θ_w::FT, ν::FT) where {FT}
+    θ_a = ν - θ_w
     return θ_a
 end
 
@@ -53,7 +51,10 @@ end
                     T_soil::FT,
                     θ_w::FT,
                     P_sfc::FT,
-                    params::SoilCO2ModelParameters{FT}
+                    θ_a100::FT,
+                    b::FT,
+                    ν::FT,
+                    params::SoilCO2ModelParameters{FT},
                     ) where {FT}
 
 Computes the diffusivity of CO₂ within the soil (D).
@@ -64,17 +65,23 @@ values of `T_ref` and `P_ref` (273 K and 101325 Pa). Here, `θ_a` is the
 volumetric air content and `θ_a100` is the volumetric air content
 at a soil water potential of
 100cm, and b is the pore size distribution of the soil.
+
+This parameterization is from Ryan et al., GMD 11, 1909-1928, 2018,
+https://doi.org/10.5194/gmd-11-1909-2018.
 """
 function co2_diffusivity(
     T_soil::FT,
     θ_w::FT,
     P_sfc::FT,
+    θ_a100::FT,
+    b::FT,
+    ν::FT,
     params::SoilCO2ModelParameters{FT},
 ) where {FT}
-    (; D_ref, θ_a100, b, ν, earth_param_set) = params
+    (; D_ref, earth_param_set) = params
     T_ref = FT(LP.T_0(earth_param_set))
     P_ref = FT(LP.P_ref(earth_param_set))
-    θ_a = volumetric_air_content(θ_w, params)
+    θ_a = volumetric_air_content(θ_w, ν)
     D0 = D_ref * (T_soil / T_ref)^FT(1.75) * (P_ref / P_sfc)
     D =
         D0 *
diff --git a/src/standalone/Soil/soil_hydrology_parameterizations.jl b/src/standalone/Soil/soil_hydrology_parameterizations.jl
index 43ccc68590..b677c19bd5 100644
--- a/src/standalone/Soil/soil_hydrology_parameterizations.jl
+++ b/src/standalone/Soil/soil_hydrology_parameterizations.jl
@@ -66,6 +66,19 @@ function inverse_matric_potential(cm::vanGenuchten{FT}, ψ::FT) where {FT}
 end
 
 
+"""
+     approximate_ψ_S_slope(cm::vanGenuchten)
+
+An estimate of the slope of the absolute value of the logψ-logS curve.
+Following Lehmann, Assouline, and Or (2008), we linearize the ψ(S) curve about the inflection point (where d²ψ/dS² = 0, at S = (1+m)^(-m)). 
+"""
+function approximate_ψ_S_slope(cm::vanGenuchten)
+    m = cm.m
+    n = cm.n
+    return (1 + m) / (n * m * m)
+end
+
+
 
 """
     pressure_head(
@@ -167,6 +180,18 @@ function inverse_matric_potential(cm::BrooksCorey{FT}, ψ::FT) where {FT}
     return S
 end
 
+
+"""
+     approximate_ψ_S_slope(cm::BrooksCorey)
+
+The slope of the logψ-logS curve for the Brooks and Corey
+model.
+"""
+function approximate_ψ_S_slope(cm::BrooksCorey)
+    return 1 / cm.c
+end
+
+
 """
    dψdϑ(cm::BrooksCorey{FT}, ϑ, ν, θ_r, S_s)
 
diff --git a/test/integrated/soil_energy_hydrology_biogeochemistry.jl b/test/integrated/soil_energy_hydrology_biogeochemistry.jl
index 9025ff6240..bb7431ef17 100644
--- a/test/integrated/soil_energy_hydrology_biogeochemistry.jl
+++ b/test/integrated/soil_energy_hydrology_biogeochemistry.jl
@@ -4,7 +4,7 @@ import ClimaComms
 using ClimaCore
 import ClimaParams
 using ClimaLand
-using ClimaLand.Domains: Column
+using ClimaLand.Domains: Column, HybridBox
 using ClimaLand.Soil
 using ClimaLand.Soil.Biogeochemistry
 using Dates
@@ -65,8 +65,7 @@ for FT in (Float32, Float64)
         # Make biogeochemistry model args
         Csom = (z, t) -> eltype(z)(5.0)
 
-        co2_parameters =
-            Soil.Biogeochemistry.SoilCO2ModelParameters(FT; ν = 0.556)
+        co2_parameters = Soil.Biogeochemistry.SoilCO2ModelParameters(FT)
         C = FT(4)
         co2_top_bc = Soil.Biogeochemistry.SoilCO2StateBC((p, t) -> C)
         co2_bot_bc = Soil.Biogeochemistry.SoilCO2StateBC((p, t) -> C)
@@ -96,25 +95,22 @@ for FT in (Float32, Float64)
             earth_param_set;
             c_co2 = TimeVaryingInput(atmos_co2),
         )
-
-        soil_drivers = Soil.Biogeochemistry.SoilDrivers(
-            Soil.Biogeochemistry.PrognosticMet{FT}(),
-            Soil.Biogeochemistry.PrescribedSOC{FT}(Csom),
-            atmos,
-        )
         soilco2_args = (;
             boundary_conditions = co2_boundary_conditions,
             sources = co2_sources,
             domain = lsm_domain,
             parameters = co2_parameters,
-            drivers = soil_drivers,
         )
 
         # Create integrated model instance
+        land_args = (; atmos = atmos, soil_organic_carbon = Csom)
         model = LandSoilBiogeochemistry{FT}(;
+            land_args = land_args,
             soil_args = soil_args,
             soilco2_args = soilco2_args,
         )
+        @test model.soilco2.drivers.met.ν == model.soil.parameters.ν
+        @test model.soilco2.drivers.met isa ClimaLand.PrognosticMet
         Y, p, coords = initialize(model)
         @test propertynames(p.drivers) ==
               (:P_liq, :P_snow, :T, :P, :u, :q, :c_co2, :thermal_state)
@@ -154,7 +150,7 @@ for FT in (Float32, Float64)
         set_initial_cache!(p, Y, t0)
 
         @test p.soil.T ≈ Soil.Biogeochemistry.soil_temperature(
-            model.soilco2.driver.met,
+            model.soilco2.drivers.met,
             p,
             Y,
             t0,
@@ -163,7 +159,7 @@ for FT in (Float32, Float64)
         @test all(
             parent(
                 Soil.Biogeochemistry.soil_SOM_C(
-                    model.soilco2.driver.soc,
+                    model.soilco2.drivers.soc,
                     p,
                     Y,
                     t0,
@@ -172,95 +168,99 @@ for FT in (Float32, Float64)
             ) .== FT(5.0),
         )
         @test p.soil.θ_l ≈ Soil.Biogeochemistry.soil_moisture(
-            model.soilco2.driver.met,
+            model.soilco2.drivers.met,
             p,
             Y,
             t0,
             z,
         )
+    end
+    @testset "PrognosticMet, FT = $FT" begin
+        earth_param_set = LP.LandParameters(FT)
+        zmax = FT(0)
+        zmin = FT(-1)
+        nelems = 10
+        xmin = ymin = zmin
+        xmax = ymax = zmax
+        col = Column(; zlim = (zmin, zmax), nelements = nelems)
+        box = HybridBox(;
+            xlim = (xmin, xmax),
+            ylim = (ymin, ymax),
+            zlim = (zmin, zmax),
+            nelements = (nelems, nelems, nelems),
+            npolynomial = 2,
+        )
 
-        try
-            co2_parameters =
-                Soil.Biogeochemistry.SoilCO2ModelParameters(FT; ν = 0.2)
-            soil_drivers = Soil.Biogeochemistry.SoilDrivers(
-                Soil.Biogeochemistry.PrognosticMet{FT}(),
-                Soil.Biogeochemistry.PrescribedSOC{FT}(Csom),
-                atmos,
-            )
-            soilco2_args = (;
-                boundary_conditions = co2_boundary_conditions,
-                sources = co2_sources,
-                domain = lsm_domain,
-                parameters = co2_parameters,
-                drivers = soil_drivers,
-            )
-
-            # Create integrated model instance
-            model = LandSoilBiogeochemistry{FT}(;
-                soil_args = soil_args,
-                soilco2_args = soilco2_args,
-            )
-            Y, p, coords = initialize(model)
-
-            function init_soil_2!(Y, z, params)
-                ν = params.ν
-                FT = eltype(Y.soil.ϑ_l)
-                Y.soil.ϑ_l .= FT(0.33)
-                Y.soil.θ_i .= FT(0.0)
-                T = FT(279.85)
-                ρc_s = Soil.volumetric_heat_capacity(
-                    FT(0.33),
-                    FT(0.0),
-                    params.ρc_ds,
-                    params.earth_params_set,
-                )
-                Y.soil.ρe_int .=
-                    Soil.volumetric_internal_energy.(
-                        FT(0.0),
-                        ρc_s,
-                        T,
-                        params.earth_param_set,
-                    )
-            end
+        ν = FT(0.556)
+        K_sat = FT(0.0443 / 3600 / 100)
+        S_s = FT(1e-3)
+        vg_n = FT(2.0)
+        vg_α = FT(2.6)
+        hydrology_cm = vanGenuchten{FT}(; α = vg_α, n = vg_n)
+        θ_r = FT(0.1)
+        ν_ss_om = FT(0.0)
+        ν_ss_quartz = FT(1.0)
+        ν_ss_gravel = FT(0.0)
 
-            function init_co2_2!(Y, C_0)
-                Y.soilco2.C .= C_0
-            end
 
-            z = coords.subsurface.z
-            init_soil_2!(Y, z, model.soil.parameters)
-            init_co2_2!(Y, C)
-            t0 = FT(0.0)
-            set_initial_cache! = make_set_initial_cache(model)
-            set_initial_cache!(p, Y, t0)
+        soil_ps_col = Soil.EnergyHydrologyParameters(
+            FT;
+            ν,
+            ν_ss_om,
+            ν_ss_quartz,
+            ν_ss_gravel,
+            K_sat,
+            hydrology_cm,
+            S_s,
+            θ_r,
+        )
 
-            @test p.soil.T ≈ Soil.Biogeochemistry.soil_temperature(
-                model.soilco2.driver.met,
-                p,
-                Y,
-                t0,
-                z,
-            )
-            soil_drivers = Soil.Biogeochemistry.SoilDrivers(
-                Soil.Biogeochemistry.PrescribedMet{FT}(Csom, Csom),
-                Soil.Biogeochemistry.PrescribedSOC{FT}(Csom),
-                atmos,
-            )
-            soilco2_args = (;
-                boundary_conditions = co2_boundary_conditions,
-                sources = co2_sources,
-                domain = lsm_domain,
-                parameters = co2_parameters,
-                drivers = soil_drivers,
+        zero_field = ClimaCore.Fields.zeros(box.space.subsurface)
+        vg_fields_to_hcm_field(α::FT, n::FT) where {FT} =
+            ClimaLand.Soil.vanGenuchten{FT}(;
+                @NamedTuple{α::FT, n::FT}((α, n))...,
             )
+        hydrology_cm_field =
+            vg_fields_to_hcm_field.(vg_α .+ zero_field, vg_n .+ zero_field)
+        ν_field = ν .+ zero_field
+        ν_ss_om_field = ν_ss_om .+ zero_field
+        ν_ss_quartz_field = ν_ss_quartz .+ zero_field
+        ν_ss_gravel_field = ν_ss_gravel .+ zero_field
+        K_sat_field = K_sat .+ zero_field
+        S_s_field = S_s .+ zero_field
+        θ_r_field = θ_r .+ zero_field
 
-            # Create integrated model instance
-            model = LandSoilBiogeochemistry{FT}(;
-                soil_args = soil_args,
-                soilco2_args = soilco2_args,
-            )
-        catch err
-            @test isa(err, AssertionError)
-        end
+        soil_ps_box = Soil.EnergyHydrologyParameters(
+            FT;
+            ν = ν_field,
+            ν_ss_om = ν_ss_om_field,
+            ν_ss_quartz = ν_ss_quartz_field,
+            ν_ss_gravel = ν_ss_gravel_field,
+            K_sat = K_sat_field,
+            hydrology_cm = hydrology_cm_field,
+            S_s = S_s_field,
+            θ_r = θ_r_field,
+        )
+
+        met_col = ClimaLand.PrognosticMet(soil_ps_col)
+        met_box = ClimaLand.PrognosticMet(soil_ps_box)
+        @test met_box.ν == soil_ps_box.ν
+        @test met_col.ν == soil_ps_col.ν
+        @test met_box.θ_a100 == @. Soil.inverse_matric_potential(
+            soil_ps_box.hydrology_cm,
+            -FT(1),
+        ) * (soil_ps_box.ν - soil_ps_box.θ_r) + soil_ps_box.θ_r
+        @test met_col.θ_a100 ==
+              Soil.inverse_matric_potential(soil_ps_col.hydrology_cm, -FT(1)) *
+              (soil_ps_col.ν - soil_ps_col.θ_r) + soil_ps_col.θ_r
+        @test met_box.b ==
+              @. Soil.approximate_ψ_S_slope(soil_ps_box.hydrology_cm)
+        @test met_col.b == Soil.approximate_ψ_S_slope(soil_ps_col.hydrology_cm)
+        vg_m = 1 - 1 / vg_n
+        @test Soil.approximate_ψ_S_slope(soil_ps_col.hydrology_cm) ==
+              (1 + vg_m) / (vg_n * vg_m^2)
+        @test Soil.approximate_ψ_S_slope(
+            BrooksCorey{FT}(; c = FT(1.0), ψb = FT(1.0)),
+        ) == FT(1)
     end
 end
diff --git a/test/standalone/Soil/Biogeochemistry/biogeochemistry_module.jl b/test/standalone/Soil/Biogeochemistry/biogeochemistry_module.jl
index 5bd4f52da8..6373c6a7c0 100644
--- a/test/standalone/Soil/Biogeochemistry/biogeochemistry_module.jl
+++ b/test/standalone/Soil/Biogeochemistry/biogeochemistry_module.jl
@@ -19,7 +19,7 @@ for FT in (Float32, Float64)
         θ_i = (z, t) -> eltype(z)(0)
         Csom = (z, t) -> eltype(z)(5.0) # 3 [kg C m-3] soil organic C content at depth z
         D_ref = FT(0.0)
-        parameters = SoilCO2ModelParameters(FT; ν = 0.556, D_ref)
+        parameters = SoilCO2ModelParameters(FT; D_ref)
 
         nelems = 50 # number of layers in the vertical
         zmin = FT(-1) # 0 to 1 m depth
@@ -54,9 +54,13 @@ for FT in (Float32, Float64)
             earth_param_set;
             c_co2 = TimeVaryingInput(atmos_co2),
         )
-
+        ν = FT(0.6)
+        θ_r = FT(0.0)
+        α = FT(0.1)
+        n = FT(2)
+        hcm = ClimaLand.Soil.vanGenuchten{FT}(; α = α, n = n)
         soil_drivers = SoilDrivers(
-            PrescribedMet{FT}(T_soil, θ_l),
+            PrescribedMet{FT}(T_soil, θ_l, ν, θ_r, hcm),
             PrescribedSOC{FT}(Csom),
             atmos,
         )
@@ -89,7 +93,7 @@ for FT in (Float32, Float64)
         θ_i = (z, t) -> eltype(z)(0.0)
         Csom = (z, t) -> eltype(z)(5.0) # 3 [kg C m-3] soil organic C content at depth z
 
-        parameters = SoilCO2ModelParameters(FT; ν = 0.556)
+        parameters = SoilCO2ModelParameters(FT)
         C = FT(4)
         nelems = 50 # number of layers in the vertical
         zmin = FT(-1) # 0 to 1 m depth
@@ -124,9 +128,13 @@ for FT in (Float32, Float64)
             earth_param_set;
             c_co2 = TimeVaryingInput(atmos_co2),
         )
-
+        ν = FT(0.6)
+        θ_r = FT(0.0)
+        α = FT(0.1)
+        n = FT(2)
+        hcm = ClimaLand.Soil.vanGenuchten{FT}(; α = α, n = n)
         soil_drivers = SoilDrivers(
-            PrescribedMet{FT}(T_soil, θ_l),
+            PrescribedMet{FT}(T_soil, θ_l, ν, θ_r, hcm),
             PrescribedSOC{FT}(Csom),
             atmos, # need to create some functions
         )
diff --git a/test/standalone/Soil/Biogeochemistry/co2_parameterizations.jl b/test/standalone/Soil/Biogeochemistry/co2_parameterizations.jl
index cee5eda5ff..38fa487577 100644
--- a/test/standalone/Soil/Biogeochemistry/co2_parameterizations.jl
+++ b/test/standalone/Soil/Biogeochemistry/co2_parameterizations.jl
@@ -22,25 +22,27 @@ for FT in (Float32, Float64)
         T_ref = FT(LP.T_0(earth_param_set))
         R = FT(LP.gas_constant(earth_param_set))
 
-        parameters = SoilCO2ModelParameters(FT; ν)
+        parameters = SoilCO2ModelParameters(FT)
 
         # Test that parameterizations functions are working properly
-        θ_a = volumetric_air_content(θ_w, parameters)
-        @test θ_a == parameters.ν - θ_w
+        θ_a = volumetric_air_content(θ_w, ν)
+        @test θ_a == ν - θ_w
         @test typeof(θ_a) == FT
+        θ_a100 = FT(0.3)
+        b = FT(0.45)
 
-        D = co2_diffusivity(T_soil, θ_w, P_sfc, parameters)
+        D = co2_diffusivity(T_soil, θ_w, P_sfc, θ_a100, b, ν, parameters)
         @test D ==
               (
                   parameters.D_ref *
                   (T_soil / T_ref)^FT(1.75) *
                   (FT(LP.P_ref(parameters.earth_param_set)) / P_sfc)
               ) *
-              (FT(2)parameters.θ_a100^FT(3) + FT(0.04)parameters.θ_a100) *
-              (θ_a / parameters.θ_a100)^(FT(2) + FT(3) / parameters.b)
+              (FT(2)θ_a100^FT(3) + FT(0.04) * θ_a100) *
+              (θ_a / θ_a100)^(FT(2) + FT(3) / b)
         @test typeof(D) == FT
 
-        ms = microbe_source(T_soil, θ_l, Csom, parameters)
+        ms = microbe_source(T_soil, θ_l, Csom, ν, parameters)
         # check that the value is correct
         (; p_sx, D_liq, kM_o2, D_oa, kM_sx, O2_a, α_sx, Ea_sx) = parameters
         MM_sx =