diff --git a/docs/developers_guide/ocean/api.md b/docs/developers_guide/ocean/api.md
index 675973bf7..43d26c0ef 100644
--- a/docs/developers_guide/ocean/api.md
+++ b/docs/developers_guide/ocean/api.md
@@ -82,6 +82,40 @@
    viz.Viz.run
 ```
 
+### sphere_transport
+
+```{eval-rst}
+.. currentmodule:: polaris.ocean.tasks.sphere_transport
+
+.. autosummary::
+   :toctree: generated/
+
+   add_sphere_transport_tasks
+
+   SphereTransport
+   SphereTransport.configure
+
+   init.Init
+   init.Init.run
+
+   forward.Forward
+
+   analysis.Analysis
+   analysis.Analysis.convergence_parameters
+
+   mixing_analysis.MixingAnalysis
+   mixing_analysis.MixingAnalysis.run
+
+   filament_analysis.FilamentAnalysis
+   filament_analysis.FilamentAnalysis.run
+
+   viz.VizMap
+   viz.VizMap.runtime_setup
+
+   viz.Viz
+   viz.Viz.run
+```
+
 ### cosine_bell
 
 ```{eval-rst}
diff --git a/docs/developers_guide/ocean/tasks/correlated_tracers_2d.md b/docs/developers_guide/ocean/tasks/correlated_tracers_2d.md
new file mode 100644
index 000000000..af67f1fcd
--- /dev/null
+++ b/docs/developers_guide/ocean/tasks/correlated_tracers_2d.md
@@ -0,0 +1,108 @@
+(dev-ocean-correlated-tracers-2d)=
+
+# correlated_tracers_2d
+
+The {py:class}`polaris.ocean.tasks.sphere_transport.SphereTransport`
+`correlated_tracers_2d` test performs a 12-day run on the sphere that has a periodic
+deforming flow which affects tracer distributions. The resolution of the
+sphere varies (by default, between 60 and 240 km). After one period, the
+tracer distributions are compared the initial condition to evaluate numerical
+errors associated with the horizontal advection scheme and determine the rate
+of convergence.
+
+## framework
+
+The config options for the `correlated_tracers_2d` test is described in 
+{ref}`ocean-correlated_tracers-2d` in the User's Guide.
+
+Additionally, the test uses a `forward.yaml` file with a few common
+model config options related to drag and default horizontal and
+vertical momentum and tracer diffusion, as well as defining `mesh`, `input`,
+`restart`, and `output` streams.  This file has Jinja templating that is
+used to update model config options based on Polaris config options, see
+{ref}`dev-ocean-spherical-convergence`.
+
+### base_mesh
+
+Sphere transport tasks use shared `base_mesh` steps for creating
+{ref}`dev-ocean-spherical-meshes` at a sequence of resolutions.
+
+### init
+
+The class {py:class}`polaris.ocean.tasks.sphere_transport.init.Init`
+defines a step for setting up the initial state at each resolution. The
+initial state is differentiated between `correlated_tracers_2d` and the other tests
+set up by the class.
+
+### forward
+
+The class {py:class}`polaris.ocean.tasks.sphere_transport.forward.Forward`
+descends from {py:class}`polaris.ocean.convergence.spherical.SphericalConvergenceForward`,
+and defines a step for running MPAS-Ocean from an initial condition produced in
+an `init` step. See {ref}`dev-ocean-convergence` for some relevant
+discussion of the parent class. The time step is determined from the resolution
+based on the `dt_per_km` config option in the `[convergence_forward]` 
+section.  Other model config options are taken from `forward.yaml`.
+
+### analysis
+
+The class {py:class}`polaris.ocean.tasks.sphere_transport.analysis.Analysis`
+descends from
+{py:class}`polaris.ocean.convergence.ConvergenceAnalysis`,
+and defines a step for computing the error norm (L2) for the results
+at each resolution for tracers and layer thickness, saving them in
+`convergence_*.csv` and plotting them in `convergence_*.png`.
+
+### mixing_analysis
+
+The class {py:class}`polaris.ocean.tasks.sphere_transport.mixing_analysis.MixingAnalysis`
+plots the relationship between two correlated tracers for each resolution in
+`triplots.png`.
+
+### viz
+
+Visualization steps are available only in the `correlated_tracers_2d/with_viz`
+tasks.  They are not included in the `correlated_tracers_2d` in order to keep regression
+as fast as possible when visualization isn't needed.
+
+The class {py:class}`polaris.ocean.tasks.sphere_transport.viz.VizMap`
+defines a step for creating a mapping file from the MPAS mesh at a given
+resolution to a lon-lat grid at a resolution and interpolation method 
+determined by config options.
+
+```cfg
+# options for visualization for the cosine bell convergence test case
+[sphere_transport_viz]
+
+# visualization latitude and longitude resolution
+dlon = 0.5
+dlat = 0.5
+
+# remapping method ('bilinear', 'neareststod', 'conserve')
+remap_method = conserve
+```
+
+The class {py:class}`polaris.ocean.tasks.sphere_transport.viz.Viz`
+is a step for plotting the initial and final states of the advection test for
+each resolution, mapped to the common lat-lon grid.  The colormap is controlled
+by these options:
+
+```cfg
+# options for visualization for the cosine bell convergence test case
+[sphere_transport_viz_*]
+
+# colormap options
+# colormap
+colormap_name = viridis
+
+# the type of norm used in the colormap
+norm_type = linear
+
+# A dictionary with keywords for the norm
+norm_args = {'vmin': 0., 'vmax': 1.}
+
+# We could provide colorbar tick marks but we'll leave the defaults
+# colorbar_ticks = np.linspace(0., 1., 9)
+```
+
+See {ref}`dev-visualization-global` for more details.
diff --git a/docs/developers_guide/ocean/tasks/divergent_2d.md b/docs/developers_guide/ocean/tasks/divergent_2d.md
new file mode 100644
index 000000000..65a27c1ee
--- /dev/null
+++ b/docs/developers_guide/ocean/tasks/divergent_2d.md
@@ -0,0 +1,102 @@
+(dev-ocean-divergent-2d)=
+
+# divergent_2d
+
+The {py:class}`polaris.ocean.tasks.sphere_transport.SphereTransport`
+`divergent_2d` test performs a 12-day run on the sphere that has a periodic
+divergent flow which affects tracer distributions. The resolution of the
+sphere varies (by default, between 60 and 240 km). After one period, the
+tracer distributions are compared the initial condition to evaluate numerical
+errors associated with the horizontal advection scheme and determine the rate
+of convergence.
+
+## framework
+
+The config options for the `divergent_2d` test is described in 
+{ref}`ocean-divergent-2d` in the User's Guide.
+
+Additionally, the test uses a `forward.yaml` file with a few common
+model config options related to drag and default horizontal and
+vertical momentum and tracer diffusion, as well as defining `mesh`, `input`,
+`restart`, and `output` streams.  This file has Jinja templating that is
+used to update model config options based on Polaris config options, see
+{ref}`dev-ocean-spherical-convergence`.
+
+### base_mesh
+
+Sphere transport tasks use shared `base_mesh` steps for creating
+{ref}`dev-ocean-spherical-meshes` at a sequence of resolutions.
+
+### init
+
+The class {py:class}`polaris.ocean.tasks.sphere_transport.init.Init`
+defines a step for setting up the initial state at each resolution. The
+initial state is differentiated between `divergent_2d` and the other tests
+set up by the class.
+
+### forward
+
+The class {py:class}`polaris.ocean.tasks.sphere_transport.forward.Forward`
+descends from {py:class}`polaris.ocean.convergence.spherical.SphericalConvergenceForward`,
+and defines a step for running MPAS-Ocean from an initial condition produced in
+an `init` step. See {ref}`dev-ocean-convergence` for some relevant
+discussion of the parent class. The time step is determined from the resolution
+based on the `dt_per_km` config option in the `[convergence_forward]` 
+section.  Other model config options are taken from `forward.yaml`.
+
+### analysis
+
+The class {py:class}`polaris.ocean.tasks.sphere_transport.analysis.Analysis`
+descends from
+{py:class}`polaris.ocean.convergence.ConvergenceAnalysis`,
+and defines a step for computing the error norm (L2) for the results
+at each resolution for tracers and layer thickness, saving them in
+`convergence_*.csv` and plotting them in `convergence_*.png`.
+
+### viz
+
+Visualization steps are available only in the `divergent_2d/with_viz`
+tasks.  They are not included in the `divergent_2d` in order to keep regression
+as fast as possible when visualization isn't needed.
+
+The class {py:class}`polaris.ocean.tasks.sphere_transport.viz.VizMap`
+defines a step for creating a mapping file from the MPAS mesh at a given
+resolution to a lon-lat grid at a resolution and interpolation method 
+determined by config options.
+
+```cfg
+# options for visualization for the cosine bell convergence test case
+[sphere_transport_viz]
+
+# visualization latitude and longitude resolution
+dlon = 0.5
+dlat = 0.5
+
+# remapping method ('bilinear', 'neareststod', 'conserve')
+remap_method = conserve
+```
+
+The class {py:class}`polaris.ocean.tasks.sphere_transport.viz.Viz`
+is a step for plotting the initial and final states of the advection test for
+each resolution, mapped to the common lat-lon grid.  The colormap is controlled
+by these options:
+
+```cfg
+# options for visualization for the cosine bell convergence test case
+[sphere_transport_viz_*]
+
+# colormap options
+# colormap
+colormap_name = viridis
+
+# the type of norm used in the colormap
+norm_type = linear
+
+# A dictionary with keywords for the norm
+norm_args = {'vmin': 0., 'vmax': 1.}
+
+# We could provide colorbar tick marks but we'll leave the defaults
+# colorbar_ticks = np.linspace(0., 1., 9)
+```
+
+See {ref}`dev-visualization-global` for more details.
diff --git a/docs/developers_guide/ocean/tasks/index.md b/docs/developers_guide/ocean/tasks/index.md
index 268e21c92..e2c45b4aa 100644
--- a/docs/developers_guide/ocean/tasks/index.md
+++ b/docs/developers_guide/ocean/tasks/index.md
@@ -6,9 +6,13 @@
 :titlesonly: true
 
 baroclinic_channel
+correlated_tracers_2d
 cosine_bell
 geostrophic
+divergent_2d
 inertial_gravity_wave
 manufactured_solution
