diff --git a/examples/ecco_inspect_temperature_salinity.jl b/examples/ecco_inspect_temperature_salinity.jl
index f055d378b..bd0fe68f3 100644
--- a/examples/ecco_inspect_temperature_salinity.jl
+++ b/examples/ecco_inspect_temperature_salinity.jl
@@ -21,11 +21,7 @@ grid = LatitudeLongitudeGrid(arch; z,
                              latitude  = (-80, 80),
                              longitude = (0, 360))
 
-bottom_height = regrid_bathymetry(grid;
-                                  minimum_depth = 10,
-                                  interpolation_passes = 5,
-                                  major_basins = 1)
-
+bottom_height = regrid_bathymetry(grid; minimum_depth = 10, major_basins = 1)
 grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height))
 
 T = CenterField(grid)
diff --git a/examples/ecco_mixed_layer_depth.jl b/examples/ecco_mixed_layer_depth.jl
index 0c6fa1ad1..8f622987a 100644
--- a/examples/ecco_mixed_layer_depth.jl
+++ b/examples/ecco_mixed_layer_depth.jl
@@ -22,11 +22,7 @@ grid = LatitudeLongitudeGrid(arch; z,
                              latitude  = (-80, 80),
                              longitude = (0, 360))
 
-bottom_height = regrid_bathymetry(grid;
-                                  minimum_depth = 10,
-                                  interpolation_passes = 5,
-                                  major_basins = 1)
-
+bottom_height = regrid_bathymetry(grid; minimum_depth = 10, major_basins = 1)
 grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height))
 
 start_date = DateTime(1993, 1, 1)
diff --git a/examples/generate_bathymetry.jl b/examples/generate_bathymetry.jl
index 5da3e4359..431d9120f 100644
--- a/examples/generate_bathymetry.jl
+++ b/examples/generate_bathymetry.jl
@@ -43,7 +43,7 @@ grid = LatitudeLongitudeGrid(size = (Nλ, Nφ, 1),
 #    and fills them with land.
 
 h_rough = regrid_bathymetry(grid)
-h_smooth = regrid_bathymetry(grid; interpolation_passes = 40)
+h_smooth = regrid_bathymetry(grid; smoothing=InterpolationPasses(40))
 h_one_basin = regrid_bathymetry(grid; major_basins = 1)
 nothing #hide
 
diff --git a/examples/mediterranean_simulation_with_ecco_restoring.jl b/examples/mediterranean_simulation_with_ecco_restoring.jl
index 53b8d8d4b..81915fcf5 100644
--- a/examples/mediterranean_simulation_with_ecco_restoring.jl
+++ b/examples/mediterranean_simulation_with_ecco_restoring.jl
@@ -54,13 +54,11 @@ grid = LatitudeLongitudeGrid(GPU();
 # ### Bathymetry Interpolation
 #
 # The script interpolates bathymetric data onto the grid, ensuring that the model accurately represents
-# the sea floor's topography. Parameters such as `minimum_depth` and `interpolation_passes`
-# are adjusted to refine the bathymetry representation.
+# the sea floor's topography. `minimum_depth` is adjusted to refine the bathymetry representation.
 
 bottom_height = regrid_bathymetry(grid,
                                   height_above_water = 1,
                                   minimum_depth = 10,
-                                  interpolation_passes = 25,
                                   connected_regions_allowed = 1)
 
 grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height))
diff --git a/examples/near_global_ocean_simulation.jl b/examples/near_global_ocean_simulation.jl
index 5877931cf..19813cbdd 100644
--- a/examples/near_global_ocean_simulation.jl
+++ b/examples/near_global_ocean_simulation.jl
@@ -55,7 +55,6 @@ grid = LatitudeLongitudeGrid(arch;
 
 bottom_height = regrid_bathymetry(grid;
                                   minimum_depth = 10meters,
-                                  interpolation_passes = 5,
                                   major_basins = 3)
 
 grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height); active_cells_map=true)
