diff --git a/docs/tutorials/Usage/domain_tutorial.jl b/docs/tutorials/Usage/domain_tutorial.jl
index 51535ecc65..2f2a21b5a9 100644
--- a/docs/tutorials/Usage/domain_tutorial.jl
+++ b/docs/tutorials/Usage/domain_tutorial.jl
@@ -16,7 +16,7 @@ system of differential equations, where at each coordinate point a value of
 `ρ`, `T`, `P`, and `u⃗` are solved for at each step. The choice of domain is a question "only"
 of geometry: you may be interested in a large eddy simulation (using a box domain), or
 in a global model (where you would need a spherical shell domain
-representing the atmosphere or ocean from some depth to a given height).
+representing the atmosphere or ocean from some depth to z_sfc = 0).
 
 For land surface
 models, each variable is not defined everywhere in space. For example,
@@ -55,7 +55,7 @@ There is a single key method which take a ClimaLSM domain as an argument.
 - [` coordinates(domain)`](https://clima.github.io/ClimaLSM.jl/dev/APIs/shared_utilities/#ClimaLSM.Domains.coordinates): under the hood, this function  uses
 the function space (domain.space) to create the coordinate field. This returns the coordinates as a ClimaCore.Fields.Field object.
 Depending on the domain, the returned coordinate field will have elements of different names and types. For example,
-the SphericalShell domain has coordinates of latitude, longitude, and height, while a Plane domain has coordinates
+the SphericalShell domain has coordinates of latitude, longitude, and depth, while a Plane domain has coordinates
 of x and y, and a Point domain only has a coordinate z_sfc.
 
 
diff --git a/experiments/integrated/ozark/ozark.jl b/experiments/integrated/ozark/ozark.jl
index b9ab014957..8a79db2162 100644
--- a/experiments/integrated/ozark/ozark.jl
+++ b/experiments/integrated/ozark/ozark.jl
@@ -425,8 +425,7 @@ Plots.plot!(
 Plots.plot(plt1, plt2, layout = (2, 1))
 Plots.savefig(joinpath(savedir, "stomatal_conductance.png"))
 
-# Current resolution is 3.333 cm per layer, our nodes are at the center of the
-# layers. The second layer is ~ 5cm
+# Current resolution has the first layer at 0.1 cm, the second at 5cm.
 plt1 = Plots.plot(size = (1500, 800))
 Plots.plot!(
     plt1,
diff --git a/experiments/integrated/ozark/ozark_domain.jl b/experiments/integrated/ozark/ozark_domain.jl
index 0eb852a043..b297f2285e 100644
--- a/experiments/integrated/ozark/ozark_domain.jl
+++ b/experiments/integrated/ozark/ozark_domain.jl
@@ -1,10 +1,16 @@
 # Domain setup
 # For soil column
-nelements = 60
+nelements = 10
 zmin = FT(-2)
 zmax = FT(0)
-land_domain =
-    LSMSingleColumnDomain(; zlim = (zmin, zmax), nelements = nelements)
+dz_bottom = FT(0.5)
+dz_top = FT(0.025)
+
+land_domain = LSMSingleColumnDomain(;
+    zlim = (zmin, zmax),
+    nelements = nelements,
+    dz_tuple = (dz_bottom, dz_top),
+)
 
 # Number of stem and leaf compartments. Leaf compartments are stacked on top of stem compartments
 n_stem = Int64(1)
diff --git a/experiments/integrated/ozark/ozark_simulation.jl b/experiments/integrated/ozark/ozark_simulation.jl
index 0fa282f6e8..9214c76483 100644
--- a/experiments/integrated/ozark/ozark_simulation.jl
+++ b/experiments/integrated/ozark/ozark_simulation.jl
@@ -1,8 +1,9 @@
 t0 = FT(120 * 3600 * 24)# start mid year
 N_days = 120
 tf = t0 + FT(3600 * 24 * N_days)
-dt = FT(60)
-n = 60
+dt = FT(240)
+# Number of timesteps between saving output
+n = 15
 saveat = Array(t0:(n * dt):tf)
 timestepper = CTS.RK4()
 # Set up timestepper
diff --git a/src/shared_utilities/Domains.jl b/src/shared_utilities/Domains.jl
index 074d1b8caf..1d9943a909 100644
--- a/src/shared_utilities/Domains.jl
+++ b/src/shared_utilities/Domains.jl
@@ -63,6 +63,8 @@ struct Column{FT, S} <: AbstractDomain{FT}
     zlim::Tuple{FT, FT}
     "Number of elements used to discretize the interval"
     nelements::Tuple{Int}
+    "Tuple for mesh stretching specifying *target* (dz_bottom, dz_top) (m). If nothing, no stretching is applied."
+    dz_tuple::Union{Tuple{FT, FT}, Nothing}
     "Boundary face identifiers"
     boundary_tags::Tuple{Symbol, Symbol}
     "The associated ClimaCore Space"
@@ -72,12 +74,28 @@ end
 """
     Column(;
         zlim::Tuple{FT, FT},
-        nelements::Int) where {FT}
-Outer constructor for the `Column` type.
+        nelements::Int,
+        dz_tuple::Union{Tuple{FT, FT}, Nothing} = nothing) where {FT}
+
+Outer constructor for the `Column` type. 
+
+Using `ClimaCore` tools, the coordinate
+mesh can be stretched such that the top of the domain has finer resolution
+than the bottom of the domain. In order to activate this, a tuple with the
+target dz_bottom and dz_top should be passed via keyword argument. The default is
+uniform spacing. Please note that in order to use this feature, ClimaCore requires that
+the elements of zlim be <=0. Additionally, the dz_tuple you supply may not be compatible
+with the domain boundaries in some cases, in which case you may need to choose
+different values.
+
 The `boundary_tags` field values are used to label the boundary faces 
 at the top and bottom of the domain.
 """
-function Column(; zlim::Tuple{FT, FT}, nelements::Int) where {FT}
+function Column(;
+    zlim::Tuple{FT, FT},
+    nelements::Int,
+    dz_tuple::Union{Tuple{FT, FT}, Nothing} = nothing,
+) where {FT}
     @assert zlim[1] < zlim[2]
     boundary_tags = (:bottom, :top)
     column = ClimaCore.Domains.IntervalDomain(
@@ -85,11 +103,26 @@ function Column(; zlim::Tuple{FT, FT}, nelements::Int) where {FT}
         ClimaCore.Geometry.ZPoint{FT}(zlim[2]);
         boundary_tags = boundary_tags,
     )
-    mesh = ClimaCore.Meshes.IntervalMesh(column; nelems = nelements)
+    if dz_tuple isa Nothing
+        mesh = ClimaCore.Meshes.IntervalMesh(column; nelems = nelements)
+    else
+        @assert zlim[2] <= 0
+        mesh = ClimaCore.Meshes.IntervalMesh(
+            column,
+            ClimaCore.Meshes.GeneralizedExponentialStretching{FT}(
+                dz_tuple[1],
+                dz_tuple[2],
+            );
+            nelems = nelements,
+            reverse_mode = true,
+        )
+    end
+
     center_space = ClimaCore.Spaces.CenterFiniteDifferenceSpace(mesh)
     return Column{FT, typeof(center_space)}(
         zlim,
         (nelements,),
+        dz_tuple,
         boundary_tags,
         center_space,
     )
@@ -174,6 +207,7 @@ end
         xlim::Tuple{FT, FT}
         ylim::Tuple{FT, FT}
         zlim::Tuple{FT, FT}
+        dz_tuple::Union{Tuple{FT, FT}, Nothing}
         nelements::Tuple{Int, Int, Int}
         npolynomial::Int
         periodic::Tuple{Bool, Bool}
@@ -194,6 +228,8 @@ struct HybridBox{FT, S} <: AbstractDomain{FT}
     ylim::Tuple{FT, FT}
     "Domain interval limits along z axis, in meters"
     zlim::Tuple{FT, FT}
+    "Tuple for mesh stretching specifying *target* (dz_bottom, dz_top) (m). If nothing, no stretching is applied."
+    dz_tuple::Union{Tuple{FT, FT}, Nothing}
     "Number of elements to discretize interval, (nx, ny,nz)"
     nelements::Tuple{Int, Int, Int}
     " Polynomial order for the horizontal directions"
@@ -211,6 +247,7 @@ end
         zlim::Tuple{FT, FT},
         nelements::Tuple{Int, Int, Int},
         npolynomial::Int,
+        dz_tuple::Union{Tuple{FT, FT}, Nothing} = nothing,
         periodic = (true, true),
     ) where {FT}
 Constructs the `HybridBox` domain
@@ -222,6 +259,15 @@ to the x-axis, the second corresponding to the y-axis, and the third
 corresponding to the z-axis. The domain is periodic at the (xy) boundaries,
 and the function space is of polynomial order `npolynomial` in the
 horizontal directions.
+
+Using `ClimaCore` tools, the coordinate
+mesh can be stretched such that the top of the domain has finer resolution
+than the bottom of the domain. In order to activate this, a tuple with the
+target dz_bottom and dz_top should be passed via keyword argument. The default is
+uniform spacing. Please note that in order to use this feature, ClimaCore requires that
+the elements of zlim be <=0. Additionally, the dz_tuple you supply may not be compatible
+with the domain boundaries in some cases, in which case you may need to choose
+different values.
 """
 function HybridBox(;
     xlim::Tuple{FT, FT},
@@ -229,6 +275,7 @@ function HybridBox(;
     zlim::Tuple{FT, FT},
     nelements::Tuple{Int, Int, Int},
     npolynomial::Int,
+    dz_tuple::Union{Tuple{FT, FT}, Nothing} = nothing,
     periodic = (true, true),
 ) where {FT}
     @assert xlim[1] < xlim[2]
@@ -240,7 +287,21 @@ function HybridBox(;
         ClimaCore.Geometry.ZPoint(zlim[2]);
         boundary_tags = (:bottom, :top),
     )
-    vertmesh = ClimaCore.Meshes.IntervalMesh(vertdomain, nelems = nelements[3])
+    if dz_tuple isa Nothing
+        vertmesh =
+            ClimaCore.Meshes.IntervalMesh(vertdomain; nelems = nelements[3])
+    else
+        @assert zlim[2] <= 0
+        vertmesh = ClimaCore.Meshes.IntervalMesh(
+            vertdomain,
+            ClimaCore.Meshes.GeneralizedExponentialStretching{FT}(
+                dz_tuple[1],
+                dz_tuple[2],
+            );
+            nelems = nelements[3],
+            reverse_mode = true,
+        )
+    end
     vert_center_space = ClimaCore.Spaces.CenterFiniteDifferenceSpace(vertmesh)
 
     horzdomain = Plane(;
@@ -260,6 +321,7 @@ function HybridBox(;
         xlim,
         ylim,
         zlim,
+        dz_tuple,
         nelements,
         npolynomial,
         periodic,
@@ -270,7 +332,8 @@ end
 """
     struct SphericalShell{FT} <: AbstractDomain{FT}
         radius::FT
-        height::FT
+        depth::FT
+        dz_tuple::Union{Tuple{FT, FT}, Nothing}
         nelements::Tuple{Int, Int}
         npolynomial::Int
     end
@@ -285,7 +348,9 @@ struct SphericalShell{FT, S} <: AbstractDomain{FT}
     "The radius of the shell"
     radius::FT
     "The radial extent of the shell"
-    height::FT
+    depth::FT
+    "Tuple for mesh stretching specifying *target* (dz_bottom, dz_top) (m). If nothing, no stretching is applied."
+    dz_tuple::Union{Tuple{FT, FT}, Nothing}
     "The number of elements to be used in the non-radial and radial directions"
     nelements::Tuple{Int, Int}
     "The polynomial order to be used in the non-radial directions"
@@ -297,31 +362,49 @@ end
 """
     SphericalShell(;
         radius::FT,
-        height::FT,
+        depth::FT,
         nelements::Tuple{Int, Int},
         npolynomial::Int,
+        dz_tuple::Union{Tuple{FT, FT}, Nothing} = nothing,
     ) where {FT}
 Outer constructor for the `SphericalShell` domain, using keyword arguments.
+
+Using `ClimaCore` tools, the coordinate
+mesh can be stretched such that the top of the domain has finer resolution
+than the bottom of the domain. In order to activate this, a tuple with the
+target dz_bottom and dz_top should be passed via keyword argument. The default is
+uniform spacing. Please note that the dz_tuple you supply may not be compatible
+with the depth/nelements chosen, in which case you may need to choose
+different values.
 """
 function SphericalShell(;
     radius::FT,
-    height::FT,
+    depth::FT,
     nelements::Tuple{Int, Int},
     npolynomial::Int,
+    dz_tuple::Union{Tuple{FT, FT}, Nothing} = nothing,
 ) where {FT}
     @assert 0 < radius
-    @assert 0 < height
+    @assert 0 < depth
     vertdomain = ClimaCore.Domains.IntervalDomain(
-        ClimaCore.Geometry.ZPoint(FT(0)),
-        ClimaCore.Geometry.ZPoint(FT(height));
+        ClimaCore.Geometry.ZPoint(FT(-depth)),
+        ClimaCore.Geometry.ZPoint(FT(0));
         boundary_tags = (:bottom, :top),
     )
-
-    vertmesh = ClimaCore.Meshes.IntervalMesh(
-        vertdomain,
-        ClimaCore.Meshes.Uniform(),
-        nelems = nelements[2],
-    )
+    if dz_tuple isa Nothing
+        vertmesh =
+            ClimaCore.Meshes.IntervalMesh(vertdomain; nelems = nelements[2])
+    else
+        vertmesh = ClimaCore.Meshes.IntervalMesh(
+            vertdomain,
+            ClimaCore.Meshes.GeneralizedExponentialStretching{FT}(
+                dz_tuple[1],
+                dz_tuple[2],
+            );
+            nelems = nelements[2],
+            reverse_mode = true,
+        )
+    end
     vert_center_space = ClimaCore.Spaces.CenterFiniteDifferenceSpace(vertmesh)
 
     horzdomain = ClimaCore.Domains.SphereDomain(radius)
@@ -336,7 +419,8 @@ function SphericalShell(;
     )
     return SphericalShell{FT, typeof(hv_center_space)}(
         radius,
-        height,
+        depth,
+        dz_tuple,
         nelements,
         npolynomial,
         hv_center_space,
@@ -424,12 +508,18 @@ end
     LSMSingleColumnDomain(;
         zlim::Tuple{FT, FT},
         nelements::Int,
+        dz_tuple::Union{Tuple{FT, FT}, Nothing} = nothing,
     ) where {FT}
 A constructor for the LSMSingleColumnDomain.
 """
-function LSMSingleColumnDomain(; zlim::Tuple{FT, FT}, nelements::Int) where {FT}
+function LSMSingleColumnDomain(;
+    zlim::Tuple{FT, FT},
+    nelements::Int,
+    dz_tuple::Union{Tuple{FT, FT}, Nothing} = nothing,
+) where {FT}
     @assert zlim[1] < zlim[2]
-    subsurface_domain = Column(; zlim = FT.(zlim), nelements = nelements)
+    subsurface_domain =
+        Column(; zlim = FT.(zlim), nelements = nelements, dz_tuple = dz_tuple)
     surface_space = obtain_surface_space(subsurface_domain.space)
     surface_domain = Point{FT, typeof(surface_space)}(zlim[2], surface_space)
     return LSMSingleColumnDomain{
@@ -469,7 +559,8 @@ end
         zlim::Tuple{FT, FT}
         nelements::Tuple{Int, Int, Int}
         npolynomial::Int
-        periodic::Tuple{Bool, Bool}
+        periodic::Tuple{Bool, Bool} = (true, true),
+        dz_tuple::Union{Tuple{FT,FT}, Nothing} = nothing,
     ) where {FT}
 A constructor for the LSMMultiColumnDomain.
 """
@@ -479,7 +570,8 @@ function LSMMultiColumnDomain(;
     zlim::Tuple{FT, FT},
     nelements::Tuple{Int, Int, Int},
     npolynomial::Int,
-    periodic::Tuple{Bool, Bool},
+    periodic::Tuple{Bool, Bool} = (true, true),
+    dz_tuple::Union{Tuple{FT, FT}, Nothing} = nothing,
 ) where {FT}
     @assert xlim[1] < xlim[2]
     @assert ylim[1] < ylim[2]
@@ -489,6 +581,7 @@ function LSMMultiColumnDomain(;
         xlim = xlim,
         ylim = ylim,
         zlim = zlim,
+        dz_tuple = dz_tuple,
         nelements = nelements,
         npolynomial = npolynomial,
         periodic = periodic,
@@ -536,23 +629,26 @@ end
 """
     LSMSphericalShellDomain(;
         radius::FT,
-        height::FT,
+        depth::FT,
         nelements::Tuple{Int, Int},
         npolynomial::Int,
+        dz_tuple::Union{Tuple{FT, FT}, Nothing} = nothing,
     ) where {FT}
 A constructor for the LSMSphericalShellDomain.
 """
 function LSMSphericalShellDomain(;
     radius::FT,
-    height::FT,
+    depth::FT,
     nelements::Tuple{Int, Int},
     npolynomial::Int,
+    dz_tuple::Union{Tuple{FT, FT}, Nothing} = nothing,
 ) where {FT}
     @assert radius > FT(0)
-    @assert height > FT(0)
+    @assert depth > FT(0)
     subsurface = SphericalShell(;
         radius = radius,
-        height = height,
+        depth = depth,
+        dz_tuple = dz_tuple,
         nelements = nelements,
         npolynomial = npolynomial,
     )
diff --git a/test/shared_utilities/domains.jl b/test/shared_utilities/domains.jl
index f2b6bc8d26..c16e30c4a7 100644
--- a/test/shared_utilities/domains.jl
+++ b/test/shared_utilities/domains.jl
@@ -20,20 +20,20 @@ TestFloatTypes = (Float32, Float64)
 
 @testset "Clima Core Domains" begin
     for FT in TestFloatTypes
-        zmin = FT(1.0)
-        zmax = FT(2.0)
+        zmin = FT(-1.0)
+        zmax = FT(0.0)
         xlim = FT.((0.0, 10.0))
         ylim = FT.((0.0, 1.0))
         zlim = FT.((zmin, zmax))
-        nelements = (1, 1, 5)
+        nelements = (1, 1, 10)
         shell = SphericalShell(;
             radius = FT(100.0),
-            height = FT(30.0),
+            depth = FT(30.0),
             nelements = (6, 20),
             npolynomial = 3,
         )
         @test shell.radius == FT(100)
-        @test shell.height == FT(30)
+        @test shell.depth == FT(30)
         @test shell.nelements == (6, 20)
         @test shell.npolynomial == 3
         shell_coords = coordinates(shell)
@@ -44,6 +44,20 @@ TestFloatTypes = (Float32, Float64)
         @test typeof(shell.space) <:
               ClimaCore.Spaces.CenterExtrudedFiniteDifferenceSpace
 
+        shell_stretch = SphericalShell(;
+            radius = FT(100.0),
+            depth = FT(1.0),
+            dz_tuple = FT.((0.3, 0.03)),
+            nelements = (6, 10),
+            npolynomial = 3,
+        )
+        shell_coords_stretch = coordinates(shell_stretch)
+        dz =
+            parent(shell_coords_stretch.z)[:, 1, 4, 1, 216][2:end] .-
+            parent(shell_coords_stretch.z)[:, 1, 4, 1, 216][1:(end - 1)]
+        @test abs(dz[1] - 0.3) < 1e-1
+        @test abs(dz[end] - 0.03) < 1e-2
+
 
         shell_surface = SphericalSurface(;
             radius = FT(100.0),
@@ -60,6 +74,13 @@ TestFloatTypes = (Float32, Float64)
         @test typeof(shell_surface.space) <:
               ClimaCore.Spaces.SpectralElementSpace2D
 
+        xy_plane = Plane(;
+            xlim = xlim,
+            ylim = ylim,
+            nelements = nelements[1:2],
+            periodic = (true, true),
+            npolynomial = 0,
+        )
 
         xyz_column_box = HybridBox(;
             xlim = xlim,
@@ -82,6 +103,20 @@ TestFloatTypes = (Float32, Float64)
         @test typeof(xyz_column_box.space) <:
               ClimaCore.Spaces.CenterExtrudedFiniteDifferenceSpace
 
+        xyz_stretch_column_box = HybridBox(;
+            xlim = xlim,
+            ylim = ylim,
+            zlim = zlim,
+            dz_tuple = FT.((0.3, 0.03)),
+            nelements = nelements,
+            npolynomial = 0,
+        )
+        box_coords_stretch = coordinates(xyz_stretch_column_box)
+        dz =
+            parent(box_coords_stretch.z)[:][2:end] .-
+            parent(box_coords_stretch.z)[:][1:(end - 1)]
+        @test abs(dz[1] - 0.3) < 1e-1
+        @test abs(dz[end] - 0.03) < 1e-2
         xy_plane = Plane(;
             xlim = xlim,
             ylim = ylim,
@@ -107,7 +142,21 @@ TestFloatTypes = (Float32, Float64)
         @test typeof(column_coords) <: ClimaCore.Fields.Field
         @test typeof(z_column.space) <:
               ClimaCore.Spaces.CenterFiniteDifferenceSpace
+        @test any(
+            parent(column_coords)[:][2:end] .-
+            parent(column_coords)[:][1:(end - 1)] .≈
+            (zmax - zmin) / nelements[3],
+        )
 
+        z_column_stretch =
+            Column(; zlim = zlim, nelements = 10, dz_tuple = FT.((0.3, 0.03)))
+        column_coords = coordinates(z_column_stretch)
+        @test z_column_stretch.zlim == FT.(zlim)
+        dz =
+            parent(column_coords)[:][2:end] .-
+            parent(column_coords)[:][1:(end - 1)]
+        @test abs(dz[1] - 0.3) < 1e-1
+        @test abs(dz[end] - 0.03) < 1e-2
     end
 end
 
@@ -130,11 +179,15 @@ end
 
 @testset "LSMSingleColumnDomain" begin
     for FT in TestFloatTypes
-        zmin = FT(1.0)
-        zmax = FT(2.0)
+        zmin = FT(-1.0)
+        zmax = FT(0.0)
         zlim = FT.((zmin, zmax))
         nelements = 5
-        domain = LSMSingleColumnDomain(; zlim = zlim, nelements = nelements)
+        domain = LSMSingleColumnDomain(;
+            zlim = zlim,
+            nelements = nelements,
+            dz_tuple = FT.((0.3, 0.03)),
+        )
         point = domain.surface
         column = domain.subsurface
 
@@ -149,8 +202,8 @@ end
 
 @testset "LSMMultiColumnDomain" begin
     for FT in TestFloatTypes
-        zmin = FT(1.0)
-        zmax = FT(2.0)
+        zmin = FT(-1.0)
+        zmax = FT(0.0)
         zlim = FT.((zmin, zmax))
         xlim = FT.((0.0, 10.0))
         ylim = FT.((0.0, 1.0))
@@ -164,6 +217,7 @@ end
             nelements = nelements,
             periodic = (true, true),
             npolynomial = npolynomial,
+            dz_tuple = FT.((0.3, 0.03)),
         )
         plane = domain.surface
         box = domain.subsurface
@@ -182,9 +236,10 @@ end
     for FT in TestFloatTypes
         domain = LSMSphericalShellDomain(;
             radius = FT(100.0),
-            height = FT(30.0),
+            depth = FT(30.0),
             nelements = (6, 20),
             npolynomial = 3,
+            dz_tuple = FT.((5.0, 0.3)),
         )
         surf = domain.surface
         shell = domain.subsurface
diff --git a/test/shared_utilities/utilities.jl b/test/shared_utilities/utilities.jl
index a49d2f2acc..0fecfa00a5 100644
--- a/test/shared_utilities/utilities.jl
+++ b/test/shared_utilities/utilities.jl
@@ -199,7 +199,7 @@ end
         )
         domain2 = ClimaLSM.Domains.SphericalShell(;
             radius = FT(2),
-            height = FT(1.0),
+            depth = FT(1.0),
             nelements = (10, 5),
             npolynomial = 3,
         )
diff --git a/test/standalone/Bucket/albedo_types.jl b/test/standalone/Bucket/albedo_types.jl
index 7c418e961b..5704d85d6e 100644
--- a/test/standalone/Bucket/albedo_types.jl
+++ b/test/standalone/Bucket/albedo_types.jl
@@ -46,7 +46,7 @@ function create_domain_2d(FT)
     np = 2
     return LSMSphericalShellDomain(;
         radius = rad,
-        height = h,
+        depth = h,
         nelements = ne,
         npolynomial = np,
     )
@@ -277,7 +277,7 @@ end
         LSMSingleColumnDomain(; zlim = (-100.0, 0.0), nelements = 10),
         LSMSphericalShellDomain(;
             radius = FT(100.0),
-            height = FT(3.5),
+            depth = FT(3.5),
             nelements = (2, 10),
             npolynomial = 2,
         ),
@@ -418,7 +418,7 @@ end
         LSMSingleColumnDomain(; zlim = (-100.0, 0.0), nelements = 10),
         LSMSphericalShellDomain(;
             radius = FT(100.0),
-            height = FT(3.5),
+            depth = FT(3.5),
             nelements = (2, 10),
             npolynomial = 2,
         ),
diff --git a/test/standalone/Bucket/snow_bucket_tests.jl b/test/standalone/Bucket/snow_bucket_tests.jl
index 77523e6842..6179e837e9 100644
--- a/test/standalone/Bucket/snow_bucket_tests.jl
+++ b/test/standalone/Bucket/snow_bucket_tests.jl
@@ -50,7 +50,7 @@ bucket_domains = [
     ),
     LSMSphericalShellDomain(;
         radius = 100.0,
-        height = 3.5,
+        depth = 3.5,
         nelements = (1, 10),
         npolynomial = 1,
     ),
diff --git a/test/standalone/Bucket/soil_bucket_tests.jl b/test/standalone/Bucket/soil_bucket_tests.jl
index 95dc225547..88b566e5b9 100644
--- a/test/standalone/Bucket/soil_bucket_tests.jl
+++ b/test/standalone/Bucket/soil_bucket_tests.jl
@@ -48,7 +48,7 @@ bucket_domains = [
     ),
     LSMSphericalShellDomain(;
         radius = 100.0,
-        height = 3.5,
+        depth = 3.5,
         nelements = (1, 10),
         npolynomial = 1,
     ),
diff --git a/test/standalone/Soil/soil_bc.jl b/test/standalone/Soil/soil_bc.jl
index 7036cf3381..2d2c898845 100644
--- a/test/standalone/Soil/soil_bc.jl
+++ b/test/standalone/Soil/soil_bc.jl
@@ -17,7 +17,7 @@ FT = Float64
 
     domain = SphericalShell(;
         radius = FT(1.0),
-        height = FT(1.0),
+        depth = FT(1.0),
         nelements = (1, 2),
         npolynomial = 3,
     )
diff --git a/test/standalone/Soil/soil_test_3d.jl b/test/standalone/Soil/soil_test_3d.jl
index 3036872307..b67e59915d 100644
--- a/test/standalone/Soil/soil_test_3d.jl
+++ b/test/standalone/Soil/soil_test_3d.jl
@@ -356,7 +356,7 @@ end
 
     soil_domain = SphericalShell(;
         radius = FT(100.0),
-        height = FT(1.0),
+        depth = FT(1.0),
         nelements = (1, 30),
         npolynomial = 3,
     )
@@ -404,13 +404,13 @@ end
 
     ## vertical change should be zero except at the boundary, where it should be ±30.
     @test (mean(unique(parent(ClimaCore.level(dY.soil.ϑ_l, 1)))) - 30.0) / ν <
-          1e-11
+          1e-10
     @test (mean(unique(parent(ClimaCore.level(dY.soil.ϑ_l, 30)))) + 30.0) / ν <
-          1e-11
+          1e-10
     @test mean([
         maximum(abs.(unique(parent(ClimaCore.level(dY.soil.ϑ_l, k))))) for
         k in 2:29
-    ]) / ν < 1e-11
+    ]) / ν < 1e-10
 
     exp_tendency!(dY, Y, p, t0)
     ClimaLSM.dss!(dY, p, t0)
@@ -419,7 +419,7 @@ end
     @test mean([
         maximum(abs.(unique(parent(ClimaCore.level(dY.soil.ϑ_l, k))))) for
         k in 1:30
-    ]) / ν < 1e-11
+    ]) / ν < 1e-10
 
 end
 
@@ -438,7 +438,7 @@ end
     )
     shell = ClimaLSM.Domains.SphericalShell(;
         radius = FT(1),
-        height = FT(1),
+        depth = FT(1),
         nelements = (2, 3),
         npolynomial = 1,
     )