diff --git a/docs/src/APIs/shared_utilities.md b/docs/src/APIs/shared_utilities.md
index faf0fb9819..9a1c473e79 100644
--- a/docs/src/APIs/shared_utilities.md
+++ b/docs/src/APIs/shared_utilities.md
@@ -19,6 +19,7 @@ ClimaLand.Domains.obtain_face_space
 ClimaLand.Domains.obtain_surface_space
 ClimaLand.Domains.obtain_surface_domain
 ClimaLand.Domains.top_center_to_surface
+ClimaLand.Domains.top_face_to_surface
 ClimaLand.Domains.linear_interpolation_to_surface!
 ClimaLand.Domains.get_Δz
 
diff --git a/docs/tutorials/standalone/Soil/freezing_front.jl b/docs/tutorials/standalone/Soil/freezing_front.jl
index 4c830f371d..4c27cb0bbd 100644
--- a/docs/tutorials/standalone/Soil/freezing_front.jl
+++ b/docs/tutorials/standalone/Soil/freezing_front.jl
@@ -122,7 +122,6 @@ params = Soil.EnergyHydrologyParameters(
 zmax = FT(0)
 zmin = FT(-0.2)
 nelems = 20
-Δz = FT(0.01)
 soil_domain = Column(; zlim = (zmin, zmax), nelements = nelems);
 
 # Set the boundary conditions:
@@ -155,7 +154,7 @@ boundary_fluxes = (;
 # Sources are added as elements of a list of sources. Here we just add freezing
 # and thawing.
 
-sources = (PhaseChange{FT}(Δz),);
+sources = (PhaseChange{FT}(),);
 
 # Now we can package this up in the
 # [`EnergyHydrology`](@ref ClimaLand.Soil.EnergyHydrology) model
diff --git a/docs/tutorials/standalone/Soil/sublimation.jl b/docs/tutorials/standalone/Soil/sublimation.jl
index 5e0b93602e..414f027a00 100644
--- a/docs/tutorials/standalone/Soil/sublimation.jl
+++ b/docs/tutorials/standalone/Soil/sublimation.jl
@@ -118,9 +118,8 @@ zmin = FT(-0.35)
 nelems = 12
 soil_domain = Column(; zlim = (zmin, zmax), nelements = nelems);
 z = ClimaCore.Fields.coordinate_field(soil_domain.space.subsurface).z;
-Δz = parent(z)[end]
 # Soil model, and create the prognostic vector Y and cache p:
-sources = (PhaseChange{FT}(Δz),);
+sources = (PhaseChange{FT}(),);
 soil = Soil.EnergyHydrology{FT}(;
     parameters = params,
     domain = soil_domain,
diff --git a/experiments/standalone/Biogeochemistry/experiment.jl b/experiments/standalone/Biogeochemistry/experiment.jl
index 2171446c22..28374bb089 100644
--- a/experiments/standalone/Biogeochemistry/experiment.jl
+++ b/experiments/standalone/Biogeochemistry/experiment.jl
@@ -46,7 +46,6 @@ for (FT, tf) in ((Float32, 2 * dt), (Float64, tf))
     zmax = FT(0)
     zmin = FT(-1)
     nelems = 20
-    Δz = abs(zmax - zmin) / nelems
 
     lsm_domain = Column(; zlim = (zmin, zmax), nelements = nelems)
 
@@ -57,7 +56,7 @@ for (FT, tf) in ((Float32, 2 * dt), (Float64, tf))
     bot_flux_bc_h = Soil.HeatFluxBC((p, t) -> 0.0)
 
 
-    sources = (PhaseChange{FT}(Δz),)
+    sources = (PhaseChange{FT}(),)
     boundary_fluxes = (;
         top = WaterHeatBC(; water = top_flux_bc_w, heat = bot_flux_bc_h),
         bottom = WaterHeatBC(; water = bot_flux_bc_w, heat = bot_flux_bc_h),
diff --git a/src/integrated/soil_canopy_model.jl b/src/integrated/soil_canopy_model.jl
index 7cf89c66fe..67cd03d063 100644
--- a/src/integrated/soil_canopy_model.jl
+++ b/src/integrated/soil_canopy_model.jl
@@ -91,12 +91,7 @@ function SoilCanopyModel{FT}(;
     # These may be passed in, or set, depending on use scenario.
     (; atmos, radiation) = land_args
     # These should always be set by the constructor.
-    Δz = minimum(
-        ClimaCore.Fields.Δz_field(
-            ClimaLand.coordinates(soil_args.domain).subsurface,
-        ),
-    )
-    sources = (RootExtraction{FT}(), Soil.PhaseChange{FT}(Δz))
+    sources = (RootExtraction{FT}(), Soil.PhaseChange{FT}())
     # add heat BC
     top_bc =
         ClimaLand.Soil.AtmosDrivenFluxBC(atmos, CanopyRadiativeFluxes{FT}())
diff --git a/src/shared_utilities/Domains.jl b/src/shared_utilities/Domains.jl
index 371eb3be59..e9f023f64b 100644
--- a/src/shared_utilities/Domains.jl
+++ b/src/shared_utilities/Domains.jl
@@ -107,6 +107,8 @@ struct Column{FT} <: AbstractDomain{FT}
     boundary_names::Tuple{Symbol, Symbol}
     "A NamedTuple of associated ClimaCore spaces: in this case, the surface space and subsurface center space"
     space::NamedTuple
+    "Fields associated with the coordinates of the domain that are useful to store"
+    fields::NamedTuple
 end
 
 """
@@ -159,9 +161,18 @@ function Column(;
     subsurface_space = ClimaCore.Spaces.CenterFiniteDifferenceSpace(mesh)
     surface_space = obtain_surface_space(subsurface_space)
     space = (; surface = surface_space, subsurface = subsurface_space)
-    return Column{FT}(zlim, (nelements,), dz_tuple, boundary_names, space)
+    fields = get_additional_domain_fields(subsurface_space)
+    return Column{FT}(
+        zlim,
+        (nelements,),
+        dz_tuple,
+        boundary_names,
+        space,
+        fields,
+    )
 end
 
+
 """
     Plane{FT} <: AbstractDomain{FT}
 A struct holding the necessary information
@@ -275,6 +286,8 @@ struct HybridBox{FT} <: AbstractDomain{FT}
     periodic::Tuple{Bool, Bool}
     "A NamedTuple of associated ClimaCore spaces: in this case, the surface space and subsurface center space"
     space::NamedTuple
+    "Fields associated with the coordinates of the domain that are useful to store"
+    fields::NamedTuple
 end
 
 """
@@ -357,6 +370,7 @@ function HybridBox(;
 
     surface_space = obtain_surface_space(subsurface_space)
     space = (; surface = surface_space, subsurface = subsurface_space)
+    fields = get_additional_domain_fields(subsurface_space)
     return HybridBox{FT}(
         xlim,
         ylim,
@@ -366,6 +380,7 @@ function HybridBox(;
         npolynomial,
         periodic,
         space,
+        fields,
     )
 end
 
@@ -401,6 +416,8 @@ struct SphericalShell{FT} <: AbstractDomain{FT}
     npolynomial::Int
     "A NamedTuple of associated ClimaCore spaces: in this case, the surface space and subsurface center space"
     space::NamedTuple
+    "Fields associated with the coordinates of the domain that are useful to store"
+    fields::NamedTuple
 end
 
 """
@@ -465,6 +482,7 @@ function SphericalShell(;
     )
     surface_space = obtain_surface_space(subsurface_space)
     space = (; surface = surface_space, subsurface = subsurface_space)
+    fields = get_additional_domain_fields(subsurface_space)
     return SphericalShell{FT}(
         radius,
         depth,
@@ -472,6 +490,7 @@ function SphericalShell(;
         nelements,
         npolynomial,
         space,
+        fields,
     )
 end
 
@@ -677,15 +696,15 @@ When `val` is a scalar (e.g. a single float or struct), returns `val`.
 top_center_to_surface(val) = val
 
 """
-    linear_interpolation_to_surface!(sfc_field, center_field, z)
+    linear_interpolation_to_surface!(sfc_field, center_field, z, Δz_top)
 
 Linearly interpolate the center field `center_field` to the surface
-defined by the top face coordinate of `z`; updates the `sfc_field`
+defined by the top face coordinate of `z` with a center to face distance
+`Δz_top` in the first layer; updates the `sfc_field`
 on the surface (face) space in place.
 """
-function linear_interpolation_to_surface!(sfc_field, center_field, z)
-    Δz_top, _ = get_Δz(z)
-    surface_space = obtain_surface_space(axes(z))
+function linear_interpolation_to_surface!(sfc_field, center_field, z, Δz_top)
+    surface_space = axes(sfc_field)
     Δz_top = ClimaCore.Fields.field_values(Δz_top)
     nz = Spaces.nlevels(axes(center_field))
     f1 = ClimaCore.Fields.field_values(ClimaCore.Fields.level(center_field, nz))
@@ -694,10 +713,8 @@ function linear_interpolation_to_surface!(sfc_field, center_field, z)
     )
     z1 = ClimaCore.Fields.field_values(ClimaCore.Fields.level(z, nz))
     z2 = ClimaCore.Fields.field_values(ClimaCore.Fields.level(z, nz - 1))
-    sfc_field .= ClimaCore.Fields.Field(
-        (@. (f1 - f2) / (z1 - z2) * (Δz_top + z1 - z2) + f2),
-        surface_space,
-    )
+    ClimaCore.Fields.field_values(sfc_field) .=
+        @. (f1 - f2) / (z1 - z2) * (Δz_top + z1 - z2) + f2
 end
 
 """
@@ -721,12 +738,61 @@ function get_Δz(z::ClimaCore.Fields.Field)
     return Δz_top ./ 2, Δz_bottom ./ 2
 end
 
+"""
+    top_face_to_surface(face_field::ClimaCore.Fields.Field, surface_space)
+
+Creates and returns a ClimaCore.Fields.Field defined on the space
+corresponding to the surface of the space on which `face_field`
+is defined, with values equal to the those at the level of the top
+face.
+
+Given a `face_field` defined on a 3D
+extruded face finite difference space, this would return a 2D field
+corresponding to the surface, with values equal to the topmost level.
+ """
+function top_face_to_surface(face_field::ClimaCore.Fields.Field, surface_space)
+    face_space = axes(face_field)
+    N = ClimaCore.Spaces.nlevels(face_space)
+    sfc_level =
+        ClimaCore.Fields.level(face_field, ClimaCore.Utilities.PlusHalf(N - 1))
+    # Project onto surface space
+    return ClimaCore.Fields.Field(
+        ClimaCore.Fields.field_values(sfc_level),
+        surface_space,
+    )
+end
+
+"""
+    get_additional_domain_fields(subsurface_space)
+
+A helper function which returns additional fields corresponding to ClimaLand
+domains which have a subsurface_space (Column, HybridBox, SphericalShell);
+these fields are the center coordinates of the subsurface space, the spacing between
+the top center and top surface and bottom center and bottom surface, as well as the
+field corresponding to the surface height z.
+
+We allocate these once, upon domain construction, so that they are accessible
+during the simulation.
+"""
+function get_additional_domain_fields(subsurface_space)
+    surface_space = obtain_surface_space(subsurface_space)
+    z = ClimaCore.Fields.coordinate_field(subsurface_space).z
+    Δz_top, Δz_bottom = get_Δz(z)
+    face_space = obtain_face_space(subsurface_space)
+    z_face = ClimaCore.Fields.coordinate_field(face_space).z
+    z_sfc = top_face_to_surface(z_face, surface_space)
+    fields = (; z = z, Δz_top = Δz_top, Δz_bottom = Δz_bottom, z_sfc = z_sfc)
+    return fields
+end
+
+
 export AbstractDomain
 export Column, Plane, HybridBox, Point, SphericalShell, SphericalSurface
 export coordinates,
     obtain_face_space,
     obtain_surface_space,
     top_center_to_surface,
+    top_face_to_surface,
     obtain_surface_domain,
     linear_interpolation_to_surface!,
     get_Δz
diff --git a/src/standalone/Soil/Biogeochemistry/Biogeochemistry.jl b/src/standalone/Soil/Biogeochemistry/Biogeochemistry.jl
index 40d1e584c8..871590745a 100644
--- a/src/standalone/Soil/Biogeochemistry/Biogeochemistry.jl
+++ b/src/standalone/Soil/Biogeochemistry/Biogeochemistry.jl
@@ -171,8 +171,8 @@ ClimaLand.auxiliary_domain_names(model::SoilCO2Model) = (
 
 function make_update_boundary_fluxes(model::SoilCO2Model)
     function update_boundary_fluxes!(p, Y, t)
-        z = ClimaCore.Fields.coordinate_field(model.domain.space.subsurface).z
-        Δz_top, Δz_bottom = ClimaLand.Domains.get_Δz(z)
+        Δz_top = model.domain.fields.Δz_top
+        Δz_bottom = model.domain.fields.Δz_bottom
         p.soilco2.top_bc .= boundary_flux(
             model.boundary_conditions.top,
             TopBoundary(),
@@ -405,7 +405,7 @@ function ClimaLand.make_update_aux(model::SoilCO2Model)
         # get FT to enforce types of variables not stored directly in `p`
         FT = eltype(Y.soilco2.C)
         params = model.parameters
-        z = ClimaCore.Fields.coordinate_field(model.domain.space.subsurface).z
+        z = model.domain.fields.z
         T_soil = FT.(soil_temperature(model.driver.met, p, Y, t, z))
         θ_l = FT.(soil_moisture(model.driver.met, p, Y, t, z))
         Csom = FT.(soil_SOM_C(model.driver.soc, p, Y, t, z))
diff --git a/src/standalone/Soil/energy_hydrology.jl b/src/standalone/Soil/energy_hydrology.jl
index 2b66e411b8..6f148f995f 100644
--- a/src/standalone/Soil/energy_hydrology.jl
+++ b/src/standalone/Soil/energy_hydrology.jl
@@ -131,8 +131,8 @@ end
 
 function make_update_boundary_fluxes(model::EnergyHydrology)
     function update_boundary_fluxes!(p, Y, t)
-        z = ClimaCore.Fields.coordinate_field(model.domain.space.subsurface).z
-        Δz_top, Δz_bottom = ClimaLand.Domains.get_Δz(z)
+        Δz_top = model.domain.fields.Δz_top
+        Δz_bottom = model.domain.fields.Δz_bottom
         soil_boundary_fluxes!(
             model.boundary_conditions.top,
             ClimaLand.TopBoundary(),
@@ -175,7 +175,7 @@ function ClimaLand.make_compute_exp_tendency(
     model::EnergyHydrology{FT},
 ) where {FT}
     function compute_exp_tendency!(dY, Y, p, t)
-        z = ClimaCore.Fields.coordinate_field(model.domain.space.subsurface).z
+        z = model.domain.fields.z
         heat_top_flux_bc = p.soil.top_bc.heat
         heat_bottom_flux_bc = p.soil.bottom_bc.heat
 
@@ -236,7 +236,7 @@ function ClimaLand.make_compute_imp_tendency(
     model::EnergyHydrology{FT},
 ) where {FT}
     function compute_imp_tendency!(dY, Y, p, t)
-        z = ClimaCore.Fields.coordinate_field(model.domain.space.subsurface).z
+        z = model.domain.fields.z
         rre_top_flux_bc = p.soil.top_bc.water
         rre_bottom_flux_bc = p.soil.bottom_bc.water
 
@@ -379,11 +379,12 @@ function horizontal_components!(
     hgrad = Operators.Gradient()
     # The flux is already covariant, from hgrad, so no need to convert.
     @. dY.soil.ϑ_l += -hdiv(-p.soil.K * hgrad(p.soil.ψ + z))
-    ρe_int_l = volumetric_internal_energy_liq.(p.soil.T, model.parameters)
     @. dY.soil.ρe_int +=
         -hdiv(
             -p.soil.κ * hgrad(p.soil.T) -
-            p.soil.K * ρe_int_l * hgrad(p.soil.ψ + z),
+            p.soil.K *
+            volumetric_internal_energy_liq(p.soil.T, model.parameters) *
+            hgrad(p.soil.ψ + z),
         )
 end
 
@@ -528,9 +529,7 @@ end
 
 PhaseChange source type.
 """
-struct PhaseChange{FT} <: AbstractSoilSource{FT}
-    Δz::FT
-end
+struct PhaseChange{FT} <: AbstractSoilSource{FT} end
 
 
 """
@@ -555,13 +554,31 @@ function ClimaLand.source!(
     (; ν, earth_param_set) = params
     _ρ_l = FT(LP.ρ_cloud_liq(earth_param_set))
     _ρ_i = FT(LP.ρ_cloud_ice(earth_param_set))
-    ρc = volumetric_heat_capacity.(p.soil.θ_l, Y.soil.θ_i, params)
-    τ = thermal_time.(ρc, src.Δz, p.soil.κ)
-
-    liquid_source =
-        phase_change_source.(p.soil.θ_l, Y.soil.θ_i, p.soil.T, τ, params)
-    @. dY.soil.ϑ_l += -liquid_source
-    @. dY.soil.θ_i += (_ρ_l / _ρ_i) * liquid_source
+    Δz_top = model.domain.fields.Δz_top # center face distance
+    @. dY.soil.ϑ_l +=
+        -phase_change_source(
+            p.soil.θ_l,
+            Y.soil.θ_i,
+            p.soil.T,
+            thermal_time(
+                volumetric_heat_capacity(p.soil.θ_l, Y.soil.θ_i, params),
+                2 * Δz_top, # the factor of 2 appears to get the face-face/layer thickness, Δz_top is center-face distance
+                p.soil.κ,
+            ),
+            params,
+        )
+    @. dY.soil.θ_i +=
+        (_ρ_l / _ρ_i) * phase_change_source(
+            p.soil.θ_l,
+            Y.soil.θ_i,
+            p.soil.T,
+            thermal_time(
+                volumetric_heat_capacity(p.soil.θ_l, Y.soil.θ_i, params),
+                2 * Δz_top, #the factor of 2 appears to get the face-face/layer thickness, Δz_top is center-face distance
+                p.soil.κ,
+            ),
+            params,
+        )
 end
 
 
@@ -594,11 +611,11 @@ function ClimaLand.source!(
 ) where {FT}
     _ρ_i = FT(LP.ρ_cloud_ice(model.parameters.earth_param_set))
     _ρ_l = FT(LP.ρ_cloud_liq(model.parameters.earth_param_set))
-    z = ClimaCore.Fields.coordinate_field(axes(Y.soil.θ_i)).z
-    Δz, _ = ClimaLand.Domains.get_Δz(z) # center to face distance, we need a factor of two
+    z = model.domain.fields.z
+    Δz_top = model.domain.fields.Δz_top # this returns the center-face distance, not layer thickness
     @. dY.soil.θ_i +=
         -p.soil.turbulent_fluxes.vapor_flux * p.soil.ice_frac * _ρ_l / _ρ_i *
-        heaviside(z + 2 * Δz) # only apply to top layer
+        heaviside(z + 2 * Δz_top) # only apply to top layer, recall that z is negative
 end
 
 ## The functions below are required to be defined
@@ -649,13 +666,13 @@ function ClimaLand.surface_resistance(
     p,
     t,
 ) where {FT}
-    cds = ClimaCore.Fields.coordinate_field(model.domain.space.subsurface)
     (; ν, θ_r, hydrology_cm) = model.parameters
     θ_l_sfc = p.soil.sfc_scratch
     ClimaLand.Domains.linear_interpolation_to_surface!(
         θ_l_sfc,
         p.soil.θ_l,
-        cds.z,
+        model.domain.fields.z,
+        model.domain.fields.Δz_top,
     )
     θ_i_sfc = ClimaLand.Domains.top_center_to_surface(Y.soil.θ_i)
     return ClimaLand.Soil.soil_resistance.(θ_l_sfc, θ_i_sfc, model.parameters)
@@ -734,9 +751,13 @@ function ClimaLand.surface_specific_humidity(
     M_w = LP.molar_mass_water(model.parameters.earth_param_set)
     thermo_params =
         LP.thermodynamic_parameters(model.parameters.earth_param_set)
-    cds = ClimaCore.Fields.coordinate_field(model.domain.space.subsurface)
     ψ_sfc = p.soil.sfc_scratch
-    ClimaLand.Domains.linear_interpolation_to_surface!(ψ_sfc, p.soil.ψ, cds.z)
+    ClimaLand.Domains.linear_interpolation_to_surface!(
+        ψ_sfc,
+        p.soil.ψ,
+        model.domain.fields.z,
+        model.domain.fields.Δz_top,
+    )
     q_over_ice =
         Thermodynamics.q_vap_saturation_generic.(
             thermo_params,
@@ -772,20 +793,7 @@ end
 Returns the surface height of the `EnergyHydrology` model.
 """
 function ClimaLand.surface_height(model::EnergyHydrology{FT}, Y, p) where {FT}
-    face_space =
-        ClimaLand.Domains.obtain_face_space(model.domain.space.subsurface)
-    N = ClimaCore.Spaces.nlevels(face_space)
-    surface_space = model.domain.space.surface
-    z_sfc = ClimaCore.Fields.Field(
-        ClimaCore.Fields.field_values(
-            ClimaCore.Fields.level(
-                ClimaCore.Fields.coordinate_field(face_space).z,
-                ClimaCore.Utilities.PlusHalf(N - 1),
-            ),
-        ),
-        surface_space,
-    )
-    return z_sfc
+    return model.domain.fields.z_sfc
 end
 
 function ClimaLand.get_drivers(model::EnergyHydrology)
diff --git a/src/standalone/Soil/rre.jl b/src/standalone/Soil/rre.jl
index 4f1c0d486a..e7eb0afc70 100644
--- a/src/standalone/Soil/rre.jl
+++ b/src/standalone/Soil/rre.jl
@@ -98,8 +98,9 @@ end
 
 function make_update_boundary_fluxes(model::RichardsModel)
     function update_boundary_fluxes!(p, Y, t)
-        z = ClimaCore.Fields.coordinate_field(model.domain.space.subsurface).z
-        Δz_top, Δz_bottom = ClimaLand.Domains.get_Δz(z)
+        z = model.domain.fields.z
+        Δz_top = model.domain.fields.Δz_top
+        Δz_bottom = model.domain.fields.Δz_bottom
         p.soil.top_bc .= boundary_flux(
             model.boundary_conditions.top,
             TopBoundary(),
@@ -145,7 +146,7 @@ with that value.
 """
 function ClimaLand.make_compute_imp_tendency(model::RichardsModel)
     function compute_imp_tendency!(dY, Y, p, t)
-        z = ClimaCore.Fields.coordinate_field(model.domain.space.subsurface).z
+        z = model.domain.fields.z
         top_flux_bc = p.soil.top_bc
         bottom_flux_bc = p.soil.bottom_bc
 
@@ -193,7 +194,7 @@ function ClimaLand.make_compute_exp_tendency(model::Soil.RichardsModel)
     function compute_exp_tendency!(dY, Y, p, t)
         # set dY before updating it
         dY.soil.ϑ_l .= eltype(dY.soil.ϑ_l)(0)
-        z = ClimaCore.Fields.coordinate_field(model.domain.space.subsurface).z
+        z = model.domain.fields.z
 
         horizontal_components!(
             dY,
diff --git a/test/shared_utilities/domains.jl b/test/shared_utilities/domains.jl
index 24dca93ac0..b4a6ef7fdd 100644
--- a/test/shared_utilities/domains.jl
+++ b/test/shared_utilities/domains.jl
@@ -14,7 +14,9 @@ using ClimaLand.Domains:
     coordinates,
     obtain_surface_space,
     obtain_face_space,
-    obtain_surface_domain
+    obtain_surface_domain,
+    get_Δz,
+    top_face_to_surface
 
 FT = Float32
 @testset "Clima Core Domains, FT = $FT" begin
@@ -35,6 +37,14 @@ FT = Float32
         nelements = n_elements_sphere,
         npolynomial = npoly_sphere,
     )
+    @test shell.fields.z ==
+          ClimaCore.Fields.coordinate_field(shell.space.subsurface).z
+    face_space = obtain_face_space(shell.space.subsurface)
+    z_face = ClimaCore.Fields.coordinate_field(face_space).z
+    @test shell.fields.z_sfc == top_face_to_surface(z_face, shell.space.surface)
+    Δz_top, Δz_bottom = get_Δz(shell.fields.z)
+    @test shell.fields.Δz_top == Δz_top
+    @test shell.fields.Δz_bottom == Δz_bottom
     @test shell.radius == radius
     @test shell.depth == depth
     @test shell.nelements == n_elements_sphere
@@ -85,41 +95,46 @@ FT = Float32
 
 
     # HybridBox
-    xyz_column_box = HybridBox(;
+    box = HybridBox(;
         xlim = xlim,
         ylim = ylim,
         zlim = zlim,
         nelements = nelements,
         npolynomial = 0,
     )
-    box_coords = coordinates(xyz_column_box).subsurface
+    @test box.fields.z ==
+          ClimaCore.Fields.coordinate_field(box.space.subsurface).z
+    face_space = obtain_face_space(box.space.subsurface)
+    z_face = ClimaCore.Fields.coordinate_field(face_space).z
+    @test box.fields.z_sfc == top_face_to_surface(z_face, box.space.surface)
+    Δz_top, Δz_bottom = get_Δz(box.fields.z)
+    @test box.fields.Δz_top == Δz_top
+    @test box.fields.Δz_bottom == Δz_bottom
+    box_coords = coordinates(box).subsurface
     @test eltype(box_coords) == ClimaCore.Geometry.XYZPoint{FT}
     @test typeof(box_coords) <: ClimaCore.Fields.Field
-    @test xyz_column_box.xlim == FT.(xlim)
-    @test xyz_column_box.ylim == FT.(ylim)
-    @test xyz_column_box.zlim == FT.(zlim)
-    @test xyz_column_box.nelements == nelements
-    @test xyz_column_box.npolynomial == 0
-    @test xyz_column_box.periodic == (true, true)
-    @test typeof(
-        ClimaCore.Spaces.horizontal_space(xyz_column_box.space.subsurface),
-    ) <: ClimaCore.Spaces.SpectralElementSpace2D
-    @test typeof(xyz_column_box.space.subsurface) <:
-          ClimaCore.Spaces.CenterExtrudedFiniteDifferenceSpace
-    @test typeof(xyz_column_box.space.surface) <:
+    @test box.xlim == FT.(xlim)
+    @test box.ylim == FT.(ylim)
+    @test box.zlim == FT.(zlim)
+    @test box.nelements == nelements
+    @test box.npolynomial == 0
+    @test box.periodic == (true, true)
+    @test typeof(ClimaCore.Spaces.horizontal_space(box.space.subsurface)) <:
           ClimaCore.Spaces.SpectralElementSpace2D
-    @test obtain_surface_space(xyz_column_box.space.subsurface) ==
-          xyz_column_box.space.surface
-    @test obtain_surface_domain(xyz_column_box) == Plane{FT}(
+    @test typeof(box.space.subsurface) <:
+          ClimaCore.Spaces.CenterExtrudedFiniteDifferenceSpace
+    @test typeof(box.space.surface) <: ClimaCore.Spaces.SpectralElementSpace2D
+    @test obtain_surface_space(box.space.subsurface) == box.space.surface
+    @test obtain_surface_domain(box) == Plane{FT}(
         xlim,
         ylim,
         nelements[1:2],
         (true, true),
         0,
-        (; surface = xyz_column_box.space.surface),
+        (; surface = box.space.surface),
     )
 
-    xyz_stretch_column_box = HybridBox(;
+    stretch_box = HybridBox(;
         xlim = xlim,
         ylim = ylim,
         zlim = zlim,
@@ -127,7 +142,7 @@ FT = Float32
         nelements = nelements,
         npolynomial = 0,
     )
-    box_coords_stretch = coordinates(xyz_stretch_column_box).subsurface
+    box_coords_stretch = coordinates(stretch_box).subsurface
     dz =
         Array(parent(box_coords_stretch.z))[:][2:end] .-
         Array(parent(box_coords_stretch.z))[:][1:(end - 1)]
@@ -157,6 +172,15 @@ FT = Float32
     # Column
 
     z_column = Column(; zlim = zlim, nelements = nelements[3])
+    @test z_column.fields.z ==
+          ClimaCore.Fields.coordinate_field(z_column.space.subsurface).z
+    face_space = obtain_face_space(z_column.space.subsurface)
+    z_face = ClimaCore.Fields.coordinate_field(face_space).z
+    @test z_column.fields.z_sfc ==
+          top_face_to_surface(z_face, z_column.space.surface)
+    Δz_top, Δz_bottom = get_Δz(z_column.fields.z)
+    @test z_column.fields.Δz_top == Δz_top
+    @test z_column.fields.Δz_bottom == Δz_bottom
     column_coords = coordinates(z_column).subsurface
     @test z_column.zlim == FT.(zlim)
     @test z_column.nelements[1] == nelements[3]
diff --git a/test/standalone/Soil/climate_drivers.jl b/test/standalone/Soil/climate_drivers.jl
index 2b0bc85e72..271828c6c4 100644
--- a/test/standalone/Soil/climate_drivers.jl
+++ b/test/standalone/Soil/climate_drivers.jl
@@ -115,6 +115,7 @@ for FT in (Float32, Float64)
             )
 
             Y, p, coords = initialize(model)
+            Δz_top = model.domain.fields.Δz_top
             @test propertynames(p.drivers) == (
                 :P_liq,
                 :P_snow,
@@ -229,6 +230,7 @@ for FT in (Float32, Float64)
                 ψ_sfc,
                 p.soil.ψ,
                 coords.subsurface.z,
+                Δz_top,
             )
             q_sfc = @. (q_sat * exp(g * ψ_sfc * M_w / (R * T_sfc)))
             @test ClimaLand.surface_specific_humidity(
@@ -274,6 +276,7 @@ for FT in (Float32, Float64)
                 S_l_sfc,
                 Soil.effective_saturation.(ν, Y.soil.ϑ_l, θ_r),
                 coords.subsurface.z,
+                Δz_top,
             )
             τ_a = ClimaLand.Domains.top_center_to_surface(
                 @. (ν - p.soil.θ_l - Y.soil.θ_i)^(FT(5 / 2)) / ν
diff --git a/test/standalone/Soil/soiltest.jl b/test/standalone/Soil/soiltest.jl
index 5e0e365358..32e146e5d3 100644
--- a/test/standalone/Soil/soiltest.jl
+++ b/test/standalone/Soil/soiltest.jl
@@ -648,7 +648,7 @@ for FT in (Float32, Float64)
             ),
         )
 
-        sources = (PhaseChange{FT}(Δz),)
+        sources = (PhaseChange{FT}(),)
 
         ###
         hyd_off_en_on = Soil.EnergyHydrologyParameters(