diff --git a/examples/one_degree_simulation.jl b/examples/one_degree_simulation.jl
index fba9da108..019b21a72 100644
--- a/examples/one_degree_simulation.jl
+++ b/examples/one_degree_simulation.jl
@@ -46,7 +46,6 @@ underlying_grid = TripolarGrid(arch;
 
 bottom_height = regrid_bathymetry(underlying_grid;
                                   minimum_depth = 10,
-                                  interpolation_passes = 75, # 75 interpolation passes smooth the bathymetry near Florida so that the Gulf Stream is able to flow
 				                  major_basins = 2)
 
 # For this bathymetry at this horizontal resolution we need to manually open the Gibraltar strait.
diff --git a/experiments/one_degree_simulation/one_degree_simulation.jl b/experiments/one_degree_simulation/one_degree_simulation.jl
index bc8bef425..ea32ecc84 100644
--- a/experiments/one_degree_simulation/one_degree_simulation.jl
+++ b/experiments/one_degree_simulation/one_degree_simulation.jl
@@ -15,7 +15,7 @@ Nz = 32
 z_faces = exponential_z_faces(; Nz, depth=6000, h=34)
 underlying_grid = TripolarGrid(arch; size=(Nx, Ny, Nz), z=z_faces)
 
-bottom_height = regrid_bathymetry(underlying_grid; minimum_depth=30, interpolation_passes=20, major_basins=1)
+bottom_height = regrid_bathymetry(underlying_grid; minimum_depth=30, major_basins=1)
 view(bottom_height, 73:78, 88:89, 1) .= -1000 # open Gibraltar strait 
 
 grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bottom_height); active_cells_map=true)
diff --git a/src/Bathymetry.jl b/src/Bathymetry.jl
index 4b6cda0e9..8a4e89f68 100644
--- a/src/Bathymetry.jl
+++ b/src/Bathymetry.jl
@@ -1,6 +1,6 @@
 module Bathymetry
 
-export regrid_bathymetry, retrieve_bathymetry, download_bathymetry
+export regrid_bathymetry, retrieve_bathymetry, download_bathymetry, InterpolationPasses
 
 using ImageMorphology
 using ..DataWrangling: download_progress
@@ -66,7 +66,7 @@ end
                       dir = download_bathymetry_cache,
                       url = ClimaOcean.Bathymetry.etopo_url,
                       filename = "ETOPO_2022_v1_60s_N90W180_surface.nc",
-                      interpolation_passes = 1,
+                      smoothing = nothing
                       major_basins = 1)
 
 Return bathymetry associated with the NetCDF file at `path = joinpath(dir, filename)` regridded onto `target_grid`.
@@ -93,9 +93,8 @@ Keyword Arguments
   2. `lon` vector of longitude nodes
   3. `z` matrix of depth values
 
-- `interpolation_passes`: regridding/interpolation passes. The bathymetry is interpolated in
-                          `interpolation_passes - 1` intermediate steps. The more the interpolation
-                          steps the smoother the final bathymetry becomes.
+- `smoothing`: smoothing method for bathymetry.
+               Default: nothing.
 
   Example
   =======
@@ -123,7 +122,7 @@ function regrid_bathymetry(target_grid;
                            dir = download_bathymetry_cache,
                            url = etopo_url,
                            filename = "ETOPO_2022_v1_60s_N90W180_surface.nc",
-                           interpolation_passes = 1,
+                           smoothing = nothing,
                            major_basins = 1) # Allow an `Inf` number of "lakes"
 
     filepath = download_bathymetry(; url, dir, filename)
