diff --git a/NEWS.md b/NEWS.md
index 66f61b030c..f735ce7231 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -3,6 +3,7 @@ ClimaCore.jl Release Notes
 
 main
 -------
+- Added hyperbolic tangent stretching. PR [#1930](https://github.com/CliMA/ClimaCore.jl/pull/1930).
 
 v0.14.11
 -------
diff --git a/docs/src/api.md b/docs/src/api.md
index e5d62ce010..d4f9df8c3e 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -105,6 +105,7 @@ Meshes.NormalizedBilinearMap
 Meshes.Uniform
 Meshes.ExponentialStretching
 Meshes.GeneralizedExponentialStretching
+Meshes.HyperbolicTangentStretching
 ```
 
 ### Mesh utilities
diff --git a/src/Meshes/interval.jl b/src/Meshes/interval.jl
index d0d799b33c..7a2ea8d755 100644
--- a/src/Meshes/interval.jl
+++ b/src/Meshes/interval.jl
@@ -16,6 +16,7 @@ Constuct a 1D mesh on `domain` with `nelems` elements, using `stretching`. Possi
 - [`Uniform()`](@ref)
 - [`ExponentialStretching(H)`](@ref)
 - [`GeneralizedExponentialStretching(dz_bottom, dz_top)`](@ref)
+- [`HyperbolicTangentStretching(dz_bottom)`](@ref)
 """
 struct IntervalMesh{I <: IntervalDomain, V <: AbstractVector} <: AbstractMesh1D
     domain::I
@@ -165,6 +166,8 @@ model configurations).
 Then, the user can define a stretched mesh via
 
     ClimaCore.Meshes.IntervalMesh(interval_domain, ExponentialStretching(H); nelems::Int, reverse_mode = false)
+
+`faces` contain reference z without any warping.
 """
 struct ExponentialStretching{FT} <: StretchingRule
     H::FT
@@ -213,6 +216,8 @@ For land configurations, use `reverse_mode` = `true` (default value `false`).
 Then, the user can define a generalized stretched mesh via
 
     ClimaCore.Meshes.IntervalMesh(interval_domain, GeneralizedExponentialStretching(dz_bottom, dz_top); nelems::Int, reverse_mode = false)
+
+`faces` contain reference z without any warping.
 """
 struct GeneralizedExponentialStretching{FT} <: StretchingRule
     dz_bottom::FT
@@ -241,9 +246,9 @@ function IntervalMesh(
         throw(ArgumentError("dz_top must be ≤ dz_bottom"))
     end
 
-    # bottom coord height value, always min, for both atmos and land, since z-axis does not change
+    # bottom coord height value is always min and top coord height value is always max
+    # since the vertical coordinate is positive upward
     z_bottom = Geometry.component(domain.coord_min, 1)
-    # top coord height value, always max, for both atmos and land, since z-axis does not change
     z_top = Geometry.component(domain.coord_max, 1)
     # but in case of reverse_mode, we temporarily swap them together with dz_bottom and dz_top
     # so that the following root solve algorithm does not need to change
@@ -256,7 +261,7 @@ function IntervalMesh(
     # define the inverse σ⁻¹ exponential stretching function
     exp_stretch(ζ, h) = ζ == 1 ? ζ : -h * log(1 - (1 - exp(-1 / h)) * ζ)
 
-    # nondimensional vertical coordinate (]0.0, 1.0])
+    # nondimensional vertical coordinate ([0.0, 1.0])
     ζ_n = LinRange(one(FT_solve), nelems, nelems) / nelems
 
     # find bottom height variation
@@ -292,7 +297,7 @@ function IntervalMesh(
         find_top,
         RootSolvers.SecantMethod(guess₋, guess₊),
         RootSolvers.CompactSolution(),
-        RootSolvers.ResidualTolerance(FT_solve(1e-3)),
+        RootSolvers.ResidualTolerance(FT_solve(tol)),
     )
     if h_top_sol.converged !== true
         error(
@@ -305,7 +310,7 @@ function IntervalMesh(
     h =
         h_bottom .+
         (ζ_n .- ζ_n[1]) * (h_top - h_bottom) / (ζ_n[end - 1] - ζ_n[1])
-    faces = (z_bottom + (z_top - z_bottom)) * exp_stretch.(ζ_n, h)
+    faces = z_bottom .+ (z_top - z_bottom) * exp_stretch.(ζ_n, h)
 
     # add the bottom level
     faces = FT_solve[z_bottom; faces...]
@@ -318,6 +323,94 @@ function IntervalMesh(
     IntervalMesh(domain, CT.(faces))
 end
 
+"""
+    HyperbolicTangentStretching(dz_surface::FT)
+
+Apply a hyperbolic tangent stretching to the domain when constructing elements.
+`dz_surface` is the target element grid spacing at the surface. In typical atmosphere
+configuration, it is the grid spacing at the bottom of the
+vertical column domain (m). On the other hand, for typical land configurations,
+it is the grid spacing at the top of the vertical column domain.
+
+For an interval ``[z_0,z_1]``, this makes the elements uniformally spaced in
+``\\zeta``, where
+```math
+\\eta = 1 - \\frac{tanh[\\gamma(1-\\zeta)]}{tanh(\\gamma)},
+```
+where ``\\eta = \\frac{z - z_0}{z_1-z_0}``. The stretching parameter ``\\gamma``
+is chosen to achieve a given resolution `dz_surface` at the surface. 
+
+Then, the user can define a stretched mesh via
+
+    ClimaCore.Meshes.IntervalMesh(interval_domain, HyperbolicTangentStretching(dz_surface); nelems::Int, reverse_mode)
+
+`reverse_mode` is default to false for atmosphere configurations. For land configurations, 
+use `reverse_mode` = `true`.
+
+`faces` contain reference z without any warping.
+"""
+struct HyperbolicTangentStretching{FT} <: StretchingRule
+    dz_surface::FT
+end
+
+function IntervalMesh(
+    domain::IntervalDomain{CT},
+    stretch::HyperbolicTangentStretching{FT};
+    nelems::Int,
+    FT_solve = Float64,
+    tol::Union{FT, Nothing} = nothing,
+    reverse_mode::Bool = false,
+) where {CT <: Geometry.Abstract1DPoint{FT}} where {FT}
+    if nelems ≤ 1
+        throw(ArgumentError("`nelems` must be ≥ 2"))
+    end
+
+    dz_surface = FT_solve(stretch.dz_surface)
+    tol === nothing && (tol = dz_surface * FT_solve(1e-6))
+
+    # bottom coord height value is always min and top coord height value is always max
+    # since the vertical coordinate is positive upward
+    z_bottom = Geometry.component(domain.coord_min, 1)
+    z_top = Geometry.component(domain.coord_max, 1)
+    # but in case of reverse_mode, we temporarily swap them
+    # so that the following root solve algorithm does not need to change
+    if reverse_mode
+        z_bottom, z_top = Geometry.component(domain.coord_max, 1),
+        -Geometry.component(domain.coord_min, 1)
+    end
+
+    # define the hyperbolic tangent stretching function
+    tanh_stretch(ζ, γ) = 1 - tanh(γ * (1 - ζ)) / tanh(γ)
+
+    # nondimensional vertical coordinate ([0.0, 1.0])
+    ζ_n = LinRange(one(FT_solve), nelems, nelems) / nelems
+
+    # find the stretching parameter given the grid spacing at the surface
+    find_surface(γ) = dz_surface - z_top * tanh_stretch(ζ_n[1], γ)
+    γ_sol = RootSolvers.find_zero(
+        find_surface,
+        RootSolvers.NewtonsMethodAD(FT_solve(1.0)),
+        RootSolvers.CompactSolution(),
+        RootSolvers.ResidualTolerance(FT_solve(tol)),
+    )
+    if γ_sol.converged !== true
+        error(
+            "gamma root failed to converge for dz_surface: $dz_surface on domain ($z_bottom, $z_top)",
+        )
+    end
+
+    faces = z_bottom .+ (z_top - z_bottom) * tanh_stretch.(ζ_n, γ_sol.root)
+
+    # add the bottom level
+    faces = FT_solve[z_bottom; faces...]
+    if reverse_mode
+        reverse!(faces)
+        faces = map(f -> eltype(faces)(-f), faces)
+        faces[end] = faces[end] == -z_bottom ? z_bottom : faces[1]
+    end
+    monotonic_check(faces)
+    IntervalMesh(domain, CT.(faces))
+end
 
 """
     truncate_mesh(
diff --git a/test/Meshes/interval.jl b/test/Meshes/interval.jl
index 1714727c94..f365b9e491 100644
--- a/test/Meshes/interval.jl
+++ b/test/Meshes/interval.jl
@@ -214,7 +214,7 @@ end
 end
 
 @testset "IntervalMesh GeneralizedExponentialStretching" begin
-    # use normalized GCM profile heights (7.5km)
+    # use normalized GCM profile heights (45km)
     @test_throws ArgumentError unit_intervalmesh(
         stretching = Meshes.GeneralizedExponentialStretching(
             0.02 / 45.0,
@@ -250,7 +250,7 @@ end
 
 
 @testset "IntervalMesh GeneralizedExponentialStretching reverse" begin
-    # use normalized GCM profile heights (7.5km)
+    # use normalized GCM profile heights (45km)
     @test_throws ArgumentError unit_intervalmesh(
         stretching = Meshes.GeneralizedExponentialStretching(
             7.0 / 45.0,
@@ -288,6 +288,55 @@ end
     @test fₑ - fₑ₋₁ ≈ 0.02 / 45.0 rtol = 1e-2
 end
 
+@testset "IntervalMesh HyperbolicTangentStretching" begin
+    # use normalized GCM profile heights (75km)
+    @test_throws ArgumentError unit_intervalmesh(
+        stretching = Meshes.HyperbolicTangentStretching(0.03 / 75.0),
+        nelems = 0,
+    )
+    @test_throws ArgumentError unit_intervalmesh(
+        stretching = Meshes.HyperbolicTangentStretching(0.03 / 75.0),
+        nelems = 1,
+    )
+    # test a gcm like configuration
+    nelems = 75
+    dom, mesh = unit_intervalmesh(
+        stretching = Meshes.HyperbolicTangentStretching(0.03 / 75.0),
+        nelems = nelems, # 76 face levels
+    )
+    # test the bottom and top of the mesh coordinates are correct
+    @test Meshes.coordinates(mesh, 1, 1) == Geometry.ZPoint(0.0)
+    @test Meshes.coordinates(mesh, 1, nelems + 1) ≈ Geometry.ZPoint(1.0)
+    # test the interval of the mesh coordinates at the surface is the same as dz_surface
+    @test Meshes.coordinates(mesh, 1, 2) ≈ Geometry.ZPoint(0.03 / 75.0)
+end
+
+@testset "IntervalMesh HyperbolicTangentStretching reverse" begin
+    # use normalized GCM profile heights (75km)
+    @test_throws ArgumentError unit_intervalmesh(
+        stretching = Meshes.HyperbolicTangentStretching(0.03 / 75.0),
+        nelems = 0,
+        reverse_mode = true,
+    )
+    @test_throws ArgumentError unit_intervalmesh(
+        stretching = Meshes.HyperbolicTangentStretching(0.03 / 75.0),
+        nelems = 1,
+        reverse_mode = true,
+    )
+    # test a gcm like configuration, for land
+    nelems = 75
+    dom, mesh = unit_intervalmesh(
+        stretching = Meshes.HyperbolicTangentStretching(0.03 / 75.0),
+        nelems = nelems, # 76 face levels
+        reverse_mode = true,
+    )
+    # test the bottom and top of the mesh coordinates are correct
+    @test Meshes.coordinates(mesh, 1, 1) == Geometry.ZPoint(-1.0)
+    @test Meshes.coordinates(mesh, 1, nelems + 1) ≈ Geometry.ZPoint(0.0)
+    # test the interval of the mesh coordinates at the surface is the same as dz_surface
+    @test Meshes.coordinates(mesh, 1, nelems) ≈ Geometry.ZPoint(-0.03 / 75.0)
+end
+
 @testset "Truncated IntervalMesh" begin
     FT = Float64
     nz = 55