diff --git a/experiments/benchmarks/land.jl b/experiments/benchmarks/land.jl
index 8d58dafc01..bf3a54805a 100644
--- a/experiments/benchmarks/land.jl
+++ b/experiments/benchmarks/land.jl
@@ -27,12 +27,9 @@ import ClimaTimeSteppers as CTS
 using Test
 using ClimaCore
 using ClimaUtilities.ClimaArtifacts
-import Interpolations
 
 import ClimaUtilities.TimeVaryingInputs:
     TimeVaryingInput, LinearInterpolation, PeriodicCalendar
-import ClimaUtilities.SpaceVaryingInputs: SpaceVaryingInput
-import ClimaUtilities.Regridders: InterpolationsRegridder
 import ClimaUtilities.ClimaArtifacts: @clima_artifact
 import ClimaParams as CP
 
@@ -49,7 +46,6 @@ import Profile, ProfileCanvas
 
 const FT = Float64;
 time_interpolation_method = LinearInterpolation(PeriodicCalendar())
-regridder_type = :InterpolationsRegridder
 context = ClimaComms.context()
 device = ClimaComms.device()
 device_suffix = device isa ClimaComms.CPUSingleThreaded ? "cpu" : "gpu"
@@ -58,7 +54,6 @@ outdir = "land_benchmark_$(device_suffix)"
 !ispath(outdir) && mkpath(outdir)
 
 function setup_prob(t0, tf, Δt; nelements = (101, 15))
-
     earth_param_set = LP.LandParameters(FT)
     radius = FT(6378.1e3)
     depth = FT(50)
@@ -85,129 +80,26 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
         earth_param_set,
         FT;
         time_interpolation_method = time_interpolation_method,
-        regridder_type = regridder_type,
-    )
-
-    soil_params_artifact_path =
-        ClimaLand.Artifacts.soil_params_artifact_folder_path(; context)
-    extrapolation_bc = (
-        Interpolations.Periodic(),
-        Interpolations.Flat(),
-        Interpolations.Flat(),
     )
