diff --git a/docs/developers_guide/ocean/api.md b/docs/developers_guide/ocean/api.md
index 57578c95a..5db67ddac 100644
--- a/docs/developers_guide/ocean/api.md
+++ b/docs/developers_guide/ocean/api.md
@@ -130,6 +130,37 @@
    viz.Viz.run
 ```
 
+### manufactured_solution
+
+```{eval-rst}
+.. currentmodule:: polaris.ocean.tests.manufactured_solution
+
+.. autosummary::
+   :toctree: generated/
+
+   ManufacturedSolution
+
+   analysis.Analysis
+   analysis.Analysis.run
+
+   exact_solution.ExactSolution
+   exact_solution.ExactSolution.ssh
+   exact_solution.ExactSolution.normalVelocity
+
+   forward.Forward
+   forward.Forward.compute_cell_count
+   forward.Forward.dynamic_model_config
+
+   initial_state.InitialState
+   initial_state.InitialState.run
+
+   viz.Viz
+   viz.Viz.run
+
+   convergence.Convergence
+   convergence.Convergence.validate
+```
+
 ### single_column
 
 ```{eval-rst}
diff --git a/docs/developers_guide/ocean/test_groups/index.md b/docs/developers_guide/ocean/test_groups/index.md
index 7d8835eb9..c8b755ee8 100644
--- a/docs/developers_guide/ocean/test_groups/index.md
+++ b/docs/developers_guide/ocean/test_groups/index.md
@@ -8,5 +8,6 @@
 baroclinic_channel
 global_convergence
 inertial_gravity_wave
+manufactured_solution
 single_column
 ```
diff --git a/docs/developers_guide/ocean/test_groups/manufactured_solution.md b/docs/developers_guide/ocean/test_groups/manufactured_solution.md
new file mode 100644
index 000000000..08c604003
--- /dev/null
+++ b/docs/developers_guide/ocean/test_groups/manufactured_solution.md
@@ -0,0 +1,82 @@
+(dev-ocean-manufactured_solution)=
+
+# manufactured_solution
+
+The `manufactured_solution` test group
+({py:class}`polaris.ocean.tests.manufactured_solution.ManufacturedSolution`)
+implements a test case according to the Method of Manufactured Solutions
+(see {ref}`ocean-manufactured-solution`) at 4 resolutions (200, 100, 50, and 25 km).  Here,
+we describe the shared framework for this test group and the 1 test case.
+
+(dev-ocean-manufactured-solution)=
+
+## framework
+
+The shared config options for the `manufactured_solution` test group
+are described in {ref}`ocean-manufactured-solution` in the User's Guide.
+
+Additionally, the test group has a shared `forward.yaml` file with
+a few common model config options related to run duration and default 
+horizontal  and vertical momentum and tracer diffusion, as well as defining 
+`mesh`, `input`, `restart`, and `output` streams.
+
+### exact_solution
+
+The class {py:class}`polaris.ocean.tests.manufactured_solution.exact_solution.ExactSolution`
+defines a class for storing attributes and methods relevant to computing the
+exact manufactured solution.  The constructor obtains the parameters from the
+config file. The
+{py:meth}`polaris.ocean.tests.manufactured_solution.exact_solution.ExactSolution.ssh(t)`
+method computes the SSH field at time `t`.  The
+{py:meth}`polaris.ocean.tests.manufactured_solution.exact_solution.ExactSolution.normalVelocity(t)`
+method computes the `normalVelocity` field at time `t`.
+
+### initial_state
+
+The class {py:class}`polaris.ocean.tests.manufactured_solution.initial_state.InitialState`
+defines a step for setting up the initial state for each test case.
+
+First, a mesh appropriate for the resolution is generated using
+{py:func}`mpas_tools.planar_hex.make_planar_hex_mesh()`.  Then, the mesh is
+culled to remove periodicity in the y direction.  A vertical grid is generated
+with 1 layer.  Finally, the initial layerThickness field is computed from the
+exact solution for the SSH field and the initial velocity is also assigned to
+the exact solution. The tracer and coriolis fields are uniform in space.
+
+### forward
+
+The class {py:class}`polaris.ocean.tests.manufactured_solution.forward.Forward`
+defines a step for running MPAS-Ocean from the initial condition produced in
+the `initial_state` step.  Namelist and streams files are updated in
+{py:meth}`polaris.ocean.tests.manufactured_solution.forward.Forward.dynamic_model_config()`
+with time steps determined algorithmically based on config options.  The
+number of cells is computed from config options in
+{py:meth}`polaris.ocean.tests.manufactured_solution.forward.Forward.compute_cell_count()`
+so that this can be used to constrain the number of MPI tasks that tests
+have as their target and minimum (if the resources are not explicitly
+prescribed).  For MPAS-Ocean, PIO namelist options are modified and a
+graph partition is generated as part of `runtime_setup()`.  Finally, the ocean 
+model is run.
+
+### analysis
+
+The class {py:class}`polaris.ocean.tests.manufactured_solution.analysis.Analysis`
+defines a step for computing the root mean-square-error from the final
+simulated field and the exact solution. It uses the config options to determine
+whether the convergence rate falls within acceptable bounds.
+
+### viz
+
+The class {py:class}`polaris.ocean.tests.manufactured_solution.viz.Viz`
+defines a step for visualization. It produces two plots: the convergence of the
+RMSE with resolution and a plan-view of the simulated, exact, and (simulated -
+exact) SSH fields.
+
+(dev-ocean-manufactured-solution-convergence)=
+
+### convergence
+
+The {py:class}`polaris.ocean.tests.manufactured_solution.convergence.Convergence`
+test performs a 10-hour run.  Then, validation of `temperature`, 
+`layerThickness` and `normalVelocity` are performed against a
+baseline if one is provided when calling {ref}`dev-polaris-setup`.
diff --git a/docs/users_guide/ocean/test_groups/index.md b/docs/users_guide/ocean/test_groups/index.md
index c0043e182..53f48960d 100644
--- a/docs/users_guide/ocean/test_groups/index.md
+++ b/docs/users_guide/ocean/test_groups/index.md
@@ -8,5 +8,6 @@
 baroclinic_channel
 global_convergence
 inertial_gravity_wave
+manufactured_solution
 single_column
 ```
