diff --git a/Artifacts.toml b/Artifacts.toml
index 6ea887b27ab..d07655e70e4 100644
--- a/Artifacts.toml
+++ b/Artifacts.toml
@@ -25,3 +25,14 @@ git-tree-sha1 = "6693f8451b447543952b3d557d733c7e5c670ed3"
     [[ozone_concentrations_lowres.download]]
     sha256 = "908384b32e90c733c2173654afbb974d238a9a516ad8965039077c6ddb29b5f0"
     url = "https://caltech.box.com/shared/static/6c16wt7i1htxq6tadpt1ng2uftjgdtpr.gz"
+
+[earth_orography_30arcseconds]
+git-tree-sha1 = "03dd08fcbf363ed055a176dd7adefb60ff1c3493"
+
+[earth_orography_60arcseconds]
+git-tree-sha1 = "fe19d8dbe7a18ff39588e1f718014b0479d9c0f7"
+lazy = true
+
+    [[earth_orography_60arcseconds.download]]
+    sha256 = "eca66c0701d1c2b9e271742314915ffbf4a0fae92709df611c323f38e019966e"
+    url = "https://caltech.box.com/shared/static/4asrxcgl6xsgenfcug9p0wkkyhtqilgk.gz"
diff --git a/NEWS.md b/NEWS.md
index f50e4d0374d..6a6053febd6 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -4,6 +4,18 @@ ClimaAtmos.jl Release Notes
 Main
 -------
 
