diff --git a/compass/ocean/inactive_top_cells.py b/compass/ocean/inactive_top_cells.py
new file mode 100644
index 0000000000..0f1ca7af62
--- /dev/null
+++ b/compass/ocean/inactive_top_cells.py
@@ -0,0 +1,45 @@
+import os
+
+import xarray
+from mpas_tools.io import write_netcdf
+
+
+def remove_inactive_top_cells_output(in_filename, out_filename=None,
+                                     mesh_filename=None):
+    """
+    Remove inactive top cells from the output netCDF file
+
+    Parameters
+    ----------
+    in_filename : str
+        Filename for the netCDF file to be cropped
+
+    out_filename : str, optional
+        Filename for the netCDF file after cropping.  Tbe default name is the
+        original file name with ``_crop`` appended before the extension
+
+    mesh_filename : str, optional
+        Filename for an MPAS mesh if not included in the file to be cropped
+
+    """
+    if not os.path.exists(in_filename):
+        raise OSError(f'File {in_filename} does not exist.')
+
+    if out_filename is None:
+        basename, ext = os.path.splitext(in_filename)
+        out_filename = f'{basename}_crop{ext}'
+
+    with xarray.open_dataset(in_filename) as ds_in:
+        if mesh_filename is None:
+            ds_mesh = ds_in
+        else:
+            ds_mesh = xarray.open_dataset(mesh_filename)
+        minLevelCell = ds_mesh.minLevelCell
+        minval = minLevelCell.min().values
+        maxval = minLevelCell.max().values
+        if minval != maxval:
+            raise ValueError('Expected minLevelCell to have a constant '
+                             'value for inactive top cell tests')
+        ds_out = ds_in.isel(nVertLevels=slice(minval - 1, None))
+
+    write_netcdf(ds_out, out_filename)
diff --git a/compass/ocean/suites/inactive_top.txt b/compass/ocean/suites/inactive_top.txt
new file mode 100644
index 0000000000..2f8b40d60c
--- /dev/null
+++ b/compass/ocean/suites/inactive_top.txt
@@ -0,0 +1,7 @@
+ocean/global_ocean/QU240/mesh
+ocean/global_ocean/QU240/WOA23/init
+ocean/global_ocean/QU240/WOA23/init_inactive_top
+ocean/global_ocean/QU240/WOA23/performance_test
+ocean/global_ocean/QU240/WOA23/performance_test_inactive_top
+ocean/global_ocean/QU240/WOA23/RK4/performance_test
+ocean/global_ocean/QU240/WOA23/RK4/performance_test_inactive_top
diff --git a/compass/ocean/tests/global_ocean/__init__.py b/compass/ocean/tests/global_ocean/__init__.py
index b68a3f6d58..d81e7b5ac9 100644
--- a/compass/ocean/tests/global_ocean/__init__.py
+++ b/compass/ocean/tests/global_ocean/__init__.py
@@ -56,7 +56,8 @@ def __init__(self, mpas_core):
                         include_rk4=True,
                         include_regression=True,
                         include_phc=True,
-                        include_en4_1900=True)
+                        include_en4_1900=True,
+                        with_inactive_top_cells=True)
 
         # for other meshes, we do fewer tests
         self._add_tests(mesh_names=['QU', 'Icos', 'QUwISC', 'IcoswISC'],
@@ -88,7 +89,7 @@ def __init__(self, mpas_core):
     def _add_tests(self, mesh_names, DynamicAdjustment,
                    high_res_topography=True, include_rk4=False,
                    include_regression=False, include_phc=False,
-                   include_en4_1900=False):
+                   include_en4_1900=False, with_inactive_top_cells=False):
         """ Add test cases for the given mesh(es) """
 
         default_ic = 'WOA23'
@@ -100,7 +101,8 @@ def _add_tests(self, mesh_names, DynamicAdjustment,
             self.add_test_case(mesh_test)
 
             init_test = Init(test_group=self, mesh=mesh_test,
-                             initial_condition=default_ic)
+                             initial_condition=default_ic,
+                             with_inactive_top_cells=False)
             self.add_test_case(init_test)
 
             time_integrator = default_time_int
@@ -163,6 +165,21 @@ def _add_tests(self, mesh_names, DynamicAdjustment,
                         test_group=self, mesh=mesh_test, init=init_test,
                         time_integrator=time_integrator))
 
