diff --git a/Artifacts.toml b/Artifacts.toml
index a4e1b6ecd3..bd9505b77f 100644
--- a/Artifacts.toml
+++ b/Artifacts.toml
@@ -1,3 +1,10 @@
+[saturated_land_ic]
+git-tree-sha1 = "cc2888b590f5158113b54c75eee3e2f797be8956"
+
+    [[saturated_land_ic.download]]
+    sha256 = "d73ec9bf8a20f1c63a0133d483a08e801db9125fda43ab51d080bc0192343f84"
+    url = "https://caltech.box.com/shared/static/99j3juuj3vnnzy16ye31dsfcx25s71z6.gz"
+
 [landsea_mask_30arcseconds]
 git-tree-sha1 = "517c925535981f9d17ac9e7a65e2c35c3ce9597d"
 
diff --git a/experiments/calibration/PBS_calibration.pbs b/experiments/calibration/PBS_calibration.pbs
index 7738cc8ede..373c137741 100644
--- a/experiments/calibration/PBS_calibration.pbs
+++ b/experiments/calibration/PBS_calibration.pbs
@@ -17,5 +17,5 @@ export CLIMACOMMS_DEVICE="CUDA"
 export CLIMACOMMS_CONTEXT="SINGLETON"
 
 julia --project=.buildkite -e 'using Pkg; Pkg.instantiate(;verbose=true)'
-julia --project=.buildkite/ experiments/calibration/generate_observations.jl
+#julia --project=.buildkite/ experiments/calibration/generate_observations.jl
 julia --project=.buildkite/ experiments/calibration/run_calibration.jl
diff --git a/experiments/calibration/api.jl b/experiments/calibration/api.jl
index c372d727b1..99f6faa22f 100644
--- a/experiments/calibration/api.jl
+++ b/experiments/calibration/api.jl
@@ -9,7 +9,7 @@ import ClimaLand
         sample_date_ranges::Vector{NTuple{2, DATE}}
         extend::EXTEND
         spinup::SPINUP
-        nelements::Tuple{Int64, Int64}
+        nelements::Union{Tuple{Int64, Int64, Int64}, Tuple{Int64, Int64}}
         output_dir::String
         rng_seed::Int64
     end
@@ -43,7 +43,7 @@ struct CalibrateConfig{SPINUP <: Dates.Period, EXTEND <: Dates.Period, MODEL}
 
     "The number of horizontal and vertical elements of the model. Used for the
     simulation and determining the ocean mask"
-    nelements::Tuple{Int64, Int64}
+    nelements::Union{Tuple{Int64, Int64, Int64}, Tuple{Int64, Int64}}
 
     "The directory to store the iterations and members of the calibration."
     output_dir::String
diff --git a/experiments/calibration/models/snowy_land.jl b/experiments/calibration/models/snowy_land.jl
index 6a579bbe31..de822354f5 100644
--- a/experiments/calibration/models/snowy_land.jl
+++ b/experiments/calibration/models/snowy_land.jl
@@ -63,8 +63,24 @@ function setup_model(
         scf,
     )
 
+    ground = ClimaLand.PrognosticGroundConditions{FT}()
+    canopy_forcing = (; atmos, radiation, ground)
+    photosynthesis = PModel{FT}(domain, toml_dict)
+    conductance = PModelConductance{FT}(toml_dict)
+    soil_moisture_stress =
+        ClimaLand.Canopy.PiecewiseMoistureStressModel{FT}(domain, toml_dict)
+    canopy = ClimaLand.Canopy.CanopyModel{FT}(
+        surface_domain,
+        canopy_forcing,
+        LAI,
+        toml_dict;
+        prognostic_land_components = (:canopy, :snow, :soil, :soilco2),
+        photosynthesis,
+        conductance,
+        soil_moisture_stress,
+    )
     # Construct the land model with all default components except for snow
-    land = LandModel{FT}(forcing, LAI, toml_dict, domain, Δt; snow)
+    land = LandModel{FT}(forcing, LAI, toml_dict, domain, Δt; snow, canopy)
     return land
 end
 