@@ -179,8 +178,7 @@ function regrid_bathymetry(target_grid;
     set!(native_z, z_data)
     fill_halo_regions!(native_z)
 
-    target_z = interpolate_bathymetry_in_passes(native_z, target_grid;
-                                                passes = interpolation_passes)
+    target_z = interpolate_bathymetry(native_z, target_grid, smoothing)
 
     if minimum_depth > 0
         zi = interior(target_z, :, :, 1)
@@ -200,67 +198,64 @@ function regrid_bathymetry(target_grid;
 end
 
 # Here we can either use `regrid!` (three dimensional version) or `interpolate!`.
-function interpolate_bathymetry_in_passes(native_z, target_grid;
-                                          passes = 10)
-    Nλt, Nφt = Nt = size(target_grid)
-    Nλn, Nφn = Nn = size(native_z)
-
-    resxt = minimum_xspacing(target_grid)
-    resyt = minimum_yspacing(target_grid)
-
-    resxn = minimum_xspacing(native_z.grid)
-    resyn = minimum_yspacing(native_z.grid)
-
-    # Check whether we are refining the grid in any directions.
-    # If so, skip interpolation passes, as they are not needed.
-    if resxt < resxn || resyt < resyn
-        target_z = Field{Center, Center, Nothing}(target_grid)
-        interpolate!(target_z, native_z)
-        @info string("Skipping passes for interpolating bathymetry of size $Nn ", '\n',
-                     "to target grid of size $Nt. Interpolation passes may only ", '\n',
-                     "be used to coarsen bathymetry and require that the bathymetry ", '\n',
-                     "is finer than the target grid in both horizontal directions.")
-        return target_z
+function interpolate_bathymetry(user_z, target_grid, smoothing)
+    smooth_z = smooth_bathymetry(user_z, target_grid, smoothing)
+    target_z = Field{Center, Center, Nothing}(target_grid)
+    interpolate!(target_z, smooth_z)
+    return target_z
+end
+
+smooth_bathymetry(z, ::Nothing) = z
+
+struct InterpolationPasses
+    passes :: Int
+
+    function InterpolationPasses(passes)
+        if passes < 1
+            throw(ArgumentError("`passes` must be larger than 0"))
+        elseif !isinteger(passes)
+            throw(ArgumentError("`passes` must be an integer"))
+        end
+
+        return new(passes)
     end
+end
 
-    # Interpolate in passes
-    latitude  = y_domain(native_z.grid)
-    longitude = x_domain(native_z.grid)
+function smooth_bathymetry(user_z, target_grid, smoothing::InterpolationPasses)
+    passes = smoothing.passes
+    Nλt, Nφt = Nt = size(target_grid)
+    Nλn, Nφn = Nn = size(user_z)
 
-    ΔNλ = floor((Nλn - Nλt) / passes)
-    ΔNφ = floor((Nφn - Nφt) / passes)
+    # Interpolate in passes
+    longitude = x_domain(user_z.grid)
+    latitude  = y_domain(user_z.grid)
 
-    Nλ = [Nλn - ΔNλ * pass for pass in 1:passes-1]
-    Nφ = [Nφn - ΔNφ * pass for pass in 1:passes-1]
+    ΔNλ = floor((Nλn - Nλt) / (passes - 1))
+    ΔNφ = floor((Nφn - Nφt) / (passes - 1))
 
-    Nλ = Int[Nλ..., Nλt]
-    Nφ = Int[Nφ..., Nφt]
+    Nλ = [Nλn - ΔNλ * pass for pass in 1:passes]
+    Nφ = [Nφn - ΔNφ * pass for pass in 1:passes]
 
-    old_z  = native_z
+    old_z = user_z
     TX, TY = topology(target_grid)
 
-    for pass = 1:passes - 1
+    for pass = 1:passes
         new_size = (Nλ[pass], Nφ[pass], 1)
-
         @debug "Bathymetry interpolation pass $pass with size $new_size"
 
         new_grid = LatitudeLongitudeGrid(architecture(target_grid), Float32,
                                          size = new_size,
-                                         latitude = (latitude[1],  latitude[2]),
+                                         latitude = (latitude[1], latitude[2]),
                                          longitude = (longitude[1], longitude[2]),
                                          z = (0, 1),
                                          topology = (TX, TY, Bounded))
 
         new_z = Field{Center, Center, Nothing}(new_grid)
-
         interpolate!(new_z, old_z)
         old_z = new_z
     end
 
-    target_z = Field{Center, Center, Nothing}(target_grid)
-    interpolate!(target_z, old_z)
-
-    return target_z
+    return new_z
 end
 
 # Regridding bathymetry for distributed grids, we handle the whole process
diff --git a/test/test_bathymetry.jl b/test/test_bathymetry.jl
index 7c5a56312..eb2d43e3d 100644
--- a/test/test_bathymetry.jl
+++ b/test/test_bathymetry.jl
@@ -48,18 +48,6 @@ using ClimaOcean.Bathymetry: remove_minor_basins!
         # of -1000 meters
         @test mean(view(bottom_height, 51:100, :, 1)) == -1000
 
-        grid = LatitudeLongitudeGrid(arch;
-                                     size = (200, 200, 10),
-                                     longitude = (0, 2),
-                                     latitude = (-0.1, 0.1),
-                                     z = (-6000, 0))
-
-        control_bottom_height = regrid_bathymetry(grid)
-        interpolated_bottom_height = regrid_bathymetry(grid; interpolation_passes=100)
-
-        # Testing that multiple passes do not change the solution when refining the grid
-        @test parent(control_bottom_height) == parent(interpolated_bottom_height)
-
         grid = LatitudeLongitudeGrid(arch;
                                      size = (200, 200, 10),
                                      longitude = (0, 100),
@@ -67,7 +55,7 @@ using ClimaOcean.Bathymetry: remove_minor_basins!
                                      z = (-6000, 0))
 
         control_bottom_height = regrid_bathymetry(grid)
-        interpolated_bottom_height = regrid_bathymetry(grid; interpolation_passes=10)
+        interpolated_bottom_height = regrid_bathymetry(grid; smoothing=InterpolationPasses(10))
 
         # Testing that multiple passes _do_ change the solution when coarsening the grid
         @test parent(control_bottom_height) != parent(interpolated_bottom_height)
diff --git a/test/test_diagnostics.jl b/test/test_diagnostics.jl
index bf4aa1ff1..b1af3fd13 100644
--- a/test/test_diagnostics.jl
+++ b/test/test_diagnostics.jl
@@ -16,7 +16,6 @@ using ClimaOcean.Diagnostics: MixedLayerDepthField, MixedLayerDepthOperand
 
         bottom_height = regrid_bathymetry(grid;
                                           minimum_depth = 10,
-                                          interpolation_passes = 5,
                                           major_basins = 1)
 
         grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height))