+            if with_inactive_top_cells:
+                init_test = Init(test_group=self, mesh=mesh_test,
+                                 initial_condition=default_ic,
+                                 with_inactive_top_cells=True)
+                self.add_test_case(init_test)
+                self.add_test_case(
+                    PerformanceTest(
+                        test_group=self, mesh=mesh_test, init=init_test,
+                        time_integrator=default_time_int))
+                if include_rk4:
+                    self.add_test_case(
+                        PerformanceTest(
+                            test_group=self, mesh=mesh_test, init=init_test,
+                            time_integrator='RK4'))
+
             initial_conditions = []
             if include_phc:
                 initial_conditions.append('PHC')
diff --git a/compass/ocean/tests/global_ocean/forward.py b/compass/ocean/tests/global_ocean/forward.py
index bcb30eb0f3..2b11836c45 100644
--- a/compass/ocean/tests/global_ocean/forward.py
+++ b/compass/ocean/tests/global_ocean/forward.py
@@ -3,6 +3,7 @@
 from importlib.resources import contents
 
 from compass.model import run_model
+from compass.ocean.inactive_top_cells import remove_inactive_top_cells_output
 from compass.ocean.tests.global_ocean.metadata import (
     add_mesh_and_init_metadata,
 )
@@ -155,6 +156,9 @@ def __init__(self, test_case, mesh, time_integrator, init=None,
                 filename='graph.info',
                 work_dir_target=f'{init.path}/initial_state/graph.info')
 
+            if init.with_inactive_top_cells:
+                self.add_output_file(filename='output_crop.nc')
+
         self.add_model_as_input()
 
     def setup(self):
@@ -214,6 +218,12 @@ def run(self):
         update_pio = self.config.getboolean('global_ocean',
                                             'forward_update_pio')
         run_model(self, update_pio=update_pio)
+
+        if self.init.with_inactive_top_cells:
+            remove_inactive_top_cells_output(in_filename='output.nc',
+                                             out_filename='output_crop.nc',
+                                             mesh_filename='init.nc')
+
         add_mesh_and_init_metadata(self.outputs, self.config,
                                    init_filename='init.nc')
 
diff --git a/compass/ocean/tests/global_ocean/init/__init__.py b/compass/ocean/tests/global_ocean/init/__init__.py
index c59841abeb..8d6081de6e 100644
--- a/compass/ocean/tests/global_ocean/init/__init__.py
+++ b/compass/ocean/tests/global_ocean/init/__init__.py
@@ -24,8 +24,14 @@ class Init(TestCase):
     init_subdir : str
         The subdirectory within the test group for all test cases with this
         initial condition
+
+    inactive_top_comp_subdir : str
+        If ``with_inactive_top_cells == True``, the subdirectory equivalent to
+        ``init_subdir`` for test cases without inactive top cells for use for
+        validation between tests with and without inactive top cells
     """
-    def __init__(self, test_group, mesh, initial_condition):
+    def __init__(self, test_group, mesh, initial_condition,
+                 with_inactive_top_cells=False):
         """
         Create the test case
 
@@ -43,17 +49,25 @@ def __init__(self, test_group, mesh, initial_condition):
         name = 'init'
         mesh_name = mesh.mesh_name
         ic_dir = initial_condition
-        self.init_subdir = os.path.join(mesh_name, ic_dir)
+        init_subdir = os.path.join(mesh_name, ic_dir)
+        if with_inactive_top_cells:
+            # Add the name of the subdir without inactive top cells whether or
+            # not is has or will be run
+            self.inactive_top_comp_subdir = init_subdir
+            name = f'{name}_inactive_top'
+        self.init_subdir = init_subdir
         subdir = os.path.join(self.init_subdir, name)
         super().__init__(test_group=test_group, name=name, subdir=subdir)
 
         self.mesh = mesh
         self.initial_condition = initial_condition
+        self.with_inactive_top_cells = with_inactive_top_cells
 
         self.add_step(
             InitialState(
                 test_case=self, mesh=mesh,
-                initial_condition=initial_condition))
+                initial_condition=initial_condition,
+                with_inactive_top_cells=with_inactive_top_cells))
 
         if mesh.with_ice_shelf_cavities:
             self.add_step(RemapIceShelfMelt(test_case=self, mesh=mesh))
@@ -84,6 +98,15 @@ def configure(self, config=None):
         config.set('global_ocean', 'init_description',
                    descriptions[initial_condition])
 