diff --git a/docs/users_guide/ocean/test_groups/manufactured_solution.md b/docs/users_guide/ocean/test_groups/manufactured_solution.md
new file mode 100644
index 000000000..2038acb21
--- /dev/null
+++ b/docs/users_guide/ocean/test_groups/manufactured_solution.md
@@ -0,0 +1,123 @@
+(ocean-manufactured-solution)=
+
+# manufactured_solution
+
+The manufactured solution test group implements configurations for surface wave
+propagation with the rotating, nonlinear shallow water equations on a doubly
+periodic domain. This test group is intended to utilize tendency terms embedded
+in the forward ocean model in order to produce the manufactured solution. This
+solution can be then used to assess the numerical accuracy and convergence of
+the discretized nonlinear momentum equation. 
+
+Currently, the test group contains one test case, which is a convergence test
+from
+[Bishnu et al.(2023)](https://doi.org/10.22541/essoar.167100170.03833124/v1)
+
+(ocean-manufactured-solution-convergence)=
+
+## convergence
+
+### description
+
+The `convergence` test case runs the manufactured solution simulation for 4
+different resolutions: 200, 100, 50, and 25 km.
+ 
+The forward step for each resolution runs the simulation for 10 hours. The
+model is configured without vertical advection and mixing. No tracers are enabled
+and the pressure gradient used is the gradient of the sea surface height.
+Horizontal mixing and bottom friction are also neglected.
+
+The analysis step computes the root mean-square-error of the difference between
+the simulated SSH field and the exact solution at the end of the simulation. It
+also computes the convergence rate with resolution
+
+The visualization step produces two plots: the convergence of the RMSE with
+resolution and a plan-view of the simulated, exact, and (simulated-exact)
+SSH fields.
+
+### mesh
+
+For each resolution, the `initial_state` step generates and planar hexagonal
+mesh that is periodic in both the x and y directions.
+
+### vertical grid
+
+Since this test case is a shallow water case, the vertical grid is set to a
+single layer configuration.
+
+```cfg
+[vertical_grid]
+
+# The type of vertical grid
+grid_type = uniform
+
+# Number of vertical levels
+vert_levels = 1 
+
+# Depth of the bottom of the ocean
+bottom_depth = 1000.0
+
+# The type of vertical coordinate (e.g. z-level, z-star)
+coord_type = z-star
+
+# Whether to use "partial" or "full", or "None" to not alter the topography
+partial_cell_type = None
+
+# The minimum fraction of a layer for partial cells
+min_pc_fraction = 0.1
+```
+
+### initial conditions
+
+The initial conditions are set to the following:
+$$
+\eta = \eta_0 \sin(k_x x + k_y y - \omega t)\\
+u = \eta_0 \cos(k_x x + k_y y - \omega t)\\
+v = u
+$$
+
+### forcing
+
+N/A
+
+### time step and run duration
+
+The time step is determined by the config option ``dt_per_km`` according to the
+mesh resolution. The run duration is 10 hours.
+
+### config options
+
+The following config options are availiable for this case:
+
+```cfg
+[manufactured_solution]
+
+# the size of the domain in km in the x and y directions
+lx = 10000.0
+
+# the coriolis parameter
+coriolis_parameter = 1.0e-4
+
+# the amplitude of the sea surface height perturbation
+ssh_amplitude = 1.0
+
+# Number of wavelengths in x direction
+n_wavelengths_x = 2
+
+# Number of wavelengths in y direction
+n_wavelengths_y = 2
+
+# Time step per resolution (s/km), since dt is proportional to resolution
+dt_per_km = 3.0
+
+# Convergence threshold below which the test fails
+conv_thresh = 1.8
+
+# Convergence rate above which a warning is issued
+conv_max = 2.2
+```
+
+### cores
+
+The number of cores is determined according to the config options
+``max_cells_per_core`` and ``goal_cells_per_core``.
diff --git a/e3sm_submodules/E3SM-Project b/e3sm_submodules/E3SM-Project
index 0bd06172c..b9c1a8c5c 160000
--- a/e3sm_submodules/E3SM-Project
+++ b/e3sm_submodules/E3SM-Project
@@ -1 +1 @@
-Subproject commit 0bd06172c98e2a216da7522564123962ed190194
+Subproject commit b9c1a8c5cfdb5045995c4af05d9ebb7c6f736d32
diff --git a/polaris/ocean/__init__.py b/polaris/ocean/__init__.py
index 01fa6bf3e..37620335f 100644
--- a/polaris/ocean/__init__.py
+++ b/polaris/ocean/__init__.py
@@ -2,6 +2,7 @@
 from polaris.ocean.tests.baroclinic_channel import BaroclinicChannel
 from polaris.ocean.tests.global_convergence import GlobalConvergence
 from polaris.ocean.tests.inertial_gravity_wave import InertialGravityWave
+from polaris.ocean.tests.manufactured_solution import ManufacturedSolution
 from polaris.ocean.tests.single_column import SingleColumn
 
 
@@ -20,6 +21,7 @@ def __init__(self):
         self.add_test_group(BaroclinicChannel(component=self))
         self.add_test_group(GlobalConvergence(component=self))
         self.add_test_group(InertialGravityWave(component=self))
+        self.add_test_group(ManufacturedSolution(component=self))
         self.add_test_group(SingleColumn(component=self))
 
     def configure(self, config):
diff --git a/polaris/ocean/suites/nightly.txt b/polaris/ocean/suites/nightly.txt
index cf28bec6c..a555a2129 100644
--- a/polaris/ocean/suites/nightly.txt
+++ b/polaris/ocean/suites/nightly.txt
@@ -2,3 +2,4 @@ ocean/baroclinic_channel/10km/threads_test
 ocean/baroclinic_channel/10km/decomp_test
 ocean/baroclinic_channel/10km/restart_test
 ocean/inertial_gravity_wave/convergence
+ocean/manufactured_solution/convergence
diff --git a/polaris/ocean/suites/pr.txt b/polaris/ocean/suites/pr.txt
index d2a0da53a..19cacd304 100644
--- a/polaris/ocean/suites/pr.txt
+++ b/polaris/ocean/suites/pr.txt
@@ -4,3 +4,4 @@ ocean/baroclinic_channel/10km/restart_test
 ocean/inertial_gravity_wave/convergence
 ocean/single_column/960km/cvmix
 ocean/single_column/960km/ideal_age
+ocean/manufactured_solution/convergence
diff --git a/polaris/ocean/tests/manufactured_solution/__init__.py b/polaris/ocean/tests/manufactured_solution/__init__.py
new file mode 100644
index 000000000..e37de57e3
--- /dev/null
+++ b/polaris/ocean/tests/manufactured_solution/__init__.py
@@ -0,0 +1,16 @@
+from polaris import TestGroup
+from polaris.ocean.tests.manufactured_solution.convergence import Convergence
+
+
+class ManufacturedSolution(TestGroup):
+    """
+    A test group for manufactured solution test cases
+    """
+    def __init__(self, component):
+        """
+        component : polaris.ocean.Ocean
+            the ocean component that this test group belongs to
+        """
+        super().__init__(component=component, name='manufactured_solution')
+
+        self.add_test_case(Convergence(test_group=self))
diff --git a/polaris/ocean/tests/manufactured_solution/analysis.py b/polaris/ocean/tests/manufactured_solution/analysis.py
new file mode 100644
index 000000000..451e61563
--- /dev/null
+++ b/polaris/ocean/tests/manufactured_solution/analysis.py
@@ -0,0 +1,81 @@
+import datetime
+import warnings
+
+import cmocean  # noqa: F401
+import numpy as np
+import xarray as xr
+
+from polaris import Step
+from polaris.ocean.tests.manufactured_solution.exact_solution import (
+    ExactSolution,
+)
+
+
+class Analysis(Step):
+    """
+    A step for analysing the output from the manufactured solution
+    test case
+
+    Attributes
+    ----------
+    resolutions : list of int
+        The resolutions of the meshes that have been run
+    """
+    def __init__(self, test_case, resolutions):
+        """
+        Create the step
+
+        Parameters
+        ----------
+        test_case : polaris.TestCase
+            The test case this step belongs to
+
+        resolutions : list of int
+            The resolutions of the meshes that have been run
+        """
+        super().__init__(test_case=test_case, name='analysis')
+        self.resolutions = resolutions
+
+        for resolution in resolutions:
+            self.add_input_file(
+                filename=f'init_{resolution}km.nc',
+                target=f'../{resolution}km/initial_state/initial_state.nc')
+            self.add_input_file(
+                filename=f'output_{resolution}km.nc',
+                target=f'../{resolution}km/forward/output.nc')
+
+    def run(self):
+        """
+        Run this step of the test case
+        """
+        config = self.config
+        resolutions = self.resolutions
+
+        section = config['manufactured_solution']
+        conv_thresh = section.getfloat('conv_thresh')
+        conv_max = section.getfloat('conv_max')
+
+        rmse = []
+        for i, res in enumerate(resolutions):
+            init = xr.open_dataset(f'init_{res}km.nc')
+            ds = xr.open_dataset(f'output_{res}km.nc')
+            exact = ExactSolution(config, init)
+
+            t0 = datetime.datetime.strptime(ds.xtime.values[0].decode(),
+                                            '%Y-%m-%d_%H:%M:%S')
+            tf = datetime.datetime.strptime(ds.xtime.values[-1].decode(),
+                                            '%Y-%m-%d_%H:%M:%S')
+            t = (tf - t0).total_seconds()
+            ssh_model = ds.ssh.values[-1, :]
+            rmse.append(np.sqrt(np.mean((ssh_model - exact.ssh(t).values)**2)))
+
+        p = np.polyfit(np.log10(resolutions), np.log10(rmse), 1)
+        conv = p[0]
+
+        if conv < conv_thresh:
+            raise ValueError(f'order of convergence '
+                             f' {conv} < min tolerence {conv_thresh}')
+
+        if conv > conv_max:
+            warnings.warn(f'order of convergence '
+                          f'{conv} > max tolerence {conv_max}')
diff --git a/polaris/ocean/tests/manufactured_solution/convergence/__init__.py b/polaris/ocean/tests/manufactured_solution/convergence/__init__.py
new file mode 100644
index 000000000..79693f11f
--- /dev/null
+++ b/polaris/ocean/tests/manufactured_solution/convergence/__init__.py
@@ -0,0 +1,52 @@
+from polaris import TestCase
+from polaris.ocean.tests.manufactured_solution.analysis import Analysis
+from polaris.ocean.tests.manufactured_solution.forward import Forward
+from polaris.ocean.tests.manufactured_solution.initial_state import (
+    InitialState,
+)
+from polaris.ocean.tests.manufactured_solution.viz import Viz
+from polaris.validate import compare_variables
+
+
+class Convergence(TestCase):
+    """
+    The convergence test case for the manufactured solution test group
+
+    Attributes
+    ----------
+    resolutions : list of floats
+        The resolutions of the test case in km
+    """
+    def __init__(self, test_group):
+        """
+        Create the test case
+
+        Parameters
+        ----------
+        test_group : polaris.ocean.tests.manufactured_solution.
+                     ManufacturedSolution
+            The test group that this test case belongs to
+        """
+        name = 'convergence'
+        super().__init__(test_group=test_group, name=name)
+
+        self.resolutions = [200, 100, 50, 25]
+        for res in self.resolutions:
+            self.add_step(InitialState(test_case=self, resolution=res))
+            self.add_step(Forward(test_case=self, resolution=res))
+
+        self.add_step(Analysis(test_case=self, resolutions=self.resolutions))
+        self.add_step(Viz(test_case=self, resolutions=self.resolutions),
+                      run_by_default=False)
+
+    def validate(self):
+        """
+        Compare ``layerThickness`` and
+        ``normalVelocity`` in the ``forward`` step with a baseline if one was
+        provided.
+        """
+        super().validate()
+        variables = ['layerThickness', 'normalVelocity']
+        for res in self.resolutions:
+            compare_variables(test_case=self, variables=variables,
+                              filename1=f'{res}km/forward/output.nc')
diff --git a/polaris/ocean/tests/manufactured_solution/exact_solution.py b/polaris/ocean/tests/manufactured_solution/exact_solution.py
new file mode 100644
index 000000000..8167be66e
--- /dev/null
+++ b/polaris/ocean/tests/manufactured_solution/exact_solution.py
@@ -0,0 +1,113 @@
+import numpy as np
+
+
+class ExactSolution():
+    """
+    Class to compute the exact solution for the manufactured solution
+    test case
+
+    Attributes
+    ----------
+    angleEdge : xr.DataArray
+        angle between edge normals and positive x direction
+
+    xCell : xr.DataArray
+        x coordinates of mesh cell centers
+
+    yCell : xr.DataArray
+        y coordinates of mesh cell centers
+
+    xEdge: xr.DataArray
+        x coordinates of mesh edges
+
+    yEdge : xr.DataArray
+        y coordinates of mesh edges
+
+    eta0 : float
+        Amplitide of sea surface height
+
+    kx : float
+        Wave number in the x direction
+
+    ky : float
+        Wave number in the y direction
+
+    omega : float
+        Angular frequency
+    """
+
+    def __init__(self, config, ds=None):
+        """
+        Create a new exact solution object
+
+        Parameters
+        ----------
+        ds : xr.DataSet
+            MPAS mesh information
+
+        config : polaris.config.PolarisConfigParser
+            Config options for test case
+        """
+        bottom_depth = config.getfloat('vertical_grid', 'bottom_depth')
+        section = config['manufactured_solution']
+        self.g = 9.80616
+        self.eta0 = section.getfloat('ssh_amplitude')
+        lx = section.getfloat('lx')
+        ly = np.sqrt(3.0) / 2.0 * lx
+        npx = section.getfloat('n_wavelengths_x')
+        npy = section.getfloat('n_wavelengths_y')
+        self.lambda_x = (lx * 1e3) / npx
+        self.lambda_y = (ly * 1e3) / npy
+        self.kx = 2.0 * np.pi / self.lambda_x
+        self.ky = 2.0 * np.pi / self.lambda_y
+        self.omega = np.sqrt(self.g * bottom_depth *
+                             (self.kx**2 + self.ky**2))
+
+        if ds is not None:
+            self.angleEdge = ds.angleEdge
+            self.xCell = ds.xCell
+            self.yCell = ds.yCell
+            self.xEdge = ds.xEdge
+            self.yEdge = ds.yEdge
+
+    def ssh(self, t):
+        """
+        Exact solution for sea surface height
+
+        Parameters
+        ----------
+        t : float
+            time at which to evaluate exact solution
+
+        Returns
+        -------
+        eta : xr.DataArray
+            the exact sea surface height solution on cells at time t
+
+        """
+
+        eta = self.eta0 * np.sin(self.kx * self.xCell +
+                                 self.ky * self.yCell -
+                                 self.omega * t)
+
+        return eta
+
+    def normal_velocity(self, t):
+        """
+        Exact solution for normal velocity
+
+        Parameters
+        ----------
+        t : float
+            time at which to evaluate exact solution
+
+        Returns
+        -------
+        normalvelocity : xr.DataArray
+            the exact normal velocity solution on edges at time t
+        """
+        u = self.eta0 * np.cos(self.kx * self.xEdge +
+                               self.ky * self.yEdge -
+                               self.omega * t)
+        v = u
+        return u * np.cos(self.angleEdge) + v * np.sin(self.angleEdge)
diff --git a/polaris/ocean/tests/manufactured_solution/forward.py b/polaris/ocean/tests/manufactured_solution/forward.py
new file mode 100644
index 000000000..feacba587
--- /dev/null
+++ b/polaris/ocean/tests/manufactured_solution/forward.py
@@ -0,0 +1,123 @@
+import time
+
+import numpy as np
+import xarray as xr
+
+from polaris.mesh.planar import compute_planar_hex_nx_ny
+from polaris.ocean.model import OceanModelStep
+from polaris.ocean.tests.manufactured_solution.exact_solution import (
+    ExactSolution,
+)
+
+
+class Forward(OceanModelStep):
+    """
+    A step for performing forward ocean component runs as part of manufactured
+    solution test cases.
+
+    Attributes
+    ----------
+    resolution : float
+        The resolution of the test case in km
+    """
+    def __init__(self, test_case, resolution,
+                 ntasks=None, min_tasks=None, openmp_threads=1):
+        """
+        Create a new test case
+
+        Parameters
+        ----------
+        test_case : polaris.TestCase
+            The test case this step belongs to
+
+        resolution : km
+            The resolution of the test case in km
+
+        ntasks : int, optional
+            the number of tasks the step would ideally use.  If fewer tasks
+            are available on the system, the step will run on all available
+            tasks as long as this is not below ``min_tasks``
+
+        min_tasks : int, optional
+            the number of tasks the step requires.  If the system has fewer
+            than this number of tasks, the step will fail
+
+        openmp_threads : int, optional
+            the number of OpenMP threads the step will use
+        """
+        self.resolution = resolution
+        super().__init__(test_case=test_case,
+                         name=f'forward_{resolution}km',
+                         subdir=f'{resolution}km/forward',
+                         ntasks=ntasks, min_tasks=min_tasks,
+                         openmp_threads=openmp_threads)
+
+        self.add_input_file(filename='initial_state.nc',
+                            target='../initial_state/initial_state.nc')
+        self.add_input_file(filename='graph.info',
+                            target='../initial_state/culled_graph.info')
+
+        self.add_output_file(filename='output.nc')
+
+        self.add_yaml_file('polaris.ocean.config',
+                           'single_layer.yaml')
+        self.add_yaml_file('polaris.ocean.tests.manufactured_solution',
+                           'forward.yaml')
+
+    def compute_cell_count(self, at_setup):
+        """
+        Compute the approximate number of cells in the mesh, used to constrain
+        resources
+
+        Parameters
+        ----------
+        at_setup : bool
+            Whether this method is being run during setup of the step, as
+            opposed to at runtime
+
+        Returns
+        -------
+        cell_count : int or None
+            The approximate number of cells in the mesh
+        """
+        if at_setup:
+            # no file to read from, so we'll compute it based on config options
+            section = self.config['manufactured_solution']
+            lx = section.getfloat('lx')
+            ly = np.sqrt(3.0) / 2.0 * lx
+            nx, ny = compute_planar_hex_nx_ny(lx, ly, self.resolution)
+            cell_count = nx * ny
+        else:
+            # get nCells from the input file
+            with xr.open_dataset('initial_state.nc') as ds:
+                cell_count = ds.sizes['nCells']
+        return cell_count
+
+    def dynamic_model_config(self, at_setup):
+        """
+        Set the model time step from config options at setup and runtime
+
+        Parameters
+        ----------
+        at_setup : bool
+            Whether this method is being run during setup of the step, as
+            opposed to at runtime
+        """
+        super().dynamic_model_config(at_setup=at_setup)
+
+        # dt is proportional to resolution
+        config = self.config
+        section = config['manufactured_solution']
+        dt_per_km = section.getfloat('dt_per_km')
+        dt = dt_per_km * self.resolution
+        # https://stackoverflow.com/a/1384565/7728169
+        dt_str = time.strftime('%H:%M:%S', time.gmtime(dt))
+        exact_solution = ExactSolution(config)
+        options = {'config_dt': dt_str,
+                   'config_manufactured_solution_amplitude':
+                   exact_solution.eta0,
+                   'config_manufactured_solution_wavelength_x':
+                   exact_solution.lambda_x,
+                   'config_manufactured_solution_wavelength_y':
+                   exact_solution.lambda_y}
+        self.add_model_config_options(options)
diff --git a/polaris/ocean/tests/manufactured_solution/forward.yaml b/polaris/ocean/tests/manufactured_solution/forward.yaml
new file mode 100644
index 000000000..6b3502fce
--- /dev/null
+++ b/polaris/ocean/tests/manufactured_solution/forward.yaml
@@ -0,0 +1,29 @@
+omega:
+  time_management:
+    config_run_duration: 10:00:00
+  time_integration:
+    config_time_integrator: RK4
+  bottom_drag:
+    config_bottom_drag_mode: implicit
+    config_implicit_bottom_drag_type: constant
+    config_implicit_constant_bottom_drag_coeff: 0.0
+  manufactured_solution:
+     config_use_manufactured_solution: true
+  debug:
+    config_disable_vel_hmix: true
+  streams:
+    mesh:
+      filename_template: initial_state.nc
+    input:
+      filename_template: initial_state.nc
+    restart: {}
+    output:
+      type: output
+      filename_template: output.nc
+      output_interval: 0000_00:20:00
+      clobber_mode: truncate
+      contents:
+      - xtime
+      - normalVelocity
+      - layerThickness
+      - ssh
diff --git a/polaris/ocean/tests/manufactured_solution/initial_state.py b/polaris/ocean/tests/manufactured_solution/initial_state.py
new file mode 100644
index 000000000..29da385cf
--- /dev/null
+++ b/polaris/ocean/tests/manufactured_solution/initial_state.py
@@ -0,0 +1,103 @@
+import numpy as np
+import xarray as xr
+from mpas_tools.io import write_netcdf
+from mpas_tools.mesh.conversion import convert, cull
+from mpas_tools.planar_hex import make_planar_hex_mesh
+
+from polaris import Step
+from polaris.mesh.planar import compute_planar_hex_nx_ny
+from polaris.ocean.tests.manufactured_solution.exact_solution import (
+    ExactSolution,
+)
+from polaris.ocean.vertical import init_vertical_coord
+
+
+class InitialState(Step):
+    """
+    A step for creating a mesh and initial condition for the
+    manufactured solution test cases
+
+    Attributes
+    ----------
+    resolution : float
+        The resolution of the test case in km
+    """
+    def __init__(self, test_case, resolution):
+        """
+        Create the step
+
+        Parameters
+        ----------
+        test_case : polaris.TestCase
+            The test case this step belongs to
+
+        resolution : float
+            The resolution of the test case in km
+        """
+        super().__init__(test_case=test_case,
+                         name=f'initial_state_{resolution}km',
+                         subdir=f'{resolution}km/initial_state')
+        self.resolution = float(resolution)
+
+    def run(self):
+        """
+        Run this step of the test case
+        """
+        logger = self.logger
+        config = self.config
+
+        section = config['manufactured_solution']
+        resolution = self.resolution
+
+        lx = section.getfloat('lx')
+        ly = np.sqrt(3.0) / 2.0 * lx
+        coriolis_parameter = section.getfloat('coriolis_parameter')
+
+        nx, ny = compute_planar_hex_nx_ny(lx, ly, resolution)
+        dc = 1e3 * resolution
+
+        ds_mesh = make_planar_hex_mesh(nx=nx, ny=ny, dc=dc,
+                                       nonperiodic_x=False,
+                                       nonperiodic_y=False)
+        write_netcdf(ds_mesh, 'base_mesh.nc')
+
+        ds_mesh = cull(ds_mesh, logger=logger)
+        ds_mesh = convert(ds_mesh, graphInfoFileName='culled_graph.info',
+                          logger=logger)
+        write_netcdf(ds_mesh, 'culled_mesh.nc')
+
+        bottom_depth = config.getfloat('vertical_grid', 'bottom_depth')
+
+        ds = ds_mesh.copy()
+
+        ds['ssh'] = xr.zeros_like(ds.xCell)
+        ds['bottomDepth'] = bottom_depth * xr.ones_like(ds.xCell)
+
+        init_vertical_coord(config, ds)
+
+        ds['fCell'] = coriolis_parameter * xr.ones_like(ds.xCell)
+        ds['fEdge'] = coriolis_parameter * xr.ones_like(ds.xEdge)
+        ds['fVertex'] = coriolis_parameter * xr.ones_like(ds.xVertex)
+
+        ds_mesh['maxLevelCell'] = ds.maxLevelCell
+
+        # Evaluate the exact solution at time=0
+        exact_solution = ExactSolution(config, ds)
+        ssh = exact_solution.ssh(0.0)
+        normal_velocity = exact_solution.normal_velocity(0.0)
+
+        ssh = ssh.expand_dims(dim='Time', axis=0)
+        ds['ssh'] = ssh
+
+        normal_velocity, _ = xr.broadcast(normal_velocity, ds.refBottomDepth)
+        normal_velocity = normal_velocity.transpose('nEdges', 'nVertLevels')
+        normal_velocity = normal_velocity.expand_dims(dim='Time', axis=0)
+        ds['normalVelocity'] = normal_velocity
+
+        layer_thickness = ssh + bottom_depth
+        layer_thickness, _ = xr.broadcast(layer_thickness, ds.refBottomDepth)
+        layer_thickness = layer_thickness.transpose('Time', 'nCells',
+                                                    'nVertLevels')
+        ds['layerThickness'] = layer_thickness
+
+        write_netcdf(ds, 'initial_state.nc')
diff --git a/polaris/ocean/tests/manufactured_solution/manufactured_solution.cfg b/polaris/ocean/tests/manufactured_solution/manufactured_solution.cfg
new file mode 100644
index 000000000..14d42ca80
--- /dev/null
+++ b/polaris/ocean/tests/manufactured_solution/manufactured_solution.cfg
@@ -0,0 +1,54 @@
+[ocean]
+# the number of cells per core to aim for
+goal_cells_per_core = 200
+
+# the approximate maximum number of cells per core (the test will fail if too
+# few cores are available)
+max_cells_per_core = 4500
+
+# config options for manufactured solution testcases
+[manufactured_solution]
+
+# the size of the domain in km in the x and y directions
+lx = 10000.0
+
+# the coriolis parameter
+coriolis_parameter = 1.0e-4
+
+# the amplitude of the sea surface height perturbation
+ssh_amplitude = 1.0
+
+# Number of wavelengths in x direction
+n_wavelengths_x = 2
+
+# Number of wavelengths in y direction
+n_wavelengths_y = 2
+
+# Time step per resolution (s/km), since dt is proportional to resolution
+dt_per_km = 3.0
+
+# Convergence threshold below which the test fails
+conv_thresh = 1.8
+
+# Convergence rate above which a warning is issued
+conv_max = 2.2
+
+[vertical_grid]
+
+# The type of vertical grid
+grid_type = uniform
+
+# Number of vertical levels
+vert_levels = 1
+
+# Depth of the bottom of the ocean
+bottom_depth = 1000.0
+
+# The type of vertical coordinate (e.g. z-level, z-star)
+coord_type = z-star
+
+# Whether to use "partial" or "full", or "None" to not alter the topography
+partial_cell_type = None
+
+# The minimum fraction of a layer for partial cells
+min_pc_fraction = 0.1
diff --git a/polaris/ocean/tests/manufactured_solution/viz.py b/polaris/ocean/tests/manufactured_solution/viz.py
new file mode 100644
index 000000000..303934ed7
--- /dev/null
+++ b/polaris/ocean/tests/manufactured_solution/viz.py
@@ -0,0 +1,125 @@
+import datetime
+
+import cmocean  # noqa: F401
+import matplotlib.pyplot as plt
+import numpy as np
+import xarray as xr
+
+from polaris import Step
+from polaris.ocean.tests.manufactured_solution.exact_solution import (
+    ExactSolution,
+)
+from polaris.viz import plot_horiz_field
+
+
+class Viz(Step):
+    """
+    A step for visualizing the output from the manufactured solution
+    test case
+
+    Attributes
+    ----------
+    resolutions : list of int
+        The resolutions of the meshes that have been run
+    """
+    def __init__(self, test_case, resolutions):
+        """
+        Create the step
+
+        Parameters
+        ----------
+        test_case : polaris.ocean.tests.manufactured_solution.convergence.Convergence # noqa: E501
+            The test case this step belongs to
+
+        resolutions : list of int
+            The resolutions of the meshes that have been run
+        """
+        super().__init__(test_case=test_case, name='viz')
+        self.resolutions = resolutions
+
+        for resolution in resolutions:
+            self.add_input_file(
+                filename=f'init_{resolution}km.nc',
+                target=f'../{resolution}km/initial_state/initial_state.nc')
+            self.add_input_file(
+                filename=f'output_{resolution}km.nc',
+                target=f'../{resolution}km/forward/output.nc')
+
+        self.add_output_file('convergence.png')
+
+    def run(self):
+        """
+        Run this step of the test case
+        """
+        plt.switch_backend('Agg')
+        config = self.config
+        resolutions = self.resolutions
+        nres = len(resolutions)
+
+        section = config['manufactured_solution']
+        eta0 = section.getfloat('ssh_amplitude')
+
+        fig, axes = plt.subplots(nrows=nres, ncols=3, figsize=(12, 2 * nres))
+        rmse = []
+        for i, res in enumerate(resolutions):
+            init = xr.open_dataset(f'init_{res}km.nc')
+            ds = xr.open_dataset(f'output_{res}km.nc')
+            exact = ExactSolution(config, init)
+
+            t0 = datetime.datetime.strptime(ds.xtime.values[0].decode(),
+                                            '%Y-%m-%d_%H:%M:%S')
+            tf = datetime.datetime.strptime(ds.xtime.values[-1].decode(),
+                                            '%Y-%m-%d_%H:%M:%S')
+            t = (tf - t0).total_seconds()
+            ssh_model = ds.ssh.values[-1, :]
+            rmse.append(np.sqrt(np.mean((ssh_model - exact.ssh(t).values)**2)))
+
+            # Comparison plots
+            ds['ssh_exact'] = exact.ssh(t)
+            ds['ssh_error'] = ssh_model - exact.ssh(t)
+            if i == 0:
+                error_range = np.max(np.abs(ds.ssh_error.values))
+
+            plot_horiz_field(ds, init, 'ssh', ax=axes[i, 0],
+                             cmap='cmo.balance', t_index=ds.sizes["Time"] - 1,
+                             vmin=-eta0, vmax=eta0, cmap_title="SSH")
+            plot_horiz_field(ds, init, 'ssh_exact', ax=axes[i, 1],
+                             cmap='cmo.balance',
+                             vmin=-eta0, vmax=eta0, cmap_title="SSH")
+            plot_horiz_field(ds, init, 'ssh_error', ax=axes[i, 2],
+                             cmap='cmo.balance', cmap_title="dSSH",
+                             vmin=-error_range, vmax=error_range)
+
+        axes[0, 0].set_title('Numerical solution')
+        axes[0, 1].set_title('Analytical solution')
+        axes[0, 2].set_title('Error (Numerical - Analytical)')
+
+        pad = 5
+        for ax, res in zip(axes[:, 0], resolutions):
+            ax.annotate(f'{res}km', xy=(0, 0.5),
+                        xytext=(-ax.yaxis.labelpad - pad, 0),
+                        xycoords=ax.yaxis.label, textcoords='offset points',
+                        size='large', ha='right', va='center')
+
+        fig.savefig('comparison.png', bbox_inches='tight', pad_inches=0.1)
+
+        # Convergence polts
+        fig = plt.figure()
+        ax = fig.add_subplot(111)
+        p = np.polyfit(np.log10(resolutions), np.log10(rmse), 1)
+        conv = np.round(p[0], 3)
+        ax.loglog(resolutions, rmse, '-ok', label=f'numerical (order={conv})')
+
+        c = rmse[0] * 1.5 / resolutions[0]
+        order1 = c * np.power(resolutions, 1)
+        c = rmse[0] * 1.5 / resolutions[0]**2
+        order2 = c * np.power(resolutions, 2)
+
+        ax.loglog(resolutions, order1, '--k', label='first order', alpha=0.3)
+        ax.loglog(resolutions, order2, 'k', label='second order', alpha=0.3)
+        ax.set_xlabel('resolution (km)')
+        ax.set_ylabel('RMS error (m)')
+        ax.set_title('Error Convergence')
+        ax.invert_xaxis()
+        ax.legend(loc='lower left')
+        fig.savefig('convergence.png', bbox_inches='tight', pad_inches=0.1)
diff --git a/polaris/run/serial.py b/polaris/run/serial.py
index 573de128b..6d50eac04 100644
--- a/polaris/run/serial.py
+++ b/polaris/run/serial.py
@@ -424,6 +424,7 @@ def _run_step(test_case, step, new_log_file, available_resources,
 
     missing_files = list()
     for input_file in step.inputs:
+        print(f'DEBUG: input file: {input_file}')
         if not os.path.exists(input_file):
             missing_files.append(input_file)
 
diff --git a/polaris/viz/polaris.mplstyle b/polaris/viz/polaris.mplstyle
index 6b533a92d..d87e4b201 100644
--- a/polaris/viz/polaris.mplstyle
+++ b/polaris/viz/polaris.mplstyle
@@ -5,7 +5,7 @@ figure.labelsize: 14
 
 legend.fontsize: 14
 
-axes.titlesize: 16
+axes.titlesize: 14
 axes.labelsize: 14
 
 patch.linewidth: 0.5