diff --git a/test/test_distributed_utils.jl b/test/test_distributed_utils.jl
index 8783d2685..2e518d113 100644
--- a/test/test_distributed_utils.jl
+++ b/test/test_distributed_utils.jl
@@ -59,8 +59,7 @@ end
 
     global_height = regrid_bathymetry(global_grid;
                                       dir = "./",
-                                      filename = "trivial_bathymetry.nc",
-                                      interpolation_passes=10)
+                                      filename = "trivial_bathymetry.nc")
 
     arch_x  = Distributed(CPU(), partition=Partition(4, 1))
     arch_y  = Distributed(CPU(), partition=Partition(1, 4))
@@ -75,8 +74,7 @@ end
 
         local_height = regrid_bathymetry(local_grid;
                                          dir = "./",
-                                         filename = "trivial_bathymetry.nc",
-                                         interpolation_passes=10)
+                                         filename = "trivial_bathymetry.nc")
 
         Nx, Ny, _ = size(local_grid)
         rx, ry, _ = arch.local_index
diff --git a/test/test_ocean_sea_ice_model.jl b/test/test_ocean_sea_ice_model.jl
index 4970826b1..a8b6d25a7 100644
--- a/test/test_ocean_sea_ice_model.jl
+++ b/test/test_ocean_sea_ice_model.jl
@@ -52,7 +52,6 @@ using ClimaSeaIce.Rheologies
 
             bottom_height = regrid_bathymetry(grid;
                                               minimum_depth = 10,
-                                              interpolation_passes = 20,
                                               major_basins = 1)
 
             grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height); active_cells_map=true)