diff --git a/docs/src/APIs/SoilBiogeochemistry.md b/docs/src/APIs/SoilBiogeochemistry.md
index 3de0d106c2..123a73be58 100644
--- a/docs/src/APIs/SoilBiogeochemistry.md
+++ b/docs/src/APIs/SoilBiogeochemistry.md
@@ -25,7 +25,6 @@ ClimaLand.Soil.Biogeochemistry.AtmosCO2StateBC
 ClimaLand.Soil.Biogeochemistry.AbstractSoilDriver
 ClimaLand.Soil.Biogeochemistry.SoilDrivers
 ClimaLand.Soil.Biogeochemistry.PrescribedMet
-ClimaLand.Soil.Biogeochemistry.PrescribedSOC
 ```
 
 ## Functions of State
diff --git a/docs/src/APIs/shared_utilities.md b/docs/src/APIs/shared_utilities.md
index 56ca3fac0d..bc75c2bb35 100644
--- a/docs/src/APIs/shared_utilities.md
+++ b/docs/src/APIs/shared_utilities.md
@@ -71,6 +71,7 @@ ClimaLand.get_drivers
 ClimaLand.PrescribedAtmosphere
 ClimaLand.PrescribedPrecipitation
 ClimaLand.PrescribedRadiativeFluxes
+ClimaLand.PrescribedSoilOrganicCarbon
 ClimaLand.CoupledAtmosphere
 ClimaLand.CoupledRadiativeFluxes
 ClimaLand.AbstractAtmosphericDrivers
diff --git a/docs/tutorials/integrated/soil_canopy_tutorial.jl b/docs/tutorials/integrated/soil_canopy_tutorial.jl
index aed1c5bb6e..e01062286b 100644
--- a/docs/tutorials/integrated/soil_canopy_tutorial.jl
+++ b/docs/tutorials/integrated/soil_canopy_tutorial.jl
@@ -177,7 +177,7 @@ soilco2_type = Soil.Biogeochemistry.SoilCO2Model{FT}
 soilco2_ps = SoilCO2ModelParameters(FT; ν = soil_ν);
 
 # soil microbes args
-Csom = (z, t) -> eltype(z)(5); # kg C m⁻³, this is a guess, not measured at the site
+Csom = ClimaLand.PrescribedSoilOrganicCarbon{FT}(TimeVaryingInput((t) -> 5))
 
 soilco2_top_bc = Soil.Biogeochemistry.AtmosCO2StateBC()
 soilco2_bot_bc = Soil.Biogeochemistry.SoilCO2StateBC((p, t) -> 0.0);
@@ -187,7 +187,7 @@ 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),
+    Csom,
     atmos,
 );
 
diff --git a/experiments/benchmarks/land.jl b/experiments/benchmarks/land.jl
index 85cd4d727f..e2d137e903 100644
--- a/experiments/benchmarks/land.jl
+++ b/experiments/benchmarks/land.jl
@@ -412,7 +412,7 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
     soilco2_type = Soil.Biogeochemistry.SoilCO2Model{FT}
 
     # soil microbes args
-    Csom = (z, t) -> eltype(z)(5.0)
+    Csom = ClimaLand.PrescribedSoilOrganicCarbon{FT}(TimeVaryingInput((t) -> 5))
 
     # Set the soil CO2 BC to being atmospheric CO2
     soilco2_top_bc = Soil.Biogeochemistry.AtmosCO2StateBC()
@@ -424,7 +424,7 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
 
     soilco2_drivers = Soil.Biogeochemistry.SoilDrivers(
         Soil.Biogeochemistry.PrognosticMet{FT}(),
-        Soil.Biogeochemistry.PrescribedSOC{FT}(Csom),
+        Csom,
         atmos,
     )
 
diff --git a/experiments/integrated/fluxnet/ozark_pft.jl b/experiments/integrated/fluxnet/ozark_pft.jl
index 22d412e0c7..3d4e8ab6fa 100644
--- a/experiments/integrated/fluxnet/ozark_pft.jl
+++ b/experiments/integrated/fluxnet/ozark_pft.jl
@@ -140,8 +140,7 @@ soilco2_ps = SoilCO2ModelParameters(
     ν = soil_ν, # same as soil
 )
 
-# soil microbes args
-Csom = (z, t) -> eltype(z)(5.0)
+Csom = ClimaLand.PrescribedSoilOrganicCarbon{FT}(TimeVaryingInput((t) -> 5))
 
 # Set the soil CO2 BC to being atmospheric CO2
 soilco2_top_bc = Soil.Biogeochemistry.AtmosCO2StateBC()
@@ -152,7 +151,7 @@ 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),
+    Csom,
     atmos,
 )
 
diff --git a/experiments/integrated/fluxnet/run_fluxnet.jl b/experiments/integrated/fluxnet/run_fluxnet.jl
index 6b89cd17e0..8974f4b876 100644
--- a/experiments/integrated/fluxnet/run_fluxnet.jl
+++ b/experiments/integrated/fluxnet/run_fluxnet.jl
@@ -99,8 +99,7 @@ soilco2_ps = SoilCO2ModelParameters(
     ν = soil_ν, # same as soil
 )
 
-# soil microbes args
-Csom = (z, t) -> eltype(z)(5.0)
+Csom = ClimaLand.PrescribedSoilOrganicCarbon{FT}(TimeVaryingInput((t) -> 5))
 
 # Set the soil CO2 BC to being atmospheric CO2
 soilco2_top_bc = Soil.Biogeochemistry.AtmosCO2StateBC()
@@ -111,7 +110,7 @@ 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),
+    Csom,
     atmos,
 )
 
diff --git a/experiments/integrated/global/global_soil_canopy.jl b/experiments/integrated/global/global_soil_canopy.jl
index adc3424b09..682e3f2431 100644
--- a/experiments/integrated/global/global_soil_canopy.jl
+++ b/experiments/integrated/global/global_soil_canopy.jl
@@ -181,7 +181,7 @@ soil_model_type = Soil.EnergyHydrology{FT}
 soilco2_type = Soil.Biogeochemistry.SoilCO2Model{FT}
 
 # soil microbes args
-Csom = (z, t) -> eltype(z)(5.0)
+Csom = ClimaLand.PrescribedSoilOrganicCarbon{FT}(TimeVaryingInput((t) -> 5))
 
 # Set the soil CO2 BC to being atmospheric CO2
 soilco2_top_bc = Soil.Biogeochemistry.AtmosCO2StateBC()
@@ -192,7 +192,7 @@ 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),
+    Csom,
     atmos,
 )
 
diff --git a/experiments/integrated/performance/conservation/ozark_conservation_setup.jl b/experiments/integrated/performance/conservation/ozark_conservation_setup.jl
index 74983ef2f8..8accb16f39 100644
--- a/experiments/integrated/performance/conservation/ozark_conservation_setup.jl
+++ b/experiments/integrated/performance/conservation/ozark_conservation_setup.jl
@@ -82,7 +82,7 @@ soilco2_type = Soil.Biogeochemistry.SoilCO2Model{FT}
 soilco2_ps = SoilCO2ModelParameters(FT; ν = soil_ν)
 
 # soil microbes args
-Csom = (z, t) -> eltype(z)(5.0)
+Csom = ClimaLand.PrescribedSoilOrganicCarbon{FT}(TimeVaryingInput((t) -> 5))
 
 # Set the soil CO2 BC to being atmospheric CO2
 soilco2_top_bc = Soil.Biogeochemistry.AtmosCO2StateBC()
@@ -93,7 +93,7 @@ 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),
+    Csom,
     atmos,
 )
 
diff --git a/experiments/integrated/performance/profile_allocations.jl b/experiments/integrated/performance/profile_allocations.jl
index 95b114da4d..7582e69055 100644
--- a/experiments/integrated/performance/profile_allocations.jl
+++ b/experiments/integrated/performance/profile_allocations.jl
@@ -229,7 +229,7 @@ soilco2_type = Soil.Biogeochemistry.SoilCO2Model{FT}
 soilco2_ps = SoilCO2ModelParameters(FT; ν = soil_ν)
 
 # soil microbes args
-Csom = (z, t) -> eltype(z)(5.0)
+Csom = ClimaLand.PrescribedSoilOrganicCarbon{FT}(TimeVaryingInput((t) -> 5))
 
 # Set the soil CO2 BC to being atmospheric CO2
 soilco2_top_bc = Soil.Biogeochemistry.AtmosCO2StateBC()
@@ -240,7 +240,7 @@ 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),
+    Csom,
     atmos,
 )
 
diff --git a/experiments/standalone/Biogeochemistry/experiment.jl b/experiments/standalone/Biogeochemistry/experiment.jl
index 44903485a7..a2544a890e 100644
--- a/experiments/standalone/Biogeochemistry/experiment.jl
+++ b/experiments/standalone/Biogeochemistry/experiment.jl
@@ -69,7 +69,7 @@ for (FT, tf) in ((Float32, 2 * dt), (Float64, tf))
     )
 
     # Make biogeochemistry model args
-    Csom = (z, t) -> eltype(z)(5.0)
+    Csom = ClimaLand.PrescribedSoilOrganicCarbon{FT}(TimeVaryingInput((t) -> 5))
 
     co2_parameters = Soil.Biogeochemistry.SoilCO2ModelParameters(FT; ν = 0.556)
     C = FT(100)
@@ -105,7 +105,7 @@ for (FT, tf) in ((Float32, 2 * dt), (Float64, tf))
 
     soil_drivers = Soil.Biogeochemistry.SoilDrivers(
         Soil.Biogeochemistry.PrognosticMet{FT}(),
-        Soil.Biogeochemistry.PrescribedSOC{FT}(Csom),
+        Csom,
         atmos,
     )
 
diff --git a/lib/ClimaLandSimulations/src/Fluxnet/run_fluxnet.jl b/lib/ClimaLandSimulations/src/Fluxnet/run_fluxnet.jl
index 72ad5333ce..55de75208a 100644
--- a/lib/ClimaLandSimulations/src/Fluxnet/run_fluxnet.jl
+++ b/lib/ClimaLandSimulations/src/Fluxnet/run_fluxnet.jl
@@ -63,7 +63,7 @@ function run_fluxnet(
     )
 
     # soil microbes args
-    Csom = (z, t) -> eltype(z)(5.0)
+    Csom = ClimaLand.PrescribedSoilOrganicCarbon{FT}(TimeVaryingInput((t) -> 5))
 
     soilco2_top_bc = Soil.Biogeochemistry.AtmosCO2StateBC()
     soilco2_bot_bc = Soil.Biogeochemistry.SoilCO2FluxBC((p, t) -> 0.0) # no flux
@@ -74,7 +74,7 @@ function run_fluxnet(
 
     soilco2_drivers = Soil.Biogeochemistry.SoilDrivers(
         Soil.Biogeochemistry.PrognosticMet{FT}(),
-        Soil.Biogeochemistry.PrescribedSOC{FT}(Csom),
+        Csom,
         drivers.atmos,
     )
 
diff --git a/src/integrated/soil_canopy_model.jl b/src/integrated/soil_canopy_model.jl
index 5c68aa3948..c5e062b82f 100644
--- a/src/integrated/soil_canopy_model.jl
+++ b/src/integrated/soil_canopy_model.jl
@@ -590,5 +590,9 @@ end
 
 
 function ClimaLand.get_drivers(model::SoilCanopyModel)
-    return (model.canopy.atmos, model.canopy.radiation)
+    return (
+        model.canopy.atmos,
+        model.canopy.radiation,
+        model.soilco2.driver.soc,
+    )
 end
diff --git a/src/integrated/soil_energy_hydrology_biogeochemistry.jl b/src/integrated/soil_energy_hydrology_biogeochemistry.jl
index 3172bd2666..b7d80005c3 100644
--- a/src/integrated/soil_energy_hydrology_biogeochemistry.jl
+++ b/src/integrated/soil_energy_hydrology_biogeochemistry.jl
@@ -94,8 +94,8 @@ function ClimaLand.get_drivers(model::LandSoilBiogeochemistry)
         <:AbstractRadiativeDrivers,
         <:Soil.AbstractRunoffModel,
     }
-        return (bc.atmos, bc.radiation)
+        return (bc.atmos, bc.radiation, model.soilco2.driver.soc)
     else
-        return (model.soilco2.driver.atmos, nothing)
+        return (model.soilco2.driver.atmos, model.soilco2.driver.soc)
     end
 end
diff --git a/src/shared_utilities/drivers.jl b/src/shared_utilities/drivers.jl
index dcaba85e06..9581e1fc1d 100644
--- a/src/shared_utilities/drivers.jl
+++ b/src/shared_utilities/drivers.jl
@@ -11,6 +11,7 @@ export AbstractAtmosphericDrivers,
     AbstractRadiativeDrivers,
     PrescribedAtmosphere,
     PrescribedPrecipitation,
+    PrescribedSoilOrganicCarbon,
     CoupledAtmosphere,
     PrescribedRadiativeFluxes,
     CoupledRadiativeFluxes,
@@ -48,6 +49,21 @@ An abstract type of radiative drivers of land models.
 abstract type AbstractRadiativeDrivers{FT} <: AbstractClimaLandDrivers{FT} end
 
 
+"""
+     PrescribedSoilOrganicCarbon{FT}
+
+A type for prescribing soil organic carbon.
+$(DocStringExtensions.FIELDS)
+"""
+struct PrescribedSoilOrganicCarbon{FT, SOC <: AbstractTimeVaryingInput} <:
+       AbstractClimaLandDrivers{FT}
+    "Soil organic carbon, function of time and space: kg C/m^3"
+    soc::SOC
+end
+
+PrescribedSoilOrganicCarbon{FT}(soc) where {FT} =
+    PrescribedSoilOrganicCarbon{FT, typeof(soc)}(soc)
+
 """
     PrescribedAtmosphere{FT, CA, DT} <: AbstractAtmosphericDrivers{FT}
 
@@ -766,6 +782,28 @@ function initialize_drivers(r::PrescribedRadiativeFluxes{FT}, coords) where {FT}
     return vars.drivers
 end
 
+"""
+    initialize_drivers(r::PrescribedSoilOrganicCarbon{FT}, coords) where {FT}
+
+Creates and returns a NamedTuple for the `PrescribedSoilOrganicCarbon` driver,
+ with variable `soc`.
+"""
+function initialize_drivers(
+    r::PrescribedSoilOrganicCarbon{FT},
+    coords,
+) where {FT}
+    keys = (:soc,)
+    types = (FT,)
+    domain_names = (:subsurface,)
+    model_name = :drivers
+    # intialize_vars packages the variables as a named tuple,
+    # as part of a named tuple with `model_name` as the key.
+    # Here we just want the variable named tuple itself
+    vars =
+        ClimaLand.initialize_vars(keys, types, domain_names, coords, model_name)
+    return vars.drivers
+end
+
 """
     initialize_drivers(d::Union{AbstractClimaLandDrivers, Nothing}, coords)
 
@@ -869,3 +907,15 @@ function make_update_drivers(r::PrescribedRadiativeFluxes{FT}) where {FT}
     end
     return update_drivers!
 end
+
+"""
+    make_update_drivers(a::PrescribedSoilOrganicCarbon{FT}) where {FT}
+Creates and returns a function which updates the driver variables
+in the case of a PrescribedSoilOrganicCarbon.
+"""
+function make_update_drivers(a::PrescribedSoilOrganicCarbon{FT}) where {FT}
+    function update_drivers!(p, t)
+        evaluate!(p.drivers.soc, a.soc, t)
+    end
+    return update_drivers!
+end
diff --git a/src/standalone/Soil/Biogeochemistry/Biogeochemistry.jl b/src/standalone/Soil/Biogeochemistry/Biogeochemistry.jl
index e149532e6f..8cb630afa5 100644
--- a/src/standalone/Soil/Biogeochemistry/Biogeochemistry.jl
+++ b/src/standalone/Soil/Biogeochemistry/Biogeochemistry.jl
@@ -27,7 +27,7 @@ import ClimaLand:
 export SoilCO2ModelParameters,
     SoilCO2Model,
     PrescribedMet,
-    PrescribedSOC,
+    PrescribedSoilOrganicCarbon,
     MicrobeProduction,
     SoilCO2FluxBC,
     AtmosCO2StateBC,
@@ -290,7 +290,7 @@ $(DocStringExtensions.FIELDS)
 struct SoilDrivers{
     FT,
     MET <: AbstractSoilDriver,
-    SOC <: AbstractSoilDriver,
+    SOC <: PrescribedSoilOrganicCarbon{FT},
     ATM <: PrescribedAtmosphere{FT},
 }
     "Soil temperature and moisture drivers - Prescribed or Prognostic"
@@ -334,24 +334,6 @@ function PrescribedMet{FT}(
 end
 
 
-"""
-    PrescribedSOC <: AbstractSoilDriver
-
-A container which holds the prescribed function for soil organic carbon
-
-This is meant for use when running the biogeochemistry model without a soil
-organic carbon model.
-
-$(DocStringExtensions.FIELDS)
-"""
-struct PrescribedSOC{FT, F <: Function} <: AbstractSoilDriver
-    "Carbon content of soil organic matter, of the form f(z::FT, t) where FT <: AbstractFloat"
-    soil_organic_carbon::F
-end
-
-function PrescribedSOC{FT}(Csom) where {FT <: AbstractFloat}
-    return PrescribedSOC{FT, typeof(Csom)}(Csom)
-end
 
 """
     soil_temperature(driver::PrescribedMet, p, Y, t, z)
@@ -373,16 +355,6 @@ function soil_moisture(driver::PrescribedMet, p, Y, t, z)
     return driver.volumetric_liquid_fraction.(z, t)
 end
 
-"""
-    soil_som_C(driver::PrescribedSOC, p, Y, t, z)
-
-Returns the carbon soil organic matter (SOM) at location (z) and time (t) for the prescribed
-soil case.
-"""
-function soil_SOM_C(driver::PrescribedSOC, p, Y, t, z)
-    return driver.soil_organic_carbon.(z, t)
-end
-
 """
     air_pressure(driver::PrescribedAtmosphere, t)
 
@@ -406,7 +378,7 @@ function ClimaLand.make_update_aux(model::SoilCO2Model)
         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)
+        Csom = p.drivers.soc
         P_sfc = air_pressure(model.driver.atmos, p, Y, t)
         θ_w = θ_l
 
@@ -581,7 +553,7 @@ function ClimaLand.boundary_flux(
 end
 
 function ClimaLand.get_drivers(model::SoilCO2Model)
-    return (model.driver.atmos, nothing)
+    return (model.driver.atmos, model.driver.soc)
 end
 
 include("./co2_parameterizations.jl")
diff --git a/src/standalone/Soil/energy_hydrology.jl b/src/standalone/Soil/energy_hydrology.jl
index c58328bbdb..24ec8575b4 100644
--- a/src/standalone/Soil/energy_hydrology.jl
+++ b/src/standalone/Soil/energy_hydrology.jl
@@ -880,6 +880,6 @@ function ClimaLand.get_drivers(model::EnergyHydrology)
     }
         return (bc.atmos, bc.radiation)
     else
-        return (nothing, nothing)
+        return (nothing,)
     end
 end
diff --git a/src/standalone/Soil/rre.jl b/src/standalone/Soil/rre.jl
index aa890f17c6..4798065340 100644
--- a/src/standalone/Soil/rre.jl
+++ b/src/standalone/Soil/rre.jl
@@ -442,8 +442,8 @@ function ClimaLand.get_drivers(model::RichardsModel)
         <:PrescribedPrecipitation,
         <:AbstractRunoffModel,
     }
-        return (bc.precip, nothing)
+        return (bc.precip,)
     else
-        return (nothing, nothing)
+        return (nothing,)
     end
 end
diff --git a/test/integrated/soil_energy_hydrology_biogeochemistry.jl b/test/integrated/soil_energy_hydrology_biogeochemistry.jl
index 9025ff6240..8507b0cea7 100644
--- a/test/integrated/soil_energy_hydrology_biogeochemistry.jl
+++ b/test/integrated/soil_energy_hydrology_biogeochemistry.jl
@@ -63,7 +63,9 @@ for FT in (Float32, Float64)
         )
 
         # Make biogeochemistry model args
-        Csom = (z, t) -> eltype(z)(5.0)
+        Csom = ClimaLand.PrescribedSoilOrganicCarbon{FT}(
+            TimeVaryingInput((t) -> 5),
+        )
 
         co2_parameters =
             Soil.Biogeochemistry.SoilCO2ModelParameters(FT; ν = 0.556)
@@ -99,7 +101,7 @@ for FT in (Float32, Float64)
 
         soil_drivers = Soil.Biogeochemistry.SoilDrivers(
             Soil.Biogeochemistry.PrognosticMet{FT}(),
-            Soil.Biogeochemistry.PrescribedSOC{FT}(Csom),
+            Csom,
             atmos,
         )
         soilco2_args = (;
@@ -179,88 +181,5 @@ for FT in (Float32, Float64)
             z,
         )
 
-        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
-
-            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)
-
-            @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,
-            )
-
-            # Create integrated model instance
-            model = LandSoilBiogeochemistry{FT}(;
-                soil_args = soil_args,
-                soilco2_args = soilco2_args,
-            )
-        catch err
-            @test isa(err, AssertionError)
-        end
     end
 end
diff --git a/test/shared_utilities/drivers.jl b/test/shared_utilities/drivers.jl
index 4c7bcc4714..42983a405b 100644
--- a/test/shared_utilities/drivers.jl
+++ b/test/shared_utilities/drivers.jl
@@ -27,9 +27,9 @@ FT = Float32
     liquid_precip = TimeVaryingInput((t) -> -1.0)
     pp = ClimaLand.PrescribedPrecipitation{FT}(liquid_precip)
     coords = (; surface = [1, 2, 3])
-    @test ClimaLand.initialize_drivers(pp, nothing, coords) ==
+    @test ClimaLand.initialize_drivers((pp,), coords) ==
           NamedTuple{(:P_liq,)}((zeros(FT, 3),))
-    @test ClimaLand.initialize_drivers(nothing, nothing, coords) == (;)
+    @test ClimaLand.initialize_drivers((nothing, nothing), coords) == (;)
     pa_keys = (:P_liq, :P_snow, :T, :P, :u, :q, :c_co2)
     zero_thermal_state = map(
         _ -> ClimaCore.RecursiveApply.rzero(Thermodynamics.PhaseEquil{FT}),
@@ -38,13 +38,13 @@ FT = Float32
     pa_vals = ([zeros(FT, 3) for k in pa_keys]...,)
     all_pa_keys = (pa_keys..., :thermal_state)
     all_pa_vals = (pa_vals..., zero_thermal_state)
-    @test ClimaLand.initialize_drivers(pa, nothing, coords) ==
+    @test ClimaLand.initialize_drivers((pa,), coords) ==
           NamedTuple{all_pa_keys}(all_pa_vals)
     pr_keys = (:SW_d, :LW_d, :θs)
     pr_vals = ([zeros(FT, 3) for k in pr_keys]...,)
     all_papr_keys = (pa_keys..., :thermal_state, pr_keys...)
     all_papr_vals = (pa_vals..., zero_thermal_state, pr_vals...)
-    @test ClimaLand.initialize_drivers(pa, pr, coords) ==
+    @test ClimaLand.initialize_drivers((pa, pr), coords) ==
           NamedTuple{all_papr_keys}(all_papr_vals)
 end
 
@@ -111,13 +111,13 @@ end
         earth_param_set,
     )
     pr = ClimaLand.PrescribedRadiativeFluxes(FT, f, f, f)
-    coords = (; surface = [1])
-    p = (; drivers = ClimaLand.initialize_drivers(nothing, nothing, coords))
-    nothing_update! = ClimaLand.make_update_drivers(nothing, nothing)
+    coords = (; surface = [1], subsurface = [1, 2])
+    p = (; drivers = ClimaLand.initialize_drivers((nothing,), coords))
+    nothing_update! = ClimaLand.make_update_drivers((nothing,))
     nothing_update!(p, 0.0)
     @test p.drivers == (;)
-    p = (; drivers = ClimaLand.initialize_drivers(pa, nothing, coords))
-    atmos_only_update! = ClimaLand.make_update_drivers(pa, nothing)
+    p = (; drivers = ClimaLand.initialize_drivers((pa,), coords))
+    atmos_only_update! = ClimaLand.make_update_drivers((pa,))
     atmos_only_update!(p, 0.0)
     @test p.drivers.P_liq == [FT(10)]
     @test p.drivers.P_snow == [FT(10)]
@@ -127,8 +127,8 @@ end
     @test p.drivers.u == [FT(10)]
     @test p.drivers.c_co2 == [FT(4.2e-4)]
 
-    p = (; drivers = ClimaLand.initialize_drivers(pa, pr, coords))
-    update! = ClimaLand.make_update_drivers(pa, pr)
+    p = (; drivers = ClimaLand.initialize_drivers((pa, pr), coords))
+    update! = ClimaLand.make_update_drivers((pa, pr))
     update!(p, 0.0)
     @test p.drivers.P_liq == [FT(10)]
     @test p.drivers.P_snow == [FT(10)]
@@ -141,8 +141,8 @@ end
     @test p.drivers.LW_d == [FT(10)]
     @test p.drivers.θs == [FT(0)]
 
-    p = (; drivers = ClimaLand.initialize_drivers(nothing, pr, coords))
-    rad_only_update! = ClimaLand.make_update_drivers(nothing, pr)
+    p = (; drivers = ClimaLand.initialize_drivers((pr,), coords))
+    rad_only_update! = ClimaLand.make_update_drivers((pr,))
     rad_only_update!(p, 0.0)
     @test p.drivers.SW_d == [FT(10)]
     @test p.drivers.LW_d == [FT(10)]
@@ -150,8 +150,16 @@ end
 
     liquid_precip = TimeVaryingInput((t) -> -1.0)
     pp = ClimaLand.PrescribedPrecipitation{FT}(liquid_precip)
-    precip_update! = ClimaLand.make_update_drivers(pp, nothing)
-    p = (; drivers = ClimaLand.initialize_drivers(pp, nothing, coords))
+    precip_update! = ClimaLand.make_update_drivers((pp,))
+    p = (; drivers = ClimaLand.initialize_drivers((pp,), coords))
     precip_update!(p, 0.0)
     @test p.drivers.P_liq == [FT(-1.0)]
+
+
+    soc = TimeVaryingInput((t) -> 1.0)
+    soc_driver = ClimaLand.PrescribedSoilOrganicCarbon{FT}(soc)
+    soc_update! = ClimaLand.make_update_drivers((soc_driver,))
+    soc_p = (; drivers = ClimaLand.initialize_drivers((soc_driver,), coords))
+    soc_update!(soc_p, 0.0)
+    @test soc_p.drivers.soc == [FT(1.0), FT(1.0)]
 end
diff --git a/test/standalone/Soil/Biogeochemistry/biogeochemistry_module.jl b/test/standalone/Soil/Biogeochemistry/biogeochemistry_module.jl
index 5bd4f52da8..686eaf65f8 100644
--- a/test/standalone/Soil/Biogeochemistry/biogeochemistry_module.jl
+++ b/test/standalone/Soil/Biogeochemistry/biogeochemistry_module.jl
@@ -17,7 +17,9 @@ for FT in (Float32, Float64)
         T_soil = (z, t) -> eltype(z)(t)
         θ_l = (z, t) -> eltype(z)(0.3)
         θ_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
+        Csom = ClimaLand.PrescribedSoilOrganicCarbon{FT}(
+            TimeVaryingInput((t) -> 5),
+        )
         D_ref = FT(0.0)
         parameters = SoilCO2ModelParameters(FT; ν = 0.556, D_ref)
 
@@ -55,11 +57,7 @@ for FT in (Float32, Float64)
             c_co2 = TimeVaryingInput(atmos_co2),
         )
 
-        soil_drivers = SoilDrivers(
-            PrescribedMet{FT}(T_soil, θ_l),
-            PrescribedSOC{FT}(Csom),
-            atmos,
-        )
+        soil_drivers = SoilDrivers(PrescribedMet{FT}(T_soil, θ_l), Csom, atmos)
 
         model = SoilCO2Model{FT}(;
             parameters = parameters,
@@ -87,7 +85,9 @@ for FT in (Float32, Float64)
         T_soil = (z, t) -> eltype(z)(303)
         θ_l = (z, t) -> eltype(z)(0.3)
         θ_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
+        Csom = ClimaLand.PrescribedSoilOrganicCarbon{FT}(
+            TimeVaryingInput((t) -> 5),
+        )
 
         parameters = SoilCO2ModelParameters(FT; ν = 0.556)
         C = FT(4)
@@ -127,7 +127,7 @@ for FT in (Float32, Float64)
 
         soil_drivers = SoilDrivers(
             PrescribedMet{FT}(T_soil, θ_l),
-            PrescribedSOC{FT}(Csom),
+            Csom,
             atmos, # need to create some functions
         )