+### Features
+
+### ETOPO2022 60arc-second topography dataset. 
+
+- Update artifacts to use 60arc-second ETOPO2022 ice-surface topography
+  dataset. Update surface smoothing functions to rely only on spectral 
+  Laplacian operations. Update raw-topo gravity wave parameterization 
+  dataset. Update interfaces in `make_hybrid_spaces` to support new 
+  inputs using `SpaceVaryingInput` utility. Include a simple example
+  to generate spectra from scalar variables. 
+  PR [3378](https://github.com/CliMA/ClimaAtmos.jl/pull/3378)
+
 v0.27.7
 -------
 
diff --git a/artifacts/artifact_funcs.jl b/artifacts/artifact_funcs.jl
index 00a1cff1d6c..491cacacea9 100644
--- a/artifacts/artifact_funcs.jl
+++ b/artifacts/artifact_funcs.jl
@@ -50,19 +50,6 @@ function mima_gwf_path()
     return AW.get_data_folder(mima_data)
 end
 
-function topo_elev_dataset_path()
-    etopo1_elev_data = AW.ArtifactWrapper(
-        @__DIR__,
-        "topo-elev-info",
-        AW.ArtifactFile[AW.ArtifactFile(
-            url = "https://caltech.box.com/shared/static/gvilybsu5avxso1wubxjthpip9skc6mf.nc",
-            filename = "ETOPO1_coarse.nc",
-        ),],
-    )
-    return AW.get_data_folder(etopo1_elev_data)
-end
-
-
 function gfdl_ogw_data_path()
     gfdl_data = AW.ArtifactWrapper(
         @__DIR__,
diff --git a/artifacts/download_artifacts.jl b/artifacts/download_artifacts.jl
index a56ec1d7d75..ac7552100e1 100644
--- a/artifacts/download_artifacts.jl
+++ b/artifacts/download_artifacts.jl
@@ -6,7 +6,6 @@ function trigger_download(lazy_download = true)
     @info "Era single column dataset path:`$(era_single_column_dataset_path())`"
     @info "topo dataset path:`$(topo_res_path())`"
     @info "MiMA convective gravity wave path:`$(mima_gwf_path())`"
-    @info "ETOPO1 arc-minute relief model:`$(topo_elev_dataset_path())`"
     @info "GFDL OGWD test data:`$(gfdl_ogw_data_path())`"
     return nothing
 end
diff --git a/config/default_configs/default_config.yml b/config/default_configs/default_config.yml
index 82c44fd0dc6..7fda1d8dcf6 100644
--- a/config/default_configs/default_config.yml
+++ b/config/default_configs/default_config.yml
@@ -49,9 +49,6 @@ sleve_s:
 topo_smoothing:
   help: "Choose whether to order-2 smoothing on the LGL mesh"
   value: false
-smoothing_order:
-  help: "Define the surface smoothing kernel factor (integer) [`3 (default)`]"
-  value: 3
 # ODE
 use_newton_rtol:
   help: "Whether to check if the current iteration of Newton's method has an error within a relative tolerance, instead of always taking the maximum number of iterations (only for ClimaTimeSteppers.jl)"
diff --git a/examples/Project.toml b/examples/Project.toml
index 1a9183986e7..1b7c3a78d65 100644
--- a/examples/Project.toml
+++ b/examples/Project.toml
@@ -1,5 +1,4 @@
 [deps]
-Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
 ArtifactWrappers = "a14bc488-3040-4b00-9dc1-f6467924858a"
 AtmosphericProfilesLibrary = "86bc3604-9858-485a-bdbe-831ec50de11d"
diff --git a/examples/topography_spectra.jl b/examples/topography_spectra.jl
new file mode 100644
index 00000000000..12b47cf0760
--- /dev/null
+++ b/examples/topography_spectra.jl
@@ -0,0 +1,129 @@
+using LazyArtifacts
+using ClimaAtmos.AtmosArtifacts
+import ClimaAtmos as CA
+import ClimaCore as CC
+import ClimaComms
+import ClimaParams as CP
+import ClimaCore.Quadratures
+using ClimaUtilities: SpaceVaryingInputs
+using ClimaCore.Hypsography: diffuse_surface_elevation!
+using ClimaCore: Remapping, Geometry, Operators, Spaces
+using CairoMakie
+
+import ClimaCoreSpectra: power_spectrum_2d
+
+const AA = AtmosArtifacts
+
+
+function mask(x::FT) where {FT}
+    return x * FT(x > 0)
+end
+
+"""
+    generate_spaces(; h_elem=60, diffiter=32)
+For a given number of elements `h_elem`, this function generates a 
+horizontal spectral space (ClimaCore.SpectralElementSpace2D) and
+generates an interpolated (using SpaceVaryingInputs) orography field
+on this space. The `diffuse_surface_elevation` function is called 
+to apply `diffiter` instances of Laplacian diffusion passes 
+(which rely on the canonical unsteady diffusion equation).
+Returns a `Field` containing surface orography. 
+"""
+function generate_spaces(; h_elem = 60, diffiter = 32, planet_radius = 6.378e6)
+    FT = Float32
+    cubed_sphere_mesh =
+        CA.cubed_sphere_mesh(; radius = FT(planet_radius), h_elem)
+    quad = Quadratures.GLL{4}()
+    comms_ctx = ClimaComms.context()
+    h_space = CA.make_horizontal_space(cubed_sphere_mesh, quad, comms_ctx, true)
+    Δh_scale = Spaces.node_horizontal_length_scale(h_space)
+    @assert h_space isa CC.Spaces.SpectralElementSpace2D
+    coords = CC.Fields.coordinate_field(h_space)
+    target_field = CC.Fields.zeros(h_space)
+    elev_from_file = SpaceVaryingInputs.SpaceVaryingInput(
+        AA.earth_orography_file_path(; context = comms_ctx),
+        "z",
+        h_space,
+    )
+    elev_from_file = @. mask(elev_from_file)
+    diffuse_surface_elevation!(
+        elev_from_file,
+        κ = FT(28e7 * (30 / h_elem)^2),
+        dt = FT(1),
+        maxiter = diffiter,
+    )
+    elev_from_file = @. mask(elev_from_file)
+    return elev_from_file
+end
+
+"""
+    remap_to_array(test_var, hcoords)
+Given an Array of `hcoords`, this function remaps a ClimaCore
+Field `test_var` onto `hcoords`. 
+Returns an `Array` containing the values which make up `test_var`.
+"""
+function remap_to_array(test_var, hcoords)
+    remapper = Remapping.Remapper(axes(test_var), hcoords)
+    orog = Array(Remapping.interpolate(remapper, test_var))
+end
+
+"""
+    gen_spectra(test_var)
+Given a Field `test_var`, this function first calls `remap_to_array`
+to remap this spectral element field onto Gaussian long-lat
+coordinates, and computes the spectra.
+Returns (`w_numbers`, `power_spectrum`) which contain the spherical
+wavenumbers and the spectral energy corresponding to this array
+of `w_numbers`. 
+"""
+function gen_spectra(test_var)
+    # Remap onto 1° resolution horizontal coordinates. 
+    longpts = range(-180.0, 180.0, 360)
+    latpts = range(-90.0, 90.0, 180)
+    hcoords =
+        [Geometry.LatLongPoint(lat, long) for long in longpts, lat in latpts]
+    test_var = remap_to_array(test_var, hcoords)
+    # Use ClimaCoreSpectra to generate power spectra for orography data. 
+    len1 = size(test_var)[1]
+    len2 = size(test_var)[2]
+    @assert len1 == 2len2
+    FT = eltype(test_var)
+    mass_weight = FT(1) # No weighting applied 
+    spectrum_data, wave_numbers, _spherical, mesh_info =
+        power_spectrum_2d(FT, test_var, mass_weight)
+    power_spectrum =
+        dropdims(sum(spectrum_data, dims = 1), dims = 1)[begin:(end - 1), :]
+    w_numbers = collect(0:1:(mesh_info.num_spherical - 1))
+    @info "Returning wavenumber array (spherical) and orography power spectrum"
+    return w_numbers, power_spectrum
+end
+
+"""
+    generate_all_spectra(;h_elem=60)
+Uses the spectral calculator and space generation tools in this example file
+to generate a series of surface orography fields with different extents of 
+Laplacian smoothing. 
+Returns a `CairoMakie.Figure`. 
+"""
+function generate_all_spectra(; h_elem = 60)
+    fig = Figure(; size = (1200, 900))
+    ax1 = Axis(
+        fig[1, 1],
+        xlabel = "spherical wavenumber",
+        ylabel = "log(elevation spectra)",
+        title = "Surface elevation spectra",
+        xticks = [2, 4, 8, 16, 32, 64, 128, 256],
+        yticks = [10^-1, 10^0, 10^1, 10^2, 10^3, 10^4],
+        xscale = log10,
+        yscale = log10,
+        xgridvisible = true,
+        ygridvisible = false,
+    )
+    for ii in (0, 4, 8, 16, 32, 64, 128, 256)
+        test_var = generate_spaces(; h_elem, diffiter = ii)
+        sph_wn, psd = gen_spectra(test_var)
+        scatterlines!(sph_wn[2:end], psd[2:end], label = "$(ii)×")
+    end
+    CairoMakie.Legend(fig[1, 2], ax1)
+    return fig
+end
diff --git a/reproducibility_tests/ref_counter.jl b/reproducibility_tests/ref_counter.jl
index 9c803a683d2..862549f2c4f 100644
--- a/reproducibility_tests/ref_counter.jl
+++ b/reproducibility_tests/ref_counter.jl
@@ -1,4 +1,4 @@
-185
+186
 
 # **README**
 #
@@ -20,6 +20,14 @@
 
 
 #=
+
+186
+- Topography dataset has been modified to the 60 arc-second ETOPO2022 dataset. 
+  This is behaviour changing for the gravity-wave (raw-topo) parameterization
+  when computing `hmax` and `T tensor`.
+ 
+185
+
 184
 - Changed default ozone profile.
   Jobs that failed:
diff --git a/src/parameterized_tendencies/gravity_wave_drag/orographic_gravity_wave.jl b/src/parameterized_tendencies/gravity_wave_drag/orographic_gravity_wave.jl
index 18a8dcb196a..b4ee306ea05 100644
--- a/src/parameterized_tendencies/gravity_wave_drag/orographic_gravity_wave.jl
+++ b/src/parameterized_tendencies/gravity_wave_drag/orographic_gravity_wave.jl
@@ -39,7 +39,8 @@ function orographic_gravity_wave_cache(Y, ogw::OrographicGravityWave)
     elseif ogw.topo_info == "raw_topo"
         # TODO: right now this option may easily crash
         # because we did not incorporate any smoothing when interpolate back to model grid
-        elevation_rll = joinpath(topo_elev_dataset_path(), "ETOPO1_coarse.nc")
+        elevation_rll =
+            AA.earth_orography_file_path(; context = ClimaComms.context(Y.c))
         radius =
             Spaces.topology(
                 Spaces.horizontal_space(axes(Y.c)),
diff --git a/src/parameterized_tendencies/gravity_wave_drag/orographic_gravity_wave_helper.jl b/src/parameterized_tendencies/gravity_wave_drag/orographic_gravity_wave_helper.jl
index ba93a126f6e..3108922d46c 100644
--- a/src/parameterized_tendencies/gravity_wave_drag/orographic_gravity_wave_helper.jl
+++ b/src/parameterized_tendencies/gravity_wave_drag/orographic_gravity_wave_helper.jl
@@ -226,10 +226,13 @@ end
 
 function compute_OGW_info(Y, elev_data, earth_radius, γ, h_frac)
     # obtain lat, lon, elevation from the elev_data
+    FT = Spaces.undertype(Spaces.axes(Y.c))
+    # downsample to elev dims (3600×1800)
+    skip_pt = 6
     nt = NCDataset(elev_data, "r") do ds
-        lon = Array(ds["longitude"])
-        lat = Array(ds["latitude"])
-        elev = Array(ds["elevation"])
+        lon = FT.(Array(ds["lon"]))[1:skip_pt:end]
+        lat = FT.(Array(ds["lat"]))[1:skip_pt:end]
+        elev = FT.(Array(ds["z"]))[1:skip_pt:end, 1:skip_pt:end]
         (; lon, lat, elev)
     end
     (; lon, lat, elev) = nt
diff --git a/src/solver/type_getters.jl b/src/solver/type_getters.jl
index b267cb5d1d5..d83152ed269 100644
--- a/src/solver/type_getters.jl
+++ b/src/solver/type_getters.jl
@@ -1,4 +1,3 @@
-using Adapt
 using Dates: DateTime, @dateformat_str
 using Interpolations
 import NCDatasets
@@ -136,39 +135,6 @@ function get_spaces(parsed_args, params, comms_ctx)
     bubble = parsed_args["bubble"]
     deep = parsed_args["deep_atmosphere"]
 
-    @assert topography in ("NoWarp", "DCMIP200", "Earth", "Agnesi", "Schar")
-    if topography == "DCMIP200"
-        warp_function = topography_dcmip200
-    elseif topography == "Agnesi"
-        warp_function = topography_agnesi
-    elseif topography == "Schar"
-        warp_function = topography_schar
-    elseif topography == "NoWarp"
-        warp_function = nothing
-    elseif topography == "Earth"
-        data_path = joinpath(topo_elev_dataset_path(), "ETOPO1_coarse.nc")
-        array_type = ClimaComms.array_type(comms_ctx.device)
-        earth_spline = NCDatasets.NCDataset(data_path) do data
-            zlevels = Array(data["elevation"])
-            lon = Array(data["longitude"])
-            lat = Array(data["latitude"])
-            # Apply Smoothing
-            smooth_degree = Int(parsed_args["smoothing_order"])
-            esmth = CA.gaussian_smooth(zlevels, smooth_degree)
-            Adapt.adapt(
-                array_type,
-                linear_interpolation(
-                    (lon, lat),
-                    esmth,
-                    extrapolation_bc = (Periodic(), Flat()),
-                ),
-            )
-        end
-        @info "Generated interpolation stencil"
-        warp_function = generate_topography_warp(earth_spline)
-    end
-    @info "Topography" topography
-
 
     h_elem = parsed_args["h_elem"]
     radius = CAP.planet_radius(params)
@@ -183,19 +149,7 @@ function get_spaces(parsed_args, params, comms_ctx)
         else
             Meshes.Uniform()
         end
-        if warp_function == nothing
-            make_hybrid_spaces(h_space, z_max, z_elem, z_stretch; deep)
-        else
-            make_hybrid_spaces(
-                h_space,
-                z_max,
-                z_elem,
-                z_stretch;
-                parsed_args = parsed_args,
-                surface_warp = warp_function,
-                deep,
-            )
-        end
+        make_hybrid_spaces(h_space, z_max, z_elem, z_stretch; deep, parsed_args)
     elseif parsed_args["config"] == "column" # single column
         @warn "perturb_initstate flag is ignored for single column configuration"
         FT = eltype(params)
@@ -240,15 +194,7 @@ function get_spaces(parsed_args, params, comms_ctx)
         else
             Meshes.Uniform()
         end
-        make_hybrid_spaces(
-            h_space,
-            z_max,
-            z_elem,
-            z_stretch;
-            parsed_args,
-            surface_warp = warp_function,
-            deep,
-        )
+        make_hybrid_spaces(h_space, z_max, z_elem, z_stretch; parsed_args, deep)
     elseif parsed_args["config"] == "plane"
         FT = eltype(params)
         nh_poly = parsed_args["nh_poly"]
@@ -264,15 +210,7 @@ function get_spaces(parsed_args, params, comms_ctx)
         else
             Meshes.Uniform()
         end
-        make_hybrid_spaces(
-            h_space,
-            z_max,
-            z_elem,
-            z_stretch;
-            parsed_args,
-            surface_warp = warp_function,
-            deep,
-        )
+        make_hybrid_spaces(h_space, z_max, z_elem, z_stretch; parsed_args, deep)
     end
     ncols = Fields.ncolumns(center_space)
     ndofs_total = ncols * z_elem
diff --git a/src/topography/topography.jl b/src/topography/topography.jl
index 57547749803..2db782847a1 100644
--- a/src/topography/topography.jl
+++ b/src/topography/topography.jl
@@ -11,11 +11,11 @@ function topography_dcmip200(coords)
     FT = eltype(λ)
     ϕₘ = FT(0) # degrees (equator)
     λₘ = FT(3 / 2 * 180)  # degrees
-    rₘ = @. FT(acos(sind(ϕₘ) * sind(ϕ) + cosd(ϕₘ) * cosd(ϕ) * cosd(λ - λₘ))) # Great circle distance (rads)
+    rₘ = FT(acos(sind(ϕₘ) * sind(ϕ) + cosd(ϕₘ) * cosd(ϕ) * cosd(λ - λₘ))) # Great circle distance (rads)
     Rₘ = FT(3π / 4) # Moutain radius
     ζₘ = FT(π / 16) # Mountain oscillation half-width
     h₀ = FT(2000)
-    zₛ = @. ifelse(
+    zₛ = ifelse(
         rₘ < Rₘ,
         FT(h₀ / 2) * (1 + cospi(rₘ / Rₘ)) * (cospi(rₘ / ζₘ))^2,
         FT(0),
@@ -23,19 +23,6 @@ function topography_dcmip200(coords)
     return zₛ
 end
 
-function generate_topography_warp(earth_spline)
-    function topography_earth(coords)
-        λ, Φ = coords.long, coords.lat
-        FT = eltype(λ)
-        @info "Load spline"
-        elevation = @. FT(earth_spline(λ, Φ))
-        zₛ = @. ifelse(elevation > FT(0), elevation, FT(0))
-        @info "Assign elevation"
-        return zₛ
-    end
-    return topography_earth
-end
-
 """
     topography_agnesi(x,z)
 x = horizontal coordinate [m]
@@ -52,7 +39,7 @@ function topography_agnesi(coords)
     h_c = FT(1)
     a_c = FT(10000)
     x_c = FT(120000)
-    zₛ = @. h_c / (1 + ((x - x_c) / a_c)^2)
+    zₛ = h_c / (1 + ((x - x_c) / a_c)^2)
     return zₛ
 end
 
@@ -74,6 +61,6 @@ function topography_schar(coords)
     λ_c = FT(4000)
     a_c = FT(5000)
     x_c = FT(60000)
-    zₛ = @. h_c * exp(-((x - x_c) / a_c)^2) * (cospi((x - x_c) / λ_c))^2
+    zₛ = h_c * exp(-((x - x_c) / a_c)^2) * (cospi((x - x_c) / λ_c))^2
     return zₛ
 end
diff --git a/src/utils/AtmosArtifacts.jl b/src/utils/AtmosArtifacts.jl
index b89e551fd18..2a3092406c4 100644
--- a/src/utils/AtmosArtifacts.jl
+++ b/src/utils/AtmosArtifacts.jl
@@ -67,4 +67,19 @@ function aerosol_concentration_file_path(; context = nothing)
     return res_file_path("aerosol_concentrations"; context)
 end
 
+"""
+    earth_orography_file_path(; context=nothing)
+
+Construct the file path for the 60arcsecond orography data NetCDF file.
+
+Downloads the 60arc-second dataset by default. 
+"""
+function earth_orography_file_path(; context = nothing)
+    filename = "ETOPO_2022_v1_60s_N90W180_surface.nc"
+    return joinpath(
+        @clima_artifact("earth_orography_60arcseconds", context),
+        filename,
+    )
+end
+
 end
diff --git a/src/utils/common_spaces.jl b/src/utils/common_spaces.jl
index 7feb5f38059..0d08bc5a41b 100644
--- a/src/utils/common_spaces.jl
+++ b/src/utils/common_spaces.jl
@@ -1,6 +1,7 @@
 using ClimaCore:
-    Geometry, Domains, Meshes, Topologies, Spaces, Grids, Hypsography
+    Geometry, Domains, Meshes, Topologies, Spaces, Grids, Hypsography, Fields
 using ClimaComms
+using ClimaUtilities: SpaceVaryingInputs.SpaceVaryingInput
 
 function periodic_line_mesh(; x_max, x_elem)
     domain = Domains.IntervalDomain(
@@ -79,13 +80,10 @@ function make_hybrid_spaces(
     z_max,
     z_elem,
     z_stretch;
-    surface_warp = nothing,
-    topo_smoothing = false,
     deep = false,
     parsed_args = nothing,
 )
     FT = eltype(z_max)
-    # TODO: change this to make_hybrid_grid
     h_grid = Spaces.grid(h_space)
     z_domain = Domains.IntervalDomain(
         Geometry.ZPoint(zero(z_max)),
@@ -100,11 +98,65 @@ function make_hybrid_spaces(
         z_mesh,
     )
     z_grid = Grids.FiniteDifferenceGrid(z_topology)
-    if isnothing(surface_warp)
+
+    topography = parsed_args["topography"]
+    @assert topography in ("NoWarp", "DCMIP200", "Earth", "Agnesi", "Schar")
+    if topography == "DCMIP200"
+        z_surface = SpaceVaryingInput(topography_dcmip200, h_space)
+        @info "Computing DCMIP200 orography on spectral horizontal space"
+    elseif topography == "Agnesi"
+        z_surface = SpaceVaryingInput(topography_agnesi, h_space)
+        @info "Computing Agnesi orography on spectral horizontal space"
+    elseif topography == "Schar"
+        z_surface = SpaceVaryingInput(topography_schar, h_space)
+        @info "Computing Schar orography on spectral horizontal space"
+    elseif topography == "NoWarp"
+        z_surface = zeros(h_space)
+        @info "No surface orography warp applied"
+    elseif topography == "Earth"
+        z_surface = SpaceVaryingInput(
+            AA.earth_orography_file_path(;
+                context = ClimaComms.context(h_space),
+            ),
+            "z",
+            h_space,
+        )
+        @info "Remapping Earth orography from ETOPO2022 data onto horizontal space"
+    end
+
+    topo_smoothing = parsed_args["topo_smoothing"]
+    if topography == "NoWarp"
         hypsography = Hypsography.Flat()
+    elseif topography == "Earth"
+        mask(x::FT) where {FT} = x * FT(x > 0)
+        z_surface = @. mask(z_surface)
+        # Coefficient for horizontal diffusion of surface orography
+        # is obtained from the V1, V2 Topography documentation from 
+        # E3SM
+        # https://acme-climate.atlassian.net/wiki/spaces/DOC/pages/1456603764/V1+Topography+GLL+grids
+        nelem = Meshes.n_elements_per_panel_direction(axes(z_surface))
+        Hypsography.diffuse_surface_elevation!(
+            z_surface;
+            κ = FT(28e7 * (30 / nelem)^2),
+            dt = FT(1),
+            maxiter = 128,
+        )
+        z_surface = @. mask(z_surface)
+        if parsed_args["mesh_warp_type"] == "SLEVE"
+            @info "SLEVE mesh warp"
+            hypsography = Hypsography.SLEVEAdaption(
+                Geometry.ZPoint.(z_surface),
+                FT(parsed_args["sleve_eta"]),
+                FT(parsed_args["sleve_s"]),
+            )
+        elseif parsed_args["mesh_warp_type"] == "Linear"
+            @info "Linear mesh warp"
+            hypsography =
+                Hypsography.LinearAdaption(Geometry.ZPoint.(z_surface))
+        else
+            @error "Undefined mesh-warping option"
+        end
     else
-        topo_smoothing = parsed_args["topo_smoothing"]
-        z_surface = surface_warp(Fields.coordinate_field(h_space))
         if topo_smoothing
             Hypsography.diffuse_surface_elevation!(z_surface)
         end
@@ -121,8 +173,8 @@ function make_hybrid_spaces(
                 Hypsography.LinearAdaption(Geometry.ZPoint.(z_surface))
         end
     end
+
     grid = Grids.ExtrudedFiniteDifferenceGrid(h_grid, z_grid, hypsography; deep)
-    # TODO: return the grid
     center_space = Spaces.CenterExtrudedFiniteDifferenceSpace(grid)
     face_space = Spaces.FaceExtrudedFiniteDifferenceSpace(grid)
     return center_space, face_space
diff --git a/src/utils/utilities.jl b/src/utils/utilities.jl
index fffbc4a7b21..9f6c8603257 100644
--- a/src/utils/utilities.jl
+++ b/src/utils/utilities.jl
@@ -352,58 +352,6 @@ function horizontal_integral_at_boundary(f::Fields.Field)
     sum(f ./ Fields.Δz_field(axes(f)) .* 2) # TODO: is there a way to ensure this is derived from face z? The 2d topology doesn't contain this info
 end
 
-"""
-    gaussian_smooth(arr, sigma = 1)
-
-Smooth the given 2D array by applying a Gaussian blur.
-
-Edges are not properly smoothed out: the edge value is extended to infinity.
-"""
-function gaussian_smooth(arr::AbstractArray, sigma::Int = 1)
-    n1, n2 = size(arr)
-
-    # We assume that the Gaussian goes to zero at 10 sigma and ignore contributions outside of that window
-    window = Int(ceil(10 * sigma))
-
-    # Prepare the 2D Gaussian kernel
-    gauss = [
-        exp.(-(x .^ 2 .+ y .^ 2) / (2 * sigma^2)) for
-        x in range(-window, window), y in range(-window, window)
-    ]
-
-    # Normalize it
-    gauss = gauss / sum(gauss)
-
-    smoothed_arr = zeros(size(arr))
-
-    # 2D convolution
-    for i in 1:n1
-        for j in 1:n2
-            # For each point, we "look left and right (up and down)" within our window
-            for wx in (-window):window
-                for wy in (-window):window
-                    # For values at the edge, we keep using the edge value
-                    k = clamp(i + wx, 1, n1)
-                    l = clamp(j + wy, 1, n2)
-
-                    # gauss has size 2window + 1, so its midpoint (when the gaussian is max)
-                    # is at 1 + window
-                    #
-                    # Eg, for window of 3, wx will go through the values -3, -2, 1, 0, 1, 2, 3,
-                    # and the midpoint is 4 (= 1 + window)
-                    mid_gauss_idx = 1 + window
-
-                    smoothed_arr[i, j] +=
-                        arr[k, l] *
-                        gauss[mid_gauss_idx + wx, mid_gauss_idx + wy]
-                end
-            end
-        end
-    end
-
-    return smoothed_arr
-end
-
 """
     isdivisible(dt_large::Dates.Period, dt_small::Dates.Period)
 
diff --git a/test/parameterized_tendencies/gravity_wave/orographic_gravity_wave/ogwd_3d.jl b/test/parameterized_tendencies/gravity_wave/orographic_gravity_wave/ogwd_3d.jl
index 512d7a3bbba..b2ded086ade 100644
--- a/test/parameterized_tendencies/gravity_wave/orographic_gravity_wave/ogwd_3d.jl
+++ b/test/parameterized_tendencies/gravity_wave/orographic_gravity_wave/ogwd_3d.jl
@@ -6,6 +6,8 @@ import ClimaAtmos as CA
 import Thermodynamics as TD
 import ClimaParams as CP
 import ClimaComms
+import ClimaAtmos: AtmosArtifacts as AA
+import ClimaUtilities: SpaceVaryingInputs.SpaceVaryingInput
 using ClimaCoreTempestRemap
 
 using Interpolations
@@ -20,6 +22,7 @@ comms_ctx = ClimaComms.SingletonCommsContext()
 (; config_file, job_id) = CA.commandline_kwargs()
 config = CA.AtmosConfig(config_file; job_id, comms_ctx)
 
+config.parsed_args["topography"] = "Earth";
 config.parsed_args["topo_smoothing"] = false;
 config.parsed_args["mesh_warp_type"] = "Linear";
 (; parsed_args) = config
@@ -77,24 +80,6 @@ for k in 1:size(bk)[1]
 end
 p_center = 0.5 * (p_half[:, :, 1:(end - 1), :] .+ p_half[:, :, 2:end, :])
 
-# earth warp
-data_path = joinpath(topo_elev_dataset_path(), "ETOPO1_coarse.nc")
-earth_spline = NCDataset(data_path) do data
-    zlevels = Array(data["elevation"])
-    lon = Array(data["longitude"])
-    lat = Array(data["latitude"])
-    # Apply Smoothing
-    smooth_degree = 15
-    esmth = CA.gaussian_smooth(zlevels, smooth_degree)
-    linear_interpolation(
-        (lon, lat),
-        zlevels,
-        extrapolation_bc = (Periodic(), Flat()),
-    )
-end
-@info "Generated interpolation stencil"
-warp_function = CA.generate_topography_warp(earth_spline)
-
 # Create meshes and spaces
 h_elem = 16
 nh_poly = 3
@@ -106,17 +91,9 @@ radius = 6.371229e6
 quad = Quadratures.GLL{nh_poly + 1}()
 horizontal_mesh = CA.cubed_sphere_mesh(; radius, h_elem)
 h_space = CA.make_horizontal_space(horizontal_mesh, quad, comms_ctx, false)
-
-
 z_stretch = Meshes.HyperbolicTangentStretching(dz_bottom)
-center_space, face_space = CA.make_hybrid_spaces(
-    h_space,
-    z_max,
-    z_elem,
-    z_stretch;
-    parsed_args,
-    surface_warp = warp_function,
-)
+center_space, face_space =
+    CA.make_hybrid_spaces(h_space, z_max, z_elem, z_stretch; parsed_args)
 
 ᶜlocal_geometry = Fields.local_geometry_field(center_space)
 ᶠlocal_geometry = Fields.local_geometry_field(face_space)
diff --git a/test/parameterized_tendencies/gravity_wave/orographic_gravity_wave/ogwd_baseflux.jl b/test/parameterized_tendencies/gravity_wave/orographic_gravity_wave/ogwd_baseflux.jl
index 4523e5352f6..84cbc85d45a 100644
--- a/test/parameterized_tendencies/gravity_wave/orographic_gravity_wave/ogwd_baseflux.jl
+++ b/test/parameterized_tendencies/gravity_wave/orographic_gravity_wave/ogwd_baseflux.jl
@@ -13,6 +13,9 @@ include(
 )
 
 comms_ctx = ClimaComms.SingletonCommsContext()
+(; config_file, job_id) = CA.commandline_kwargs()
+config = CA.AtmosConfig(config_file; job_id, comms_ctx)
+config.parsed_args["topography"] = "NoWarp"
 
 # Create meshes and spaces
 h_elem = 6
@@ -25,8 +28,13 @@ quad = Quadratures.GLL{nh_poly + 1}()
 horizontal_mesh = CA.cubed_sphere_mesh(; radius, h_elem)
 h_space = CA.make_horizontal_space(horizontal_mesh, quad, comms_ctx, false)
 z_stretch = Meshes.Uniform()
-center_space, face_space =
-    CA.make_hybrid_spaces(h_space, z_max, z_elem, z_stretch)
+center_space, face_space = CA.make_hybrid_spaces(
+    h_space,
+    z_max,
+    z_elem,
+    z_stretch;
+    parsed_args = config.parsed_args,
+)
 
 ᶜlocal_geometry = Fields.local_geometry_field(center_space)
 ᶠlocal_geometry = Fields.local_geometry_field(face_space)
diff --git a/test/topography.jl b/test/topography.jl
index a99e83dc411..66743e5e5af 100644
--- a/test/topography.jl
+++ b/test/topography.jl
@@ -30,8 +30,8 @@ include("test_helpers.jl")
 
 @testset "test topography functions" begin
     (; FT, coords) = get_spherical_spaces()
-    @test extrema(CA.topography_dcmip200(coords)) == (FT(0), FT(2000))
-    loc = findmax(parent(CA.topography_dcmip200(coords)))
+    @test extrema(CA.topography_dcmip200.(coords)) == (FT(0), FT(2000))
+    loc = findmax(parent(CA.topography_dcmip200.(coords)))
     @test parent(coords.lat)[loc[2]] == FT(0)
     @test parent(coords.long)[loc[2]] == FT(-90)
 end
diff --git a/test/utilities.jl b/test/utilities.jl
index c7247ce5a39..d934ff095ae 100644
--- a/test/utilities.jl
+++ b/test/utilities.jl
@@ -20,17 +20,6 @@ include("test_helpers.jl")
     @test CA.sort_files_by_time(fns) == filenames.(t_sorted)
 end
 
-@testset "gaussian_smooth" begin
-    # No smooth on constant
-    @test CA.gaussian_smooth(3.0 * ones(132, 157)) ≈ 3.0 * ones(132, 157)
-    randy = rand(123, 145)
-    smoothed = CA.gaussian_smooth(randy)
-    # min
-    @test extrema(randy)[1] <= smoothed[1]
-    # max
-    @test extrema(randy)[2] >= smoothed[2]
-end
-
 @testset "isdivisible" begin
     @test CA.isdivisible(Dates.Month(1), Dates.Day(1))
     @test !CA.isdivisible(Dates.Month(1), Dates.Day(25))
@@ -243,6 +232,9 @@ end
 @testset "make hybrid spaces" begin
     (; cent_space, face_space, xlim, zlim, velem, helem, npoly, quad) =
         get_cartesian_spaces()
+    config = CA.AtmosConfig(
+        Dict("topography" => "NoWarp", "topo_smoothing" => false),
+    )
     device = ClimaComms.CPUSingleThreaded()
     comms_ctx = ClimaComms.context(device)
     z_stretch = Meshes.Uniform()
@@ -259,9 +251,8 @@ end
         zlim[2],
         velem,
         z_stretch;
-        surface_warp = nothing,
-        topo_smoothing = false,
         deep = false,
+        parsed_args = config.parsed_args,
     )
     @test test_cent_space == cent_space
     @test test_face_space == face_space