+        if self.with_inactive_top_cells:
+            # Since we start at minLevelCell = 2, we need to increase the
+            # number of vertical levels in the cfg file to end up with the
+            # intended number in the initial state
+            vert_levels = config.getint('vertical_grid', 'vert_levels')
+            config.set('vertical_grid', 'vert_levels', f'{vert_levels + 1}',
+                       comment='active vertical levels + inactive_top_cells')
+            config.set('vertical_grid', 'inactive_top_cells', '1')
+
     def validate(self):
         """
         Test cases can override this method to perform validation of variables
@@ -93,6 +116,27 @@ def validate(self):
         compare_variables(test_case=self, variables=variables,
                           filename1='initial_state/initial_state.nc')
 
+        if self.with_inactive_top_cells:
+            # construct the work directory for the other test
+            filename2 = os.path.join(self.base_work_dir, self.mpas_core.name,
+                                     self.test_group.name,
+                                     self.inactive_top_comp_subdir,
+                                     'init/initial_state/initial_state.nc')
+            if os.path.exists(filename2):
+                variables = ['temperature', 'salinity', 'layerThickness',
+                             'normalVelocity']
+                compare_variables(
+                    test_case=self, variables=variables,
+                    filename1='initial_state/initial_state_crop.nc',
+                    filename2=filename2,
+                    quiet=False, check_outputs=False,
+                    skip_if_step_not_run=False)
+
+            else:
+                self.logger.warn('The version of "init" without inactive top '
+                                 'cells was not run.  Skipping\n'
+                                 'validation.')
+
         if self.mesh.with_ice_shelf_cavities:
             variables = ['ssh', 'landIcePressure']
             compare_variables(test_case=self, variables=variables,
diff --git a/compass/ocean/tests/global_ocean/init/init.cfg b/compass/ocean/tests/global_ocean/init/init.cfg
new file mode 100644
index 0000000000..3688d8510d
--- /dev/null
+++ b/compass/ocean/tests/global_ocean/init/init.cfg
@@ -0,0 +1,5 @@
+# Options related to the vertical grid
+[vertical_grid]
+
+# no inactive top cells by default
+inactive_top_cells = 0
diff --git a/compass/ocean/tests/global_ocean/init/initial_state.py b/compass/ocean/tests/global_ocean/init/initial_state.py
index d9ea96eedf..a47b0d3e05 100644
--- a/compass/ocean/tests/global_ocean/init/initial_state.py
+++ b/compass/ocean/tests/global_ocean/init/initial_state.py
@@ -1,7 +1,11 @@
 import os
 from importlib.resources import contents
 
+import xarray
+from mpas_tools.io import write_netcdf
+
 from compass.model import run_model
+from compass.ocean.inactive_top_cells import remove_inactive_top_cells_output
 from compass.ocean.plot import plot_initial_state, plot_vertical_grid
 from compass.ocean.tests.global_ocean.metadata import (
     add_mesh_and_init_metadata,
@@ -23,7 +27,8 @@ class InitialState(Step):
     initial_condition : {'WOA23', 'PHC', 'EN4_1900'}
         The initial condition dataset to use
     """
-    def __init__(self, test_case, mesh, initial_condition):
+    def __init__(self, test_case, mesh, initial_condition,
+                 with_inactive_top_cells):
         """
         Create the step
 
@@ -44,6 +49,7 @@ def __init__(self, test_case, mesh, initial_condition):
         super().__init__(test_case=test_case, name='initial_state')
         self.mesh = mesh
         self.initial_condition = initial_condition
+        self.with_inactive_top_cells = with_inactive_top_cells
 
         package = 'compass.ocean.tests.global_ocean.init'
 
@@ -134,6 +140,9 @@ def __init__(self, test_case, mesh, initial_condition):
                      'graph.info']:
             self.add_output_file(filename=file)
 
+        if with_inactive_top_cells:
+            self.add_output_file(filename='initial_state_crop.nc')
+
     def setup(self):
         """
         Get resources at setup from config options