-    soil_params_mask = SpaceVaryingInput(
-        joinpath(
-            soil_params_artifact_path,
-            "vGalpha_map_gupta_etal2020_1.0x1.0x4.nc",
-        ),
-        "α",
-        subsurface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-        file_reader_kwargs = (; preprocess_func = (data) -> data > 0,),
-    )
-    oceans_to_value(field, mask, value) =
-        mask == 1.0 ? field : eltype(field)(value)
-
-    vg_α = SpaceVaryingInput(
-        joinpath(
-            soil_params_artifact_path,
-            "vGalpha_map_gupta_etal2020_1.0x1.0x4.nc",
-        ),
-        "α",
-        subsurface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
-    vg_n = SpaceVaryingInput(
-        joinpath(
-            soil_params_artifact_path,
-            "vGn_map_gupta_etal2020_1.0x1.0x4.nc",
-        ),
-        "n",
-        subsurface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
-    x = parent(vg_α)
-    μ = mean(log10.(x[x .> 0]))
-    vg_α .= oceans_to_value.(vg_α, soil_params_mask, 10.0^μ)
-
-    x = parent(vg_n)
-    μ = mean(x[x .> 0])
-    vg_n .= oceans_to_value.(vg_n, soil_params_mask, μ)
-
-    vg_fields_to_hcm_field(α::FT, n::FT) where {FT} =
-        ClimaLand.Soil.vanGenuchten{FT}(; @NamedTuple{α::FT, n::FT}((α, n))...)
-    hydrology_cm = vg_fields_to_hcm_field.(vg_α, vg_n)
-
-    θ_r = SpaceVaryingInput(
-        joinpath(
-            soil_params_artifact_path,
-            "residual_map_gupta_etal2020_1.0x1.0x4.nc",
-        ),
-        "θ_r",
-        subsurface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
-
-    ν = SpaceVaryingInput(
-        joinpath(
-            soil_params_artifact_path,
-            "porosity_map_gupta_etal2020_1.0x1.0x4.nc",
-        ),
-        "ν",
-        subsurface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
-    K_sat = SpaceVaryingInput(
-        joinpath(
-            soil_params_artifact_path,
-            "ksat_map_gupta_etal2020_1.0x1.0x4.nc",
-        ),
-        "Ksat",
-        subsurface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
-
-    x = parent(K_sat)
-    μ = mean(log10.(x[x .> 0]))
-    K_sat .= oceans_to_value.(K_sat, soil_params_mask, 10.0^μ)
-
-    ν .= oceans_to_value.(ν, soil_params_mask, 1)
-
-    θ_r .= oceans_to_value.(θ_r, soil_params_mask, 0)
-
-
-    S_s =
-        oceans_to_value.(
-            ClimaCore.Fields.zeros(subsurface_space) .+ FT(1e-3),
-            soil_params_mask,
-            1,
-        )
-    ν_ss_om =
-        oceans_to_value.(
-            ClimaCore.Fields.zeros(subsurface_space) .+ FT(0.1),
-            soil_params_mask,
-            0,
-        )
-    ν_ss_quartz =
-        oceans_to_value.(
-            ClimaCore.Fields.zeros(subsurface_space) .+ FT(0.1),
-            soil_params_mask,
-            0,
-        )
-    ν_ss_gravel =
-        oceans_to_value.(
-            ClimaCore.Fields.zeros(subsurface_space) .+ FT(0.1),
-            soil_params_mask,
-            0,
+    spatially_varying_soil_params =
+        ClimaLand.default_spatially_varying_soil_parameters(
+            subsurface_space,
+            surface_space,
+            FT,
         )
-    PAR_albedo = ClimaCore.Fields.zeros(surface_space) .+ FT(0.2)
-    NIR_albedo = ClimaCore.Fields.zeros(surface_space) .+ FT(0.2)
+    (;
+        ν,
+        ν_ss_om,
+        ν_ss_quartz,
+        ν_ss_gravel,
+        hydrology_cm,
+        K_sat,
+        S_s,
+        θ_r,
+        PAR_albedo,
+        NIR_albedo,
+        f_max,
+    ) = spatially_varying_soil_params
     soil_params = Soil.EnergyHydrologyParameters(
         FT;
         ν,
@@ -221,23 +113,6 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
         PAR_albedo = PAR_albedo,
         NIR_albedo = NIR_albedo,
     )
-
-    soil_params_mask_sfc =
-        ClimaLand.Domains.top_center_to_surface(soil_params_mask)
-
-    # Read in f_max data and land sea mask
-    infile_path = ClimaLand.Artifacts.topmodel_data_path()
-    f_max =
-        SpaceVaryingInput(infile_path, "fmax", surface_space; regridder_type)
-    mask = SpaceVaryingInput(
-        infile_path,
-        "landsea_mask",
-        surface_space;
-        regridder_type,
-    )
-    # Unsure how to handle two masks
-    f_max = oceans_to_value.(f_max, mask, FT(0.0))
-    f_max = oceans_to_value.(f_max, soil_params_mask_sfc, FT(0.0))
     f_over = FT(3.28) # 1/m
     R_sb = FT(1.484e-4 / 1000) # m/s
     runoff_model = ClimaLand.Soil.Runoff.TOPMODELRunoff{FT}(;
@@ -245,94 +120,23 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
         f_max = f_max,
         R_sb = R_sb,
     )
+    # Spatially varying canopy parameters from CLM
+    clm_parameters = ClimaLand.clm_canopy_parameters(surface_space)
+    (;
+        Ω,
+        rooting_depth,
+        is_c3,
+        Vcmax25,
+        g1,
+        G_Function,
+        α_PAR_leaf,
+        τ_PAR_leaf,
+        α_NIR_leaf,
+        τ_NIR_leaf,
+    ) = clm_parameters
 
     # Energy Balance model
     ac_canopy = FT(2.5e3)
-
-    # clm_data is used for g1, vcmax, rooting, and two_stream param maps
-    clm_artifact_path = ClimaLand.Artifacts.clm_data_folder_path(; context)
-
-    # Foliage clumping index data derived from MODIS
-    modis_ci_artifact_path =
-        ClimaLand.Artifacts.modis_ci_data_folder_path(; context)
-
-    # TwoStreamModel parameters
-    nans_to_one(x) = isnan(x) ? eltype(x)(1) : x
-    Ω = SpaceVaryingInput(
-        joinpath(modis_ci_artifact_path, "He_et_al_2012_1x1.nc"),
-        "ci",
-        surface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-        file_reader_kwargs = (; preprocess_func = nans_to_one,),
-    )
-    χl = SpaceVaryingInput(
-        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
-        "xl",
-        surface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
-    G_Function = CLMGFunction(χl)
-    α_PAR_leaf = SpaceVaryingInput(
-        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
-        "rholvis",
-        surface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
-    τ_PAR_leaf = SpaceVaryingInput(
-        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
-        "taulvis",
-        surface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
-    α_NIR_leaf = SpaceVaryingInput(
-        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
-        "rholnir",
-        surface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
-    τ_NIR_leaf = SpaceVaryingInput(
-        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
-        "taulnir",
-        surface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
-
-    # 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),),
-    )
-
-    #Photosynthesis model
-    # 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"),
-        "vcmx25",
-        surface_space;
-        regridder_type,
-        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, "vegetation_properties_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)
@@ -346,13 +150,6 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
     retention_model = Canopy.PlantHydraulics.LinearRetentionCurve{FT}(a)
     plant_ν = FT(1.44e-4)
     plant_S_s = FT(1e-2 * 0.0098) # m3/m3/MPa to m3/m3/m
-    rooting_depth = SpaceVaryingInput(
-        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
-        "rooting_depth",
-        surface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
     n_stem = 0
     n_leaf = 1
     h_stem = FT(0.0)
@@ -432,7 +229,6 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
         surface_space,
         start_date;
         time_interpolation_method = time_interpolation_method,
-        regridder_type = regridder_type,
     )
     ai_parameterization =
         Canopy.PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI)
diff --git a/experiments/benchmarks/richards.jl b/experiments/benchmarks/richards.jl
index 41c91101e0..fc070ff068 100644
--- a/experiments/benchmarks/richards.jl
+++ b/experiments/benchmarks/richards.jl
@@ -26,12 +26,10 @@ import SciMLBase
 using Dates
 using Test
 import NCDatasets
-import Interpolations
 
 import ClimaComms
 ClimaComms.@import_required_backends
 import ClimaUtilities.TimeVaryingInputs: TimeVaryingInput
-import ClimaUtilities.SpaceVaryingInputs: SpaceVaryingInput
 import ClimaUtilities.Regridders: InterpolationsRegridder
 import ClimaTimeSteppers as CTS
 using ClimaCore
@@ -55,14 +53,10 @@ regridder_type = :InterpolationsRegridder
 context = ClimaComms.context()
 device = ClimaComms.device()
 device_suffix = device isa ClimaComms.CPUSingleThreaded ? "cpu" : "gpu"
-
 outdir = "richards_benchmark_$(device_suffix)"
 !ispath(outdir) && mkpath(outdir)
 
 function setup_prob(t0, tf, Δt; nelements = (101, 15))
-    # We set up the problem in a function so that we can make multiple copies (for profiling)
-
-    # Set up simulation domain
     soil_depth = FT(50)
     domain = ClimaLand.Domains.SphericalShell(;
         radius = FT(6.3781e6),
@@ -75,19 +69,13 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
     subsurface_space = domain.space.subsurface
 
     start_date = DateTime(2021)
-
-    # Read in f_max data and land sea mask
-    infile_path = ClimaLand.Artifacts.topmodel_data_path()
-    f_max =
-        SpaceVaryingInput(infile_path, "fmax", surface_space; regridder_type)
-    mask = SpaceVaryingInput(
-        infile_path,
-        "landsea_mask",
-        surface_space;
-        regridder_type,
-    )
-
-    oceans_to_zero(field, mask) = mask > 0.5 ? field : eltype(field)(0)
+    spatially_varying_soil_params =
+        ClimaLand.default_spatially_varying_soil_parameters(
+            subsurface_space,
+            surface_space,
+            FT,
+        )
+    (; ν, hydrology_cm, K_sat, S_s, θ_r, f_max) = spatially_varying_soil_params
     f_over = FT(3.28) # 1/m
     R_sb = FT(1.484e-4 / 1000) # m/s
     runoff_model = ClimaLand.Soil.Runoff.TOPMODELRunoff{FT}(;
@@ -95,92 +83,6 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
         f_max = f_max,
         R_sb = R_sb,
     )
-    soil_params_artifact_path =
-        ClimaLand.Artifacts.soil_params_artifact_folder_path(; context)
-
-    extrapolation_bc = (
-        Interpolations.Periodic(),
-        Interpolations.Flat(),
-        Interpolations.Flat(),
-    )
-    # We use this mask to set values of these parameters over the ocean, in order
-    # to keep them in the physical range
-    function mask_vg(var, value)
-        if var < 1e-8
-            return value
-        else
-            return var
-        end
-    end
-    vg_α = SpaceVaryingInput(
-        joinpath(
-            soil_params_artifact_path,
-            "vGalpha_map_gupta_etal2020_1.0x1.0x4.nc",
-        ),
-        "α",
-        subsurface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
-    vg_α .= mask_vg.(vg_α, 1e-3)
-    # We use this mask to set values of this parameter over the ocean, in order
-    # to keep it in the physical range
-    function mask_vg_n(var, value)
-        if var < 1
-            return value
-        else
-            return var
-        end
-    end
-    vg_n = SpaceVaryingInput(
-        joinpath(
-            soil_params_artifact_path,
-            "vGn_map_gupta_etal2020_1.0x1.0x4.nc",
-        ),
-        "n",
-        subsurface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
-    vg_n .= mask_vg_n.(vg_n, 1.001)
-    vg_fields_to_hcm_field(α::FT, n::FT) where {FT} =
-        ClimaLand.Soil.vanGenuchten{FT}(; @NamedTuple{α::FT, n::FT}((α, n))...)
-    hydrology_cm = vg_fields_to_hcm_field.(vg_α, vg_n)
-
-    θ_r = SpaceVaryingInput(
-        joinpath(
-            soil_params_artifact_path,
-            "residual_map_gupta_etal2020_1.0x1.0x4.nc",
-        ),
-        "θ_r",
-        subsurface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
-
-    ν = SpaceVaryingInput(
-        joinpath(
-            soil_params_artifact_path,
-            "porosity_map_gupta_etal2020_1.0x1.0x4.nc",
-        ),
-        "ν",
-        subsurface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
-    ν .= mask_vg.(ν, 1.0)
-    K_sat = SpaceVaryingInput(
-        joinpath(
-            soil_params_artifact_path,
-            "ksat_map_gupta_etal2020_1.0x1.0x4.nc",
-        ),
-        "Ksat",
-        subsurface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
-    S_s = ClimaCore.Fields.zeros(subsurface_space) .+ FT(1e-3)
-
     soil_params = ClimaLand.Soil.RichardsParameters(;
         hydrology_cm = hydrology_cm,
         ν = ν,
@@ -219,8 +121,8 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
         lateral_flow = false,
     )
 
-    Y, p, t = initialize(model)
-    z = ClimaCore.Fields.coordinate_field(domain.space.subsurface).z
+    Y, p, cds = initialize(model)
+    z = ClimaCore.Fields.coordinate_field(cds.subsurface).z
     lat = ClimaCore.Fields.coordinate_field(domain.space.subsurface).lat
     function hydrostatic_profile(
         lat::FT,
@@ -247,9 +149,9 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
     end
 
     # Set initial state values
+    vg_α = hydrology_cm.α
+    vg_n = hydrology_cm.n
     Y.soil.ϑ_l .= hydrostatic_profile.(lat, z, ν, θ_r, vg_α, vg_n, S_s, f_max)
-    @. Y.soil.ϑ_l = oceans_to_zero(Y.soil.ϑ_l, mask)
-
     # Create model update functions
     set_initial_cache! = make_set_initial_cache(model)
     exp_tendency! = make_exp_tendency(model)
diff --git a/experiments/integrated/global/global_parameters.jl b/experiments/integrated/global/global_parameters.jl
index e9941022bb..479a56ba0a 100644
--- a/experiments/integrated/global/global_parameters.jl
+++ b/experiments/integrated/global/global_parameters.jl
@@ -1,113 +1,22 @@
-soil_params_artifact_path =
-    ClimaLand.Artifacts.soil_params_artifact_folder_path(; context)
-extrapolation_bc =
-    (Interpolations.Periodic(), Interpolations.Flat(), Interpolations.Flat())
-soil_params_mask = SpaceVaryingInput(
-    joinpath(
-        soil_params_artifact_path,
-        "vGalpha_map_gupta_etal2020_1.0x1.0x4.nc",
-    ),
-    "α",
-    subsurface_space;
-    regridder_type,
-    regridder_kwargs = (; extrapolation_bc,),
-    file_reader_kwargs = (; preprocess_func = (data) -> data > 0,),
-)
-oceans_to_value(field, mask, value) = mask == 1.0 ? field : eltype(field)(value)
-
-vg_α = SpaceVaryingInput(
-    joinpath(
-        soil_params_artifact_path,
-        "vGalpha_map_gupta_etal2020_1.0x1.0x4.nc",
-    ),
-    "α",
-    subsurface_space;
-    regridder_type,
-    regridder_kwargs = (; extrapolation_bc,),
-)
-vg_n = SpaceVaryingInput(
-    joinpath(soil_params_artifact_path, "vGn_map_gupta_etal2020_1.0x1.0x4.nc"),
-    "n",
-    subsurface_space;
-    regridder_type,
-    regridder_kwargs = (; extrapolation_bc,),
-)
-x = parent(vg_α)
-μ = mean(log10.(x[x .> 0]))
-vg_α .= oceans_to_value.(vg_α, soil_params_mask, 10.0^μ)
-
-x = parent(vg_n)
-μ = mean(x[x .> 0])
-vg_n .= oceans_to_value.(vg_n, soil_params_mask, μ)
-
-vg_fields_to_hcm_field(α::FT, n::FT) where {FT} =
-    ClimaLand.Soil.vanGenuchten{FT}(; @NamedTuple{α::FT, n::FT}((α, n))...)
-hydrology_cm = vg_fields_to_hcm_field.(vg_α, vg_n)
-
-θ_r = SpaceVaryingInput(
-    joinpath(
-        soil_params_artifact_path,
-        "residual_map_gupta_etal2020_1.0x1.0x4.nc",
-    ),
-    "θ_r",
-    subsurface_space;
-    regridder_type,
-    regridder_kwargs = (; extrapolation_bc,),
-)
-
-ν = SpaceVaryingInput(
-    joinpath(
-        soil_params_artifact_path,
-        "porosity_map_gupta_etal2020_1.0x1.0x4.nc",
-    ),
-    "ν",
-    subsurface_space;
-    regridder_type,
-    regridder_kwargs = (; extrapolation_bc,),
-)
-K_sat = SpaceVaryingInput(
-    joinpath(soil_params_artifact_path, "ksat_map_gupta_etal2020_1.0x1.0x4.nc"),
-    "Ksat",
-    subsurface_space;
-    regridder_type,
-    regridder_kwargs = (; extrapolation_bc,),
-)
-
-x = parent(K_sat)
-μ = mean(log10.(x[x .> 0]))
-K_sat .= oceans_to_value.(K_sat, soil_params_mask, 10.0^μ)
-
-ν .= oceans_to_value.(ν, soil_params_mask, 1)
-
-θ_r .= oceans_to_value.(θ_r, soil_params_mask, 0)
-
-
-S_s =
-    oceans_to_value.(
-        ClimaCore.Fields.zeros(subsurface_space) .+ FT(1e-3),
-        soil_params_mask,
-        1,
-    )
-ν_ss_om =
-    oceans_to_value.(
-        ClimaCore.Fields.zeros(subsurface_space) .+ FT(0.1),
-        soil_params_mask,
-        0,
+spatially_varying_soil_params =
+    ClimaLand.default_spatially_varying_soil_parameters(
+        subsurface_space,
+        surface_space,
+        FT,
     )
-ν_ss_quartz =
-    oceans_to_value.(
-        ClimaCore.Fields.zeros(subsurface_space) .+ FT(0.1),
-        soil_params_mask,
-        0,
-    )
-ν_ss_gravel =
-    oceans_to_value.(
-        ClimaCore.Fields.zeros(subsurface_space) .+ FT(0.1),
-        soil_params_mask,
-        0,
-    )
-PAR_albedo = ClimaCore.Fields.zeros(surface_space) .+ FT(0.2)
-NIR_albedo = ClimaCore.Fields.zeros(surface_space) .+ FT(0.2)
+(;
+    ν,
+    ν_ss_om,
+    ν_ss_quartz,
+    ν_ss_gravel,
+    hydrology_cm,
+    K_sat,
+    S_s,
+    θ_r,
+    PAR_albedo,
+    NIR_albedo,
+    f_max,
+) = spatially_varying_soil_params
 soil_params = Soil.EnergyHydrologyParameters(
     FT;
     ν,
@@ -121,20 +30,6 @@ soil_params = Soil.EnergyHydrologyParameters(
     PAR_albedo = PAR_albedo,
     NIR_albedo = NIR_albedo,
 );
-soil_params_mask_sfc = ClimaLand.Domains.top_center_to_surface(soil_params_mask)
-
-# Read in f_max data and land sea mask
-infile_path = ClimaLand.Artifacts.topmodel_data_path()
-f_max = SpaceVaryingInput(infile_path, "fmax", surface_space; regridder_type)
-mask = SpaceVaryingInput(
-    infile_path,
-    "landsea_mask",
-    surface_space;
-    regridder_type,
-)
-# Unsure how to handle two masks
-f_max = oceans_to_value.(f_max, mask, FT(0.0))
-f_max = oceans_to_value.(f_max, soil_params_mask_sfc, FT(0.0))
 f_over = FT(3.28) # 1/m
 R_sb = FT(1.484e-4 / 1000) # m/s
 runoff_model = ClimaLand.Soil.Runoff.TOPMODELRunoff{FT}(;
@@ -143,23 +38,23 @@ runoff_model = ClimaLand.Soil.Runoff.TOPMODELRunoff{FT}(;
     R_sb = R_sb,
 )
 
-# TwoStreamModel parameters
-Ω = FT(0.69)
-ld = FT(0.5)
-α_PAR_leaf = FT(0.1)
-τ_PAR_leaf = FT(0.05)
-α_NIR_leaf = FT(0.45)
-τ_NIR_leaf = FT(0.25)
 
+# Spatially varying canopy parameters from CLM
+clm_parameters = ClimaLand.clm_canopy_parameters(surface_space)
+(;
+    Ω,
+    rooting_depth,
+    is_c3,
+    Vcmax25,
+    g1,
+    G_Function,
+    α_PAR_leaf,
+    τ_PAR_leaf,
+    α_NIR_leaf,
+    τ_NIR_leaf,
+) = clm_parameters
 # Energy Balance model
-ac_canopy = FT(2.5e4)
-
-# Conductance Model
-g1 = FT(141) # Wang et al: 141 sqrt(Pa) for Medlyn model; Natan used 300.
-
-#Photosynthesis model
-Vcmax25 = FT(9e-5) # from Yujie's paper 4.5e-5 , Natan used 9e-5
-
+ac_canopy = FT(2.5e3)
 # Plant Hydraulics and general plant parameters
 SAI = FT(0.0) # m2/m2
 f_root_to_shoot = FT(3.5)
@@ -173,7 +68,6 @@ conductivity_model =
 retention_model = Canopy.PlantHydraulics.LinearRetentionCurve{FT}(a)
 plant_ν = FT(1.44e-4)
 plant_S_s = FT(1e-2 * 0.0098) # m3/m3/MPa to m3/m3/m
-rooting_depth = FT(0.5) # from Natan
 n_stem = 0
 n_leaf = 1
 h_stem = FT(0.0)
@@ -187,5 +81,4 @@ compartment_surfaces = n_stem > 0 ? [zmax, h_stem, h_canopy] : [zmax, h_leaf]
 z0_m = FT(0.13) * h_canopy
 z0_b = FT(0.1) * z0_m
 
-
 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 83ec410fcf..82405e7291 100644
--- a/experiments/integrated/global/global_soil_canopy.jl
+++ b/experiments/integrated/global/global_soil_canopy.jl
@@ -4,11 +4,8 @@ ClimaComms.@import_required_backends
 import ClimaTimeSteppers as CTS
 using ClimaCore
 using ClimaUtilities.ClimaArtifacts
-import Interpolations
 import ClimaUtilities.TimeVaryingInputs:
     TimeVaryingInput, LinearInterpolation, PeriodicCalendar
-import ClimaUtilities.SpaceVaryingInputs: SpaceVaryingInput
-import ClimaUtilities.Regridders: InterpolationsRegridder
 import ClimaUtilities.OutputPathGenerator: generate_output_path
 import ClimaUtilities.ClimaArtifacts: @clima_artifact
 import ClimaParams as CP
@@ -30,9 +27,6 @@ import ClimaAnalysis
 import ClimaAnalysis.Visualize as viz
 import ClimaUtilities
 time_interpolation_method = LinearInterpolation(PeriodicCalendar())
-regridder_type = :InterpolationsRegridder
-extrapolation_bc =
-    (Interpolations.Periodic(), Interpolations.Flat(), Interpolations.Flat())
 context = ClimaComms.context()
 outdir = generate_output_path("experiments/integrated/global")
 
@@ -68,7 +62,6 @@ atmos, radiation = ClimaLand.prescribed_forcing_era5(
     earth_param_set,
     FT;
     time_interpolation_method = time_interpolation_method,
-    regridder_type = regridder_type,
 )
 
 include(
@@ -128,7 +121,6 @@ radiative_transfer_args = (;
 # Set up conductance
 conductance_args = (; parameters = Canopy.MedlynConductanceParameters(FT; g1))
 # Set up photosynthesis
-is_c3 = FT(1) # set the photosynthesis mechanism to C3
 photosynthesis_args =
     (; parameters = Canopy.FarquharParameters(FT, is_c3; Vcmax25 = Vcmax25))
 # Set up plant hydraulics
@@ -137,7 +129,6 @@ LAIfunction = ClimaLand.prescribed_lai_era5(
     surface_space,
     start_date;
     time_interpolation_method = time_interpolation_method,
-    regridder_type = regridder_type,
 )
 ai_parameterization = Canopy.PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI)
 
@@ -202,7 +193,7 @@ land = SoilCanopyModel{FT}(;
 Y, p, cds = initialize(land)
 
 t0 = 0.0
-dt = 900.0
+dt = 450.0
 tf = 3600
 
 init_soil(ν, θ_r) = θ_r + (ν - θ_r) / 2
diff --git a/experiments/long_runs/land.jl b/experiments/long_runs/land.jl
index d2afce92e3..c92dac07bd 100644
--- a/experiments/long_runs/land.jl
+++ b/experiments/long_runs/land.jl
@@ -20,7 +20,6 @@ ClimaComms.@import_required_backends
 import ClimaTimeSteppers as CTS
 using ClimaCore
 using ClimaUtilities.ClimaArtifacts
-import Interpolations
 
 import ClimaDiagnostics
 import ClimaAnalysis
@@ -29,8 +28,6 @@ import ClimaUtilities
 
 import ClimaUtilities.TimeVaryingInputs:
     TimeVaryingInput, LinearInterpolation, PeriodicCalendar
-import ClimaUtilities.SpaceVaryingInputs: SpaceVaryingInput
-import ClimaUtilities.Regridders: InterpolationsRegridder
 import ClimaUtilities.ClimaArtifacts: @clima_artifact
 import ClimaParams as CP
 
@@ -48,7 +45,6 @@ import NCDatasets
 
 const FT = Float64;
 time_interpolation_method = LinearInterpolation(PeriodicCalendar())
-regridder_type = :InterpolationsRegridder
 context = ClimaComms.context()
 device = ClimaComms.device()
 device_suffix = device isa ClimaComms.CPUSingleThreaded ? "cpu" : "gpu"
@@ -58,7 +54,6 @@ outdir =
     ClimaUtilities.OutputPathGenerator.generate_output_path(diagnostics_outdir)
 
 function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
-
     earth_param_set = LP.LandParameters(FT)
     radius = FT(6378.1e3)
     depth = FT(50)
@@ -82,131 +77,28 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         surface_space,
         start_date,
         earth_param_set,
-        FT;
-        time_interpolation_method = time_interpolation_method,
-        regridder_type = regridder_type,
-    )
-
-    soil_params_artifact_path =
-        ClimaLand.Artifacts.soil_params_artifact_folder_path(; context)
-    extrapolation_bc = (
-        Interpolations.Periodic(),
-        Interpolations.Flat(),
-        Interpolations.Flat(),
-    )
-    soil_params_mask = SpaceVaryingInput(
-        joinpath(
-            soil_params_artifact_path,
-            "vGalpha_map_gupta_etal2020_1.0x1.0x4.nc",
-        ),
-        "α",
-        subsurface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-        file_reader_kwargs = (; preprocess_func = (data) -> data > 0,),
+        FT,
     )
-    oceans_to_value(field, mask, value) =
-        mask == 1.0 ? field : eltype(field)(value)
 
-    vg_α = SpaceVaryingInput(
-        joinpath(
-            soil_params_artifact_path,
-            "vGalpha_map_gupta_etal2020_1.0x1.0x4.nc",
-        ),
-        "α",
-        subsurface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
-    vg_n = SpaceVaryingInput(
-        joinpath(
-            soil_params_artifact_path,
-            "vGn_map_gupta_etal2020_1.0x1.0x4.nc",
-        ),
-        "n",
-        subsurface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
-    x = parent(vg_α)
-    μ = mean(log10.(x[x .> 0]))
-    vg_α .= oceans_to_value.(vg_α, soil_params_mask, 10.0^μ)
-
-    x = parent(vg_n)
-    μ = mean(x[x .> 0])
-    vg_n .= oceans_to_value.(vg_n, soil_params_mask, μ)
-
-    vg_fields_to_hcm_field(α::FT, n::FT) where {FT} =
-        ClimaLand.Soil.vanGenuchten{FT}(; @NamedTuple{α::FT, n::FT}((α, n))...)
-    hydrology_cm = vg_fields_to_hcm_field.(vg_α, vg_n)
-
-    θ_r = SpaceVaryingInput(
-        joinpath(
-            soil_params_artifact_path,
-            "residual_map_gupta_etal2020_1.0x1.0x4.nc",
-        ),
-        "θ_r",
-        subsurface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
-
-    ν = SpaceVaryingInput(
-        joinpath(
-            soil_params_artifact_path,
-            "porosity_map_gupta_etal2020_1.0x1.0x4.nc",
-        ),
-        "ν",
-        subsurface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
-    K_sat = SpaceVaryingInput(
-        joinpath(
-            soil_params_artifact_path,
-            "ksat_map_gupta_etal2020_1.0x1.0x4.nc",
-        ),
-        "Ksat",
-        subsurface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
-
-    x = parent(K_sat)
-    μ = mean(log10.(x[x .> 0]))
-    K_sat .= oceans_to_value.(K_sat, soil_params_mask, 10.0^μ)
-
-    ν .= oceans_to_value.(ν, soil_params_mask, 1)
-
-    θ_r .= oceans_to_value.(θ_r, soil_params_mask, 0)
-
-
-    S_s =
-        oceans_to_value.(
-            ClimaCore.Fields.zeros(subsurface_space) .+ FT(1e-3),
-            soil_params_mask,
-            1,
-        )
-    ν_ss_om =
-        oceans_to_value.(
-            ClimaCore.Fields.zeros(subsurface_space) .+ FT(0.1),
-            soil_params_mask,
-            0,
-        )
-    ν_ss_quartz =
-        oceans_to_value.(
-            ClimaCore.Fields.zeros(subsurface_space) .+ FT(0.1),
-            soil_params_mask,
-            0,
+    spatially_varying_soil_params =
+        ClimaLand.default_spatially_varying_soil_parameters(
+            subsurface_space,
+            surface_space,
+            FT,
         )
-    ν_ss_gravel =
-        oceans_to_value.(
-            ClimaCore.Fields.zeros(subsurface_space) .+ FT(0.1),
-            soil_params_mask,
-            0,
-        )
-    PAR_albedo = ClimaCore.Fields.zeros(surface_space) .+ FT(0.2)
-    NIR_albedo = ClimaCore.Fields.zeros(surface_space) .+ FT(0.2)
+    (;
+        ν,
+        ν_ss_om,
+        ν_ss_quartz,
+        ν_ss_gravel,
+        hydrology_cm,
+        K_sat,
+        S_s,
+        θ_r,
+        PAR_albedo,
+        NIR_albedo,
+        f_max,
+    ) = spatially_varying_soil_params
     soil_params = Soil.EnergyHydrologyParameters(
         FT;
         ν,
@@ -220,23 +112,6 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         PAR_albedo = PAR_albedo,
         NIR_albedo = NIR_albedo,
     )
-
-    soil_params_mask_sfc =
-        ClimaLand.Domains.top_center_to_surface(soil_params_mask)
-
-    # Read in f_max data and land sea mask
-    infile_path = ClimaLand.Artifacts.topmodel_data_path()
-    f_max =
-        SpaceVaryingInput(infile_path, "fmax", surface_space; regridder_type)
-    mask = SpaceVaryingInput(
-        infile_path,
-        "landsea_mask",
-        surface_space;
-        regridder_type,
-    )
-    # Unsure how to handle two masks
-    f_max = oceans_to_value.(f_max, mask, FT(0.0))
-    f_max = oceans_to_value.(f_max, soil_params_mask_sfc, FT(0.0))
     f_over = FT(3.28) # 1/m
     R_sb = FT(1.484e-4 / 1000) # m/s
     runoff_model = ClimaLand.Soil.Runoff.TOPMODELRunoff{FT}(;
@@ -245,90 +120,23 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         R_sb = R_sb,
     )
 
+    # Spatially varying canopy parameters from CLM
+    clm_parameters = ClimaLand.clm_canopy_parameters(surface_space)
+    (;
+        Ω,
+        rooting_depth,
+        is_c3,
+        Vcmax25,
+        g1,
+        G_Function,
+        α_PAR_leaf,
+        τ_PAR_leaf,
+        α_NIR_leaf,
+        τ_NIR_leaf,
+    ) = clm_parameters
+
     # Energy Balance model
     ac_canopy = FT(2.5e3)
-
-    #clm_data is used for g1, vcmax, rooting, and two_stream param maps
-    clm_artifact_path = ClimaLand.Artifacts.clm_data_folder_path(; context)
-    # Foliage clumping index data derived from MODIS
-    modis_ci_artifact_path =
-        ClimaLand.Artifacts.modis_ci_data_folder_path(; context)
-
-    # TwoStreamModel parameters
-    nans_to_one(x) = isnan(x) ? eltype(x)(1) : x
-    Ω = SpaceVaryingInput(
-        joinpath(modis_ci_artifact_path, "He_et_al_2012_1x1.nc"),
-        "ci",
-        surface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-        file_reader_kwargs = (; preprocess_func = nans_to_one,),
-    )
-    χl = SpaceVaryingInput(
-        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
-        "xl",
-        surface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
-    G_Function = CLMGFunction(χl)
-    α_PAR_leaf = SpaceVaryingInput(
-        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
-        "rholvis",
-        surface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
-    τ_PAR_leaf = SpaceVaryingInput(
-        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
-        "taulvis",
-        surface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
-    α_NIR_leaf = SpaceVaryingInput(
-        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
-        "rholnir",
-        surface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
-    τ_NIR_leaf = SpaceVaryingInput(
-        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
-        "taulnir",
-        surface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
-    # 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),),
-    )
-
-    #Photosynthesis model
-    # 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"),
-        "vcmx25",
-        surface_space;
-        regridder_type,
-        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, "vegetation_properties_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)
@@ -342,13 +150,6 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
     retention_model = Canopy.PlantHydraulics.LinearRetentionCurve{FT}(a)
     plant_ν = FT(1.44e-4)
     plant_S_s = FT(1e-2 * 0.0098) # m3/m3/MPa to m3/m3/m
-    rooting_depth = SpaceVaryingInput(
-        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
-        "rooting_depth",
-        surface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
     n_stem = 0
     n_leaf = 1
     h_stem = FT(0.0)
@@ -429,7 +230,6 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         surface_space,
         start_date;
         time_interpolation_method = time_interpolation_method,
-        regridder_type = regridder_type,
     )
     ai_parameterization =
         Canopy.PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI)
diff --git a/experiments/long_runs/land_region.jl b/experiments/long_runs/land_region.jl
index 8307b8a0da..2b4ccba819 100644
--- a/experiments/long_runs/land_region.jl
+++ b/experiments/long_runs/land_region.jl
@@ -20,7 +20,6 @@ ClimaComms.@import_required_backends
 import ClimaTimeSteppers as CTS
 using ClimaCore
 using ClimaUtilities.ClimaArtifacts
-import Interpolations
 
 import ClimaDiagnostics
 import ClimaAnalysis
@@ -29,8 +28,6 @@ import ClimaUtilities
 
 import ClimaUtilities.TimeVaryingInputs:
     TimeVaryingInput, LinearInterpolation, PeriodicCalendar
-import ClimaUtilities.SpaceVaryingInputs: SpaceVaryingInput
-import ClimaUtilities.Regridders: InterpolationsRegridder
 import ClimaUtilities.ClimaArtifacts: @clima_artifact
 import ClimaParams as CP
 
@@ -47,7 +44,6 @@ import NCDatasets
 
 const FT = Float64;
 time_interpolation_method = LinearInterpolation(PeriodicCalendar())
-regridder_type = :InterpolationsRegridder
 context = ClimaComms.context()
 device = ClimaComms.device()
 device_suffix = device isa ClimaComms.CPUSingleThreaded ? "cpu" : "gpu"
@@ -57,7 +53,6 @@ outdir =
     ClimaUtilities.OutputPathGenerator.generate_output_path(diagnostics_outdir)
 
 function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (10, 10, 15))
-
     earth_param_set = LP.LandParameters(FT)
     radius = FT(6378.1e3)
     depth = FT(50)
@@ -85,130 +80,27 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (10, 10, 15))
         surface_space,
         start_date,
         earth_param_set,
-        FT;
-        time_interpolation_method = time_interpolation_method,
-        regridder_type = regridder_type,
-    )
-    soil_params_artifact_path =
-        ClimaLand.Artifacts.soil_params_artifact_folder_path(; context)
-    extrapolation_bc = (
-        Interpolations.Periodic(),
-        Interpolations.Flat(),
-        Interpolations.Flat(),
-    )
-    soil_params_mask = SpaceVaryingInput(
-        joinpath(
-            soil_params_artifact_path,
-            "vGalpha_map_gupta_etal2020_1.0x1.0x4.nc",
-        ),
-        "α",
-        subsurface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-        file_reader_kwargs = (; preprocess_func = (data) -> data > 0,),
-    )
-    oceans_to_value(field, mask, value) =
-        mask == 1.0 ? field : eltype(field)(value)
-
-    vg_α = SpaceVaryingInput(
-        joinpath(
-            soil_params_artifact_path,
-            "vGalpha_map_gupta_etal2020_1.0x1.0x4.nc",
-        ),
-        "α",
-        subsurface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
-    vg_n = SpaceVaryingInput(
-        joinpath(
-            soil_params_artifact_path,
-            "vGn_map_gupta_etal2020_1.0x1.0x4.nc",
-        ),
-        "n",
-        subsurface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
-    x = parent(vg_α)
-    μ = mean(log10.(x[x .> 0]))
-    vg_α .= oceans_to_value.(vg_α, soil_params_mask, 10.0^μ)
-
-    x = parent(vg_n)
-    μ = mean(x[x .> 0])
-    vg_n .= oceans_to_value.(vg_n, soil_params_mask, μ)
-
-    vg_fields_to_hcm_field(α::FT, n::FT) where {FT} =
-        ClimaLand.Soil.vanGenuchten{FT}(; @NamedTuple{α::FT, n::FT}((α, n))...)
-    hydrology_cm = vg_fields_to_hcm_field.(vg_α, vg_n)
-
-    θ_r = SpaceVaryingInput(
-        joinpath(
-            soil_params_artifact_path,
-            "residual_map_gupta_etal2020_1.0x1.0x4.nc",
-        ),
-        "θ_r",
-        subsurface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
+        FT,
     )
-
-    ν = SpaceVaryingInput(
-        joinpath(
-            soil_params_artifact_path,
-            "porosity_map_gupta_etal2020_1.0x1.0x4.nc",
-        ),
-        "ν",
-        subsurface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
-    K_sat = SpaceVaryingInput(
-        joinpath(
-            soil_params_artifact_path,
-            "ksat_map_gupta_etal2020_1.0x1.0x4.nc",
-        ),
-        "Ksat",
-        subsurface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
-
-    x = parent(K_sat)
-    μ = mean(log10.(x[x .> 0]))
-    K_sat .= oceans_to_value.(K_sat, soil_params_mask, 10.0^μ)
-
-    ν .= oceans_to_value.(ν, soil_params_mask, 1)
-
-    θ_r .= oceans_to_value.(θ_r, soil_params_mask, 0)
-
-
-    S_s =
-        oceans_to_value.(
-            ClimaCore.Fields.zeros(subsurface_space) .+ FT(1e-3),
-            soil_params_mask,
-            1,
-        )
-    ν_ss_om =
-        oceans_to_value.(
-            ClimaCore.Fields.zeros(subsurface_space) .+ FT(0.1),
-            soil_params_mask,
-            0,
-        )
-    ν_ss_quartz =
-        oceans_to_value.(
-            ClimaCore.Fields.zeros(subsurface_space) .+ FT(0.1),
-            soil_params_mask,
-            0,
-        )
-    ν_ss_gravel =
-        oceans_to_value.(
-            ClimaCore.Fields.zeros(subsurface_space) .+ FT(0.1),
-            soil_params_mask,
-            0,
+    spatially_varying_soil_params =
+        ClimaLand.default_spatially_varying_soil_parameters(
+            subsurface_space,
+            surface_space,
+            FT,
         )
-    PAR_albedo = ClimaCore.Fields.zeros(surface_space) .+ FT(0.2)
-    NIR_albedo = ClimaCore.Fields.zeros(surface_space) .+ FT(0.2)
+    (;
+        ν,
+        ν_ss_om,
+        ν_ss_quartz,
+        ν_ss_gravel,
+        hydrology_cm,
+        K_sat,
+        S_s,
+        θ_r,
+        PAR_albedo,
+        NIR_albedo,
+        f_max,
+    ) = spatially_varying_soil_params
     soil_params = Soil.EnergyHydrologyParameters(
         FT;
         ν,
@@ -222,23 +114,6 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (10, 10, 15))
         PAR_albedo = PAR_albedo,
         NIR_albedo = NIR_albedo,
     )
-
-    soil_params_mask_sfc =
-        ClimaLand.Domains.top_center_to_surface(soil_params_mask)
-
-    # Read in f_max data and land sea mask
-    infile_path = ClimaLand.Artifacts.topmodel_data_path()
-    f_max =
-        SpaceVaryingInput(infile_path, "fmax", surface_space; regridder_type)
-    mask = SpaceVaryingInput(
-        infile_path,
-        "landsea_mask",
-        surface_space;
-        regridder_type,
-    )
-    # Unsure how to handle two masks
-    f_max = oceans_to_value.(f_max, mask, FT(0.0))
-    f_max = oceans_to_value.(f_max, soil_params_mask_sfc, FT(0.0))
     f_over = FT(3.28) # 1/m
     R_sb = FT(1.484e-4 / 1000) # m/s
     runoff_model = ClimaLand.Soil.Runoff.TOPMODELRunoff{FT}(;
@@ -247,93 +122,22 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (10, 10, 15))
         R_sb = R_sb,
     )
 
+    # Spatially varying canopy parameters from CLM
+    clm_parameters = ClimaLand.clm_canopy_parameters(surface_space)
+    (;
+        Ω,
+        rooting_depth,
+        is_c3,
+        Vcmax25,
+        g1,
+        G_Function,
+        α_PAR_leaf,
+        τ_PAR_leaf,
+        α_NIR_leaf,
+        τ_NIR_leaf,
+    ) = clm_parameters
     # Energy Balance model
-    ac_canopy = FT(2.5e4) # this will likely be 10x smaller!
-    #clm_data is used for g1, vcmax, rooting, and two_stream param maps
-    clm_artifact_path = ClimaLand.Artifacts.clm_data_folder_path(; context)
-
-    # Foliage clumping index data derived from MODIS
-    modis_ci_artifact_path =
-        ClimaLand.Artifacts.modis_ci_data_folder_path(; context)
-
-    # TwoStreamModel parameters
-    nans_to_one(x) = isnan(x) ? eltype(x)(1) : x
-    Ω = SpaceVaryingInput(
-        joinpath(modis_ci_artifact_path, "He_et_al_2012_1x1.nc"),
-        "ci",
-        surface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-        file_reader_kwargs = (; preprocess_func = nans_to_one,),
-    )
-
-    χl = SpaceVaryingInput(
-        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
-        "xl",
-        surface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
-    G_Function = CLMGFunction(χl)
-    α_PAR_leaf = SpaceVaryingInput(
-        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
-        "rholvis",
-        surface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
-    τ_PAR_leaf = SpaceVaryingInput(
-        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
-        "taulvis",
-        surface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
-    α_NIR_leaf = SpaceVaryingInput(
-        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
-        "rholnir",
-        surface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
-    τ_NIR_leaf = SpaceVaryingInput(
-        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
-        "taulnir",
-        surface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
-
-    # 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"),
-        "vcmx25",
-        surface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-        file_reader_kwargs = (; preprocess_func = (data) -> data / 1_000_000,),
-    )
-
-    # c3 is read in as a float
-    is_c3 = SpaceVaryingInput(
-        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
-        "c3_dominant",
-        surface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
-
+    ac_canopy = FT(2.5e3)
     # Plant Hydraulics and general plant parameters
     SAI = FT(0.0) # m2/m2
     f_root_to_shoot = FT(3.5)
@@ -347,13 +151,6 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (10, 10, 15))
     retention_model = Canopy.PlantHydraulics.LinearRetentionCurve{FT}(a)
     plant_ν = FT(1.44e-4)
     plant_S_s = FT(1e-2 * 0.0098) # m3/m3/MPa to m3/m3/m
-    rooting_depth = SpaceVaryingInput(
-        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
-        "rooting_depth",
-        surface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
     n_stem = 0
     n_leaf = 1
     h_stem = FT(0.0)
@@ -433,7 +230,6 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (10, 10, 15))
         surface_space,
         start_date;
         time_interpolation_method = time_interpolation_method,
-        regridder_type = regridder_type,
     )
     ai_parameterization =
         Canopy.PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI)
diff --git a/experiments/long_runs/soil.jl b/experiments/long_runs/soil.jl
index ce667e0d81..14936f5269 100644
--- a/experiments/long_runs/soil.jl
+++ b/experiments/long_runs/soil.jl
@@ -20,7 +20,6 @@ ClimaComms.@import_required_backends
 import ClimaTimeSteppers as CTS
 using ClimaCore
 using ClimaUtilities.ClimaArtifacts
-import Interpolations
 
 using ClimaDiagnostics
 using ClimaAnalysis
@@ -28,8 +27,6 @@ import ClimaAnalysis.Visualize as viz
 using ClimaUtilities
 
 import ClimaUtilities.TimeVaryingInputs: LinearInterpolation, PeriodicCalendar
-import ClimaUtilities.SpaceVaryingInputs: SpaceVaryingInput
-import ClimaUtilities.Regridders: InterpolationsRegridder
 import ClimaUtilities.ClimaArtifacts: @clima_artifact
 import ClimaParams as CP
 
@@ -46,7 +43,6 @@ import NCDatasets
 
 const FT = Float64;
 time_interpolation_method = LinearInterpolation(PeriodicCalendar())
-regridder_type = :InterpolationsRegridder
 context = ClimaComms.context()
 device = ClimaComms.device()
 device_suffix = device isa ClimaComms.CPUSingleThreaded ? "cpu" : "gpu"
@@ -56,7 +52,6 @@ outdir =
     ClimaUtilities.OutputPathGenerator.generate_output_path(diagnostics_outdir)
 
 function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
-
     earth_param_set = LP.LandParameters(FT)
     radius = FT(6378.1e3)
     depth = FT(50)
@@ -81,129 +76,26 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         earth_param_set,
         FT;
         time_interpolation_method = time_interpolation_method,
-        regridder_type = regridder_type,
-    )
-
-    soil_params_artifact_path =
-        ClimaLand.Artifacts.soil_params_artifact_folder_path(; context)
-    extrapolation_bc = (
-        Interpolations.Periodic(),
-        Interpolations.Flat(),
-        Interpolations.Flat(),
-    )
-    soil_params_mask = SpaceVaryingInput(
-        joinpath(
-            soil_params_artifact_path,
-            "vGalpha_map_gupta_etal2020_1.0x1.0x4.nc",
-        ),
-        "α",
-        subsurface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-        file_reader_kwargs = (; preprocess_func = (data) -> data > 0,),
-    )
-    oceans_to_value(field, mask, value) =
-        mask == 1.0 ? field : eltype(field)(value)
-
-    vg_α = SpaceVaryingInput(
-        joinpath(
-            soil_params_artifact_path,
-            "vGalpha_map_gupta_etal2020_1.0x1.0x4.nc",
-        ),
-        "α",
-        subsurface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
-    vg_n = SpaceVaryingInput(
-        joinpath(
-            soil_params_artifact_path,
-            "vGn_map_gupta_etal2020_1.0x1.0x4.nc",
-        ),
-        "n",
-        subsurface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
-    x = parent(vg_α)
-    μ = mean(log10.(x[x .> 0]))
-    vg_α .= oceans_to_value.(vg_α, soil_params_mask, 10.0^μ)
-
-    x = parent(vg_n)
-    μ = mean(x[x .> 0])
-    vg_n .= oceans_to_value.(vg_n, soil_params_mask, μ)
-
-    vg_fields_to_hcm_field(α::FT, n::FT) where {FT} =
-        ClimaLand.Soil.vanGenuchten{FT}(; @NamedTuple{α::FT, n::FT}((α, n))...)
-    hydrology_cm = vg_fields_to_hcm_field.(vg_α, vg_n)
-
-    θ_r = SpaceVaryingInput(
-        joinpath(
-            soil_params_artifact_path,
-            "residual_map_gupta_etal2020_1.0x1.0x4.nc",
-        ),
-        "θ_r",
-        subsurface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
-
-    ν = SpaceVaryingInput(
-        joinpath(
-            soil_params_artifact_path,
-            "porosity_map_gupta_etal2020_1.0x1.0x4.nc",
-        ),
-        "ν",
-        subsurface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
-    )
-    K_sat = SpaceVaryingInput(
-        joinpath(
-            soil_params_artifact_path,
-            "ksat_map_gupta_etal2020_1.0x1.0x4.nc",
-        ),
-        "Ksat",
-        subsurface_space;
-        regridder_type,
-        regridder_kwargs = (; extrapolation_bc,),
     )
-
-    x = parent(K_sat)
-    μ = mean(log10.(x[x .> 0]))
-    K_sat .= oceans_to_value.(K_sat, soil_params_mask, 10.0^μ)
-
-    ν .= oceans_to_value.(ν, soil_params_mask, 1)
-
-    θ_r .= oceans_to_value.(θ_r, soil_params_mask, 0)
-
-
-    S_s =
-        oceans_to_value.(
-            ClimaCore.Fields.zeros(subsurface_space) .+ FT(1e-3),
-            soil_params_mask,
-            1,
+    spatially_varying_soil_params =
+        ClimaLand.default_spatially_varying_soil_parameters(
+            subsurface_space,
+            surface_space,
+            FT,
         )
-    ν_ss_om =
-        oceans_to_value.(
-            ClimaCore.Fields.zeros(subsurface_space) .+ FT(0.1),
-            soil_params_mask,
-            0,
-        )
-    ν_ss_quartz =
-        oceans_to_value.(
-            ClimaCore.Fields.zeros(subsurface_space) .+ FT(0.1),
-            soil_params_mask,
-            0,
-        )
-    ν_ss_gravel =
-        oceans_to_value.(
-            ClimaCore.Fields.zeros(subsurface_space) .+ FT(0.1),
-            soil_params_mask,
-            0,
-        )
-    PAR_albedo = ClimaCore.Fields.zeros(surface_space) .+ FT(0.2)
-    NIR_albedo = ClimaCore.Fields.zeros(surface_space) .+ FT(0.2)
+    (;
+        ν,
+        ν_ss_om,
+        ν_ss_quartz,
+        ν_ss_gravel,
+        hydrology_cm,
+        K_sat,
+        S_s,
+        θ_r,
+        PAR_albedo,
+        NIR_albedo,
+        f_max,
+    ) = spatially_varying_soil_params
     soil_params = Soil.EnergyHydrologyParameters(
         FT;
         ν,
@@ -218,22 +110,6 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         NIR_albedo = NIR_albedo,
     )
 
-    soil_params_mask_sfc =
-        ClimaLand.Domains.top_center_to_surface(soil_params_mask)
-
-    # Read in f_max data and land sea mask
-    infile_path = ClimaLand.Artifacts.topmodel_data_path()
-    f_max =
-        SpaceVaryingInput(infile_path, "fmax", surface_space; regridder_type)
-    mask = SpaceVaryingInput(
-        infile_path,
-        "landsea_mask",
-        surface_space;
-        regridder_type,
-    )
-    # Unsure how to handle two masks
-    f_max = oceans_to_value.(f_max, mask, FT(0.0))
-    f_max = oceans_to_value.(f_max, soil_params_mask_sfc, FT(0.0))
     f_over = FT(3.28) # 1/m
     R_sb = FT(1.484e-4 / 1000) # m/s
     runoff_model = ClimaLand.Soil.Runoff.TOPMODELRunoff{FT}(;
@@ -267,9 +143,41 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
 
 
     Y, p, cds = initialize(soil)
-
-    init_soil(ν, θ_r) = θ_r + (ν - θ_r) / 2
-    Y.soil.ϑ_l .= init_soil.(ν, θ_r)
+    z = ClimaCore.Fields.coordinate_field(cds.subsurface).z
+    # This function approximates the hydrostatic equilibrium solution in
+    # the vadose and unsaturated regimes by solving for ∂(ψ+z)/∂z = 0,
+    # assuming the transition between the two is at a coordinate of z_∇.
+
+    # The approximation arises because the porosity, residual water fraction,
+    # and van Genuchtem parameters are spatially varying but treated constant
+    # in solving for equilibrium. Finally, we make a plausible but total guess
+    # for the water table depth using the TOPMODEL parameters. 
+    function hydrostatic_profile(
+        lat::FT,
+        z::FT,
+        ν::FT,
+        θ_r::FT,
+        α::FT,
+        n::FT,
+        S_s::FT,
+        fmax,
+    ) where {FT}
+        m = 1 - 1 / n
+        zmin = FT(-50.0)
+        zmax = FT(0.0)
+
+        z_∇ = FT(zmin / 5.0 + (zmax - zmin) / 2.5 * (fmax - 0.35) / 0.7)
+        if z > z_∇
+            S = FT((FT(1) + (α * (z - z_∇))^n)^(-m))
+            ϑ_l = S * (ν - θ_r) + θ_r
+        else
+            ϑ_l = -S_s * (z - z_∇) + ν
+        end
+        return FT(ϑ_l)
+    end
+    vg_α = hydrology_cm.α
+    vg_n = hydrology_cm.n
+    Y.soil.ϑ_l .= hydrostatic_profile.(lat, z, ν, θ_r, vg_α, vg_n, S_s, f_max)
     Y.soil.θ_i .= FT(0.0)
     T = FT(276.85)
     ρc_s =
diff --git a/experiments/standalone/Soil/richards_runoff.jl b/experiments/standalone/Soil/richards_runoff.jl
index e9b2674daa..558023e667 100644
--- a/experiments/standalone/Soil/richards_runoff.jl
+++ b/experiments/standalone/Soil/richards_runoff.jl
@@ -10,7 +10,6 @@ using ClimaUtilities.ClimaArtifacts
 import Interpolations
 
 import ClimaUtilities.TimeVaryingInputs: TimeVaryingInput
-import ClimaUtilities.SpaceVaryingInputs: SpaceVaryingInput
 import ClimaUtilities.Regridders: InterpolationsRegridder
 import ClimaUtilities.OutputPathGenerator: generate_output_path
 import NCDatasets
@@ -43,18 +42,15 @@ domain = ClimaLand.Domains.SphericalShell(;
 );
 surface_space = domain.space.surface
 subsurface_space = domain.space.subsurface
+spatially_varying_soil_params =
+    ClimaLand.default_spatially_varying_soil_parameters(
+        subsurface_space,
+        surface_space,
+        FT,
+    )
+(; ν, hydrology_cm, K_sat, S_s, θ_r, f_max, mask) =
+    spatially_varying_soil_params
 
-# Read in f_max data and land sea mask
-infile_path = ClimaLand.Artifacts.topmodel_data_path()
-f_max = SpaceVaryingInput(infile_path, "fmax", surface_space; regridder_type)
-mask = SpaceVaryingInput(
-    infile_path,
-    "landsea_mask",
-    surface_space;
-    regridder_type,
-)
-
-oceans_to_zero(field, mask) = mask > 0.5 ? field : eltype(field)(0)
 f_over = FT(3.28) # 1/m
 R_sb = FT(1.484e-4 / 1000) # m/s
 runoff_model = ClimaLand.Soil.Runoff.TOPMODELRunoff{FT}(;
@@ -62,83 +58,6 @@ runoff_model = ClimaLand.Soil.Runoff.TOPMODELRunoff{FT}(;
     f_max = f_max,
     R_sb = R_sb,
 )
-soil_params_artifact_path =
-    ClimaLand.Artifacts.soil_params_artifact_folder_path(; context)
-
-extrapolation_bc =
-    (Interpolations.Periodic(), Interpolations.Flat(), Interpolations.Flat())
-# We use this mask to set values of these parameters over the ocean, in order
-# to keep them in the physical range
-function mask_vg(var, value)
-    if var < 1e-8
-        return value
-    else
-        return var
-    end
-end
-vg_α = SpaceVaryingInput(
-    joinpath(
-        soil_params_artifact_path,
-        "vGalpha_map_gupta_etal2020_1.0x1.0x4.nc",
-    ),
-    "α",
-    subsurface_space;
-    regridder_type,
-    regridder_kwargs = (; extrapolation_bc,),
-)
-vg_α .= mask_vg.(vg_α, 1e-3)
-# We use this mask to set values of this parameter over the ocean, in order
-# to keep it in the physical range
-function mask_vg_n(var, value)
-    if var < 1
-        return value
-    else
-        return var
-    end
-end
-vg_n = SpaceVaryingInput(
-    joinpath(soil_params_artifact_path, "vGn_map_gupta_etal2020_1.0x1.0x4.nc"),
-    "n",
-    subsurface_space;
-    regridder_type,
-    regridder_kwargs = (; extrapolation_bc,),
-)
-vg_n .= mask_vg_n.(vg_n, 1.001)
-vg_fields_to_hcm_field(α::FT, n::FT) where {FT} =
-    ClimaLand.Soil.vanGenuchten{FT}(; @NamedTuple{α::FT, n::FT}((α, n))...)
-hydrology_cm = vg_fields_to_hcm_field.(vg_α, vg_n)
-
-θ_r = SpaceVaryingInput(
-    joinpath(
-        soil_params_artifact_path,
-        "residual_map_gupta_etal2020_1.0x1.0x4.nc",
-    ),
-    "θ_r",
-    subsurface_space;
-    regridder_type,
-    regridder_kwargs = (; extrapolation_bc,),
-)
-
-ν = SpaceVaryingInput(
-    joinpath(
-        soil_params_artifact_path,
-        "porosity_map_gupta_etal2020_1.0x1.0x4.nc",
-    ),
-    "ν",
-    subsurface_space;
-    regridder_type,
-    regridder_kwargs = (; extrapolation_bc,),
-)
-ν .= mask_vg.(ν, 1.0)
-K_sat = SpaceVaryingInput(
-    joinpath(soil_params_artifact_path, "ksat_map_gupta_etal2020_1.0x1.0x4.nc"),
-    "Ksat",
-    subsurface_space;
-    regridder_type,
-    regridder_kwargs = (; extrapolation_bc,),
-)
-S_s = ClimaCore.Fields.zeros(subsurface_space) .+ FT(1e-3)
-
 soil_params = ClimaLand.Soil.RichardsParameters(;
     hydrology_cm = hydrology_cm,
     ν = ν,
@@ -205,8 +124,9 @@ end
 t0 = 0.0
 tf = 3600.0 * 24 * 2
 dt = 1800.0
+vg_α = hydrology_cm.α
+vg_n = hydrology_cm.n
 Y.soil.ϑ_l .= hydrostatic_profile.(lat, z, ν, θ_r, vg_α, vg_n, S_s, f_max)
-@. Y.soil.ϑ_l = oceans_to_zero(Y.soil.ϑ_l, mask)
 set_initial_cache! = make_set_initial_cache(model)
 exp_tendency! = make_exp_tendency(model);
 imp_tendency! = ClimaLand.make_imp_tendency(model);
@@ -250,6 +170,7 @@ cb = SciMLBase.CallbackSet(driver_cb, saving_cb)
 sol = @time SciMLBase.solve(prob, ode_algo; dt = dt, saveat = dt, callback = cb)
 
 # Make plots on CPU
+oceans_to_zero(x, mask) = mask == 1 ? x : eltype(x)(0)
 if context.device isa ClimaComms.CPUSingleThreaded
     longpts = range(-180.0, 180.0, 101)
     latpts = range(-90.0, 90.0, 101)
diff --git a/src/ClimaLand.jl b/src/ClimaLand.jl
index dc8f5dbbd2..2c9f927121 100644
--- a/src/ClimaLand.jl
+++ b/src/ClimaLand.jl
@@ -326,4 +326,7 @@ include("integrated/soil_snow_model.jl")
 include(joinpath("diagnostics", "Diagnostics.jl"))
 import .Diagnostics: default_diagnostics
 
+# Simulations
+include(joinpath("simulations", "spatial_parameters.jl"))
+
 end
diff --git a/src/simulations/spatial_parameters.jl b/src/simulations/spatial_parameters.jl
new file mode 100644
index 0000000000..35627ff6fe
--- /dev/null
+++ b/src/simulations/spatial_parameters.jl
@@ -0,0 +1,346 @@
+using ClimaComms
+using ClimaUtilities.ClimaArtifacts
+import Interpolations
+import ClimaUtilities.SpaceVaryingInputs: SpaceVaryingInput
+import ClimaUtilities.Regridders: InterpolationsRegridder
+export clm_canopy_parameters
+function mean(x::AbstractArray{T}) where {T}
+    sum(x) / length(x)
+end
+
+"""
+    clm_canopy_parameters(
+        surface_space;
+        regridder_type = :InterpolationsRegridder,
+        extrapolation_bc = (
+            Interpolations.Periodic(),
+            Interpolations.Flat(),
+            Interpolations.Flat(),
+        ),
+    )
+
+Reads spatially varying parameters for the canopy, from NetCDF files
+based on CLM and MODIS data, and regrids them to the grid defined by the
+`surface_space` of the Clima simulation. Returns a NamedTuple of ClimaCore
+Fields. 
+
+In particular, this file returns a field for 
+- clumping index Ω
+- albedo and transmissitivy in PAR and NIR bands
+- leaf angle distribution G function parameter χl
+- Medlyn g1
+- C3 flag 
+- VCmax25
+
+The values correspond to the value of the dominant PFT at each point.
+
+The NetCDF files are stored in ClimaArtifacts and more detail on their origin
+is provided there. The keyword arguments `regridder_type` and `extrapolation_bc` 
+affect the regridding by (1) changing how we interpolate to ClimaCore points which 
+are not in the data, and (2) changing how extrapolate to points beyond the range of the 
+data.
+"""
+function clm_canopy_parameters(
+    surface_space;
+    regridder_type = :InterpolationsRegridder,
+    extrapolation_bc = (
+        Interpolations.Periodic(),
+        Interpolations.Flat(),
+        Interpolations.Flat(),
+    ),
+)
+    context = ClimaComms.context(surface_space)
+    clm_artifact_path = ClimaLand.Artifacts.clm_data_folder_path(; context)
+    # Foliage clumping index data derived from MODIS
+    modis_ci_artifact_path =
+        ClimaLand.Artifacts.modis_ci_data_folder_path(; context)
+
+    # TwoStreamModel parameters
+    nans_to_one(x) = isnan(x) ? eltype(x)(1) : x
+    Ω = SpaceVaryingInput(
+        joinpath(modis_ci_artifact_path, "He_et_al_2012_1x1.nc"),
+        "ci",
+        surface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+        file_reader_kwargs = (; preprocess_func = nans_to_one,),
+    )
+    χl = SpaceVaryingInput(
+        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
+        "xl",
+        surface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+    )
+    G_Function = CLMGFunction(χl)
+    α_PAR_leaf = SpaceVaryingInput(
+        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
+        "rholvis",
+        surface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+    )
+    τ_PAR_leaf = SpaceVaryingInput(
+        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
+        "taulvis",
+        surface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+    )
+    α_NIR_leaf = SpaceVaryingInput(
+        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
+        "rholnir",
+        surface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+    )
+    τ_NIR_leaf = SpaceVaryingInput(
+        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
+        "taulnir",
+        surface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+    )
+    # 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),),
+    )
+
+    #Photosynthesis model
+    # 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"),
+        "vcmx25",
+        surface_space;
+        regridder_type,
+        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, "vegetation_properties_map.nc"),
+        "c3_dominant",
+        surface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+    )
+    rooting_depth = SpaceVaryingInput(
+        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
+        "rooting_depth",
+        surface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+    )
+    return (;
+        Ω = Ω,
+        rooting_depth = rooting_depth,
+        is_c3 = is_c3,
+        Vcmax25 = Vcmax25,
+        g1 = g1,
+        G_Function = G_Function,
+        α_PAR_leaf = α_PAR_leaf,
+        τ_PAR_leaf = τ_PAR_leaf,
+        α_NIR_leaf = α_PAR_leaf,
+        τ_NIR_leaf = τ_NIR_leaf,
+    )
+end
+
+"""
+    default_spatially_varying_soil_parameters(
+        surface_space,
+        subsurface_space,
+        regridder_type = :InterpolationsRegridder,
+        extrapolation_bc = (
+            Interpolations.Periodic(),
+            Interpolations.Flat(),
+            Interpolations.Flat(),
+        ),
+    )
+
+Reads spatially varying parameters for the soil model, from NetCDF files
+based on SoilGrids and the van Genuchten data product from Gupta et al 2020,
+ and regrids them to the grid defined by the
+`subsurface_space` and `surface_space` of the Clima simulation, as appropriate. 
+Returns a NamedTuple of ClimaCore Fields. 
+
+In particular, this file returns a field for 
+- (α, n, m) (van Genuchten parameters)
+- Ksat
+- porosity
+- residual water content
+- various texture variables: volumetric fractions of organic matter, coarse fragments, and quartz.
+- TOPMODEL parameters
+- soil albedo parameters
+- specific storativity.
+
+The NetCDF files are stored in ClimaArtifacts and more detail on their origin
+is provided there. The keyword arguments `regridder_type` and `extrapolation_bc` 
+affect the regridding by (1) changing how we interpolate to ClimaCore points which 
+are not in the data, and (2) changing how extrapolate to points beyond the range of the 
+data.
+"""
+function default_spatially_varying_soil_parameters(
+    subsurface_space,
+    surface_space,
+    FT;
+    regridder_type = :InterpolationsRegridder,
+    extrapolation_bc = (
+        Interpolations.Periodic(),
+        Interpolations.Flat(),
+        Interpolations.Flat(),
+    ),
+)
+    context = ClimaComms.context(surface_space)
+    soil_params_artifact_path =
+        ClimaLand.Artifacts.soil_params_artifact_folder_path(; context)
+    soil_params_mask = SpaceVaryingInput(
+        joinpath(
+            soil_params_artifact_path,
+            "vGalpha_map_gupta_etal2020_1.0x1.0x4.nc",
+        ),
+        "α",
+        subsurface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+        file_reader_kwargs = (; preprocess_func = (data) -> data > 0,),
+    )
+    oceans_to_value(field, mask, value) =
+        mask == 1.0 ? field : eltype(field)(value)
+
+    vg_α = SpaceVaryingInput(
+        joinpath(
+            soil_params_artifact_path,
+            "vGalpha_map_gupta_etal2020_1.0x1.0x4.nc",
+        ),
+        "α",
+        subsurface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+    )
+    vg_n = SpaceVaryingInput(
+        joinpath(
+            soil_params_artifact_path,
+            "vGn_map_gupta_etal2020_1.0x1.0x4.nc",
+        ),
+        "n",
+        subsurface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+    )
+    x = parent(vg_α)
+    μ = mean(log10.(x[x .> 0]))
+    vg_α .= oceans_to_value.(vg_α, soil_params_mask, 10.0^μ)
+
+    x = parent(vg_n)
+    μ = mean(x[x .> 0])
+    vg_n .= oceans_to_value.(vg_n, soil_params_mask, μ)
+
+    vg_fields_to_hcm_field(α::FT, n::FT) where {FT} =
+        ClimaLand.Soil.vanGenuchten{FT}(; @NamedTuple{α::FT, n::FT}((α, n))...)
+    hydrology_cm = vg_fields_to_hcm_field.(vg_α, vg_n)
+
+    θ_r = SpaceVaryingInput(
+        joinpath(
+            soil_params_artifact_path,
+            "residual_map_gupta_etal2020_1.0x1.0x4.nc",
+        ),
+        "θ_r",
+        subsurface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+    )
+
+    ν = SpaceVaryingInput(
+        joinpath(
+            soil_params_artifact_path,
+            "porosity_map_gupta_etal2020_1.0x1.0x4.nc",
+        ),
+        "ν",
+        subsurface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+    )
+    K_sat = SpaceVaryingInput(
+        joinpath(
+            soil_params_artifact_path,
+            "ksat_map_gupta_etal2020_1.0x1.0x4.nc",
+        ),
+        "Ksat",
+        subsurface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+    )
+
+    x = parent(K_sat)
+    μ = mean(log10.(x[x .> 0]))
+    K_sat .= oceans_to_value.(K_sat, soil_params_mask, 10.0^μ)
+
+    ν .= oceans_to_value.(ν, soil_params_mask, 1)
+
+    θ_r .= oceans_to_value.(θ_r, soil_params_mask, 0)
+
+
+    S_s =
+        oceans_to_value.(
+            ClimaCore.Fields.zeros(subsurface_space) .+ FT(1e-3),
+            soil_params_mask,
+            1,
+        )
+    ν_ss_om =
+        oceans_to_value.(
+            ClimaCore.Fields.zeros(subsurface_space) .+ FT(0.1),
+            soil_params_mask,
+            0,
+        )
+    ν_ss_quartz =
+        oceans_to_value.(
+            ClimaCore.Fields.zeros(subsurface_space) .+ FT(0.1),
+            soil_params_mask,
+            0,
+        )
+    ν_ss_gravel =
+        oceans_to_value.(
+            ClimaCore.Fields.zeros(subsurface_space) .+ FT(0.1),
+            soil_params_mask,
+            0,
+        )
+    PAR_albedo = ClimaCore.Fields.zeros(surface_space) .+ FT(0.2)
+    NIR_albedo = ClimaCore.Fields.zeros(surface_space) .+ FT(0.2)
+
+    soil_params_mask_sfc =
+        ClimaLand.Domains.top_center_to_surface(soil_params_mask)
+
+    # Read in f_max data and land sea mask
+    infile_path = ClimaLand.Artifacts.topmodel_data_path()
+    f_max =
+        SpaceVaryingInput(infile_path, "fmax", surface_space; regridder_type)
+    mask = SpaceVaryingInput(
+        infile_path,
+        "landsea_mask",
+        surface_space;
+        regridder_type,
+    )
+    # Unsure how to handle two masks
+    f_max = oceans_to_value.(f_max, mask, FT(0.0))
+    f_max = oceans_to_value.(f_max, soil_params_mask_sfc, FT(0.0))
+    return (;
+        ν = ν,
+        ν_ss_om = ν_ss_om,
+        ν_ss_quartz = ν_ss_quartz,
+        ν_ss_gravel = ν_ss_gravel,
+        hydrology_cm = hydrology_cm,
+        K_sat = K_sat,
+        S_s = S_s,
+        θ_r = θ_r,
+        PAR_albedo = PAR_albedo,
+        NIR_albedo = NIR_albedo,
+        f_max = f_max,
+        mask = mask,
+    )
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 46f331e71c..c4573fe752 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -111,3 +111,8 @@ end
 @safetestset "Integrated soil and snow" begin
     include("integrated/soil_snow.jl")
 end
+
+# Simulations
+@safetestset "Default parameters" begin
+    include("simulations/spatial_parameters.jl")
+end
diff --git a/test/simulations/spatial_parameters.jl b/test/simulations/spatial_parameters.jl
new file mode 100644
index 0000000000..8606ae33f3
--- /dev/null
+++ b/test/simulations/spatial_parameters.jl
@@ -0,0 +1,81 @@
+using Test
+
+import ClimaComms
+ClimaComms.@import_required_backends
+using ClimaCore
+using ClimaUtilities.ClimaArtifacts
+import Interpolations
+import ClimaUtilities
+import ClimaUtilities.Regridders: InterpolationsRegridder
+import ClimaUtilities.ClimaArtifacts: @clima_artifact
+import ClimaParams as CP
+
+using ClimaLand
+import ClimaLand
+import ClimaLand.Parameters as LP
+const FT = Float64;
+regridder_type = :InterpolationsRegridder
+extrapolation_bc =
+    (Interpolations.Periodic(), Interpolations.Flat(), Interpolations.Flat())
+context = ClimaComms.context()
+
+earth_param_set = LP.LandParameters(FT)
+radius = FT(6378.1e3)
+depth = FT(50)
+domain = ClimaLand.Domains.SphericalShell(;
+    radius = radius,
+    depth = depth,
+    nelements = (101, 15),
+    npolynomial = 1,
+    dz_tuple = FT.((10.0, 0.05)),
+)
+surface_space = domain.space.surface
+subsurface_space = domain.space.subsurface
+spatially_varying_soil_params =
+    ClimaLand.default_spatially_varying_soil_parameters(
+        subsurface_space,
+        surface_space,
+        FT;
+        regridder_type = regridder_type,
+        extrapolation_bc = extrapolation_bc,
+    )
+param_names3d = (:ν, :ν_ss_om, :ν_ss_quartz, :ν_ss_gravel, :K_sat, :S_s, :θ_r)
+param_names2d = (:PAR_albedo, :NIR_albedo, :f_max, :mask)
+for p in param_names3d
+    @test p ∈ propertynames(spatially_varying_soil_params)
+    @test axes(getproperty(spatially_varying_soil_params, p)) ==
+          subsurface_space
+end
+
+for p in param_names2d
+    @test p ∈ propertynames(spatially_varying_soil_params)
+    @test axes(getproperty(spatially_varying_soil_params, p)) == surface_space
+end
+@test :hydrology_cm ∈ propertynames(spatially_varying_soil_params)
+hcm = spatially_varying_soil_params.hydrology_cm
+@test axes(hcm) == subsurface_space
+@test eltype(hcm) == ClimaLand.Soil.vanGenuchten{FT}
+
+clm_parameters = ClimaLand.clm_canopy_parameters(
+    surface_space;
+    regridder_type = regridder_type,
+    extrapolation_bc = extrapolation_bc,
+)
+param_names = (
+    :Ω,
+    :rooting_depth,
+    :is_c3,
+    :Vcmax25,
+    :g1,
+    :α_PAR_leaf,
+    :τ_PAR_leaf,
+    :α_NIR_leaf,
+    :τ_NIR_leaf,
+)
+for p in param_names
+    @test p ∈ propertynames(clm_parameters)
+    @test axes(getproperty(clm_parameters, p)) == surface_space
+end
+
+@test :G_Function ∈ propertynames(clm_parameters)
+@test axes(getproperty(clm_parameters, :G_Function).χl) == surface_space