diff --git a/experiments/calibration/run_calibration.jl b/experiments/calibration/run_calibration.jl
index 749a770c46..2c0b3ce768 100644
--- a/experiments/calibration/run_calibration.jl
+++ b/experiments/calibration/run_calibration.jl
@@ -11,14 +11,16 @@ import JLD2
 include(joinpath(pkgdir(ClimaLand), "experiments/calibration/api.jl"))
 
 const CALIBRATE_CONFIG = CalibrateConfig(;
-    short_names = ["lwu"],
+    short_names = ["lwu", "shf", "lhf"],
     minibatch_size = 1,
-    n_iterations = 1,
-    sample_date_ranges = [("2007-12-1", "2007-12-1")],
+    n_iterations = 10,
+    sample_date_ranges = [
+        ("$(2000 + 2*i)-12-1", "$(2002 + 2*i)-9-1") for i in 0:9
+    ], # 2000 to 2020
     extend = Dates.Month(3),
-    spinup = Dates.Month(0),
+    spinup = Dates.Month(3),
     nelements = (180, 360, 15),
-    output_dir = "experiments/calibration/land_model",
+    output_dir = "/glade/derecho/scratch/kdeck/recalibrate_saturated",
     rng_seed = 42,
     obs_vec_filepath = "experiments/calibration/land_observation_vector.jld2",
     model_type = ClimaLand.LandModel,
@@ -26,9 +28,13 @@ const CALIBRATE_CONFIG = CalibrateConfig(;
 
 
 if abspath(PROGRAM_FILE) == @__FILE__
-    # true solution is at 0.96
-    priors =
-        [EKP.constrained_gaussian("emissivity_bare_soil", 0.82, 0.12, 0.0, 2.0)]
+    priors = [
+        EKP.constrained_gaussian("moisture_stress_c", 1.0, 0.5, 0, 2),
+        EKP.constrained_gaussian("pmodel_cstar", 0.41, 0.11, 0, Inf),
+        EKP.constrained_gaussian("pmodel_β", 146, 10, 0, Inf),
+        EKP.constrained_gaussian("leaf_Cd", 0.01, 0.05, 0, Inf),
+        EKP.constrained_gaussian("ac_canopy", 3.5e4, 2e4, 2.5e3, 2.5e5),
+    ]
     prior = EKP.combine_distributions(priors)
 
     observation_vector = JLD2.load_object(CALIBRATE_CONFIG.obs_vec_filepath)
diff --git a/src/Artifacts.jl b/src/Artifacts.jl
index cb09fb6cbf..fd9839ecd4 100644
--- a/src/Artifacts.jl
+++ b/src/Artifacts.jl
@@ -6,6 +6,11 @@ import ClimaUtilities.ClimaArtifacts: @clima_artifact
 
 import LazyArtifacts
 
+function saturated_land_ic_path(; context = nothing)
+    dir = @clima_artifact("saturated_land_ic", context)
+    return joinpath(dir, "saturated_land_ic.nc")
+end
+
 """
     soil_ic_2008_50m_path(; context)
 
diff --git a/src/shared_utilities/Domains.jl b/src/shared_utilities/Domains.jl
index 017c7432e3..2b434b82a3 100644
--- a/src/shared_utilities/Domains.jl
+++ b/src/shared_utilities/Domains.jl
@@ -1247,8 +1247,8 @@ end
         apply_mask = true,
         mask_threshold = 0.5,
         nelements = (101, 15),
-        dz_tuple = (10.0, 0.05),
-        depth = 50.0,
+        dz_tuple = (3.0, 0.05),
+        depth = 15.0,
         npolynomial = 0,
         context = ClimaComms.context(),
         filepath = landseamask_file_path(;context),
@@ -1401,7 +1401,7 @@ function global_box_domain(
     ylim = FT.((π / 2 * radius_earth, π / 2 * radius_earth))
     longlat = FT.((-0.5, -0.5))
     dz_tuple = FT.(dz_tuple)
-    device = ClimaComms.device()
+    device = ClimaComms.device(context)
     domain = HybridBox(;
         xlim,
         ylim,
diff --git a/src/simulations/Simulations.jl b/src/simulations/Simulations.jl
index 4d2c08ba68..422f3580a7 100644
--- a/src/simulations/Simulations.jl
+++ b/src/simulations/Simulations.jl
@@ -70,7 +70,7 @@ end
         model;
         outdir = ".",
         set_ic! = make_set_initial_state_from_file(
-            ClimaLand.Artifacts.soil_ic_2008_50m_path(;
+            ClimaLand.Artifacts.saturated_land_ic_path(;
                 context = ClimaComms.context(model),
             ),
             model,
@@ -116,7 +116,7 @@ function LandSimulation(
     model;
     outdir = ".",
     set_ic! = make_set_initial_state_from_file(
-        ClimaLand.Artifacts.soil_ic_2008_50m_path(;
+        ClimaLand.Artifacts.saturated_land_ic_path(;
             context = ClimaComms.context(model),
         ),
         model,
diff --git a/src/simulations/initial_conditions.jl b/src/simulations/initial_conditions.jl
index 384ced4296..7101e950d2 100644
--- a/src/simulations/initial_conditions.jl
+++ b/src/simulations/initial_conditions.jl
@@ -49,11 +49,11 @@ function set_soil_initial_conditions!(
         regridder_kwargs = (; extrapolation_bc, interpolation_method),
     )
 
-    Y.soil.ϑ_l .= enforce_residual_constraint.(Y.soil.ϑ_l, θ_r)
-    Y.soil.ϑ_l .= enforce_porosity_constraint.(Y.soil.ϑ_l, ν)
-    Y.soil.θ_i .=
-        enforce_residual_constraint.(Y.soil.θ_i, eltype(Y.soil.θ_i)(0))
-    Y.soil.θ_i .= enforce_porosity_constraint.(Y.soil.ϑ_l, Y.soil.θ_i, ν)
+#    Y.soil.ϑ_l .= enforce_residual_constraint.(Y.soil.ϑ_l, θ_r)
+#    Y.soil.ϑ_l .= enforce_porosity_constraint.(Y.soil.ϑ_l, ν)
+#    Y.soil.θ_i .=
+#        enforce_residual_constraint.(Y.soil.θ_i, eltype(Y.soil.θ_i)(0))
+#    Y.soil.θ_i .= enforce_porosity_constraint.(Y.soil.ϑ_l, Y.soil.θ_i, ν)
     ρc_s =
         ClimaLand.Soil.volumetric_heat_capacity.(
             Y.soil.ϑ_l,
@@ -68,21 +68,21 @@ function set_soil_initial_conditions!(
         regridder_type,
         regridder_kwargs = (; extrapolation_bc, interpolation_method),
     )
-    T =
-        ClimaLand.Soil.temperature_from_ρe_int.(
-            Y.soil.ρe_int,
-            Y.soil.θ_i,
-            ρc_s,
-            soil.parameters.earth_param_set,
-        )
-    T .= clip_to_bounds.(T, T_bounds[1], T_bounds[2])
-    Y.soil.ρe_int .=
-        ClimaLand.Soil.volumetric_internal_energy.(
-            Y.soil.θ_i,
-            ρc_s,
-            T,
-            soil.parameters.earth_param_set,
-        )
+#    T =
+#        ClimaLand.Soil.temperature_from_ρe_int.(
+#            Y.soil.ρe_int,
+#            Y.soil.θ_i,
+#            ρc_s,
+#            soil.parameters.earth_param_set,
+#        )
+#    T .= clip_to_bounds.(T, T_bounds[1], T_bounds[2])
+#    Y.soil.ρe_int .=
+#        ClimaLand.Soil.volumetric_internal_energy.(
+#            Y.soil.θ_i,
+#            ρc_s,
+#            T,
+#            soil.parameters.earth_param_set,
+#        )
     return nothing
 end
 
@@ -247,20 +247,12 @@ function make_set_initial_state_from_file(
         # to soil potential (soil moisture), averaged over the soil layers,
         # which would correspond to approximate steady state
         if land.canopy.hydraulics isa ClimaLand.Canopy.PlantHydraulicsModel
-            @. p.soil.ψ = ClimaLand.Soil.pressure_head(
-                land.soil.parameters.hydrology_cm,
-                land.soil.parameters.θ_r,
-                Y.soil.ϑ_l,
-                land.soil.parameters.ν - Y.soil.θ_i,
-                land.soil.parameters.S_s,
-            )
-            ψ_roots = ClimaCore.Fields.zeros(axes(Y.canopy.hydraulics.ϑ_l.:1))
-            z = land.soil.domain.fields.z
-            tmp = @. ClimaLand.Canopy.root_distribution(
-                z,
-                land.canopy.biomass.rooting_depth,
-            ) * p.soil.ψ / land.soil.domain.fields.depth
-            ClimaCore.Operators.column_integral_definite!(ψ_roots, tmp)
+            ψ_roots = SpaceVaryingInput(
+                ic_path,
+                "lwp",
+                land.snow.domain.space.surface;
+                regridder_type,
+                regridder_kwargs = (; extrapolation_bc, interpolation_method))
             Y.canopy.hydraulics.ϑ_l.:1 .=
                 ClimaLand.Canopy.PlantHydraulics.inverse_water_retention_curve.(
                     land.canopy.hydraulics.parameters.retention_model,
diff --git a/src/standalone/Soil/Soil.jl b/src/standalone/Soil/Soil.jl
index 17ac24f9cd..f50fb92136 100644
--- a/src/standalone/Soil/Soil.jl
+++ b/src/standalone/Soil/Soil.jl
@@ -278,7 +278,10 @@ function EnergyHydrology{FT}(
         runoff;
         prognostic_land_components,
     )
-    bottom_bc = EnergyWaterFreeDrainage()
+    bottom_water_flux = WaterFluxBC((p, t) -> 0.0)
+    bottom_heat_flux = HeatFluxBC((p, t) -> 0.0)
+    bottom_bc =
+        WaterHeatBC(; water = bottom_water_flux, heat = bottom_heat_flux)
     boundary_conditions = (; top = top_bc, bottom = bottom_bc)
     # sublimation and subsurface runoff are added automatically
     sources = (additional_sources..., PhaseChange{FT}())
diff --git a/toml/default_parameters.toml b/toml/default_parameters.toml
index fb89e9e839..ed23c7ba84 100644
--- a/toml/default_parameters.toml
+++ b/toml/default_parameters.toml
@@ -1,6 +1,6 @@
 # Moisture stress models
 [moisture_stress_c]
-value = 0.47
+value = 1.0
 type = "float"
 description = "unitless"
 tag = "SoilMoistureStress"
@@ -33,7 +33,7 @@ tag = "TOPMODELRunoff"
 # Big leaf energy parameters
 
 ["ac_canopy"]
-value = 2.5e3
+value = 2.7e3
 type = "float"
 description = "Specific heat per emitting area [J/m^2/K]"
 tag = "BigLeafEnergyModel"
@@ -165,13 +165,13 @@ tag = "WuWuSnowCoverFractionModel"
 
 # P model parameters
 ["pmodel_cstar"]
-value = 0.30
+value = 0.41
 type = "float"
 description = "4 * dA/dJmax, assumed to be a constant marginal cost (Wang 2017, Stocker 2020)"
 tag = "PModelParameters"
 
 ["pmodel_β"]
-value = 192
+value = 146
 type = "float"
 description = "Unit cost ratio of Vcmax to transpiration (Stocker 2020)"
 tag = "PModelParameters"
@@ -481,7 +481,7 @@ description = "Scalar multipling canopy height to get displacement height"
 tag = "MoninObukhovCanopyFluxes"
 
 ["leaf_Cd"]
-value = 0.01
+value = 0.05
 type = "float"
 description = ""
 tag = "MoninObukhovCanopyFluxes"