@@ -162,6 +171,7 @@ def run(self):
         Run this step of the testcase
         """
         config = self.config
+        logger = self.logger
         interfaces = generate_1d_grid(config=config)
 
         write_1d_grid(interfaces=interfaces, out_filename='vertical_grid.nc')
@@ -171,6 +181,39 @@ def run(self):
         update_pio = config.getboolean('global_ocean', 'init_update_pio')
         run_model(self, update_pio=update_pio)
 
+        if self.with_inactive_top_cells:
+
+            logger.info("   * Updating minLevelCell for inactive top cells")
+
+            in_filename = 'initial_state.nc'
+            out_filename = in_filename
+
+            with xarray.open_dataset(in_filename) as ds:
+                ds.load()
+
+                # keep the data set with Time for output
+                ds_out = ds
+
+                ds = ds.isel(Time=0)
+
+                offset = config.getint('vertical_grid',
+                                       'inactive_top_cells')
+
+                if 'minLevelCell' in ds:
+                    minLevelCell = ds.minLevelCell + offset
+                    ds_out['minLevelCell'] = minLevelCell
+                else:
+                    logger.info("   - Variable minLevelCell, needed for "
+                                "inactive top cells, is missing from the "
+                                "initial condition")
+
+            write_netcdf(ds_out, out_filename)
+
+            remove_inactive_top_cells_output(
+                in_filename=in_filename, out_filename='initial_state_crop.nc')
+
+            logger.info("   - Complete")
+
         add_mesh_and_init_metadata(self.outputs, config,
                                    init_filename='initial_state.nc')
 
diff --git a/compass/ocean/tests/global_ocean/init/streams.init b/compass/ocean/tests/global_ocean/init/streams.init
index c105b6041e..aff8c8bc29 100644
--- a/compass/ocean/tests/global_ocean/init/streams.init
+++ b/compass/ocean/tests/global_ocean/init/streams.init
@@ -23,6 +23,7 @@
     <var name="bottomDepth"/>
     <var name="bottomDepthObserved"/>
     <var name="oceanFracObserved"/>
+    <var name="minLevelCell"/>
     <var name="maxLevelCell"/>
     <var name="vertCoordMovementWeights"/>
     <var name="edgeMask"/>
diff --git a/compass/ocean/tests/global_ocean/mesh/qu/icos.cfg b/compass/ocean/tests/global_ocean/mesh/qu/icos.cfg
index a1a448fea0..36774a8afc 100644
--- a/compass/ocean/tests/global_ocean/mesh/qu/icos.cfg
+++ b/compass/ocean/tests/global_ocean/mesh/qu/icos.cfg
@@ -9,3 +9,9 @@ prefix = Icos
 mesh_description = MPAS subdivided icosahedral mesh for E3SM version
                    ${e3sm_version} at ${min_res}-km global resolution with
                    <<<levels>>> vertical level
+
+# Options related to the vertical grid
+[vertical_grid]
+
+# Number of inactive top cell layers
+inactive_top_cells = 0
diff --git a/compass/ocean/tests/global_ocean/mesh/qu/qu.cfg b/compass/ocean/tests/global_ocean/mesh/qu/qu.cfg
index eb4a478d49..7b20022a59 100644
--- a/compass/ocean/tests/global_ocean/mesh/qu/qu.cfg
+++ b/compass/ocean/tests/global_ocean/mesh/qu/qu.cfg
@@ -20,6 +20,8 @@ max_layer_thickness = 250.0
 # the min and max occurs
 transition_levels = 28
 
+# Number of inactive top cell layers
+inactive_top_cells = 0
 
 # options for global ocean testcases
 [global_ocean]
diff --git a/compass/ocean/tests/global_ocean/mesh/qu240/qu240.cfg b/compass/ocean/tests/global_ocean/mesh/qu240/qu240.cfg
index de9bba14f9..4176455b76 100644
--- a/compass/ocean/tests/global_ocean/mesh/qu240/qu240.cfg
+++ b/compass/ocean/tests/global_ocean/mesh/qu240/qu240.cfg
@@ -16,6 +16,9 @@ min_layer_thickness = 3.0
 # The maximum layer thickness
 max_layer_thickness = 500.0
 
+# Number of inactive top cell layers
+inactive_top_cells = 0
+
 
 # options for spherical meshes
 [spherical_mesh]
diff --git a/compass/ocean/tests/global_ocean/performance_test/__init__.py b/compass/ocean/tests/global_ocean/performance_test/__init__.py
index 0487f4f9f3..990fe397d7 100644
--- a/compass/ocean/tests/global_ocean/performance_test/__init__.py
+++ b/compass/ocean/tests/global_ocean/performance_test/__init__.py
@@ -1,3 +1,6 @@
+import os
+
+from compass.ocean.inactive_top_cells import remove_inactive_top_cells_output
 from compass.ocean.tests.global_ocean.forward import (
     ForwardStep,
     ForwardTestCase,
@@ -31,9 +34,18 @@ def __init__(self, test_group, mesh, init, time_integrator):
         time_integrator : {'split_explicit', 'RK4'}
             The time integrator to use for the forward run
         """