+nondivergent_2d
+rotation_2d
 single_column
 ```
diff --git a/docs/developers_guide/ocean/tasks/nondivergent_2d.md b/docs/developers_guide/ocean/tasks/nondivergent_2d.md
new file mode 100644
index 000000000..9eb2fbc3f
--- /dev/null
+++ b/docs/developers_guide/ocean/tasks/nondivergent_2d.md
@@ -0,0 +1,108 @@
+(dev-ocean-nondivergent-2d)=
+
+# nondivergent_2d
+
+The {py:class}`polaris.ocean.tasks.sphere_transport.SphereTransport`
+`nondivergent_2d` test performs a 12-day run on the sphere that has a periodic
+deforming flow which affects tracer distributions. The resolution of the
+sphere varies (by default, between 60 and 240 km). After one period, the
+tracer distributions are compared the initial condition to evaluate numerical
+errors associated with the horizontal advection scheme and determine the rate
+of convergence.
+
+## framework
+
+The config options for the `nondivergent_2d` test is described in 
+{ref}`ocean-nondivergent-2d` in the User's Guide.
+
+Additionally, the test uses a `forward.yaml` file with a few common
+model config options related to drag and default horizontal and
+vertical momentum and tracer diffusion, as well as defining `mesh`, `input`,
+`restart`, and `output` streams.  This file has Jinja templating that is
+used to update model config options based on Polaris config options, see
+{ref}`dev-ocean-spherical-convergence`.
+
+### base_mesh
+
+Sphere transport tasks use shared `base_mesh` steps for creating
+{ref}`dev-ocean-spherical-meshes` at a sequence of resolutions.
+
+### init
+
+The class {py:class}`polaris.ocean.tasks.sphere_transport.init.Init`
+defines a step for setting up the initial state at each resolution. The
+initial state is differentiated between `nondivergent_2d` and the other tests
+set up by the class.
+
+### forward
+
+The class {py:class}`polaris.ocean.tasks.sphere_transport.forward.Forward`
+descends from {py:class}`polaris.ocean.convergence.spherical.SphericalConvergenceForward`,
+and defines a step for running MPAS-Ocean from an initial condition produced in
+an `init` step. See {ref}`dev-ocean-convergence` for some relevant
+discussion of the parent class. The time step is determined from the resolution
+based on the `dt_per_km` config option in the `[convergence_forward]` 
+section.  Other model config options are taken from `forward.yaml`.
+
+### analysis
+
+The class {py:class}`polaris.ocean.tasks.sphere_transport.analysis.Analysis`
+descends from
+{py:class}`polaris.ocean.convergence.ConvergenceAnalysis`,
+and defines a step for computing the error norm (L2) for the results
+at each resolution for tracers and layer thickness, saving them in
+`convergence_*.csv` and plotting them in `convergence_*.png`.
+
+### filament_analysis
+
+The class {py:class}`polaris.ocean.tasks.sphere_transport.filament_analysis.FilamentAnalysis`
+computes a filament diagnostic for each threshold value and each resolution
+and plots that diagnostic in `filament.png`.
+
+### viz
+
+Visualization steps are available only in the `nondivergent_2d/with_viz`
+tasks.  They are not included in the `nondivergent_2d` in order to keep regression
+as fast as possible when visualization isn't needed.
+
+The class {py:class}`polaris.ocean.tasks.sphere_transport.viz.VizMap`
+defines a step for creating a mapping file from the MPAS mesh at a given
+resolution to a lon-lat grid at a resolution and interpolation method 
+determined by config options.
+
+```cfg
+# options for visualization for the cosine bell convergence test case
+[sphere_transport_viz]
+
+# visualization latitude and longitude resolution
+dlon = 0.5
+dlat = 0.5
+
+# remapping method ('bilinear', 'neareststod', 'conserve')
+remap_method = conserve
+```
+
+The class {py:class}`polaris.ocean.tasks.sphere_transport.viz.Viz`
+is a step for plotting the initial and final states of the advection test for
+each resolution, mapped to the common lat-lon grid.  The colormap is controlled
+by these options:
+
+```cfg
+# options for visualization for the cosine bell convergence test case
+[sphere_transport_viz_*]
+
+# colormap options
+# colormap
+colormap_name = viridis
+
+# the type of norm used in the colormap
+norm_type = linear
+
+# A dictionary with keywords for the norm
+norm_args = {'vmin': 0., 'vmax': 1.}
+
+# We could provide colorbar tick marks but we'll leave the defaults
+# colorbar_ticks = np.linspace(0., 1., 9)
+```
+
+See {ref}`dev-visualization-global` for more details.
diff --git a/docs/developers_guide/ocean/tasks/rotation_2d.md b/docs/developers_guide/ocean/tasks/rotation_2d.md
new file mode 100644
index 000000000..7be67239e
--- /dev/null
+++ b/docs/developers_guide/ocean/tasks/rotation_2d.md
@@ -0,0 +1,102 @@
+(dev-ocean-rotation-2d)=
+
+# rotation_2d
+
+The {py:class}`polaris.ocean.tasks.sphere_transport.SphereTransport`
+`rotation_2d` test performs a 12-day run on the sphere that has a flow with
+a constant rotation rate which affects tracer distributions. The resolution of
+the sphere varies (by default, between 60 and 240 km). After one period, the
+tracer distributions are compared the initial condition to evaluate numerical
+errors associated with the horizontal advection scheme and determine the rate
+of convergence.
+
+## framework
+
+The config options for the `rotation_2d` test is described in 
+{ref}`ocean-rotation-2d` in the User's Guide.
+
+Additionally, the test uses a `forward.yaml` file with a few common
+model config options related to drag and default horizontal and
+vertical momentum and tracer diffusion, as well as defining `mesh`, `input`,
+`restart`, and `output` streams.  This file has Jinja templating that is
+used to update model config options based on Polaris config options, see
+{ref}`dev-ocean-spherical-convergence`.
+
+### base_mesh
+
+Sphere transport tasks use shared `base_mesh` steps for creating
+{ref}`dev-ocean-spherical-meshes` at a sequence of resolutions.
+
+### init
+
+The class {py:class}`polaris.ocean.tasks.sphere_transport.init.Init`
+defines a step for setting up the initial state at each resolution. The
+initial state is differentiated between `rotation_2d` and the other tests
+set up by the class.
+
+### forward
+
+The class {py:class}`polaris.ocean.tasks.sphere_transport.forward.Forward`
+descends from {py:class}`polaris.ocean.convergence.spherical.SphericalConvergenceForward`,
+and defines a step for running MPAS-Ocean from an initial condition produced in
+an `init` step. See {ref}`dev-ocean-convergence` for some relevant
+discussion of the parent class. The time step is determined from the resolution
+based on the `dt_per_km` config option in the `[convergence_forward]` 
+section.  Other model config options are taken from `forward.yaml`.
+
+### analysis
+
+The class {py:class}`polaris.ocean.tasks.sphere_transport.analysis.Analysis`
+descends from
+{py:class}`polaris.ocean.convergence.ConvergenceAnalysis`,
+and defines a step for computing the error norm (L2) for the results
+at each resolution for tracers and layer thickness, saving them in
+`convergence_*.csv` and plotting them in `convergence_*.png`.
+
+### viz
+
+Visualization steps are available only in the `rotation_2d/with_viz`
+tasks.  They are not included in the `rotation_2d` in order to keep regression
+as fast as possible when visualization isn't needed.
+
+The class {py:class}`polaris.ocean.tasks.sphere_transport.viz.VizMap`
+defines a step for creating a mapping file from the MPAS mesh at a given
+resolution to a lon-lat grid at a resolution and interpolation method 
+determined by config options.
+
+```cfg
+# options for visualization for the cosine bell convergence test case
+[sphere_transport_viz]
+
+# visualization latitude and longitude resolution
+dlon = 0.5
+dlat = 0.5
+
+# remapping method ('bilinear', 'neareststod', 'conserve')
+remap_method = conserve
+```
+
+The class {py:class}`polaris.ocean.tasks.sphere_transport.viz.Viz`
+is a step for plotting the initial and final states of the advection test for
+each resolution, mapped to the common lat-lon grid.  The colormap is controlled
+by these options:
+
+```cfg
+# options for visualization for the cosine bell convergence test case
+[sphere_transport_viz_*]
+
+# colormap options
+# colormap
+colormap_name = viridis
+
+# the type of norm used in the colormap
+norm_type = linear
+
+# A dictionary with keywords for the norm
+norm_args = {'vmin': 0., 'vmax': 1.}
+
+# We could provide colorbar tick marks but we'll leave the defaults
+# colorbar_ticks = np.linspace(0., 1., 9)
+```
+
+See {ref}`dev-visualization-global` for more details.
diff --git a/docs/users_guide/ocean/tasks/correlated_tracers_2d.md b/docs/users_guide/ocean/tasks/correlated_tracers_2d.md
new file mode 100644
index 000000000..aaadcab5d
--- /dev/null
+++ b/docs/users_guide/ocean/tasks/correlated_tracers_2d.md
@@ -0,0 +1,336 @@
+(ocean-correlated-tracers-2d)=
+
+# correlated tracers 2-d
+
+## description
+
+The `correlated_tracers_2d` and `correlated_tracers_2d/with_viz` tasks implement the
+non-divergent flow field test of (1) numerical order of convergence and (2)
+mixing as described in 
+[Lauritzen et al. 2012](<https://gmd.copernicus.org/articles/5/887/2012/>). 
+
+The numerical order of convergence is analyzed in the `analysis` step and
+produces a figure similar to the following showing L2 error norm as a function
+of horizontal resolution:
+
+```{image} images/correlated_tracers_2d_convergence.png
+:align: center
+:width: 500 px
+```
+
+The ability of the horizontal advection scheme to maintain nonlinear
+relationships between tracers is addressed in the `mixing_analysis` step. This
+step produces a visualization of the relationship between two debug tracers at
+the end of the simulation:
+
+```{image} images/correlated_tracers_2d_mixing.png
+:align: center
+:width: 500 px
+```
+
+## mesh
+
+Two global mesh variants are tested, quasi-uniform (QU) and icosohydral. Thus,
+there are 4 variants of the task:
+```
+ocean/spherical/icos/correlated_tracers_2d
+ocean/spherical/icos/correlated_tracers_2d/with_viz
+ocean/spherical/qu/correlated_tracers_2d
+ocean/spherical/qu/correlated_tracers_2d/with_viz
+```
+The default resolutions used in the task depends on the mesh type.
+
+For the `icos` mesh type, the defaults are:
+
+```cfg
+# config options for spherical convergence tests
+[spherical_convergence]
+
+# a list of icosahedral mesh resolutions (km) to test
+icos_resolutions = 60, 120, 240, 480
+```
+
+for the `qu` mesh type, they are:
+
+```cfg
+# config options for spherical convergence tests
+[spherical_convergence]
+
+# a list of quasi-uniform mesh resolutions (km) to test
+qu_resolutions = 60, 90, 120, 150, 180, 210, 240
+```
+
+To alter the resolutions used in this task, you will need to create your own
+config file (or add a `spherical_convergence` section to a config file if 
+you're already using one).  The resolutions are a comma-separated list of the
+resolution of the mesh in km.  If you specify a different list
+before setting up `correlated_tracers_2d`, steps will be generated with the requested
+resolutions.  (If you alter `icos_resolutions` or `qu_resolutions`) in the 
+task's config file in the work directory, nothing will happen.)  For `icos` 
+meshes, make sure you use a resolution close to those listed in 
+{ref}`dev-spherical-meshes`.  Each resolution will be rounded to the nearest 
+allowed icosahedral resolution.
+
+The `base_mesh` steps are shared with other tasks so they are not housed in
+the `correlated_tracers_2d` work directory.  Instead, they are in work directories
+like:
+
+```
+ocean/spherical/icos/base_mesh/60km
+ocean/spherical/qu/base_mesh/60km
+```
+
+For convenience, there are symlinks inside of the `correlated_tracers_2d` and
+`correlated_tracers_2d/with_viz` work directories, e.g.:
+```
+ocean/spherical/icos/correlated_tracers_2d/base_mesh/60km
+ocean/spherical/qu/correlated_tracers_2d/base_mesh/60km
+ocean/spherical/icos/correlated_tracers_2d/with_viz/base_mesh/60km
+ocean/spherical/qu/correlated_tracers_2d/with_viz/base_mesh/60km
+```
+
+## vertical grid
+
+This task only exercises the shallow water dynamics. As such, a single
+vertical level may be used. The bottom depth is constant and the
+results should be insensitive to the choice of `bottom_depth`.
+
+```cfg
+# Options related to the vertical grid
+[vertical_grid]
+
+# the type of vertical grid
+grid_type = uniform
+
+# Number of vertical levels
+vert_levels = 3
+
+# Depth of the bottom of the ocean
+bottom_depth = 300.0
+
+# The type of vertical coordinate (e.g. z-level, z-star)
+coord_type = z-level
+
+# 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 condition is characterized by three separate tracer distributions
+stored in three `debugTracers`:
+
+- `tracer1`: A c-infinity function used for convergence analysis
+- `tracer2`: A pair of c-2 cosine bells
+- `tracer3`: A second pair of cosine bells that are nonlinearly correlated with
+tracer2.
+
+```{image} images/correlated_tracers_2d_init_tracer1.png
+:align: center
+:width: 500 px
+```
+
+```{image} images/correlated_tracers_2d_init_tracer2.png
+:align: center
+:width: 500 px
+```
+
+```{image} images/correlated_tracers_2d_init_tracer3.png
+:align: center
+:width: 500 px
+```
+
+The velocity is
+
+$$
+u(t) = \frac{10 R}{\tau} \sin^2(\lambda - \frac{2 \pi t}{\tau}) \sin(2\theta) \cos(\frac{\pi t}{\tau}) + \frac{2 \pi R}{\tau}\cos(\theta)
+$$
+
+$$
+v(t) = \frac{10 R}{\tau} \sin(2\lambda - \frac{4 \pi t}{\tau}) \cos(\theta) \cos(\frac{\pi t}{\tau})
+$$
+
+Where $\lambda$ is longitude, $\theta$ is latitude and $R$ is the radius of the
+sphere. $\tau$ is the time it takes to transit the equator. The default is 12
+days and is given by the cfg option `vel_pd`.
+
+Temperature and salinity are not evolved in this task and are given
+constant values determined by config options `temperature` and `salinity`.
+
+The Coriolis parameters `fCell`, `fEdge`, and `fVertex` do not need to be
+specified for a global mesh and are initialized as zeros.
+
+## forcing
+
+This case is forced to follow $u(t)$ and $v(t)$ given above.
+
+## time step and run duration
+
+This task uses the Runge-Kutta 4th-order (RK4) time integrator. The time step 
+for forward integration is determined by multiplying the resolution by a config
+option, `rk4_dt_per_km`, so that coarser meshes have longer time steps. You can
+alter this before setup (in a user config file) or before running the task (in 
+the config file in the work directory). 
+
+```cfg
+# config options for spherical convergence tests
+[spherical_convergence_forward]
+
+# time integrator: {'split_explicit', 'RK4'}
+time_integrator = RK4
+
+# RK4 time step per resolution (s/km), since dt is proportional to resolution
+rk4_dt_per_km = 3.0
+```
+
+The `convergence_eval_time`, `run_duration` and `output_interval` are the
+period for advection to make a full rotation around the globe, 12 days:
+
+```cfg
+# config options for spherical convergence tests
+[spherical_convergence_forward]
+
+# Run duration in days
+run_duration = ${sphere_transport:vel_pd}
+
+# Output interval in days
+output_interval = ${sphere_transport:vel_pd}
+```
+
+Here, `${sphere_transport:vel_pd}` means that the same value is used as in the 
+option `vel_pd` in section `[sphere_transport]`, see below.
+
+## config options
+
+The `correlated_tracer_2d` config options include:
+
+```cfg
+# options for all sphere transport test cases
+[sphere_transport]
+
+# temperature
+temperature = 15.
+
+# salinity
+salinity = 35.
+
+# time (days) for bell to transit equator once
+vel_pd = 12.0
+
+
+# options for tracer visualization for the sphere transport test case
+[sphere_transport_viz]
+
+# visualization latitude and longitude resolution
+dlon = 0.5
+dlat = 0.5
+
+# remapping method ('bilinear', 'neareststod', 'conserve')
+remap_method = conserve
+
+
+# options for tracer visualization for the sphere transport test case
+[sphere_transport_viz_tracer]
+
+# colormap options
+# colormap
+colormap_name = viridis
+
+# the type of norm used in the colormap
+norm_type = linear
+
+# A dictionary with keywords for the norm
+norm_args = {'vmin': 0., 'vmax': 1.}
+
+# We could provide colorbar tick marks but we'll leave the defaults
+# colorbar_ticks = np.linspace(0., 1., 9)
+
+# options for plotting tracer differences from sphere transport tests
+[sphere_transport_viz_tracer_diff]
+
+# colormap options
+# colormap
+colormap_name = cmo.balance
+
+# the type of norm used in the colormap
+norm_type = linear
+
+# A dictionary with keywords for the norm
+norm_args = {'vmin': -0.25, 'vmax': 0.25}
+
+
+# options for thickness visualization for the sphere transport test case
+[sphere_transport_viz_h]
+
+# colormap options
+# colormap
+colormap_name = viridis
+
+# the type of norm used in the colormap
+norm_type = linear
+
+# A dictionary with keywords for the norm
+norm_args = {'vmin': 99., 'vmax': 101.}
+
+# We could provide colorbar tick marks but we'll leave the defaults
+# colorbar_ticks = np.linspace(0., 1., 9)
+
+# options for plotting tracer differences from sphere transport tests
+[sphere_transport_viz_h_diff]
+
+# colormap options
+# colormap
+colormap_name = cmo.balance
+
+# the type of norm used in the colormap
+norm_type = linear
+
+# A dictionary with keywords for the norm
+norm_args = {'vmin': -0.25, 'vmax': 0.25}
+
+# options for correlated tracers 2-d test case
+[correlated_tracers_2d]
+
+# velocity amplitude in meters per second
+vel_amp = 10.
+
+# convergence threshold below which the test fails
+convergence_thresh_tracer1 = 1.7
+convergence_thresh_tracer2 = 1.66
+convergence_thresh_tracer3 = 1.28
+
+# time in days at which to evaluate mixing
+mixing_evaluation_time = 6.0
+```
+
+The options in section `sphere_transport` are used by all 4 test cases based
+on Lauritzen et al. (2012) and control the initial condition. The options in
+section `correlated_tracers_2d` control the initial condition of that case
+only, the convergence rate threshold, and the time at which to evaluate mixing
+diagnostics.
+
+The options in sections `sphere_transport_viz*` control properties of the `viz`
+step of the test case. 
+
+The default options for the convergence analysis step can be changed here:
+
+```cfg
+# config options for spherical convergence tests
+[spherical_convergence]
+
+# Evaluation time for convergence analysis (in days)
+convergence_eval_time = ${sphere_transport:vel_pd}
+
+# Type of error to compute
+error_type = l2
+```
+
+## cores
+
+The target and minimum number of cores are determined by `goal_cells_per_core`
+and `max_cells_per_core` from the `ocean` section of the config file,
+respectively. This ensures that the number of cells per core is roughly
+constant across the different resolutions in the convergence study.
diff --git a/docs/users_guide/ocean/tasks/divergent_2d.md b/docs/users_guide/ocean/tasks/divergent_2d.md
new file mode 100644
index 000000000..024d8c0f3
--- /dev/null
+++ b/docs/users_guide/ocean/tasks/divergent_2d.md
@@ -0,0 +1,321 @@
+(ocean-divergent-2d)=
+
+# divergent 2-d
+
+## description
+
+The `divergent_2d` and `divergent_2d/with_viz` tasks implement the
+divergent flow field test of numerical order of convergence as described in 
+[Lauritzen et al. 2012](<https://gmd.copernicus.org/articles/5/887/2012/>).
+
+The numerical order of convergence is analyzed in the `analysis` step and
+produces a figure similar to the following showing L2 error norm as a function
+of horizontal resolution:
+
+```{image} images/divergent_2d_convergence.png
+:align: center
+:width: 500 px
+```
+
+## mesh
+
+Two global mesh variants are tested, quasi-uniform (QU) and icosohydral. Thus,
+there are 4 variants of the task:
+```
+ocean/spherical/icos/divergent_2d
+ocean/spherical/icos/divergent_2d/with_viz
+ocean/spherical/qu/divergent_2d
+ocean/spherical/qu/divergent_2d/with_viz
+```
+The default resolutions used in the task depends on the mesh type.
+
+For the `icos` mesh type, the defaults are:
+
+```cfg
+# config options for spherical convergence tests
+[spherical_convergence]
+
+# a list of icosahedral mesh resolutions (km) to test
+icos_resolutions = 60, 120, 240, 480
+```
+
+for the `qu` mesh type, they are:
+
+```cfg
+# config options for spherical convergence tests
+[spherical_convergence]
+
+# a list of quasi-uniform mesh resolutions (km) to test
+qu_resolutions = 60, 90, 120, 150, 180, 210, 240
+```
+
+To alter the resolutions used in this task, you will need to create your own
+config file (or add a `spherical_convergence` section to a config file if 
+you're already using one).  The resolutions are a comma-separated list of the
+resolution of the mesh in km.  If you specify a different list
+before setting up `divergent_2d`, steps will be generated with the requested
+resolutions.  (If you alter `icos_resolutions` or `qu_resolutions`) in the 
+task's config file in the work directory, nothing will happen.)  For `icos` 
+meshes, make sure you use a resolution close to those listed in 
+{ref}`dev-spherical-meshes`.  Each resolution will be rounded to the nearest 
+allowed icosahedral resolution.
+
+The `base_mesh` steps are shared with other tasks so they are not housed in
+the `divergent_2d` work directory.  Instead, they are in work directories
+like:
+
+```
+ocean/spherical/icos/base_mesh/60km
+ocean/spherical/qu/base_mesh/60km
+```
+
+For convenience, there are symlinks inside of the `divergent_2d` and
+`divergent_2d/with_viz` work directories, e.g.:
+```
+ocean/spherical/icos/divergent_2d/base_mesh/60km
+ocean/spherical/qu/divergent_2d/base_mesh/60km
+ocean/spherical/icos/divergent_2d/with_viz/base_mesh/60km
+ocean/spherical/qu/divergent_2d/with_viz/base_mesh/60km
+```
+
+## vertical grid
+
+This task only exercises the shallow water dynamics. As such, a single
+vertical level may be used. The bottom depth is constant and the
+results should be insensitive to the choice of `bottom_depth`.
+
+```cfg
+# Options related to the vertical grid
+[vertical_grid]
+
+# the type of vertical grid
+grid_type = uniform
+
+# Number of vertical levels
+vert_levels = 3
+
+# Depth of the bottom of the ocean
+bottom_depth = 300.0
+
+# The type of vertical coordinate (e.g. z-level, z-star)
+coord_type = z-level
+
+# 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 condition is characterized by three separate tracer distributions
+stored in three `debugTracers`:
+
+- `tracer1`: A c-infinity function used for convergence analysis
+- `tracer2`: A pair of c-2 cosine bells
+- `tracer3`: A discontinuous pair of slotted cylinders
+
+```{image} images/divergent_2d_init_tracer1.png
+:align: center
+:width: 500 px
+```
+
+```{image} images/divergent_2d_init_tracer2.png
+:align: center
+:width: 500 px
+```
+
+```{image} images/divergent_2d_init_tracer3.png
+:align: center
+:width: 500 px
+```
+
+The velocity is
+
+$$
+u(t) = \frac{u_0 R}{\tau} \sin^2(\frac{\lambda}{2} - \frac{\pi t}{\tau}) \sin(2\theta) \cos^2(\theta) \cos(\frac{\pi t}{\tau}) + \frac{2 \pi R}{\tau}\cos(\theta)
+$$
+
+$$
+v(t) = \frac{u_0 R}{2 \tau} \sin(\lambda - \frac{2 \pi t}{\tau}) \cos^3(\theta) \cos(\frac{\pi t}{\tau})
+$$
+
+Where $\lambda$ is longitude, $\theta$ is latitude and $R$ is the radius of the
+sphere. $\tau$ is the time it takes to transit the equator. The default is 12
+days and is given by the cfg option `vel_pd`. The default velocity amplitude,
+$u_0$ is 10 m/s and is given by the cfg option `vel_amp`.
+
+Temperature and salinity are not evolved in this task and are given
+constant values determined by config options `temperature` and `salinity`.
+
+The Coriolis parameters `fCell`, `fEdge`, and `fVertex` do not need to be
+specified for a global mesh and are initialized as zeros.
+
+## forcing
+
+This case is forced to follow $u(t)$ and $v(t)$ given above.
+
+## time step and run duration
+
+This task uses the Runge-Kutta 4th-order (RK4) time integrator. The time step 
+for forward integration is determined by multiplying the resolution by a config
+option, `rk4_dt_per_km`, so that coarser meshes have longer time steps. You can
+alter this before setup (in a user config file) or before running the task (in 
+the config file in the work directory). 
+
+```cfg
+# config options for spherical convergence tests
+[spherical_convergence_forward]
+
+# time integrator: {'split_explicit', 'RK4'}
+time_integrator = RK4
+
+# RK4 time step per resolution (s/km), since dt is proportional to resolution
+rk4_dt_per_km = 3.0
+```
+
+The `convergence_eval_time`, `run_duration` and `output_interval` are the
+period for advection to make a full rotation around the globe, 12 days:
+
+```cfg
+# config options for spherical convergence tests
+[spherical_convergence_forward]
+
+# Run duration in days
+run_duration = ${sphere_transport:vel_pd}
+
+# Output interval in days
+output_interval = ${sphere_transport:vel_pd}
+```
+
+Here, `${sphere_transport:vel_pd}` means that the same value is used as in the 
+option `vel_pd` in section `[sphere_transport]`, see below.
+
+## config options
+
+The `divergent_2d` config options include:
+
+```cfg
+# options for all sphere transport test cases
+[sphere_transport]
+
+# temperature
+temperature = 15.
+
+# salinity
+salinity = 35.
+
+# time (days) for bell to transit equator once
+vel_pd = 12.0
+
+
+# options for tracer visualization for the sphere transport test case
+[sphere_transport_viz]
+
+# visualization latitude and longitude resolution
+dlon = 0.5
+dlat = 0.5
+
+# remapping method ('bilinear', 'neareststod', 'conserve')
+remap_method = conserve
+
+
+# options for tracer visualization for the sphere transport test case
+[sphere_transport_viz_tracer]
+
+# colormap options
+# colormap
+colormap_name = viridis
+
+# the type of norm used in the colormap
+norm_type = linear
+
+# A dictionary with keywords for the norm
+norm_args = {'vmin': 0., 'vmax': 1.}
+
+# We could provide colorbar tick marks but we'll leave the defaults
+# colorbar_ticks = np.linspace(0., 1., 9)
+
+# options for plotting tracer differences from sphere transport tests
+[sphere_transport_viz_tracer_diff]
+
+# colormap options
+# colormap
+colormap_name = cmo.balance
+
+# the type of norm used in the colormap
+norm_type = linear
+
+# A dictionary with keywords for the norm
+norm_args = {'vmin': -0.25, 'vmax': 0.25}
+
+
+# options for thickness visualization for the sphere transport test case
+[sphere_transport_viz_h]
+
+# colormap options
+# colormap
+colormap_name = viridis
+
+# the type of norm used in the colormap
+norm_type = linear
+
+# A dictionary with keywords for the norm
+norm_args = {'vmin': 99., 'vmax': 101.}
+
+# We could provide colorbar tick marks but we'll leave the defaults
+# colorbar_ticks = np.linspace(0., 1., 9)
+
+# options for plotting tracer differences from sphere transport tests
+[sphere_transport_viz_h_diff]
+
+# colormap options
+# colormap
+colormap_name = cmo.balance
+
+# the type of norm used in the colormap
+norm_type = linear
+
+# A dictionary with keywords for the norm
+norm_args = {'vmin': -0.25, 'vmax': 0.25}
+
+# options for divergent 2-d test case
+[divergent_2d]
+
+# velocity amplitude in meters per second
+vel_amp = 5.
+
+# convergence threshold below which the test fails
+convergence_thresh_tracer1 = 1.48
+convergence_thresh_tracer2 = 1.86
+convergence_thresh_tracer3 = 0.4
+```
+
+The options in section `sphere_transport` are used by all 4 test cases based
+on Lauritzen et al. (2012) and control the initial condition. The options in
+section `divergent_2d` control the initial condition of that case
+only and the convergence rate threshold.
+
+The options in sections `sphere_transport_viz*` control properties of the `viz`
+step of the test case. 
+
+The default options for the convergence analysis step can be changed here:
+
+```cfg
+# config options for spherical convergence tests
+[spherical_convergence]
+
+# Evaluation time for convergence analysis (in days)
+convergence_eval_time = ${sphere_transport:vel_pd}
+
+# Type of error to compute
+error_type = l2
+```
+
+## cores
+
+The target and minimum number of cores are determined by `goal_cells_per_core`
+and `max_cells_per_core` from the `ocean` section of the config file,
+respectively. This ensures that the number of cells per core is roughly
+constant across the different resolutions in the convergence study.
diff --git a/docs/users_guide/ocean/tasks/images/correlated_tracers_2d_convergence.png b/docs/users_guide/ocean/tasks/images/correlated_tracers_2d_convergence.png
new file mode 100644
index 000000000..e410687f2
Binary files /dev/null and b/docs/users_guide/ocean/tasks/images/correlated_tracers_2d_convergence.png differ
diff --git a/docs/users_guide/ocean/tasks/images/correlated_tracers_2d_init_tracer1.png b/docs/users_guide/ocean/tasks/images/correlated_tracers_2d_init_tracer1.png
new file mode 100644
index 000000000..acfece6ed
Binary files /dev/null and b/docs/users_guide/ocean/tasks/images/correlated_tracers_2d_init_tracer1.png differ
diff --git a/docs/users_guide/ocean/tasks/images/correlated_tracers_2d_init_tracer2.png b/docs/users_guide/ocean/tasks/images/correlated_tracers_2d_init_tracer2.png
new file mode 100644
index 000000000..e1e8b157c
Binary files /dev/null and b/docs/users_guide/ocean/tasks/images/correlated_tracers_2d_init_tracer2.png differ
diff --git a/docs/users_guide/ocean/tasks/images/correlated_tracers_2d_init_tracer3.png b/docs/users_guide/ocean/tasks/images/correlated_tracers_2d_init_tracer3.png
new file mode 100644
index 000000000..c0c098b84
Binary files /dev/null and b/docs/users_guide/ocean/tasks/images/correlated_tracers_2d_init_tracer3.png differ
diff --git a/docs/users_guide/ocean/tasks/images/correlated_tracers_2d_mixing.png b/docs/users_guide/ocean/tasks/images/correlated_tracers_2d_mixing.png
new file mode 100644
index 000000000..c2a24582f
Binary files /dev/null and b/docs/users_guide/ocean/tasks/images/correlated_tracers_2d_mixing.png differ
diff --git a/docs/users_guide/ocean/tasks/images/divergent_2d_convergence.png b/docs/users_guide/ocean/tasks/images/divergent_2d_convergence.png
new file mode 100644
index 000000000..9e02fa489
Binary files /dev/null and b/docs/users_guide/ocean/tasks/images/divergent_2d_convergence.png differ
diff --git a/docs/users_guide/ocean/tasks/images/divergent_2d_init_tracer1.png b/docs/users_guide/ocean/tasks/images/divergent_2d_init_tracer1.png
new file mode 100644
index 000000000..acfece6ed
Binary files /dev/null and b/docs/users_guide/ocean/tasks/images/divergent_2d_init_tracer1.png differ
diff --git a/docs/users_guide/ocean/tasks/images/divergent_2d_init_tracer2.png b/docs/users_guide/ocean/tasks/images/divergent_2d_init_tracer2.png
new file mode 100644
index 000000000..e1e8b157c
Binary files /dev/null and b/docs/users_guide/ocean/tasks/images/divergent_2d_init_tracer2.png differ
diff --git a/docs/users_guide/ocean/tasks/images/divergent_2d_init_tracer3.png b/docs/users_guide/ocean/tasks/images/divergent_2d_init_tracer3.png
new file mode 100644
index 000000000..179f42ede
Binary files /dev/null and b/docs/users_guide/ocean/tasks/images/divergent_2d_init_tracer3.png differ
diff --git a/docs/users_guide/ocean/tasks/images/nondivergent_2d_convergence.png b/docs/users_guide/ocean/tasks/images/nondivergent_2d_convergence.png
new file mode 100644
index 000000000..e410687f2
Binary files /dev/null and b/docs/users_guide/ocean/tasks/images/nondivergent_2d_convergence.png differ
diff --git a/docs/users_guide/ocean/tasks/images/nondivergent_2d_filament.png b/docs/users_guide/ocean/tasks/images/nondivergent_2d_filament.png
new file mode 100644
index 000000000..b30297f98
Binary files /dev/null and b/docs/users_guide/ocean/tasks/images/nondivergent_2d_filament.png differ
diff --git a/docs/users_guide/ocean/tasks/images/nondivergent_2d_init_tracer1.png b/docs/users_guide/ocean/tasks/images/nondivergent_2d_init_tracer1.png
new file mode 100644
index 000000000..acfece6ed
Binary files /dev/null and b/docs/users_guide/ocean/tasks/images/nondivergent_2d_init_tracer1.png differ
diff --git a/docs/users_guide/ocean/tasks/images/nondivergent_2d_init_tracer2.png b/docs/users_guide/ocean/tasks/images/nondivergent_2d_init_tracer2.png
new file mode 100644
index 000000000..e1e8b157c
Binary files /dev/null and b/docs/users_guide/ocean/tasks/images/nondivergent_2d_init_tracer2.png differ
diff --git a/docs/users_guide/ocean/tasks/images/nondivergent_2d_init_tracer3.png b/docs/users_guide/ocean/tasks/images/nondivergent_2d_init_tracer3.png
new file mode 100644
index 000000000..179f42ede
Binary files /dev/null and b/docs/users_guide/ocean/tasks/images/nondivergent_2d_init_tracer3.png differ
diff --git a/docs/users_guide/ocean/tasks/images/rotation_2d_convergence.png b/docs/users_guide/ocean/tasks/images/rotation_2d_convergence.png
new file mode 100644
index 000000000..6733f7962
Binary files /dev/null and b/docs/users_guide/ocean/tasks/images/rotation_2d_convergence.png differ
diff --git a/docs/users_guide/ocean/tasks/images/rotation_2d_init_tracer1.png b/docs/users_guide/ocean/tasks/images/rotation_2d_init_tracer1.png
new file mode 100644
index 000000000..acfece6ed
Binary files /dev/null and b/docs/users_guide/ocean/tasks/images/rotation_2d_init_tracer1.png differ
diff --git a/docs/users_guide/ocean/tasks/images/rotation_2d_init_tracer2.png b/docs/users_guide/ocean/tasks/images/rotation_2d_init_tracer2.png
new file mode 100644
index 000000000..e1e8b157c
Binary files /dev/null and b/docs/users_guide/ocean/tasks/images/rotation_2d_init_tracer2.png differ
diff --git a/docs/users_guide/ocean/tasks/images/rotation_2d_init_tracer3.png b/docs/users_guide/ocean/tasks/images/rotation_2d_init_tracer3.png
new file mode 100644
index 000000000..179f42ede
Binary files /dev/null and b/docs/users_guide/ocean/tasks/images/rotation_2d_init_tracer3.png differ
diff --git a/docs/users_guide/ocean/tasks/index.md b/docs/users_guide/ocean/tasks/index.md
index 49c1cac45..ff7993372 100644
--- a/docs/users_guide/ocean/tasks/index.md
+++ b/docs/users_guide/ocean/tasks/index.md
@@ -6,9 +6,13 @@
 :titlesonly: true
 
 baroclinic_channel
+correlated_tracers_2d
 cosine_bell
 geostrophic
+divergent_2d
 inertial_gravity_wave
 manufactured_solution
+nondivergent_2d
+rotation_2d
 single_column
 ```
diff --git a/docs/users_guide/ocean/tasks/nondivergent_2d.md b/docs/users_guide/ocean/tasks/nondivergent_2d.md
new file mode 100644
index 000000000..15d1ef71e
--- /dev/null
+++ b/docs/users_guide/ocean/tasks/nondivergent_2d.md
@@ -0,0 +1,344 @@
+(ocean-nondivergent-2d)=
+
+# nondivergent 2-d
+
+## description
+
+The `nondivergent_2d` and `nondivergent_2d/with_viz` tasks implement the
+non-divergent flow field test of (1) numerical order of convergence, (2)
+filament preservation and (3) rough distribution as described in 
+[Lauritzen et al. 2012](<https://gmd.copernicus.org/articles/5/887/2012/>).
+The reversing deformational flow in this test case explores large-scale to small-scale transport.
+
+The numerical order of convergence is analyzed in the `analysis` step and
+produces a figure similar to the following showing L2 error norm as a function
+of horizontal resolution:
+
+```{image} images/nondivergent_2d_convergence.png
+:align: center
+:width: 500 px
+```
+
+The ability of the horizontal advection scheme to preserve filaments is
+addressed in the `filament_analysis` step. It produces a figure similar to the
+following, which provides a comparison across resolutions, similar to Figure 6
+from Lauritzen et al. 2012:
+
+```{image} images/nondivergent_2d_filament.png
+:align: center
+:width: 500 px
+```
+
+The ability of the horizontal advection scheme to preserve rough distributions
+is addressed in the `viz` steps for each resolution of the
+`nondivergent_2d/with_viz` test. These steps produce visualizations of three
+debug tracers at initial time, thhe mid-point and final time. The rough
+distributions are represented as slotted cylinders in the `tracer3`
+field.
+
+## mesh
+
+Two global mesh variants are tested, quasi-uniform (QU) and icosohydral. Thus,
+there are 4 variants of the task:
+```
+ocean/spherical/icos/nondivergent_2d
+ocean/spherical/icos/nondivergent_2d/with_viz
+ocean/spherical/qu/nondivergent_2d
+ocean/spherical/qu/nondivergent_2d/with_viz
+```
+The default resolutions used in the task depends on the mesh type.
+
+For the `icos` mesh type, the defaults are:
+
+```cfg
+# config options for spherical convergence tests
+[spherical_convergence]
+
+# a list of icosahedral mesh resolutions (km) to test
+icos_resolutions = 60, 120, 240, 480
+```
+
+for the `qu` mesh type, they are:
+
+```cfg
+# config options for spherical convergence tests
+[spherical_convergence]
+
+# a list of quasi-uniform mesh resolutions (km) to test
+qu_resolutions = 60, 90, 120, 150, 180, 210, 240
+```
+
+To alter the resolutions used in this task, you will need to create your own
+config file (or add a `spherical_convergence` section to a config file if 
+you're already using one).  The resolutions are a comma-separated list of the
+resolution of the mesh in km.  If you specify a different list
+before setting up `nondivergent_2d`, steps will be generated with the requested
+resolutions.  (If you alter `icos_resolutions` or `qu_resolutions`) in the 
+task's config file in the work directory, nothing will happen.)  For `icos` 
+meshes, make sure you use a resolution close to those listed in 
+{ref}`dev-spherical-meshes`.  Each resolution will be rounded to the nearest 
+allowed icosahedral resolution.
+
+The `base_mesh` steps are shared with other tasks so they are not housed in
+the `nondivergent_2d` work directory.  Instead, they are in work directories
+like:
+
+```
+ocean/spherical/icos/base_mesh/60km
+ocean/spherical/qu/base_mesh/60km
+```
+
+For convenience, there are symlinks inside of the `nondivergent_2d` and
+`nondivergent_2d/with_viz` work directories, e.g.:
+```
+ocean/spherical/icos/nondivergent_2d/base_mesh/60km
+ocean/spherical/qu/nondivergent_2d/base_mesh/60km
+ocean/spherical/icos/nondivergent_2d/with_viz/base_mesh/60km
+ocean/spherical/qu/nondivergent_2d/with_viz/base_mesh/60km
+```
+
+## vertical grid
+
+This task only exercises the shallow water dynamics. As such, a single
+vertical level may be used. The bottom depth is constant and the
+results should be insensitive to the choice of `bottom_depth`.
+
+```cfg
+# Options related to the vertical grid
+[vertical_grid]
+
+# the type of vertical grid
+grid_type = uniform
+
+# Number of vertical levels
+vert_levels = 3
+
+# Depth of the bottom of the ocean
+bottom_depth = 300.0
+
+# The type of vertical coordinate (e.g. z-level, z-star)
+coord_type = z-level
+
+# 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 condition is characterized by three separate tracer distributions
+stored in three `debugTracers`:
+
+- `tracer1`: A c-infinity function used for convergence analysis
+- `tracer2`: A pair of c-2 cosine bells
+- `tracer3`: A discontinuous pair of slotted cylinders
+
+```{image} images/nondivergent_2d_init_tracer1.png
+:align: center
+:width: 500 px
+```
+
+```{image} images/nondivergent_2d_init_tracer2.png
+:align: center
+:width: 500 px
+```
+
+```{image} images/nondivergent_2d_init_tracer3.png
+:align: center
+:width: 500 px
+```
+
+The velocity is
+
+$$
+u(t) = \frac{u_0 R}{\tau} \sin^2(\lambda - \frac{2 \pi t}{\tau}) \sin(2\theta) \cos(\frac{\pi t}{\tau}) + \frac{2 \pi R}{\tau}\cos(\theta)
+$$
+
+$$
+v(t) = \frac{u_0 R}{\tau} \sin(2\lambda - \frac{4 \pi t}{\tau}) \cos(\theta) \cos(\frac{\pi t}{\tau})
+$$
+
+Where $\lambda$ is longitude, $\theta$ is latitude and $R$ is the radius of the
+sphere. $\tau$ is the time it takes to transit the equator. The default is 12
+days and is given by the cfg option `vel_pd`. The default velocity amplitude,
+$u_0$ is 10 m/s and is given by the cfg option `vel_amp`.
+
+Temperature and salinity are not evolved in this task and are given
+constant values determined by config options `temperature` and `salinity`.
+
+The Coriolis parameters `fCell`, `fEdge`, and `fVertex` do not need to be
+specified for a global mesh and are initialized as zeros.
+
+## forcing
+
+This case is forced to follow $u(t)$ and $v(t)$ given above.
+
+## time step and run duration
+
+This task uses the Runge-Kutta 4th-order (RK4) time integrator. The time step 
+for forward integration is determined by multiplying the resolution by a config
+option, `rk4_dt_per_km`, so that coarser meshes have longer time steps. You can
+alter this before setup (in a user config file) or before running the task (in 
+the config file in the work directory). 
+
+```cfg
+# config options for spherical convergence tests
+[spherical_convergence_forward]
+
+# time integrator: {'split_explicit', 'RK4'}
+time_integrator = RK4
+
+# RK4 time step per resolution (s/km), since dt is proportional to resolution
+rk4_dt_per_km = 3.0
+```
+
+The `convergence_eval_time`, `run_duration` and `output_interval` are the
+period for advection to make a full rotation around the globe, 12 days:
+
+```cfg
+# config options for spherical convergence tests
+[spherical_convergence_forward]
+
+# Run duration in days
+run_duration = ${sphere_transport:vel_pd}
+
+# Output interval in days
+output_interval = ${sphere_transport:vel_pd}
+```
+
+Here, `${sphere_transport:vel_pd}` means that the same value is used as in the 
+option `vel_pd` in section `[sphere_transport]`, see below.
+
+## config options
+
+The `nondivergent_2d` config options include:
+
+```cfg
+# options for all sphere transport test cases
+[sphere_transport]
+
+# temperature
+temperature = 15.
+
+# salinity
+salinity = 35.
+
+# time (days) for bell to transit equator once
+vel_pd = 12.0
+
+
+# options for tracer visualization for the sphere transport test case
+[sphere_transport_viz]
+
+# visualization latitude and longitude resolution
+dlon = 0.5
+dlat = 0.5
+
+# remapping method ('bilinear', 'neareststod', 'conserve')
+remap_method = conserve
+
+
+# options for tracer visualization for the sphere transport test case
+[sphere_transport_viz_tracer]
+
+# colormap options
+# colormap
+colormap_name = viridis
+
+# the type of norm used in the colormap
+norm_type = linear
+
+# A dictionary with keywords for the norm
+norm_args = {'vmin': 0., 'vmax': 1.}
+
+# We could provide colorbar tick marks but we'll leave the defaults
+# colorbar_ticks = np.linspace(0., 1., 9)
+
+# options for plotting tracer differences from sphere transport tests
+[sphere_transport_viz_tracer_diff]
+
+# colormap options
+# colormap
+colormap_name = cmo.balance
+
+# the type of norm used in the colormap
+norm_type = linear
+
+# A dictionary with keywords for the norm
+norm_args = {'vmin': -0.25, 'vmax': 0.25}
+
+
+# options for thickness visualization for the sphere transport test case
+[sphere_transport_viz_h]
+
+# colormap options
+# colormap
+colormap_name = viridis
+
+# the type of norm used in the colormap
+norm_type = linear
+
+# A dictionary with keywords for the norm
+norm_args = {'vmin': 99., 'vmax': 101.}
+
+# We could provide colorbar tick marks but we'll leave the defaults
+# colorbar_ticks = np.linspace(0., 1., 9)
+
+# options for plotting tracer differences from sphere transport tests
+[sphere_transport_viz_h_diff]
+
+# colormap options
+# colormap
+colormap_name = cmo.balance
+
+# the type of norm used in the colormap
+norm_type = linear
+
+# A dictionary with keywords for the norm
+norm_args = {'vmin': -0.25, 'vmax': 0.25}
+
+# options for nondivergent 2-d test case
+[nondivergent_2d]
+
+# velocity amplitude in meters per second
+vel_amp = 10.
+
+# convergence threshold below which the test fails
+convergence_thresh_tracer1 = 1.5
+convergence_thresh_tracer2 = 1.3
+convergence_thresh_tracer3 = 0.3
+
+# time in days at which to evaluate filament preservation
+filament_evaluation_time = 6.0
+```
+
+The options in section `sphere_transport` are used by all 4 test cases based
+on Lauritzen et al. (2012) and control the initial condition. The options in
+section `nondivergent_2d` control the initial condition of that case
+only, the convergence rate threshold, and the time at which to evaluate
+filament preservation.
+
+The options in sections `sphere_transport_viz*` control properties of the `viz`
+step of the test case. 
+
+The default options for the convergence analysis step can be changed here:
+
+```cfg
+# config options for spherical convergence tests
+[spherical_convergence]
+
+# Evaluation time for convergence analysis (in days)
+convergence_eval_time = ${sphere_transport:vel_pd}
+
+# Type of error to compute
+error_type = l2
+```
+
+## cores
+
+The target and minimum number of cores are determined by `goal_cells_per_core`
+and `max_cells_per_core` from the `ocean` section of the config file,
+respectively. This ensures that the number of cells per core is roughly
+constant across the different resolutions in the convergence study.
diff --git a/docs/users_guide/ocean/tasks/rotation_2d.md b/docs/users_guide/ocean/tasks/rotation_2d.md
new file mode 100644
index 000000000..4fd90417e
--- /dev/null
+++ b/docs/users_guide/ocean/tasks/rotation_2d.md
@@ -0,0 +1,310 @@
+(ocean-rotation-2d)=
+
+# rotation 2-d
+
+## description
+
+The `rotation_2d` and `rotation_2d/with_viz` tasks implement the
+rotational flow field test of numerical order of convergence.
+This test is similar to `cosine_bell` except the axis of rotation is a config
+option and can be offset from the z-axis.
+
+The numerical order of convergence is analyzed in the `analysis` step and
+produces a figure similar to the following showing L2 error norm as a function
+of horizontal resolution:
+
+```{image} images/rotation_2d_convergence.png
+:align: center
+:width: 500 px
+```
+
+## mesh
+
+Two global mesh variants are tested, quasi-uniform (QU) and icosohydral. Thus,
+there are 4 variants of the task:
+```
+ocean/spherical/icos/rotation_2d
+ocean/spherical/icos/rotation_2d/with_viz
+ocean/spherical/qu/rotation_2d
+ocean/spherical/qu/rotation_2d/with_viz
+```
+The default resolutions used in the task depends on the mesh type.
+
+For the `icos` mesh type, the defaults are:
+
+```cfg
+# config options for spherical convergence tests
+[spherical_convergence]
+
+# a list of icosahedral mesh resolutions (km) to test
+icos_resolutions = 60, 120, 240, 480
+```
+
+for the `qu` mesh type, they are:
+
+```cfg
+# config options for spherical convergence tests
+[spherical_convergence]
+
+# a list of quasi-uniform mesh resolutions (km) to test
+qu_resolutions = 60, 90, 120, 150, 180, 210, 240
+```
+
+To alter the resolutions used in this task, you will need to create your own
+config file (or add a `spherical_convergence` section to a config file if 
+you're already using one).  The resolutions are a comma-separated list of the
+resolution of the mesh in km.  If you specify a different list
+before setting up `rotation_2d`, steps will be generated with the requested
+resolutions.  (If you alter `icos_resolutions` or `qu_resolutions`) in the 
+task's config file in the work directory, nothing will happen.)  For `icos` 
+meshes, make sure you use a resolution close to those listed in 
+{ref}`dev-spherical-meshes`.  Each resolution will be rounded to the nearest 
+allowed icosahedral resolution.
+
+The `base_mesh` steps are shared with other tasks so they are not housed in
+the `rotation_2d` work directory.  Instead, they are in work directories
+like:
+
+```
+ocean/spherical/icos/base_mesh/60km
+ocean/spherical/qu/base_mesh/60km
+```
+
+For convenience, there are symlinks inside of the `rotation_2d` and
+`rotation_2d/with_viz` work directories, e.g.:
+```
+ocean/spherical/icos/rotation_2d/base_mesh/60km
+ocean/spherical/qu/rotation_2d/base_mesh/60km
+ocean/spherical/icos/rotation_2d/with_viz/base_mesh/60km
+ocean/spherical/qu/rotation_2d/with_viz/base_mesh/60km
+```
+
+## vertical grid
+
+This task only exercises the shallow water dynamics. As such, a single
+vertical level may be used. The bottom depth is constant and the
+results should be insensitive to the choice of `bottom_depth`.
+
+```cfg
+# Options related to the vertical grid
+[vertical_grid]
+
+# the type of vertical grid
+grid_type = uniform
+
+# Number of vertical levels
+vert_levels = 3
+
+# Depth of the bottom of the ocean
+bottom_depth = 300.0
+
+# The type of vertical coordinate (e.g. z-level, z-star)
+coord_type = z-level
+
+# 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 condition is characterized by three separate tracer distributions
+stored in three `debugTracers`:
+
+- `tracer1`: A c-infinity function used for convergence analysis
+- `tracer2`: A pair of c-2 cosine bells
+- `tracer3`: A discontinuous pair of slotted cylinders
+
+```{image} images/rotation_2d_init_tracer1.png
+:align: center
+:width: 500 px
+```
+
+```{image} images/rotation_2d_init_tracer2.png
+:align: center
+:width: 500 px
+```
+
+```{image} images/rotation_2d_init_tracer3.png
+:align: center
+:width: 500 px
+```
+
+The velocity is that of rigid rotation about an axis offset from the z-axis of
+the sphere. It is not given in Lauritzen et al. The axis of rotation is defined by a vector given by the cfg option `rotation_vector`.
+
+Temperature and salinity are not evolved in this task and are given
+constant values determined by config options `temperature` and `salinity`.
+
+The Coriolis parameters `fCell`, `fEdge`, and `fVertex` do not need to be
+specified for a global mesh and are initialized as zeros.
+
+## forcing
+
+This flow velocity case is forced to follow the constant rotation rate given
+in the config options.
+
+## time step and run duration
+
+This task uses the Runge-Kutta 4th-order (RK4) time integrator. The time step 
+for forward integration is determined by multiplying the resolution by a config
+option, `rk4_dt_per_km`, so that coarser meshes have longer time steps. You can
+alter this before setup (in a user config file) or before running the task (in 
+the config file in the work directory). 
+
+```cfg
+# config options for spherical convergence tests
+[spherical_convergence_forward]
+
+# time integrator: {'split_explicit', 'RK4'}
+time_integrator = RK4
+
+# RK4 time step per resolution (s/km), since dt is proportional to resolution
+rk4_dt_per_km = 3.0
+```
+
+The `convergence_eval_time`, `run_duration` and `output_interval` are the
+period for advection to make a full rotation around the globe, 12 days:
+
+```cfg
+# config options for spherical convergence tests
+[spherical_convergence_forward]
+
+# Run duration in days
+run_duration = ${sphere_transport:vel_pd}
+
+# Output interval in days
+output_interval = ${sphere_transport:vel_pd}
+```
+
+Here, `${sphere_transport:vel_pd}` means that the same value is used as in the 
+option `vel_pd` in section `[sphere_transport]`, see below.
+
+## config options
+
+The `rotation_2d` config options include:
+
+```cfg
+# options for all sphere transport test cases
+[sphere_transport]
+
+# temperature
+temperature = 15.
+
+# salinity
+salinity = 35.
+
+# time (days) for bell to transit equator once
+vel_pd = 12.0
+
+
+# options for tracer visualization for the sphere transport test case
+[sphere_transport_viz]
+
+# visualization latitude and longitude resolution
+dlon = 0.5
+dlat = 0.5
+
+# remapping method ('bilinear', 'neareststod', 'conserve')
+remap_method = conserve
+
+
+# options for tracer visualization for the sphere transport test case
+[sphere_transport_viz_tracer]
+
+# colormap options
+# colormap
+colormap_name = viridis
+
+# the type of norm used in the colormap
+norm_type = linear
+
+# A dictionary with keywords for the norm
+norm_args = {'vmin': 0., 'vmax': 1.}
+
+# We could provide colorbar tick marks but we'll leave the defaults
+# colorbar_ticks = np.linspace(0., 1., 9)
+
+# options for plotting tracer differences from sphere transport tests
+[sphere_transport_viz_tracer_diff]
+
+# colormap options
+# colormap
+colormap_name = cmo.balance
+
+# the type of norm used in the colormap
+norm_type = linear
+
+# A dictionary with keywords for the norm
+norm_args = {'vmin': -0.25, 'vmax': 0.25}
+
+
+# options for thickness visualization for the sphere transport test case
+[sphere_transport_viz_h]
+
+# colormap options
+# colormap
+colormap_name = viridis
+
+# the type of norm used in the colormap
+norm_type = linear
+
+# A dictionary with keywords for the norm
+norm_args = {'vmin': 99., 'vmax': 101.}
+
+# We could provide colorbar tick marks but we'll leave the defaults
+# colorbar_ticks = np.linspace(0., 1., 9)
+
+# options for plotting tracer differences from sphere transport tests
+[sphere_transport_viz_h_diff]
+
+# colormap options
+# colormap
+colormap_name = cmo.balance
+
+# the type of norm used in the colormap
+norm_type = linear
+
+# A dictionary with keywords for the norm
+norm_args = {'vmin': -0.25, 'vmax': 0.25}
+
+# options for rotation 2-d test case
+[rotation_2d]
+
+# rotation vector in cartesian coordinates
+rotation_vector = 0.2, 0.7, 1.0
+
+# convergence threshold below which the test fails
+convergence_thresh_tracer1 = 1.4
+convergence_thresh_tracer2 = 1.8
+convergence_thresh_tracer3 = 0.4
+```
+
+The options in section `sphere_transport` are used by all 4 test cases based
+on Lauritzen et al. (2012) and control the initial condition. The options in
+section `rotation_2d` control the convergence rate threshold.
+
+The options in sections `sphere_transport_viz*` control properties of the `viz`
+step of the test case. 
+
+The default options for the convergence analysis step can be changed here:
+
+```cfg
+# config options for spherical convergence tests
+[spherical_convergence]
+
+# Evaluation time for convergence analysis (in days)
+convergence_eval_time = ${sphere_transport:vel_pd}
+
+# Type of error to compute
+error_type = l2
+```
+
+## cores
+
+The target and minimum number of cores are determined by `goal_cells_per_core`
+and `max_cells_per_core` from the `ocean` section of the config file,
+respectively. This ensures that the number of cells per core is roughly
+constant across the different resolutions in the convergence study.
diff --git a/polaris/ocean/__init__.py b/polaris/ocean/__init__.py
index 6de0e5811..90d5f07e5 100644
--- a/polaris/ocean/__init__.py
+++ b/polaris/ocean/__init__.py
@@ -9,6 +9,7 @@
     add_manufactured_solution_tasks,
 )
 from polaris.ocean.tasks.single_column import add_single_column_tasks
+from polaris.ocean.tasks.sphere_transport import add_sphere_transport_tasks
 
 
 class Ocean(Component):
@@ -29,6 +30,7 @@ def __init__(self):
 
         # single column
         add_single_column_tasks(component=self)
+        add_sphere_transport_tasks(component=self)
 
         # spherical: please keep these in alphabetical order
         add_cosine_bell_tasks(component=self)
diff --git a/polaris/ocean/suites/convergence.txt b/polaris/ocean/suites/convergence.txt
index 97a07c5e7..56a8e63c7 100644
--- a/polaris/ocean/suites/convergence.txt
+++ b/polaris/ocean/suites/convergence.txt
@@ -1,5 +1,10 @@
 ocean/planar/inertial_gravity_wave
 ocean/planar/manufactured_solution
+ocean/spherical/icos/correlated_tracers_2d
+  cached: icos_base_mesh_60km icos_base_mesh_120km icos_base_mesh_240km icos_base_mesh_480km
+ocean/spherical/qu/correlated_tracers_2d
+  cached: qu_base_mesh_60km qu_base_mesh_90km qu_base_mesh_120km qu_base_mesh_150km
+  cached: qu_base_mesh_180km qu_base_mesh_210km qu_base_mesh_240km
 ocean/spherical/icos/cosine_bell
   cached: icos_base_mesh_60km icos_init_60km icos_base_mesh_120km icos_init_120km
   cached: icos_base_mesh_240km icos_init_240km icos_base_mesh_480km icos_init_480km
@@ -16,3 +21,18 @@ ocean/spherical/qu/geostrophic
   cached: qu_base_mesh_120km qu_init_120km qu_base_mesh_150km qu_init_150km
   cached: qu_base_mesh_180km qu_init_180km qu_base_mesh_210km qu_init_210km
   cached: qu_base_mesh_240km qu_init_240km
+ocean/spherical/icos/divergent_2d
+  cached: icos_base_mesh_60km icos_base_mesh_120km icos_base_mesh_240km icos_base_mesh_480km
+ocean/spherical/qu/divergent_2d
+  cached: qu_base_mesh_60km qu_base_mesh_90km qu_base_mesh_120km qu_base_mesh_150km
+  cached: qu_base_mesh_180km qu_base_mesh_210km qu_base_mesh_240km
+ocean/spherical/icos/nondivergent_2d
+  cached: icos_base_mesh_60km icos_base_mesh_120km icos_base_mesh_240km icos_base_mesh_480km
+ocean/spherical/qu/nondivergent_2d
+  cached: qu_base_mesh_60km qu_base_mesh_90km qu_base_mesh_120km qu_base_mesh_150km
+  cached: qu_base_mesh_180km qu_base_mesh_210km qu_base_mesh_240km
+ocean/spherical/icos/rotation_2d
+  cached: icos_base_mesh_60km icos_base_mesh_120km icos_base_mesh_240km icos_base_mesh_480km
+ocean/spherical/qu/rotation_2d
+  cached: qu_base_mesh_60km qu_base_mesh_90km qu_base_mesh_120km qu_base_mesh_150km
+  cached: qu_base_mesh_180km qu_base_mesh_210km qu_base_mesh_240km
diff --git a/polaris/ocean/suites/sphere_transport.txt b/polaris/ocean/suites/sphere_transport.txt
new file mode 100644
index 000000000..a93d53fb8
--- /dev/null
+++ b/polaris/ocean/suites/sphere_transport.txt
@@ -0,0 +1,4 @@
+ocean/spherical/icos/rotation_2d
+ocean/spherical/icos/nondivergent_2d
+ocean/spherical/icos/divergent_2d
+ocean/spherical/icos/correlated_tracers_2d
diff --git a/polaris/ocean/suites/sphere_transport_with_viz.txt b/polaris/ocean/suites/sphere_transport_with_viz.txt
new file mode 100644
index 000000000..b120cc1f5
--- /dev/null
+++ b/polaris/ocean/suites/sphere_transport_with_viz.txt
@@ -0,0 +1,4 @@
+ocean/spherical/icos/rotation_2d/with_viz
+ocean/spherical/icos/nondivergent_2d/with_viz
+ocean/spherical/icos/divergent_2d/with_viz
+ocean/spherical/icos/correlated_tracers_2d/with_viz
diff --git a/polaris/ocean/tasks/cosine_bell/init.py b/polaris/ocean/tasks/cosine_bell/init.py
index 4c3b26ef9..788f0b6cf 100644
--- a/polaris/ocean/tasks/cosine_bell/init.py
+++ b/polaris/ocean/tasks/cosine_bell/init.py
@@ -124,5 +124,10 @@ def cosine_bell(max_value, ri, r):
 
     r : float
         Radius of the cosine bell in meters
+
+    Returns
+    -------
+    f : np.ndarray of type float
+        Cosine bell tracer values
     """
     return max_value / 2.0 * (1.0 + np.cos(np.pi * np.divide(ri, r)))
diff --git a/polaris/ocean/tasks/sphere_transport/__init__.py b/polaris/ocean/tasks/sphere_transport/__init__.py
new file mode 100644
index 000000000..893e05ee8
--- /dev/null
+++ b/polaris/ocean/tasks/sphere_transport/__init__.py
@@ -0,0 +1,253 @@
+from typing import Dict
+
+from polaris import Step, Task
+from polaris.config import PolarisConfigParser
+from polaris.ocean.mesh.spherical import add_spherical_base_mesh_step
+from polaris.ocean.tasks.sphere_transport.analysis import Analysis
+from polaris.ocean.tasks.sphere_transport.filament_analysis import (
+    FilamentAnalysis,
+)
+from polaris.ocean.tasks.sphere_transport.forward import Forward
+from polaris.ocean.tasks.sphere_transport.init import Init
+from polaris.ocean.tasks.sphere_transport.mixing_analysis import MixingAnalysis
+from polaris.ocean.tasks.sphere_transport.viz import Viz, VizMap
+
+
+def add_sphere_transport_tasks(component):
+    """
+    Add tasks that define variants of sphere transport test cases:
+    nondivergent_2d, divergent_2d, correlated_tracers_2d, rotation_2d
+
+    component : polaris.ocean.Ocean
+        the ocean component that the tasks will be added to
+    """
+
+    for icosahedral, prefix in [(True, 'icos'), (False, 'qu')]:
+        for case_name in ['rotation_2d', 'nondivergent_2d', 'divergent_2d',
+                          'correlated_tracers_2d']:
+            filepath = f'spherical/{prefix}/{case_name}/{case_name}.cfg'
+            config = PolarisConfigParser(filepath=filepath)
+            config.add_from_package('polaris.ocean.convergence',
+                                    'convergence.cfg')
+            config.add_from_package('polaris.ocean.convergence.spherical',
+                                    'spherical.cfg')
+            package = 'polaris.ocean.tasks.sphere_transport'
+            config.add_from_package(package, 'sphere_transport.cfg')
+            config.add_from_package(package, f'{case_name}.cfg')
+            for include_viz in [False, True]:
+                component.add_task(SphereTransport(
+                    component=component, case_name=case_name,
+                    icosahedral=icosahedral, config=config,
+                    include_viz=include_viz))
+
+
+class SphereTransport(Task):
+    """
+    A test case for testing properties of tracer advection
+
+    Attributes
+    ----------
+    resolutions : list of float
+        A list of mesh resolutions
+
+    icosahedral : bool
+        Whether to use icosahedral, as opposed to less regular, JIGSAW meshes
+
+    include_viz : bool
+        Include VizMap and Viz steps for each resolution
+    """
+    def __init__(self, component, config, case_name, icosahedral, include_viz):
+        """
+        Create test case for creating a global MPAS-Ocean mesh
+
+        Parameters
+        ----------
+        component : polaris.ocean.Ocean
+            The ocean component that this task belongs to
+
+        config : polaris.config.PolarisConfigParser
+            A shared config parser
+
+        case_name: string
+            The name of the case which determines what variant of the
+            configuration to use
+
+        icosahedral : bool
+            Whether to use icosahedral, as opposed to less regular, JIGSAW
+            meshes
+
+        include_viz : bool
+            Include VizMap and Viz steps for each resolution
+        """
+        if icosahedral:
+            prefix = 'icos'
+        else:
+            prefix = 'qu'
+
+        subdir = f'spherical/{prefix}/{case_name}'
+        name = f'{prefix}_{case_name}'
+        if include_viz:
+            subdir = f'{subdir}/with_viz'
+            name = f'{name}_with_viz'
+            link = f'{case_name}.cfg'
+        else:
+            # config options live in the task already so no need for a symlink
+            link = None
+        super().__init__(component=component, name=name, subdir=subdir)
+        self.resolutions = list()
+        self.icosahedral = icosahedral
+        self.case_name = case_name
+        self.include_viz = include_viz
+
+        self.set_shared_config(config, link=link)
+
+        # add the steps with default resolutions so they can be listed
+        self._setup_steps()
+
+    def configure(self):
+        """
+        Set config options for the test case
+        """
+        super().configure()
+
+        # set up the steps again in case a user has provided new resolutions
+        self._setup_steps()
+
+    def _setup_steps(self):
+        """ setup steps given resolutions """
+        case_name = self.case_name
+        icosahedral = self.icosahedral
+        config = self.config
+        config_filename = self.config_filename
+
+        if icosahedral:
+            prefix = 'icos'
+        else:
+            prefix = 'qu'
+
+        resolutions = config.getlist('spherical_convergence',
+                                     f'{prefix}_resolutions', dtype=float)
+
+        if self.resolutions == resolutions:
+            return
+
+        # start fresh with no steps
+        for step in list(self.steps.values()):
+            self.remove_step(step)
+
+        self.resolutions = resolutions
+
+        component = self.component
+
+        analysis_dependencies: Dict[str, Dict[str, Step]] = (
+            dict(mesh=dict(), init=dict(), forward=dict()))
+        for resolution in resolutions:
+            base_mesh_step, mesh_name = add_spherical_base_mesh_step(
+                component, resolution, icosahedral)
+            self.add_step(base_mesh_step, symlink=f'base_mesh/{mesh_name}')
+            analysis_dependencies['mesh'][resolution] = base_mesh_step
+
+            sph_trans_dir = f'spherical/{prefix}/{case_name}'
+
+            name = f'{prefix}_init_{mesh_name}'
+            subdir = f'{sph_trans_dir}/init/{mesh_name}'
+            if self.include_viz:
+                symlink = f'init/{mesh_name}'
+            else:
+                symlink = None
+            if subdir in component.steps:
+                init_step = component.steps[subdir]
+            else:
+                init_step = Init(component=component, name=name, subdir=subdir,
+                                 base_mesh=base_mesh_step, case_name=case_name)
+                init_step.set_shared_config(config, link=config_filename)
+            self.add_step(init_step, symlink=symlink)
+            analysis_dependencies['init'][resolution] = init_step
+
+            name = f'{prefix}_forward_{mesh_name}'
+            subdir = f'{sph_trans_dir}/forward/{mesh_name}'
+            if self.include_viz:
+                symlink = f'forward/{mesh_name}'
+            else:
+                symlink = None
+            if subdir in component.steps:
+                forward_step = component.steps[subdir]
+            else:
+                forward_step = Forward(component=component, name=name,
+                                       subdir=subdir, resolution=resolution,
+                                       base_mesh=base_mesh_step,
+                                       init=init_step,
+                                       case_name=case_name)
+                forward_step.set_shared_config(config, link=config_filename)
+            self.add_step(forward_step, symlink=symlink)
+            analysis_dependencies['forward'][resolution] = forward_step
+
+            if self.include_viz:
+                with_viz_dir = f'{sph_trans_dir}/with_viz'
+
+                name = f'{prefix}_map_{mesh_name}'
+                subdir = f'{with_viz_dir}/map/{mesh_name}'
+                viz_map = VizMap(component=component, name=name,
+                                 subdir=subdir, base_mesh=base_mesh_step,
+                                 mesh_name=mesh_name)
+                viz_map.set_shared_config(config, link=config_filename)
+                self.add_step(viz_map)
+
+                name = f'{prefix}_viz_{mesh_name}'
+                subdir = f'{with_viz_dir}/viz/{mesh_name}'
+                step = Viz(component=component, name=name,
+                           subdir=subdir, base_mesh=base_mesh_step,
+                           init=init_step, forward=forward_step,
+                           viz_map=viz_map, mesh_name=mesh_name)
+                step.set_shared_config(config, link=config_filename)
+                self.add_step(step)
+
+        subdir = f'{sph_trans_dir}/analysis'
+        if self.include_viz:
+            symlink = 'analysis'
+        else:
+            symlink = None
+        if subdir in component.steps:
+            step = component.steps[subdir]
+            step.resolutions = resolutions
+            step.dependencies_dict = analysis_dependencies
+        else:
+            step = Analysis(component=component, resolutions=resolutions,
+                            subdir=subdir, case_name=case_name,
+                            dependencies=analysis_dependencies)
+            step.set_shared_config(config, link=config_filename)
+        self.add_step(step, symlink=symlink)
+
+        if case_name == 'correlated_tracers_2d':
+            subdir = f'{sph_trans_dir}/mixing_analysis'
+            if self.include_viz:
+                symlink = 'mixing_analysis'
+            else:
+                symlink = None
+            if subdir in component.steps:
+                step = component.steps[subdir]
+            else:
+                step = MixingAnalysis(component=component,
+                                      resolutions=resolutions,
+                                      icosahedral=icosahedral, subdir=subdir,
+                                      case_name=case_name,
+                                      dependencies=analysis_dependencies)
+            step.set_shared_config(config, link=config_filename)
+            self.add_step(step, symlink=symlink)
+
+        if case_name == 'nondivergent_2d':
+            subdir = f'{sph_trans_dir}/filament_analysis'
+            if self.include_viz:
+                symlink = 'filament_analysis'
+            else:
+                symlink = None
+            if subdir in component.steps:
+                step = component.steps[subdir]
+            else:
+                step = FilamentAnalysis(component=component,
+                                        resolutions=resolutions,
+                                        icosahedral=icosahedral, subdir=subdir,
+                                        case_name=case_name,
+                                        dependencies=analysis_dependencies)
+            step.set_shared_config(config, link=config_filename)
+            self.add_step(step, symlink=symlink)
diff --git a/polaris/ocean/tasks/sphere_transport/analysis.py b/polaris/ocean/tasks/sphere_transport/analysis.py
new file mode 100644
index 000000000..620a80651
--- /dev/null
+++ b/polaris/ocean/tasks/sphere_transport/analysis.py
@@ -0,0 +1,79 @@
+from polaris.ocean.convergence import ConvergenceAnalysis
+
+
+class Analysis(ConvergenceAnalysis):
+    """
+    A step for analyzing the output from sphere transport test cases
+
+    Attributes
+    ----------
+    resolutions : list of float
+        The resolutions of the meshes that have been run
+
+    case_name : str
+        The name of the test case
+    """
+    def __init__(self, component, resolutions, subdir, case_name,
+                 dependencies):
+        """
+        Create the step
+
+        Parameters
+        ----------
+        component : polaris.Component
+            The component the step belongs to
+
+        resolutions : list of float
+            The resolutions of the meshes that have been run
+
+        subdir : str
+            The subdirectory that the step resides in
+
+        case_name: str
+            The name of the test case
+
+        dependencies : dict of dict of polaris.Steps
+            The dependencies of this step
+        """
+        self.case_name = case_name
+        convergence_vars = [{'name': 'tracer1',
+                             'title': 'tracer1',
+                             'zidx': 1},
+                            {'name': 'tracer2',
+                             'title': 'tracer2',
+                             'zidx': 1},
+                            {'name': 'tracer3',
+                             'title': 'tracer3',
+                             'zidx': 1}]
+        super().__init__(component=component, subdir=subdir,
+                         resolutions=resolutions,
+                         dependencies=dependencies,
+                         convergence_vars=convergence_vars)
+        # Note: there is no need to overwrite the default method exact_solution
+        # which uses the initial condition
+
+    def convergence_parameters(self, field_name=None):
+        """
+        Get convergence parameters
+
+        Parameters
+        ----------
+        field_name : str
+            The name of the variable of which we evaluate convergence
+            For cosine_bell, we use the same convergence rate for all fields
+        Returns
+        -------
+        conv_thresh: float
+            The minimum convergence rate
+
+        conv_thresh: float
+            The maximum convergence rate
+        """
+        config = self.config
+        section = config[self.case_name]
+        conv_thresh = section.getfloat(f'convergence_thresh_{field_name}')
+
+        section = config['convergence']
+        error_type = section.get('error_type')
+
+        return conv_thresh, error_type
diff --git a/polaris/ocean/tasks/sphere_transport/correlated_tracers_2d.cfg b/polaris/ocean/tasks/sphere_transport/correlated_tracers_2d.cfg
new file mode 100644
index 000000000..35679737a
--- /dev/null
+++ b/polaris/ocean/tasks/sphere_transport/correlated_tracers_2d.cfg
@@ -0,0 +1,17 @@
+# options for correlated_tracers_2d test case
+[correlated_tracers_2d]
+
+# velocity amplitude in meters per second
+vel_amp = 10.
+
+# Quadratic coefficients for the correlated tracer function in decreasing
+# order of the terms
+correlation_coefficients = -0.8, 0.0, 0.9
+
+# convergence threshold below which the test fails
+convergence_thresh_tracer1 = 1.5
+convergence_thresh_tracer2 = 1.3
+convergence_thresh_tracer3 = 1.0
+
+# time in days at which to evaluate mixing
+mixing_evaluation_time = 6.0
diff --git a/polaris/ocean/tasks/sphere_transport/divergent_2d.cfg b/polaris/ocean/tasks/sphere_transport/divergent_2d.cfg
new file mode 100644
index 000000000..e442e313a
--- /dev/null
+++ b/polaris/ocean/tasks/sphere_transport/divergent_2d.cfg
@@ -0,0 +1,10 @@
+# options for divergent_2d test case
+[divergent_2d]
+
+# velocity amplitude in meters per second
+vel_amp = 5.
+
+# convergence threshold below which the test fails
+convergence_thresh_tracer1 = 1.8
+convergence_thresh_tracer2 = 1.5
+convergence_thresh_tracer3 = 0.3
diff --git a/polaris/ocean/tasks/sphere_transport/filament_analysis.py b/polaris/ocean/tasks/sphere_transport/filament_analysis.py
new file mode 100644
index 000000000..a3f196439
--- /dev/null
+++ b/polaris/ocean/tasks/sphere_transport/filament_analysis.py
@@ -0,0 +1,123 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+import xarray as xr
+
+from polaris import Step
+from polaris.mpas import time_index_from_xtime
+from polaris.ocean.resolution import resolution_to_subdir
+from polaris.viz import use_mplstyle
+
+
+class FilamentAnalysis(Step):
+    """
+    A step for analyzing the output from sphere transport test cases
+
+    Attributes
+    ----------
+    resolutions : list of float
+        The resolutions of the meshes that have been run
+
+    icosahedral : bool
+        Whether to use icosahedral, as opposed to less regular, JIGSAW
+        meshes
+
+    case_name : str
+        The name of the test case
+    """
+    def __init__(self, component, resolutions, icosahedral, subdir,
+                 case_name, dependencies):
+        """
+        Create the step
+
+        Parameters
+        ----------
+        component : polaris.Component
+            The component the step belongs to
+
+        resolutions : list of float
+            The resolutions of the meshes that have been run
+
+        icosahedral : bool
+            Whether to use icosahedral, as opposed to less regular, JIGSAW
+            meshes
+
+        subdir : str
+            The subdirectory that the step resides in
+
+        case_name: str
+            The name of the test case
+
+        dependencies : dict of dict of polaris.Steps
+            The dependencies of this step
+        """
+        super().__init__(component=component, name='filament_analysis',
+                         subdir=subdir)
+        self.resolutions = resolutions
+        self.case_name = case_name
+
+        for resolution in resolutions:
+            mesh_name = resolution_to_subdir(resolution)
+            base_mesh = dependencies['mesh'][resolution]
+            init = dependencies['init'][resolution]
+            forward = dependencies['forward'][resolution]
+            self.add_input_file(
+                filename=f'{mesh_name}_mesh.nc',
+                work_dir_target=f'{base_mesh.path}/base_mesh.nc')
+            self.add_input_file(
+                filename=f'{mesh_name}_init.nc',
+                work_dir_target=f'{init.path}/initial_state.nc')
+            self.add_input_file(
+                filename=f'{mesh_name}_output.nc',
+                work_dir_target=f'{forward.path}/output.nc')
+        self.add_output_file('filament.png')
+
+    def run(self):
+        """
+        Run this step of the test case
+        """
+        plt.switch_backend('Agg')
+        resolutions = self.resolutions
+        config = self.config
+        section = config[self.case_name]
+        eval_time = section.getfloat('filament_evaluation_time')
+        s_per_day = 86400.0
+        zidx = 1
+        variable_name = 'tracer2'
+        num_tau = 21
+        filament_tau = np.linspace(0, 1, num_tau)
+        filament_norm = np.zeros((len(resolutions), num_tau))
+        use_mplstyle()
+        fig, ax = plt.subplots()
+        for i, resolution in enumerate(resolutions):
+            mesh_name = resolution_to_subdir(resolution)
+            ds = xr.open_dataset(f'{mesh_name}_output.nc')
+            tidx = time_index_from_xtime(ds.xtime.values,
+                                         eval_time * s_per_day)
+            tracer = ds[variable_name]
+            area_cell = ds["areaCell"]
+            for j, tau in enumerate(filament_tau):
+                cells_above_tau = tracer[tidx, :, zidx] >= tau
+                cells_above_tau0 = tracer[0, :, zidx] >= tau
+                if np.sum(cells_above_tau0 * area_cell) == 0.:
+                    filament_norm[i, j] = np.nan
+                else:
+                    filament_norm[i, j] = np.divide(
+                        np.sum(area_cell * cells_above_tau),
+                        np.sum(cells_above_tau0 * area_cell))
+            plt.plot(filament_tau, filament_norm[i, :], '.-', label=mesh_name)
+        plt.plot([filament_tau[0], filament_tau[-1]], [1., 1.], 'k--')
+        ax.set_xlim([filament_tau[0], filament_tau[-1]])
+        ax.set_xlabel(r'$\tau$')
+        ax.set_ylabel(r'$l_f$')
+        plt.title(f'Filament preservation diagnostic for {variable_name}')
+        plt.legend()
+        fig.savefig('filament.png', bbox_inches='tight')
+
+        res_array = np.array(resolutions, dtype=float)
+        data = np.column_stack((res_array, filament_norm))
+        col_headers = ['resolution']
+        for tau in filament_tau:
+            col_headers.append(f'{tau:g}')
+        df = pd.DataFrame(data, columns=col_headers)
+        df.to_csv('filament.csv', index=False)
diff --git a/polaris/ocean/tasks/sphere_transport/forward.py b/polaris/ocean/tasks/sphere_transport/forward.py
new file mode 100644
index 000000000..98007a458
--- /dev/null
+++ b/polaris/ocean/tasks/sphere_transport/forward.py
@@ -0,0 +1,52 @@
+from polaris.ocean.convergence.spherical import SphericalConvergenceForward
+
+
+class Forward(SphericalConvergenceForward):
+    """
+    A step for performing forward ocean component runs as part of the sphere
+    transport test case
+    """
+
+    def __init__(self, component, name, subdir, resolution, base_mesh, init,
+                 case_name):
+        """
+        Create a new step
+
+        Parameters
+        ----------
+        component : polaris.Component
+            The component the step belongs to
+
+        name : str
+            The name of the step
+
+        subdir : str
+            The subdirectory for the step
+
+        resolution : float
+            The resolution of the (uniform) mesh in km
+
+        base_mesh : polaris.Step
+            The base mesh step
+
+        init : polaris.Step
+            The init step
+
+        case_name: str
+            The name of the test case
+        """
+        package = 'polaris.ocean.tasks.sphere_transport'
+        flow_id = {'rotation_2d': 1,
+                   'nondivergent_2d': 2,
+                   'divergent_2d': 3,
+                   'correlated_tracers_2d': 4}
+        namelist_options = dict(
+            config_transport_tests_flow_id=flow_id[case_name])
+        validate_vars = ['normalVelocity', 'tracer1', 'tracer2', 'tracer3']
+        super().__init__(component=component, name=name, subdir=subdir,
+                         resolution=resolution, mesh=base_mesh,
+                         init=init, package=package,
+                         yaml_filename='forward.yaml',
+                         output_filename='output.nc',
+                         validate_vars=validate_vars,
+                         options=namelist_options)
diff --git a/polaris/ocean/tasks/sphere_transport/forward.yaml b/polaris/ocean/tasks/sphere_transport/forward.yaml
new file mode 100644
index 000000000..c870f49aa
--- /dev/null
+++ b/polaris/ocean/tasks/sphere_transport/forward.yaml
@@ -0,0 +1,56 @@
+omega:
+  run_modes:
+    config_ocean_run_mode: forward
+  time_management:
+    config_run_duration: {{ run_duration }}
+  decomposition:
+    config_block_decomp_file_prefix: graph.info.part.
+  advection:
+    config_vert_coord_movement: impermeable_interfaces
+  time_integration:
+    config_dt: {{ dt }}
+    config_time_integrator: {{ time_integrator }}
+  debug:
+    config_disable_thick_sflux: true
+    config_disable_vel_all_tend: true
+    config_disable_vel_coriolis: true
+    config_disable_vel_pgrad: true
+    config_disable_vel_hmix: true
+    config_disable_vel_surface_stress: true
+    config_disable_vel_explicit_bottom_drag: true
+    config_disable_vel_vmix: true
+    config_disable_vel_vadv: true
+    config_disable_tr_hmix: true
+    config_disable_tr_vmix: true
+    config_disable_tr_sflux: true
+    config_disable_tr_nonlocalflux: true
+    config_check_ssh_consistency: false
+  cvmix:
+    config_use_cvmix: false
+  eos:
+    config_eos_type: linear
+  tracer_forcing_debugTracers:
+    config_use_debugTracers: true
+  streams:
+    mesh:
+      filename_template: init.nc
+    input:
+      filename_template: init.nc
+    restart:
+      output_interval: 0030_00:00:00
+    output:
+      type: output
+      filename_template: output.nc
+      output_interval: {{ output_interval }}
+      clobber_mode: truncate
+      reference_time: 0001-01-01_00:00:00
+      contents:
+      - tracers
+      - mesh
+      - xtime
+      - normalVelocity
+      - layerThickness
+      - refZMid
+      - refLayerThickness
+      - kineticEnergyCell
+      - relativeVorticityCell
diff --git a/polaris/ocean/tasks/sphere_transport/init.py b/polaris/ocean/tasks/sphere_transport/init.py
new file mode 100644
index 000000000..815921f6c
--- /dev/null
+++ b/polaris/ocean/tasks/sphere_transport/init.py
@@ -0,0 +1,159 @@
+import numpy as np
+import xarray as xr
+from mpas_tools.io import write_netcdf
+
+from polaris import Step
+from polaris.ocean.tasks.sphere_transport.resources.flow_types import (
+    flow_divergent,
+    flow_nondivergent,
+    flow_rotation,
+)
+from polaris.ocean.tasks.sphere_transport.resources.tracer_distributions import (  # noqa: E501
+    correlation_fn,
+    cosine_bells,
+    slotted_cylinders,
+    xyztrig,
+)
+from polaris.ocean.vertical import init_vertical_coord
+
+
+class Init(Step):
+    """
+    A step for an initial condition for for the cosine bell test case
+    """
+    def __init__(self, component, name, subdir, base_mesh, case_name):
+        """
+        Create the step
+
+        Parameters
+        ----------
+        component : polaris.Component
+            The component the step belongs to
+
+        name : str
+            The name of the step
+
+        subdir : str
+            The subdirectory for the step
+
+        base_mesh : polaris.Step
+            The base mesh step
+
+        case_name: str
+            The name of the test case
+        """
+        super().__init__(component=component, name=name, subdir=subdir)
+
+        self.case_name = case_name
+        self.add_input_file(
+            filename='mesh.nc',
+            work_dir_target=f'{base_mesh.path}/base_mesh.nc')
+
+        self.add_input_file(
+            filename='graph.info',
+            work_dir_target=f'{base_mesh.path}/graph.info')
+
+        self.add_output_file(filename='initial_state.nc')
+
+    def run(self):
+        """
+        Run this step of the task
+        """
+        config = self.config
+        case_name = self.case_name
+
+        section = config['sphere_transport']
+        temperature = section.getfloat('temperature')
+        salinity = section.getfloat('salinity')
+        vel_pd = section.getfloat('vel_pd')
+
+        section = config['vertical_grid']
+        bottom_depth = section.getfloat('bottom_depth')
+
+        ds_mesh = xr.open_dataset('mesh.nc')
+        angleEdge = ds_mesh.angleEdge
+        latCell = ds_mesh.latCell
+        latEdge = ds_mesh.latEdge
+        lonCell = ds_mesh.lonCell
+        lonEdge = ds_mesh.lonEdge
+        sphere_radius = ds_mesh.sphere_radius
+
+        ds = ds_mesh.copy()
+
+        ds['bottomDepth'] = bottom_depth * xr.ones_like(latCell)
+        ds['ssh'] = xr.zeros_like(latCell)
+
+        init_vertical_coord(config, ds)
+
+        temperature_array = temperature * xr.ones_like(latCell)
+        temperature_array, _ = xr.broadcast(temperature_array, ds.refZMid)
+        ds['temperature'] = temperature_array.expand_dims(dim='Time', axis=0)
+        ds['salinity'] = salinity * xr.ones_like(ds.temperature)
+
+        # tracer1
+        tracer1 = xyztrig(lonCell, latCell, sphere_radius)
+
+        # tracer2
+        section = config['sphere_transport']
+        radius = section.getfloat('cosine_bells_radius')
+        background_value = section.getfloat('cosine_bells_background')
+        amplitude = section.getfloat('cosine_bells_amplitude')
+        tracer2 = cosine_bells(lonCell, latCell, radius, background_value,
+                               amplitude, sphere_radius)
+
+        # tracer3
+        if case_name == 'correlated_tracers_2d':
+            coeff = config.getlist(case_name, 'correlation_coefficients',
+                                   dtype=float)
+            tracer3 = correlation_fn(tracer2, coeff[0], coeff[1], coeff[2])
+        else:
+            section = config['sphere_transport']
+            radius = section.getfloat('slotted_cylinders_radius')
+            background_value = section.getfloat('slotted_cylinders_background')
+            amplitude = section.getfloat('slotted_cylinders_amplitude')
+            tracer3 = slotted_cylinders(lonCell, latCell, radius,
+                                        background_value, amplitude,
+                                        sphere_radius)
+        _, tracer1_array = np.meshgrid(ds.refZMid.values, tracer1)
+        _, tracer2_array = np.meshgrid(ds.refZMid.values, tracer2)
+        _, tracer3_array = np.meshgrid(ds.refZMid.values, tracer3)
+
+        ds['tracer1'] = (('nCells', 'nVertLevels',), tracer1_array)
+        ds['tracer1'] = ds.tracer1.expand_dims(dim='Time', axis=0)
+        ds['tracer2'] = (('nCells', 'nVertLevels',), tracer2_array)
+        ds['tracer2'] = ds.tracer2.expand_dims(dim='Time', axis=0)
+        ds['tracer3'] = (('nCells', 'nVertLevels',), tracer3_array)
+        ds['tracer3'] = ds.tracer3.expand_dims(dim='Time', axis=0)
+
+        # Initialize velocity
+        s_per_hour = 3600.
+        if case_name == 'rotation_2d':
+            rotation_vector = config.getlist(case_name, 'rotation_vector',
+                                             dtype=float)
+            vector = np.array(rotation_vector)
+            u, v = flow_rotation(lonEdge, latEdge, vector,
+                                 vel_pd * s_per_hour, sphere_radius)
+        elif case_name == 'divergent_2d':
+            section = config[case_name]
+            vel_amp = section.getfloat('vel_amp')
+            u, v = flow_divergent(0., lonEdge, latEdge,
+                                  vel_amp, vel_pd * s_per_hour)
+        elif (case_name == 'nondivergent_2d' or
+              case_name == 'correlated_tracers_2d'):
+            section = config[case_name]
+            vel_amp = section.getfloat('vel_amp')
+            u, v = flow_nondivergent(0., lonEdge, latEdge,
+                                     vel_amp, vel_pd * s_per_hour)
+        else:
+            raise ValueError(f'Unexpected test case name {case_name}')
+
+        normalVelocity = sphere_radius * (u * np.cos(angleEdge) +
+                                          v * np.sin(angleEdge))
+        normalVelocity, _ = xr.broadcast(normalVelocity, ds.refZMid)
+        ds['normalVelocity'] = normalVelocity.expand_dims(dim='Time', axis=0)
+
+        ds['fCell'] = xr.zeros_like(ds_mesh.xCell)
+        ds['fEdge'] = xr.zeros_like(ds_mesh.xEdge)
+        ds['fVertex'] = xr.zeros_like(ds_mesh.xVertex)
+
+        write_netcdf(ds, 'initial_state.nc')
diff --git a/polaris/ocean/tasks/sphere_transport/mixing_analysis.py b/polaris/ocean/tasks/sphere_transport/mixing_analysis.py
new file mode 100644
index 000000000..750518a7f
--- /dev/null
+++ b/polaris/ocean/tasks/sphere_transport/mixing_analysis.py
@@ -0,0 +1,149 @@
+from math import ceil
+
+import matplotlib.pyplot as plt
+import numpy as np
+import xarray as xr
+from matplotlib.lines import Line2D
+
+from polaris import Step
+from polaris.mpas import time_index_from_xtime
+from polaris.ocean.resolution import resolution_to_subdir
+from polaris.viz import use_mplstyle
+
+
+class MixingAnalysis(Step):
+    """
+    A step for analyzing the output from sphere transport test cases
+
+    Attributes
+    ----------
+    resolutions : list of float
+        The resolutions of the meshes that have been run
+
+    icosahedral : bool
+        Whether to use icosahedral, as opposed to less regular, JIGSAW
+        meshes
+
+    case_name : str
+        The name of the test case
+    """
+    def __init__(self, component, resolutions, icosahedral, subdir,
+                 case_name, dependencies):
+        """
+        Create the step
+
+        Parameters
+        ----------
+        component : polaris.Component
+            The component the step belongs to
+
+        resolutions : list of float
+            The resolutions of the meshes that have been run
+
+        icosahedral : bool
+            Whether to use icosahedral, as opposed to less regular, JIGSAW
+            meshes
+
+        subdir : str
+            The subdirectory that the step resides in
+
+        case_name: str
+            The name of the test case
+
+        dependencies : dict of dict of polaris.Steps
+            The dependencies of this step
+        """
+        super().__init__(component=component, name='mixing_analysis',
+                         subdir=subdir)
+        self.resolutions = resolutions
+        self.case_name = case_name
+
+        for resolution in resolutions:
+            mesh_name = resolution_to_subdir(resolution)
+            base_mesh = dependencies['mesh'][resolution]
+            init = dependencies['init'][resolution]
+            forward = dependencies['forward'][resolution]
+            self.add_input_file(
+                filename=f'{mesh_name}_mesh.nc',
+                work_dir_target=f'{base_mesh.path}/base_mesh.nc')
+            self.add_input_file(
+                filename=f'{mesh_name}_init.nc',
+                work_dir_target=f'{init.path}/initial_state.nc')
+            self.add_input_file(
+                filename=f'{mesh_name}_output.nc',
+                work_dir_target=f'{forward.path}/output.nc')
+        self.add_output_file('triplots.png')
+
+    def run(self):
+        """
+        Run this step of the test case
+        """
+        plt.switch_backend('Agg')
+        resolutions = self.resolutions
+        config = self.config
+        section = config[self.case_name]
+        eval_time = section.getfloat('mixing_evaluation_time')
+        s_per_day = 86400.0
+        zidx = 1
+        nrows = int(ceil(len(resolutions) / 2))
+        use_mplstyle()
+        fig, axes = plt.subplots(nrows=nrows, ncols=2, sharex=True,
+                                 sharey=True, figsize=(5.5, 7))
+        for i, resolution in enumerate(resolutions):
+            ax = axes[int(i / 2), i % 2]
+            _init_triplot_axes(ax)
+            mesh_name = resolution_to_subdir(resolution)
+            ax.set(title=mesh_name)
+            ds = xr.open_dataset(f'{mesh_name}_output.nc')
+            if i % 2 == 0:
+                ax.set_ylabel("tracer3")
+            if int(i / 2) == nrows - 1:
+                ax.set_xlabel("tracer2")
+            tidx = time_index_from_xtime(ds.xtime.values,
+                                         eval_time * s_per_day)
+            ds = ds.isel(Time=tidx)
+            ds = ds.isel(nVertLevels=zidx)
+            tracer2 = ds["tracer2"].values
+            tracer3 = ds["tracer3"].values
+            ax.plot(tracer2, tracer3, '.', markersize=1)
+            ax.set_aspect('equal')
+        if i % 2 < 1:
+            ax = axes[int(i / 2), 1]
+            ax.set_axis_off()
+        plt.subplots_adjust(wspace=0.1, hspace=0.1)
+        fig.suptitle('Correlated tracers 2-d')
+        fig.savefig('triplots.png', bbox_inches='tight')
+
+
+def _init_triplot_axes(ax):
+    lw = 0.4
+    topline = Line2D([0.1, 1.0], [0.9, 0.9], color='k',
+                     linestyle='-', linewidth=lw)
+    midline = Line2D([0.1, 1.0], [0.9, 0.1], color='k',
+                     linestyle='-', linewidth=lw)
+    rightline = Line2D([1, 1], [0.1, 0.9], color='k',
+                       linestyle='-', linewidth=lw)
+    leftline = Line2D([0.1, 0.1], [0.1, 0.9], color='k',
+                      linestyle='-', linewidth=lw)
+    botline = Line2D([0.1, 1.0], [0.1, 0.1], color='k',
+                     linestyle='-', linewidth=lw)
+    crvx = np.linspace(0.1, 1)
+    crvy = -0.8 * np.square(crvx) + 0.9
+    ticks = np.array(range(6)) / 5
+    ax.plot(crvx, crvy, 'k-', linewidth=1.25 * lw)
+    ax.set_xticks(ticks)
+    ax.set_yticks(ticks)
+    ax.add_artist(topline)
+    ax.add_artist(midline)
+    ax.add_artist(botline)
+    ax.add_artist(rightline)
+    ax.add_artist(leftline)
+    ax.set_xlim([0, 1.1])
+    ax.set_ylim([0, 1.0])
+    ax.text(0.98, 0.87, 'Range-preserving\n unmixing', fontsize=8,
+            horizontalalignment='right', verticalalignment='top')
+    ax.text(0.12, 0.12, 'Range-preserving\n unmixing', fontsize=8,
+            horizontalalignment='left', verticalalignment='bottom')
+    ax.text(0.5, 0.27, 'Real mixing', rotation=-40., fontsize=8)
+    ax.text(0.02, 0.1, 'Overshooting', rotation=90., fontsize=8)
+    ax.grid(color='lightgray')
diff --git a/polaris/ocean/tasks/sphere_transport/nondivergent_2d.cfg b/polaris/ocean/tasks/sphere_transport/nondivergent_2d.cfg
new file mode 100644
index 000000000..d78be5b13
--- /dev/null
+++ b/polaris/ocean/tasks/sphere_transport/nondivergent_2d.cfg
@@ -0,0 +1,13 @@
+# options for nondivergent_2d test case
+[nondivergent_2d]
+
+# velocity amplitude in meters per second
+vel_amp = 10.
+
+# convergence threshold below which the test fails
+convergence_thresh_tracer1 = 1.5
+convergence_thresh_tracer2 = 1.1
+convergence_thresh_tracer3 = 0.3
+
+# time in days at which to evaluate filament preservation
+filament_evaluation_time = 6.0
diff --git a/polaris/ocean/tasks/sphere_transport/resources/__init__.py b/polaris/ocean/tasks/sphere_transport/resources/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/polaris/ocean/tasks/sphere_transport/resources/flow_types.py b/polaris/ocean/tasks/sphere_transport/resources/flow_types.py
new file mode 100644
index 000000000..78a9adf95
--- /dev/null
+++ b/polaris/ocean/tasks/sphere_transport/resources/flow_types.py
@@ -0,0 +1,129 @@
+import numpy as np
+from mpas_tools.transects import lon_lat_to_cartesian
+from numpy import cos, pi, sin
+
+
+def flow_nondivergent(t, lon, lat, u_0, tau):
+    """
+    Compute a nondivergent velocity field
+
+    Parameters
+    ----------
+    t : np.ndarray of type float
+        times in seconds at which to compute the velocity field
+
+    lon : np.ndarray of type float
+        longitude
+
+    lat : np.ndarray of type float
+        latitude
+
+    u_0 : float
+        velocity amplitude in meters per second
+
+    tau : float
+        time in seconds for the flow to circumnavigate the sphere
+
+    Returns
+    -------
+    u : np.ndarray of type float
+       zonal velocity
+
+    v : np.ndarray of type float
+       meridional velocity
+    """
+    lon_p = lon - 2. * pi * t / tau
+    coslat = cos(lat)
+    cost = cos(pi * t / tau)
+    u = (1 / tau) * (u_0 * (sin(lon_p)**2) * sin(2 * lat) * cost +
+                     2. * pi * coslat)
+    v = (u_0 / tau) * sin(2 * lon_p) * coslat * cost
+    return u, v
+
+
+def flow_divergent(t, lon, lat, u_0, tau):
+    """
+    Compute a nondivergent velocity field
+
+    Parameters
+    ----------
+    t : np.ndarray of type float
+        times in seconds at which to compute the velocity field
+
+    lon : np.ndarray of type float
+        longitude
+
+    lat : np.ndarray of type float
+        latitude
+
+    u_0 : float
+        velocity amplitude in meters per second
+
+    tau : float
+        time in seconds for the flow to circumnavigate the sphere
+
+    Returns
+    -------
+    u : np.ndarray of type float
+       zonal velocity
+
+    v : np.ndarray of type float
+       meridional velocity
+    """
+    lon_p = lon - 2. * pi * t / tau
+    coslat = cos(lat)
+    cost = cos(pi * t / tau)
+    u = (1 / tau) * (-u_0 * (sin(lon_p / 2)**2) * sin(2 * lat) *
+                     (coslat**2) * cost + 2. * pi * coslat)
+    v = (u_0 / (2 * tau)) * sin(lon_p) * (coslat**3) * cost
+    return u, v
+
+
+def flow_rotation(lon, lat, omega, tau, sphere_radius):
+    """
+    Compute a nondivergent velocity field
+
+    Parameters
+    ----------
+    lon : np.ndarray of type float
+        longitude
+
+    lat : np.ndarray of type float
+        latitude
+
+    omega : np.ndarray of type float
+        vector defining the axis of rotation of the flow in cartesian
+        coordinates
+
+    tau : float
+        time in seconds for the flow to circumnavigate the sphere
+
+    sphere_radius : float
+        radius of the sphere
+
+    Returns
+    -------
+    u : np.ndarray of type float
+       zonal velocity
+
+    v : np.ndarray of type float
+       meridional velocity
+    """
+    omega = (2. * pi / tau) * (omega / np.linalg.norm(omega))
+    x, y, z = lon_lat_to_cartesian(lon, lat, sphere_radius, degrees=False)
+    xyz = np.stack((x, y, z), axis=1)
+    vel = np.cross(omega, np.divide(np.transpose(xyz), sphere_radius), axis=0)
+    east, north = calc_local_east_north(x, y, z)
+    u = np.sum(vel * east, axis=0)
+    v = np.sum(vel * north, axis=0)
+    return u, v
+
+
+def calc_local_east_north(x, y, z):
+    axis = [0, 0, 1]
+    xyz = np.stack((x, y, z), axis=1)
+    east = np.cross(axis, np.transpose(xyz), axis=0)
+    north = np.cross(np.transpose(xyz), east, axis=0)
+    east = east / np.linalg.norm(east, axis=0)
+    north = north / np.linalg.norm(north, axis=0)
+    return east, north
diff --git a/polaris/ocean/tasks/sphere_transport/resources/tracer_distributions.py b/polaris/ocean/tasks/sphere_transport/resources/tracer_distributions.py
new file mode 100644
index 000000000..6eef52e71
--- /dev/null
+++ b/polaris/ocean/tasks/sphere_transport/resources/tracer_distributions.py
@@ -0,0 +1,213 @@
+import numpy as np
+from mpas_tools.transects import lon_lat_to_cartesian
+from mpas_tools.vector import Vector
+
+
+def slotted_cylinders(lon, lat, r, b, c, sphere_radius):
+    """
+    Compute two slotted cylinders on the sphere
+    Lauritzen et al. 2012 eqn. (12)
+
+    Parameters
+    ----------
+    lon : np.ndarray of type float
+        longitude of cells on the sphere
+
+    lat : np.ndarray of type float
+        latitude of cells on the sphere
+
+    r : float
+        radius of each slotted cylinder
+
+    b : float
+        background value of the tracer
+
+    c : float
+        value of each slotted cylinder
+
+    sphere_radius : float
+        radius of the sphere
+
+    Returns
+    -------
+    scs : np.ndarray
+        slotted cylinder tracer values
+    """
+    lon_thr = 1 / (6 * 2)
+    lat_thr = 5 / (12 * 2)
+    lon1 = 5 * (np.pi / 6)
+    lat1 = 0
+    lon2 = -5 * (np.pi / 6)
+    lat2 = 0
+    lon0 = np.where(lon > np.pi,
+                    lon - 2 * np.pi,
+                    lon)
+    x, y, z = lon_lat_to_cartesian(lon, lat, sphere_radius, degrees=False)
+    x1, y1, z1 = lon_lat_to_cartesian(lon1, lat1, sphere_radius, degrees=False)
+    x2, y2, z2 = lon_lat_to_cartesian(lon2, lat2, sphere_radius, degrees=False)
+    xyz = Vector(x, y, z)
+    xyz1 = Vector(x1, y1, z1)
+    xyz2 = Vector(x2, y2, z2)
+    r1 = xyz.angular_distance(xyz1)
+    r2 = xyz.angular_distance(xyz2)
+    scs = np.where(r1 <= r,
+                   np.where(np.logical_or((abs(lon0 - lon1) >= lon_thr),
+                                          np.logical_and(
+                                              abs(lon0 - lon1) < lon_thr,
+                                              lat - lat1 < -lat_thr)),
+                            c,
+                            b),
+                   np.where(np.logical_and(r2 <= r,
+                                           np.logical_or(
+                                               (abs(lon0 - lon2) >= lon_thr),
+                                               np.logical_and(
+                                                   abs(lon0 - lon2) < lon_thr,
+                                                   lat - lat2 > lat_thr))),
+                            c,
+                            b))
+    return scs
+
+
+def cosine_bells(lon, lat, r, b, c, sphere_radius):
+    """
+    Compute two cosine bells on the sphere
+    Lauritzen et al. 2012 eqn. (11)
+
+    Parameters
+    ----------
+    lon : np.ndarray of type float
+        longitude of cells on the sphere
+
+    lat : np.ndarray of type float
+        latitude of cells on the sphere
+
+    r : float
+        radius of each cosine bell
+
+    b : float
+        background value of the tracer
+
+    c : float
+        maximum value of each cosine bell
+
+    sphere_radius : float
+        radius of the sphere
+
+    Returns
+    -------
+    cbs : np.ndarray
+        cosine bell tracer values
+    """
+    # Location of the center of the first cosine bell
+    lon1 = 5 * (np.pi / 6)
+    lat1 = 0
+
+    # Location of the center of the second cosine bell
+    lon2 = -5 * (np.pi / 6)
+    lat2 = 0
+
+    x, y, z = lon_lat_to_cartesian(lon, lat, sphere_radius, degrees=False)
+    x1, y1, z1 = lon_lat_to_cartesian(lon1, lat1, sphere_radius, degrees=False)
+    x2, y2, z2 = lon_lat_to_cartesian(lon2, lat2, sphere_radius, degrees=False)
+    xyz = Vector(x, y, z)
+    xyz1 = Vector(x1, y1, z1)
+    xyz2 = Vector(x2, y2, z2)
+    # Distance of each cell from the center of the first cosine bell
+    r1 = xyz.angular_distance(xyz1)
+    # Distance of each cell from the center of the second cosine bell
+    r2 = xyz.angular_distance(xyz2)
+
+    cbs = np.where(r1 < r,
+                   cosine_bell(1.0, r1, r),
+                   np.where(r2 < r,
+                            cosine_bell(1.0, r2, r),
+                            0.))
+    return b + c * cbs
+
+
+def xyztrig(lon, lat, sphere_radius):
+    """
+    Compute C-infinity tracer (not included in Lauritzen et al. 2012)
+
+    Parameters
+    ----------
+    lon : np.ndarray of type float
+        longitude of cells on the sphere
+
+    lat : np.ndarray of type float
+        latitude of cells on the sphere
+
+    r : float
+        radius of each cosine bell
+
+    b : float
+        background value of the tracer
+
+    c : float
+        maximum value of each cosine bell
+
+    sphere_radius : float
+        radius of the sphere
+
+    Returns
+    -------
+    f : np.ndarray
+        C-infinity tracer values
+    """
+    x, y, z = lon_lat_to_cartesian(lon, lat, sphere_radius, degrees=False)
+    x = np.divide(x, sphere_radius)
+    y = np.divide(y, sphere_radius)
+    z = np.divide(z, sphere_radius)
+    f = 0.5 * (1 + np.sin(np.pi * x) * np.sin(np.pi * y) * np.sin(np.pi * z))
+    return f
+
+
+def cosine_bell(max_value, ri, r):
+    """
+    Compute values according to cosine bell function
+    Lauritzen et al. 2012 eqn. (10)
+
+    Parameters
+    ----------
+    max_value : float
+        Maximum value of the cosine bell function
+
+    ri : np.ndarray of type float
+        Distance from the center of the cosine bell in meters
+
+    r : float
+        Radius of the cosine bell in meters
+
+    Returns
+    -------
+    f : np.ndarray of type float
+        Cosine bell tracer values
+    """
+    return max_value / 2.0 * (1.0 + np.cos(np.pi * np.divide(ri, r)))
+
+
+def correlation_fn(q1, a, b, c):
+    """
+    Compute a quadratic function for nonlinear tracer correlation following
+    Lauritzen et al. 2012 eqns. (14) and (15)
+
+    Parameters
+    ----------
+    q1 : np.ndarray
+        tracer values
+
+    a : float
+       quadratic coefficient
+
+    b : float
+       linear coefficient
+
+    c : float
+       offset
+
+    Returns
+    -------
+    q2 : np.ndarray
+        correlated tracer values
+    """
+    return a * q1**2. + b * q1 + c
diff --git a/polaris/ocean/tasks/sphere_transport/resources/wh-bl-gr-ye-re.rgb b/polaris/ocean/tasks/sphere_transport/resources/wh-bl-gr-ye-re.rgb
new file mode 100644
index 000000000..d3d0950a1
--- /dev/null
+++ b/polaris/ocean/tasks/sphere_transport/resources/wh-bl-gr-ye-re.rgb
@@ -0,0 +1,202 @@
+ncolors= 199
+
+#  r   g   b
+ 255 255 255
+ 250 250 255
+ 245 245 255
+ 240 240 255
+ 235 235 255
+ 230 230 255
+ 224 224 255
+ 219 219 255
+ 214 214 255
+ 209 209 255
+ 204 204 255
+ 199 199 255
+ 194 194 255
+ 189 189 255
+ 184 184 255
+ 179 179 255
+ 173 173 255
+ 168 168 255
+ 163 163 255
+ 158 158 255
+ 153 153 255
+ 148 148 255
+ 143 143 255
+ 138 138 255
+ 133 133 255
+ 128 128 255
+ 122 122 255
+ 117 117 255
+ 112 112 255
+ 107 107 255
+ 102 102 255
+  97  97 255
+  92  92 255
+  87  87 255
+  82  82 255
+  77  77 255
+  71  71 255
+  66  66 255
+  61  61 255
+  56  56 255
+  51  51 255
+  46  46 255
+  41  41 255
+  36  36 255
+  31  31 255
+  26  26 255
+  20  20 255
+  15  15 255
+  10  10 255
+   5   5 255
+   0   0 255
+   0   4 250
+   0   7 245
+   0  11 240
+   1  15 235
+   1  18 231
+   1  22 226
+   1  26 221
+   1  29 216
+   1  33 211
+   2  37 206
+   2  40 201
+   2  44 196
+   2  47 191
+   2  51 186
+   2  55 182
+   3  58 177
+   3  62 172
+   3  66 167
+   3  69 162
+   3  73 157
+   3  77 152
+   4  80 147
+   4  84 142
+   4  88 137
+   4  91 133
+   4  95 128
+   4  99 123
+   5 102 118
+   5 106 113
+   5 110 108
+   5 113 103
+   5 117  98
+   5 121  93
+   6 124  88
+   6 128  84
+   6 132  79
+   6 135  74
+   6 139  69
+   6 142  64
+   7 146  59
+   7 150  54
+   7 153  49
+   7 157  44
+   7 161  39
+   7 164  35
+   8 168  30
+   8 172  25
+   8 175  20
+   8 179  15
+  13 181  15
+  18 182  14
+  23 184  14
+  28 185  14
+  33 187  14
+  38 188  13
+  43 190  13
+  48 191  13
+  52 193  12
+  57 194  12
+  62 196  12
+  67 197  11
+  72 199  11
+  77 200  11
+  82 202  11
+  87 203  10
+  92 205  10
+  97 206  10
+ 102 208   9
+ 107 209   9
+ 112 211   9
+ 117 212   8
+ 122 214   8
+ 127 215   8
+ 132 217   8
+ 136 219   7
+ 141 220   7
+ 146 222   7
+ 151 223   6
+ 156 225   6
+ 161 226   6
+ 166 228   5
+ 171 229   5
+ 176 231   5
+ 181 232   5
+ 186 234   4
+ 191 235   4
+ 196 237   4
+ 201 238   3
+ 206 240   3
+ 211 241   3
+ 215 243   2
+ 220 244   2
+ 225 246   2
+ 230 247   2
+ 235 249   1
+ 240 250   1
+ 245 252   1
+ 250 253   0
+ 255 255   0
+ 255 250   0
+ 255 245   0
+ 255 239   0
+ 255 234   0
+ 255 229   0
+ 255 224   0
+ 255 219   0
+ 255 213   0
+ 255 208   0
+ 255 203   0
+ 255 198   0
+ 255 193   0
+ 255 187   0
+ 255 182   0
+ 255 177   0
+ 255 172   0
+ 255 167   0
+ 255 161   0
+ 255 156   0
+ 255 151   0
+ 255 146   0
+ 255 141   0
+ 255 135   0
+ 255 130   0
+ 255 125   0
+ 255 120   0
+ 255 114   0
+ 255 109   0
+ 255 104   0
+ 255  99   0
+ 255  94   0
+ 255  88   0
+ 255  83   0
+ 255  78   0
+ 255  73   0
+ 255  68   0
+ 255  62   0
+ 255  57   0
+ 255  52   0
+ 255  47   0
+ 255  42   0
+ 255  36   0
+ 255  31   0
+ 255  26   0
+ 255  21   0
+ 255  16   0
+ 255  10   0
+ 255   5   0
+ 255   0   0
diff --git a/polaris/ocean/tasks/sphere_transport/rotation_2d.cfg b/polaris/ocean/tasks/sphere_transport/rotation_2d.cfg
new file mode 100644
index 000000000..075270de7
--- /dev/null
+++ b/polaris/ocean/tasks/sphere_transport/rotation_2d.cfg
@@ -0,0 +1,10 @@
+# options for rotation_2d test case
+[rotation_2d]
+
+# rotation vector in cartesian coordinates
+rotation_vector = 0.2, 0.7, 1.0
+
+# convergence threshold below which the test fails
+convergence_thresh_tracer1 = 1.8
+convergence_thresh_tracer2 = 2.0
+convergence_thresh_tracer3 = 0.4
diff --git a/polaris/ocean/tasks/sphere_transport/sphere_transport.cfg b/polaris/ocean/tasks/sphere_transport/sphere_transport.cfg
new file mode 100644
index 000000000..30b425993
--- /dev/null
+++ b/polaris/ocean/tasks/sphere_transport/sphere_transport.cfg
@@ -0,0 +1,147 @@
+# Options related to the vertical grid
+[vertical_grid]
+
+# the type of vertical grid
+grid_type = uniform
+
+# Number of vertical levels
+vert_levels = 3
+
+# Depth of the bottom of the ocean
+bottom_depth = 300.0
+
+# The type of vertical coordinate (e.g. z-level, z-star)
+coord_type = z-level
+
+# 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
+
+
+# config options for spherical convergence tests
+[convergence]
+
+# Evaluation time for convergence analysis (in hours)
+convergence_eval_time = ${sphere_transport:vel_pd}
+
+# Error type
+error_type = l2
+
+
+# config options for spherical convergence tests
+[convergence_forward]
+
+# time integrator: {'split_explicit', 'RK4'}
+time_integrator = RK4
+
+# RK4 time step per resolution (s/km), since dt is proportional to resolution
+rk4_dt_per_km = 8.0
+
+# Run duration in hours
+run_duration = ${sphere_transport:vel_pd}
+
+# Output interval in hours
+output_interval = 24.0
+
+
+# options for all sphere transport test cases
+[sphere_transport]
+
+# temperature
+temperature = 15.
+
+# salinity
+salinity = 35.
+
+# time (hours) for bell to transit equator once
+vel_pd = 288.0
+
+# radius of cosine bells tracer distributions
+cosine_bells_radius = 0.5
+
+# background value of cosine bells tracer distribution
+cosine_bells_background = 0.1
+
+# amplitude of cosine bells tracer distribution
+cosine_bells_amplitude = 0.9
+
+# radius of slotted cylinders tracer distributions
+slotted_cylinders_radius = 0.5
+
+# background value of slotted cylinders tracer distribution
+slotted_cylinders_background = 0.1
+
+# amplitude of slotted cylinders tracer distribution
+slotted_cylinders_amplitude = 1.0
+
+# options for tracer visualization for the sphere transport test case
+[sphere_transport_viz]
+
+# visualization latitude and longitude resolution
+dlon = 0.5
+dlat = 0.5
+
+# remapping method ('bilinear', 'neareststod', 'conserve')
+remap_method = conserve
+
+
+# options for tracer visualization for the sphere transport test case
+[sphere_transport_viz_tracer]
+
+# colormap options
+# colormap
+colormap_name = viridis
+
+# the type of norm used in the colormap
+norm_type = linear
+
+# A dictionary with keywords for the norm
+norm_args = {'vmin': 0., 'vmax': 1.}
+
+# We could provide colorbar tick marks but we'll leave the defaults
+# colorbar_ticks = np.linspace(0., 1., 9)
+
+# options for plotting tracer differences from sphere transport tests
+[sphere_transport_viz_tracer_diff]
+
+# colormap options
+# colormap
+colormap_name = cmo.balance
+
+# the type of norm used in the colormap
+norm_type = linear
+
+# A dictionary with keywords for the norm
+norm_args = {'vmin': -0.25, 'vmax': 0.25}
+
+
+# options for thickness visualization for the sphere transport test case
+[sphere_transport_viz_h]
+
+# colormap options
+# colormap
+colormap_name = viridis
+
+# the type of norm used in the colormap
+norm_type = linear
+
+# A dictionary with keywords for the norm
+norm_args = {'vmin': 99., 'vmax': 101.}
+
+# We could provide colorbar tick marks but we'll leave the defaults
+# colorbar_ticks = np.linspace(0., 1., 9)
+
+# options for plotting tracer differences from sphere transport tests
+[sphere_transport_viz_h_diff]
+
+# colormap options
+# colormap
+colormap_name = cmo.balance
+
+# the type of norm used in the colormap
+norm_type = linear
+
+# A dictionary with keywords for the norm
+norm_args = {'vmin': -0.25, 'vmax': 0.25}
diff --git a/polaris/ocean/tasks/sphere_transport/viz.py b/polaris/ocean/tasks/sphere_transport/viz.py
new file mode 100644
index 000000000..03dee7b55
--- /dev/null
+++ b/polaris/ocean/tasks/sphere_transport/viz.py
@@ -0,0 +1,192 @@
+import cmocean  # noqa: F401
+import xarray as xr
+
+from polaris import Step
+from polaris.mpas import time_index_from_xtime
+from polaris.remap import MappingFileStep
+from polaris.viz import plot_global_field
+
+
+class VizMap(MappingFileStep):
+    """
+    A step for making a mapping file for cosine bell viz
+
+    Attributes
+    ----------
+    mesh_name : str
+        The name of the mesh
+    """
+    def __init__(self, component, name, subdir, base_mesh, mesh_name):
+        """
+        Create the step
+
+        Parameters
+        ----------
+        component : polaris.Component
+            The component the step belongs to
+
+        name : str
+            The name of the step
+
+        subdir : str
+            The subdirectory for the step
+
+        base_mesh : polaris.Step
+            The base mesh step
+
+        mesh_name : str
+            The name of the mesh
+        """
+        super().__init__(component=component, name=name, subdir=subdir,
+                         ntasks=128, min_tasks=1)
+        self.mesh_name = mesh_name
+        self.add_input_file(
+            filename='mesh.nc',
+            work_dir_target=f'{base_mesh.path}/base_mesh.nc')
+
+    def runtime_setup(self):
+        """
+        Set up the source and destination grids for this step
+        """
+        config = self.config
+        section = config['sphere_transport_viz']
+        dlon = section.getfloat('dlon')
+        dlat = section.getfloat('dlat')
+        method = section.get('remap_method')
+        self.src_from_mpas(filename='mesh.nc', mesh_name=self.mesh_name)
+        self.dst_global_lon_lat(dlon=dlon, dlat=dlat, lon_min=0.)
+        self.method = method
+
+        super().runtime_setup()
+
+
+class Viz(Step):
+    """
+    A step for plotting fields from the cosine bell output
+
+    Attributes
+    ----------
+    mesh_name : str
+        The name of the mesh
+    """
+    def __init__(self, component, name, subdir, base_mesh, init, forward,
+                 viz_map, mesh_name):
+        """
+        Create the step
+
+        Parameters
+        ----------
+        component : polaris.Component
+            The component the step belongs to
+
+        name : str
+            The name of the step
+
+        subdir : str
+            The subdirectory in the test case's work directory for the step
+
+        base_mesh : polaris.Step
+            The base mesh step
+
+        init : polaris.Step
+            The init step
+
+        forward : polaris.Step
+            The init step
+
+        viz_map : polaris.ocean.tasks.sphere_transport.viz.VizMap
+            The step for creating a mapping files, also used to remap data
+            from the MPAS mesh to a lon-lat grid
+
+        mesh_name : str
+            The name of the mesh
+        """
+        super().__init__(component=component, name=name, subdir=subdir)
+        self.add_input_file(
+            filename='mesh.nc',
+            work_dir_target=f'{base_mesh.path}/base_mesh.nc')
+        self.add_input_file(
+            filename='initial_state.nc',
+            work_dir_target=f'{init.path}/initial_state.nc')
+        self.add_input_file(
+            filename='output.nc',
+            work_dir_target=f'{forward.path}/output.nc')
+        self.add_dependency(viz_map, name='viz_map')
+        self.mesh_name = mesh_name
+        variables_to_plot = dict({'tracer1': 'tracer',
+                                  'tracer2': 'tracer',
+                                  'tracer3': 'tracer',
+                                  'layerThickness': 'h'})
+        self.variables_to_plot = variables_to_plot
+        for var in variables_to_plot.keys():
+            self.add_output_file(f'{var}_init.png')
+            self.add_output_file(f'{var}_final.png')
+            self.add_output_file(f'{var}_diff.png')
+
+    def run(self):
+        """
+        Run this step of the test case
+        """
+        config = self.config
+        mesh_name = self.mesh_name
+        run_duration = config.getfloat('convergence_forward',
+                                       'run_duration')
+
+        viz_map = self.dependencies['viz_map']
+
+        remapper = viz_map.get_remapper()
+
+        variables_to_plot = self.variables_to_plot
+        ds_init = xr.open_dataset('initial_state.nc')
+        ds_init = ds_init[variables_to_plot.keys()].isel(Time=0, nVertLevels=0)
+        ds_init = remapper.remap(ds_init)
+        ds_init.to_netcdf('remapped_init.nc')
+
+        ds_out = xr.open_dataset('output.nc')
+        s_per_hour = 3600.0
+
+        # Visualization at halfway around the globe (provided run duration is
+        # set to the time needed to circumnavigate the globe)
+        tidx = time_index_from_xtime(ds_out.xtime.values,
+                                     run_duration * s_per_hour / 2.)
+        ds_mid = ds_out[variables_to_plot.keys()].isel(Time=tidx,
+                                                       nVertLevels=0)
+        ds_mid = remapper.remap(ds_mid)
+        ds_mid.to_netcdf('remapped_mid.nc')
+
+        # Visualization at all the way around the globe
+        tidx = time_index_from_xtime(ds_out.xtime.values,
+                                     run_duration * s_per_hour)
+        ds_final = ds_out[variables_to_plot.keys()].isel(Time=tidx,
+                                                         nVertLevels=0)
+        ds_final = remapper.remap(ds_final)
+        ds_final.to_netcdf('remapped_final.nc')
+
+        for var, section_name in variables_to_plot.items():
+            colormap_section = f'sphere_transport_viz_{section_name}'
+            plot_global_field(ds_init.lon.values, ds_init.lat.values,
+                              ds_init[var].values,
+                              out_filename=f'{var}_init.png', config=config,
+                              colormap_section=colormap_section,
+                              title=f'{mesh_name} {var} at init',
+                              plot_land=False)
+            plot_global_field(ds_init.lon.values, ds_init.lat.values,
+                              ds_mid[var].values,
+                              out_filename=f'{var}_mid.png', config=config,
+                              colormap_section=colormap_section,
+                              title=f'{mesh_name} {var} after '
+                              f'{run_duration / 48.:g} days',
+                              plot_land=False)
+            plot_global_field(ds_init.lon.values, ds_init.lat.values,
+                              ds_final[var].values,
+                              out_filename=f'{var}_final.png', config=config,
+                              colormap_section=colormap_section,
+                              title=f'{mesh_name} {var} after '
+                              f'{run_duration / 24.:g} days', plot_land=False)
+            plot_global_field(ds_init.lon.values, ds_init.lat.values,
+                              ds_final[var].values - ds_init[var].values,
+                              out_filename=f'{var}_diff.png', config=config,
+                              colormap_section=f'{colormap_section}_diff',
+                              title=f'Difference in {mesh_name} {var} from '
+                              f'initial condition after {run_duration / 24.:g}'
+                              ' days', plot_land=False)