+        name = 'performance_test'
+        if init.with_inactive_top_cells:
+            if time_integrator == 'RK4':
+                self.inactive_top_comp_subdir = os.path.join(
+                    mesh.mesh_name, init.initial_condition, time_integrator,
+                    name)
+            else:
+                self.inactive_top_comp_subdir = os.path.join(
+                    mesh.mesh_name, init.initial_condition, name)
+            name = f'{name}_inactive_top'
         super().__init__(test_group=test_group, mesh=mesh, init=init,
-                         time_integrator=time_integrator,
-                         name='performance_test')
+                         time_integrator=time_integrator, name=name)
 
         if mesh.with_ice_shelf_cavities:
             this_module = self.__module__
@@ -71,6 +83,24 @@ def validate(self):
             compare_variables(test_case=self, variables=variables,
                               filename1=f'{step_subdir}/output.nc')
 
+            if self.init.with_inactive_top_cells:
+                filename2 = os.path.join(self.base_work_dir,
+                                         self.mpas_core.name,
+                                         self.test_group.name,
+                                         self.inactive_top_comp_subdir,
+                                         'forward/output.nc')
+                if os.path.exists(filename2):
+                    compare_variables(
+                        test_case=self, variables=variables,
+                        filename1=f'{step_subdir}/output_crop.nc',
+                        filename2=filename2,
+                        quiet=False, check_outputs=False,
+                        skip_if_step_not_run=False)
+                else:
+                    self.logger.warn('The version of "performance_test" '
+                                     'without inactive top cells was not run. '
+                                     'Skipping validation.')
+
             if self.mesh.with_ice_shelf_cavities:
                 variables = [
                     'ssh', 'landIcePressure', 'landIceDraft',
diff --git a/compass/ocean/vertical/grid_1d/__init__.py b/compass/ocean/vertical/grid_1d/__init__.py
index fab57b5f23..5731da0c90 100644
--- a/compass/ocean/vertical/grid_1d/__init__.py
+++ b/compass/ocean/vertical/grid_1d/__init__.py
@@ -28,16 +28,21 @@ def generate_1d_grid(config):
         A 1D array of positive depths for layer interfaces in meters
     """
     section = config['vertical_grid']
+    offset = 0
+    if config.has_option('vertical_grid', 'inactive_top_cells'):
+        offset = section.getint('inactive_top_cells')
+        print(f'offset = {offset}')
+
     grid_type = section.get('grid_type')
     if grid_type == 'uniform':
         vert_levels = section.getint('vert_levels')
-        interfaces = _generate_uniform(vert_levels)
+        interfaces = _generate_uniform(vert_levels - offset)
     elif grid_type == 'tanh_dz':
         vert_levels = section.getint('vert_levels')
         min_layer_thickness = section.getfloat('min_layer_thickness')
         max_layer_thickness = section.getfloat('max_layer_thickness')
         bottom_depth = section.getfloat('bottom_depth')
-        interfaces = create_tanh_dz_grid(vert_levels,
+        interfaces = create_tanh_dz_grid(vert_levels - offset,
                                          bottom_depth,
                                          min_layer_thickness,
                                          max_layer_thickness)
@@ -66,6 +71,9 @@ def generate_1d_grid(config):
         # renormalize to the requested range
         interfaces = (bottom_depth / interfaces[-1]) * interfaces
 
+    if config.has_option('vertical_grid', 'inactive_top_cells'):
+        interfaces = np.append(np.zeros((offset)), interfaces)
+
     return interfaces
 
 
diff --git a/docs/developers_guide/ocean/test_groups/global_ocean.rst b/docs/developers_guide/ocean/test_groups/global_ocean.rst
index 91e28e8af4..5cc2646880 100644
--- a/docs/developers_guide/ocean/test_groups/global_ocean.rst
+++ b/docs/developers_guide/ocean/test_groups/global_ocean.rst
@@ -1013,6 +1013,11 @@ If a baseline is provided, prognostic variables and ice-shelf melt fluxes (if
 ice-shelf cavities are included in the mesh) are compared with a baseline, and
 the ``time integration`` timer is compared with that of the baseline.
 
+This test case can also serve as a verification of the inactive top cell
+implementation. If ``inactive_top_cells`` is given in the config file with
+values greater than zero, the vertical domain is shifted downward by the given
+number of layers.
+
 .. _dev_ocean_global_ocean_restart_test:
 
 restart_test test case