diff --git a/docs/users_guide/ocean/tasks/single_column.md b/docs/users_guide/ocean/tasks/single_column.md
index 9d656dacd3..9f6b2939f6 100644
--- a/docs/users_guide/ocean/tasks/single_column.md
+++ b/docs/users_guide/ocean/tasks/single_column.md
@@ -14,7 +14,7 @@ the vertical dynamics of the ocean model only. The test cases are:
 
 ## suppported models
 
-These tasks support only MPAS-Ocean.
+These tasks support MPAS-Ocean and Omega.
 
 ## mesh
 
@@ -58,30 +58,8 @@ min_pc_fraction = 0.1
 
 ## initial conditions
 
-The temperature and salinity profiles are defined using the following equations:
-
-$$
-\Phi(z) = \begin{cases}
-    \Phi_0 &\text{ if } z = z[0]\\
-    \Phi_0 + {d\Phi/dz}_{ML} z &
-    \text{ if } z > z_{MLD}\\
-    (\Phi_0 + {\Delta\Phi}_{ML}) + {d\Phi/dz}_{int} (z - z_{MLD}) &
-    \text{ if } z \le z_{MLD}
-\end{cases}
-$$
-
-where $\Phi_0 = $`surface_X`, ${d\Phi/dz}_{ML} = $`X_gradient_mixed_layer`,
-$z_{MLD} = -$`mixed_layer_depth_X`, ${\Delta\Phi}_{ML} = $
-`X_difference_across_mixed_layer`, and ${d\Phi/dz}_{int} = $
-`X_gradient_interior`. `X` in the config options above is either `temperature`
-or `salinity`.
-
-The initial velocity is vertically uniform and given by
-`single_column:zonal_velocity` and `single_column:meridional_velocity`, which
-are 0 by default (at rest).
-
-The Coriolis parameter is spatially constant and set equal to
-`coriolis_parameter`.
+The initial conditions are either stably stratified or uniform. See each task
+below.
 
 ### forcing
 
@@ -245,8 +223,35 @@ See {ref}`ocean-single-column`.
 
 See {ref}`ocean-single-column`.
 
+(ocean-single-column-stable)=
+
 ### initial conditions
 
+The temperature and salinity profiles are defined using the following equations:
+
+$$
+\Phi(z) = \begin{cases}
+    \Phi_0 &\text{ if } z = z[0]\\
+    \Phi_0 + {d\Phi/dz}_{ML} z &
+    \text{ if } z > z_{MLD}\\
+    (\Phi_0 + {\Delta\Phi}_{ML}) + {d\Phi/dz}_{int} (z - z_{MLD}) &
+    \text{ if } z \le z_{MLD}
+\end{cases}
+$$
+
+where $\Phi_0 = $`surface_X`, ${d\Phi/dz}_{ML} = $`X_gradient_mixed_layer`,
+$z_{MLD} = -$`mixed_layer_depth_X`, ${\Delta\Phi}_{ML} = $
+`X_difference_across_mixed_layer`, and ${d\Phi/dz}_{int} = $
+`X_gradient_interior`. `X` in the config options above is either `temperature`
+or `salinity`.
+
+The initial velocity is vertically uniform and given by
+`single_column:zonal_velocity` and `single_column:meridional_velocity`, which
+are 0 by default (at rest).
+
+The Coriolis parameter is spatially constant and set equal to
+`coriolis_parameter`.
+
 These config options overwrite those in {ref}`ocean-single-column`:
 
 ```cfg
@@ -263,7 +268,7 @@ mixed_layer_depth_temperature =  25.0
 salinity_difference_across_mixed_layer = 1.0
 ```
 
-The rest of the config options are given in {ref}`ocean-single-column`.
+(ocean-single-column-wind-evap)=
 
 ### forcing
 
@@ -301,8 +306,7 @@ days.
 
 ### config options
 
-See {ref}`ocean-single-column`. Currently, config options are only given in the
-shared framework.
+See {ref}`ocean-single-column-wind-evap` and {ref}`ocean-single-column-stable`.
 
 ### cores
 
@@ -345,7 +349,8 @@ The rest of the vertical grid features follow {ref}`ocean-single-column`.
 ### initial conditions
 
 The temperature and salinity are constant and the flow is at rest.
-See {ref}`ocean-single-column`.
+
+(ocean-single-column-wind)=
 
 ### forcing
 
@@ -378,7 +383,8 @@ viscosity:
 vertical_viscosity = 1.e-3
 ```
 
-All other config options derive from {ref}`ocean-single-column`.
+All other config options derive from {ref}`ocean-single-column` and
+{ref}`ocean-single-column-wind`.
 
 ### cores
 
@@ -412,6 +418,7 @@ See {ref}`ocean-single-column`.
 ### initial conditions
 
 `idealAgeTracers` is initialized as zero seconds throughout the water column.
+See {ref}`ocean-single-column-stable`.
 
 ### forcing
 
diff --git a/e3sm_submodules/Omega b/e3sm_submodules/Omega
index f2e951ac69..1b4057ed7c 160000
--- a/e3sm_submodules/Omega
+++ b/e3sm_submodules/Omega
@@ -1 +1 @@
-Subproject commit f2e951ac69e0b5667900906ff22ebb1587040d91
+Subproject commit 1b4057ed7c8898d4022b5d10754c76e0dcec780b
diff --git a/polaris/ocean/model/mpaso_to_omega.yaml b/polaris/ocean/model/mpaso_to_omega.yaml
index 15b5882835..9c7ac493a8 100644
--- a/polaris/ocean/model/mpaso_to_omega.yaml
+++ b/polaris/ocean/model/mpaso_to_omega.yaml
@@ -22,6 +22,9 @@ variables:
   cellsOnVertex: CellsOnVertex
   edgesOnVertex: EdgesOnVertex
 
+  # vertical
+  zMid: ZMid
+
   # tracers
   temperature: Temperature
   salinity: Salinity
@@ -91,6 +94,11 @@ config:
   options:
     config_disable_vel_explicit_bottom_drag: BottomDragTendencyEnable
 
+- section:
+    debug: Tendencies
+  options:
+    config_disable_vel_vmix: VelVertMixTendencyEnable
+
 - section:
     bottom_drag: Tendencies
   options:
@@ -114,3 +122,24 @@ config:
     config_eos_linear_alpha: DRhoDT
     config_eos_linear_beta: DRhoDS
     config_eos_linear_densityref: RhoT0S0
+
+- section:
+    cvmix: [VertMix, Background]
+  options:
+    config_cvmix_background_diffusion: Diffusivity
+    config_cvmix_background_viscosity: Viscosity
+
+- section:
+    cvmix: [VertMix, Convective]
+  options:
+    config_use_cvmix_convection: Enable
+    config_cvmix_convective_diffusion: Diffusivity
+    config_cvmix_convective_triggerBVF: TriggerBVF
+
+- section:
+    cvmix: [VertMix, Shear]
+  options:
+    config_use_cvmix_shear: Enable
+    config_cvmix_shear_PP_nu_zero: NuZero
+    config_cvmix_shear_PP_alpha: Alpha
+    config_cvmix_shear_PP_exp: Exponent
diff --git a/polaris/ocean/time.py b/polaris/ocean/time.py
new file mode 100644
index 0000000000..8f2bf9e5a7
--- /dev/null
+++ b/polaris/ocean/time.py
@@ -0,0 +1,34 @@
+from datetime import datetime
+
+import numpy as np
+import pandas as pd
+
+
+def get_days_since_start(ds):
+    """
+    Ocean model output may or may not include 'daysSinceStartOfSim'. This
+    routine uses 'daysSinceStartOfSim' if available, otherwise it uses 'Time'
+    """
+    if 'daysSinceStartOfSim' in ds.keys():
+        t_arr = ds.daysSinceStartOfSim.values.astype(float)
+    elif 'xtime' in ds.keys():
+        timestamps = []
+        for time_str in ds.xtime.values.astype(str):
+            try:
+                timestamp = datetime.strptime(time_str, '%Y-%m-%d_%H:%M:%S.%f')
+            except ValueError:
+                timestamp = datetime.strptime(time_str, '%Y-%m-%d_%H:%M:%S')
+            timestamps.append(timestamp)
+        # Calculate seconds since the first timestamp
+        seconds_since_start = [
+            (ts - timestamps[0]).total_seconds() for ts in timestamps
+        ]
+        t_arr = np.array(seconds_since_start, dtype=float) / 86400.0
+    elif 'Time' in ds.keys():
+        t_vals = ds['Time'].values
+        t_pd = pd.to_datetime(t_vals)
+        t_arr = 1.0e9 * (t_pd - t_pd[0]) / np.timedelta64(1, 's')
+        t_arr = t_arr.astype(float) / 86400.0
+    else:
+        raise ValueError('Could not find a time variable in dataset')
+    return t_arr
diff --git a/polaris/suites/ocean/mpaso_pr.txt b/polaris/suites/ocean/mpaso_pr.txt
index 0520e8903c..f7821c9504 100644
--- a/polaris/suites/ocean/mpaso_pr.txt
+++ b/polaris/suites/ocean/mpaso_pr.txt
@@ -11,6 +11,7 @@ ocean/planar/internal_wave/vlr/default
 # ocean/planar/merry_go_round/default
 ocean/planar/overflow/default
 ocean/single_column/cvmix
+ocean/single_column/cvmix_no_vadv
 ocean/single_column/ekman
 ocean/single_column/ideal_age
 ocean/single_column/inertial
diff --git a/polaris/tasks/ocean/single_column/__init__.py b/polaris/tasks/ocean/single_column/__init__.py
index 9999db60f4..808f791baa 100644
--- a/polaris/tasks/ocean/single_column/__init__.py
+++ b/polaris/tasks/ocean/single_column/__init__.py
@@ -1,7 +1,9 @@
-from polaris.tasks.ocean.single_column.cvmix import CVMix as CVMix
+from polaris.config import PolarisConfigParser as PolarisConfigParser
 from polaris.tasks.ocean.single_column.ekman import Ekman as Ekman
 from polaris.tasks.ocean.single_column.ideal_age import IdealAge as IdealAge
 from polaris.tasks.ocean.single_column.inertial import Inertial as Inertial
+from polaris.tasks.ocean.single_column.init import Init
+from polaris.tasks.ocean.single_column.vmix import VMix as VMix
 
 
 def add_single_column_tasks(component):
@@ -11,7 +13,152 @@ def add_single_column_tasks(component):
     component : polaris.tasks.ocean.Ocean
         the ocean component that the tasks will be added to
     """
-    component.add_task(Ekman(component=component))
-    component.add_task(CVMix(component=component))
-    component.add_task(IdealAge(component=component))
-    component.add_task(Inertial(component=component))
+    group_name = 'single_column'
+
+    name = 'vmix_stable'
+    forcing = ['wind']
+    forcing_dir = '_'.join(forcing) if forcing else 'no_forcing'
+    filepath = f'{component.name}/{group_name}/{name}/{name}.cfg'
+    config = PolarisConfigParser(filepath=filepath)
+    config.add_from_package(
+        'polaris.tasks.ocean.single_column', f'{group_name}.cfg'
+    )
+    for forcing_name in forcing:
+        config.add_from_package(
+            'polaris.tasks.ocean.single_column', f'{forcing_name}.cfg'
+        )
+    config.add_from_package(
+        'polaris.tasks.ocean.single_column', 'stable_stratification.cfg'
+    )
+    init_step = component.get_or_create_shared_step(
+        step_cls=Init,
+        subdir=f'{group_name}/{forcing_dir}/init_stable',
+        config=config,
+        config_filename=f'{name}.cfg',
+    )
+    component.add_task(
+        VMix(
+            component=component,
+            config=config,
+            name=name,
+            init=init_step,
+            indir=group_name,
+        )
+    )
+
+    name = 'vmix_unstable'
+    forcing = []
+    forcing_dir = '_'.join(forcing) if forcing else 'no_forcing'
+    filepath = f'{component.name}/{group_name}/{name}/{name}.cfg'
+    config = PolarisConfigParser(filepath=filepath)
+    config.add_from_package(
+        'polaris.tasks.ocean.single_column', f'{group_name}.cfg'
+    )
+    for forcing_name in forcing:
+        config.add_from_package(
+            'polaris.tasks.ocean.single_column', f'{forcing_name}.cfg'
+        )
+    config.add_from_package(
+        'polaris.tasks.ocean.single_column', 'unstable_stratification.cfg'
+    )
+    init_step = component.get_or_create_shared_step(
+        step_cls=Init,
+        subdir=f'{group_name}/{forcing_dir}/init_unstable',
+        config=config,
+        config_filename=f'{name}.cfg',
+    )
+    component.add_task(
+        VMix(
+            component=component,
+            config=config,
+            name=name,
+            init=init_step,
+            indir=group_name,
+        )
+    )
+
+    forcing = ['wind']
+    name = 'ekman'
+    filepath = f'{component.name}/{group_name}/{name}/{name}.cfg'
+    config = PolarisConfigParser(filepath=filepath)
+    config.add_from_package(
+        'polaris.tasks.ocean.single_column', f'{group_name}.cfg'
+    )
+    for forcing_name in forcing:
+        config.add_from_package(
+            'polaris.tasks.ocean.single_column', f'{forcing_name}.cfg'
+        )
+    init_step = component.get_or_create_shared_step(
+        step_cls=Init,
+        subdir=f'{group_name}/{forcing_name}/init',
+        config=config,
+        config_filename=f'{name}.cfg',
+    )
+    component.add_task(
+        Ekman(
+            component=component,
+            config=config,
+            init=init_step,
+            indir=group_name,
+        )
+    )
+
+    name = 'ideal_age'
+    forcing = ['evap']
+    forcing_dir = '_'.join(forcing) if forcing else 'no_forcing'
+    filepath = f'{component.name}/{group_name}/{name}/{name}.cfg'
+    config = PolarisConfigParser(filepath=filepath)
+    config.add_from_package(
+        'polaris.tasks.ocean.single_column', f'{group_name}.cfg'
+    )
+    for forcing_name in forcing:
+        config.add_from_package(
+            'polaris.tasks.ocean.single_column', f'{forcing_name}.cfg'
+        )
+    config.add_from_package(
+        'polaris.tasks.ocean.single_column', 'stable_stratification.cfg'
+    )
+    init_step = component.get_or_create_shared_step(
+        step_cls=Init,
+        subdir=f'{group_name}/{forcing_dir}/init_stable',
+        config=config,
+        config_filename=f'{name}.cfg',
+    )
+    component.add_task(
+        IdealAge(
+            component=component,
+            init=init_step,
+            config=config,
+            indir=group_name,
+        )
+    )
+
+    name = 'inertial'
+    forcing = []
+    forcing_dir = '_'.join(forcing) if forcing else 'no_forcing'
+    filepath = f'{component.name}/{group_name}/{name}/{name}.cfg'
+    config = PolarisConfigParser(filepath=filepath)
+    config.add_from_package(
+        'polaris.tasks.ocean.single_column', f'{group_name}.cfg'
+    )
+    for forcing_name in forcing:
+        config.add_from_package(
+            'polaris.tasks.ocean.single_column', f'{forcing_name}.cfg'
+        )
+    config.add_from_package(
+        'polaris.tasks.ocean.single_column', 'stable_stratification.cfg'
+    )
+    init_step = component.get_or_create_shared_step(
+        step_cls=Init,
+        subdir=f'{group_name}/{forcing_dir}/init_stable',
+        config=config,
+        config_filename=f'{name}.cfg',
+    )
+    component.add_task(
+        Inertial(
+            component=component,
+            init=init_step,
+            config=config,
+            indir=group_name,
+        )
+    )
diff --git a/polaris/tasks/ocean/single_column/cvmix/__init__.py b/polaris/tasks/ocean/single_column/cvmix/__init__.py
deleted file mode 100644
index 38349a7113..0000000000
--- a/polaris/tasks/ocean/single_column/cvmix/__init__.py
+++ /dev/null
@@ -1,55 +0,0 @@
-import os
-
-from polaris import Task
-from polaris.tasks.ocean.single_column.forward import Forward
-from polaris.tasks.ocean.single_column.init import Init
-from polaris.tasks.ocean.single_column.viz import Viz
-
-
-class CVMix(Task):
-    """
-    The CVMix single-column test case creates the mesh and initial condition,
-    then performs a short forward run testing vertical mixing on 1 core.
-    """
-
-    def __init__(self, component):
-        """
-        Create the test case
-        Parameters
-        ----------
-        component : polaris.tasks.ocean.Ocean
-            The ocean component that this task belongs to
-        """
-        name = 'cvmix'
-        subdir = os.path.join('single_column', name)
-        super().__init__(component=component, name=name, subdir=subdir)
-        self.config.add_from_package(
-            'polaris.tasks.ocean.single_column', 'single_column.cfg'
-        )
-        self.config.add_from_package(
-            'polaris.tasks.ocean.single_column.cvmix', 'cvmix.cfg'
-        )
-        self.add_step(Init(component=component, indir=self.subdir))
-
-        validate_vars = [
-            'temperature',
-            'salinity',
-            'layerThickness',
-            'normalVelocity',
-        ]
-        self.add_step(
-            Forward(
-                component=component,
-                indir=self.subdir,
-                ntasks=1,
-                min_tasks=1,
-                openmp_threads=1,
-                validate_vars=validate_vars,
-                task_name=name,
-            )
-        )
-
-        self.add_step(
-            Viz(component=component, indir=self.subdir),
-            run_by_default=False,
-        )
diff --git a/polaris/tasks/ocean/single_column/cvmix/forward.yaml b/polaris/tasks/ocean/single_column/cvmix/forward.yaml
deleted file mode 100644
index f2a7bd5282..0000000000
--- a/polaris/tasks/ocean/single_column/cvmix/forward.yaml
+++ /dev/null
@@ -1,21 +0,0 @@
-mpas-ocean:
-  time_management:
-    config_run_duration: 0010_00:00:00
-  hmix_del4:
-    config_use_mom_del4: true
-    config_mom_del4: 2.0e14
-  frazil_ice:
-    config_use_frazil_ice_formation: true
-  cvmix:
-    config_use_cvmix_convection: true
-    config_use_cvmix_shear: true
-    config_cvmix_shear_mixing_scheme: KPP
-    config_use_cvmix_kpp: true
-    config_cvmix_kpp_matching: MatchBoth
-    config_cvmix_kpp_interpolationOMLType: cubic
-  eos:
-    config_eos_type: linear
-  eos_linear:
-    config_eos_linear_beta: 0.01
-  streams:
-    KPP_testing: {}
diff --git a/polaris/tasks/ocean/single_column/ekman/__init__.py b/polaris/tasks/ocean/single_column/ekman/__init__.py
index 822642265f..85c43af9ee 100644
--- a/polaris/tasks/ocean/single_column/ekman/__init__.py
+++ b/polaris/tasks/ocean/single_column/ekman/__init__.py
@@ -3,7 +3,6 @@
 from polaris import Task
 from polaris.tasks.ocean.single_column.ekman.analysis import Analysis
 from polaris.tasks.ocean.single_column.forward import Forward
-from polaris.tasks.ocean.single_column.init import Init
 from polaris.tasks.ocean.single_column.viz import Viz
 
 
@@ -13,7 +12,7 @@ class Ekman(Task):
     then performs a short forward run spinning up an ekman layer on 1 core.
     """
 
-    def __init__(self, component):
+    def __init__(self, component, config, init, indir):
         """
         Create the test case
         Parameters
@@ -22,17 +21,18 @@ def __init__(self, component):
             The ocean component that this task belongs to
         """
         name = 'ekman'
-        subdir = os.path.join('single_column', name)
+        subdir = os.path.join(indir, name)
         super().__init__(component=component, name=name, subdir=subdir)
 
         self.config.add_from_package(
             'polaris.tasks.ocean.single_column', 'single_column.cfg'
         )
+        self.set_shared_config(config, link='ekman.cfg')
         self.config.add_from_package(
             'polaris.tasks.ocean.single_column.ekman', 'ekman.cfg'
         )
 
-        self.add_step(Init(component=component, indir=self.subdir))
+        self.add_step(init, symlink='init')
 
         validate_vars = [
             'temperature',
@@ -40,21 +40,29 @@ def __init__(self, component):
             'layerThickness',
             'normalVelocity',
         ]
+        forward_step = Forward(
+            component=component,
+            indir=self.subdir,
+            ntasks=1,
+            min_tasks=1,
+            openmp_threads=1,
+            validate_vars=validate_vars,
+            task_name=name,
+            constant_diff=True,
+        )
+        self.add_step(forward_step)
+
         self.add_step(
-            Forward(
-                component=component,
-                indir=self.subdir,
-                ntasks=1,
-                min_tasks=1,
-                openmp_threads=1,
-                validate_vars=validate_vars,
-                task_name=name,
+            Analysis(
+                component=component, indir=self.subdir, forward=forward_step
             )
         )
 
-        self.add_step(Analysis(component=component, indir=self.subdir))
-
         self.add_step(
-            Viz(component=component, indir=self.subdir),
+            Viz(
+                component=component,
+                indir=self.subdir,
+                comparisons={'forward': '../forward_constant'},
+            ),
             run_by_default=False,
         )
diff --git a/polaris/tasks/ocean/single_column/ekman/analysis.py b/polaris/tasks/ocean/single_column/ekman/analysis.py
index 299d6ed24e..93498bc4cb 100644
--- a/polaris/tasks/ocean/single_column/ekman/analysis.py
+++ b/polaris/tasks/ocean/single_column/ekman/analysis.py
@@ -14,7 +14,7 @@ class Analysis(Step):
     A step for comparing the velocity profile to an analytic solution
     """
 
-    def __init__(self, component, indir):
+    def __init__(self, component, indir, forward):
         """
         Create the step
 
@@ -27,13 +27,20 @@ def __init__(self, component, indir):
             The subdirectory that the task belongs to, that this step will
             go into a subdirectory of
 
+        forward : polaris.Step
+            The forward step for this test case
         """
         super().__init__(component=component, name='analysis', indir=indir)
+        self.forward = forward
+
+    def setup(self):
         self.add_input_file(
-            filename='init.nc', target='../forward/initial_state.nc'
+            filename='init.nc',
+            target=f'{self.base_work_dir}/{self.forward.path}/init.nc',
         )
         self.add_input_file(
-            filename='output.nc', target='../forward/output.nc'
+            filename='output.nc',
+            target=f'{self.base_work_dir}/{self.forward.path}/output.nc',
         )
 
     def run(self):
diff --git a/polaris/tasks/ocean/single_column/ekman/ekman.cfg b/polaris/tasks/ocean/single_column/ekman/ekman.cfg
index 29aabf77f2..ac837d8541 100644
--- a/polaris/tasks/ocean/single_column/ekman/ekman.cfg
+++ b/polaris/tasks/ocean/single_column/ekman/ekman.cfg
@@ -3,98 +3,10 @@
 # Bottom depth
 bottom_depth = 100.
 
-# config options for single column testcases
 [single_column]
 
-# Surface temperature
-surface_temperature = 20.0
-
-# Temperature gradient in the mixed layer in degC/m
-temperature_gradient_mixed_layer = 0.0
-
-# The temperature below the mixed layer
-temperature_difference_across_mixed_layer = 0.0
-
-# Temperature gradient below the mixed layer
-temperature_gradient_interior = 0.0
-
-# Depth of the temperature mixed layer
-mixed_layer_depth_temperature =  0.0
-
-# Surface salinity
-surface_salinity = 35.0
-
-# Salinity gradient in the mixed layer in PSU/m
-salinity_gradient_mixed_layer = 0.0
-
-# The salinity below the mixed layer
-salinity_difference_across_mixed_layer = 0.0
-
-# Salinity gradient below the mixed layer
-salinity_gradient_interior = 0.0
-
-# Depth of the salinity mixed layer
-mixed_layer_depth_salinity = 0.0
-
-# coriolis parameter
-coriolis_parameter = 1.0e-4
-
-# config options for forcing single column testcases
-[single_column_forcing]
-
-# Piston velocity to control rate of restoring toward temperature_surface_restoring_value
-temperature_piston_velocity = 0.0
-
-# Piston velocity to control rate of restoring toward salinity_surface_restoring_value
-salinity_piston_velocity = 0.0
-
-# Rate at which temperature is restored toward the initial condition
-temperature_interior_restoring_rate = 0.0
-
-# Rate at which salinity is restored toward the initial condition
-salinity_interior_restoring_rate = 0.0
-
-# Net latent heat flux applied when bulk forcing is used. Positive values indicate a net
-# input of heat to ocean
-latent_heat_flux = 0.0
-
-# Net sensible heat flux applied when bulk forcing is used. Positive values indicate a
-# net input of heat to ocean
-sensible_heat_flux = 0.0
-
-# Net solar shortwave heat flux applied when bulk forcing is used. Positive values
-# indicate a net input of heat to ocean
-shortwave_heat_flux = 0.0
-
-# Net surface evaporation when bulk forcing is used. Positive values indicate a net
-# input of water to ocean
-evaporation_flux = 0.0
-
-# Net surface rain flux when bulk forcing is used. Positive values indicate a net input
-# of water to ocean
-rain_flux = 0.0
-
-# Net surface river runoff flux when bulk forcing is used. Positive values indicate a net
-#flux of water to ocean
-river_runoff_flux = 0.0
-
-# Net surface subglacial runoff flux when bulk forcing is used. Positive values indicate a net
-#flux of water to ocean
-subglacial_runoff_flux = 0.0
-
-# Net surface ice runoff flux when bulk forcing is used. Positive values indicate a net
-#flux of water to ocean
-ice_runoff_flux = 0.0
-
-# Net iceberg freshwater flux when bulk forcing is used. Positive values indicate a net
-#flux of water to ocean
-iceberg_flux = 0.0
-
-# Zonal surface wind stress over the domain
-wind_stress_zonal = 0.1
-
-# Meridional surface wind stress over the domain
-wind_stress_meridional = 0.0
+# Run duration in days
+run_duration = 5.
 
 [single_column_ekman]
 
diff --git a/polaris/tasks/ocean/single_column/ekman/forward.yaml b/polaris/tasks/ocean/single_column/ekman/forward.yaml
index 85740192a9..bcec302721 100644
--- a/polaris/tasks/ocean/single_column/ekman/forward.yaml
+++ b/polaris/tasks/ocean/single_column/ekman/forward.yaml
@@ -1,33 +1,36 @@
-mpas-ocean:
+ocean:
   time_management:
     config_run_duration: 0005_00:00:00
   hmix_del4:
     config_use_mom_del4: false
-  frazil_ice:
-    config_use_frazil_ice_formation: false
-  bottom_drag:
-    config_implicit_constant_bottom_drag_coeff: 0.0
   cvmix:
-    config_cvmix_background_scheme: constant
     config_cvmix_background_diffusion: 0.0
     config_use_cvmix_convection: false
     config_use_cvmix_shear: false
-    config_use_cvmix_kpp: false
   eos:
     config_eos_type: linear
   eos_linear:
     config_eos_linear_beta: 0.01
   forcing:
     config_use_bulk_wind_stress: true
+
+mpas-ocean:
+  bottom_drag:
+    config_implicit_constant_bottom_drag_coeff: 0.0
+  frazil_ice:
+    config_use_frazil_ice_formation: false
+  cvmix:
+    config_cvmix_background_scheme: constant
+    config_use_cvmix_kpp: false
   tracer_forcing_activeTracers:
     config_use_activeTracers_surface_bulk_forcing: false
     config_use_activeTracers_surface_restoring: false
     config_use_activeTracers_interior_restoring: false
   streams:
     mesh:
-      filename_template: initial_state.nc
+      filename_template: mesh.nc
     input:
-      filename_template: initial_state.nc
+      filename_template: init.nc
     restart: {}
     KPP_testing: {}
     mixedLayerDepthsOutput:
diff --git a/polaris/tasks/ocean/single_column/cvmix/cvmix.cfg b/polaris/tasks/ocean/single_column/evap.cfg
similarity index 63%
rename from polaris/tasks/ocean/single_column/cvmix/cvmix.cfg
rename to polaris/tasks/ocean/single_column/evap.cfg
index 0252aced50..688bae9293 100644
--- a/polaris/tasks/ocean/single_column/cvmix/cvmix.cfg
+++ b/polaris/tasks/ocean/single_column/evap.cfg
@@ -1,15 +1,3 @@
-# config options for single column testcases
-[single_column]
-
-# Temperature gradient below the mixed layer
-temperature_gradient_interior = 0.01
-
-# Depth of the temperature mixed layer
-mixed_layer_depth_temperature =  25.0
-
-# The salinity below the mixed layer
-salinity_difference_across_mixed_layer = 1.0
-
 # config options for forcing single column testcases
 [single_column_forcing]
 
@@ -28,6 +16,3 @@ shortwave_heat_flux = 200.0
 # Net surface evaporation when bulk forcing is used. Positive values indicate a net
 # input of water to ocean
 evaporation_flux = 6.5E-4
-
-# Zonal surface wind stress over the domain
-wind_stress_zonal = 0.1
diff --git a/polaris/tasks/ocean/single_column/forward.py b/polaris/tasks/ocean/single_column/forward.py
index ed9c992aea..6431035ccc 100644
--- a/polaris/tasks/ocean/single_column/forward.py
+++ b/polaris/tasks/ocean/single_column/forward.py
@@ -24,6 +24,8 @@ def __init__(
         openmp_threads=1,
         validate_vars=None,
         task_name='',
+        enable_vadv=True,
+        constant_diff=False,
     ):
         """
         Create a new test case
@@ -62,6 +64,10 @@ def __init__(
         task_name : str, optional
             the name of the test case
         """
+        if not enable_vadv:
+            name = f'{name}_no_vadv'
+        if constant_diff:
+            name = f'{name}_constant'
         super().__init__(
             component=component,
             name=name,
@@ -75,12 +81,19 @@ def __init__(
         self.add_yaml_file('polaris.ocean.config', 'output.yaml')
 
         self.add_input_file(
-            filename='initial_state.nc', target='../init/initial_state.nc'
+            filename='mesh.nc',
+            target='../init/initial_state.nc',
+            # filename='mesh.nc', target='../init/culled_mesh.nc'
+        )
+        self.add_input_file(
+            filename='init.nc', target='../init/initial_state.nc'
         )
-        self.add_input_file(filename='forcing.nc', target='../init/forcing.nc')
         self.add_input_file(
             filename='graph.info', target='../init/culled_graph.info'
         )
+        self.add_input_file(
+            filename='forcing.nc', target='../init/initial_state.nc'
+        )
 
         self.add_yaml_file('polaris.tasks.ocean.single_column', 'forward.yaml')
         self.add_yaml_file(
@@ -93,6 +106,20 @@ def __init__(
 
         self.task_name = task_name
 
+        self.enable_vadv = enable_vadv
+
+        self.constant_diff = constant_diff
+
+    def setup(self):
+        """
+        TEMP: symlink initial condition to name hard-coded in Omega
+        """
+        super().setup()
+        model = self.config.get('ocean', 'model')
+        # TODO: remove as soon as Omega no longer hard-codes this file
+        if model == 'omega':
+            self.add_input_file(filename='OmegaMesh.nc', target='init.nc')
+
     def dynamic_model_config(self, at_setup):
         if self.task_name == 'ekman':
             nu = self.config.getfloat(
@@ -102,3 +129,30 @@ def dynamic_model_config(self, at_setup):
                 options={'config_cvmix_background_viscosity': nu},
                 config_model='mpas-ocean',
             )
+        if not self.enable_vadv:
+            self.add_model_config_options(
+                options={
+                    'config_vert_coord_movement': 'impermeable_interfaces',
+                    'config_disable_thick_vadv': True,
+                    'config_disable_vel_vadv': True,
+                    'config_disable_tr_adv': True,
+                },
+                config_model='mpas-ocean',
+            )
+
+        if self.constant_diff:
+            self.add_model_config_options(
+                options={
+                    'config_use_cvmix_convection': False,
+                    'config_use_cvmix_shear': False,
+                },
+                config_model='ocean',
+            )
+        else:
+            self.add_model_config_options(
+                options={
+                    'config_use_cvmix_convection': True,
+                    'config_use_cvmix_shear': True,
+                },
+                config_model='ocean',
+            )
diff --git a/polaris/tasks/ocean/single_column/forward.yaml b/polaris/tasks/ocean/single_column/forward.yaml
index bdaab85fb9..1a995a5fc0 100644
--- a/polaris/tasks/ocean/single_column/forward.yaml
+++ b/polaris/tasks/ocean/single_column/forward.yaml
@@ -1,23 +1,24 @@
-mpas-ocean:
-  run_modes:
-    config_ocean_run_mode: forward
+ocean:
   time_management:
     config_stop_time: none
   time_integration:
-    config_dt: 00:10:00
+    config_dt: 0000_00:10:00
+  forcing:
+    config_use_bulk_wind_stress: true
+
+mpas-ocean:
+  run_modes:
+    config_ocean_run_mode: forward
   split_explicit_ts:
     config_btr_dt: 0000_00:00:30
   forcing:
-    config_use_bulk_wind_stress: true
-  tracer_forcing_activeTracers:
-    config_use_activeTracers_surface_bulk_forcing: true
-    config_use_activeTracers_surface_restoring: true
-    config_use_activeTracers_interior_restoring: true
+    config_bulk_wind_stress_interp_isotropic: true
+    config_use_bulk_thickness_flux: true
   streams:
     mesh:
-      filename_template: initial_state.nc
+      filename_template: mesh.nc
     input:
-      filename_template: initial_state.nc
+      filename_template: init.nc
     restart: {}
     mixedLayerDepthsOutput:
       contents:
@@ -48,7 +49,7 @@ mpas-ocean:
       - frazilSalinityTendency
     forcing_data:
       type: input
-      filename_template: forcing.nc
+      filename_template: init.nc
       input_interval: initial_only
       contents:
       - tracersSurfaceRestoringFields
@@ -69,3 +70,54 @@ mpas-ocean:
       - icebergFreshWaterFlux
       - riverRunoffFlux
       - iceRunoffFlux
+
+Omega:
+  IOStreams:
+    InitialVertCoord:
+      Filename: init.nc
+    InitialState:
+      UsePointerFile: false
+      Filename: init.nc
+      Contents:
+      - State
+      - Tracers
+      - WindStressZonal
+      - WindStressMeridional
+    History:
+      Filename: output.nc
+      Freq: 1
+      FreqUnits: Seconds
+      IfExists: append
+      # effectively never
+      FileFreq: 9999
+      FileFreqUnits: years
+      Contents:
+      - AuxiliaryState
+      - State
+      - ZMid
+      - Tracers
+      - SshCell
+    RestartRead:
+      UsePointerFile: true
+      PointerFilename: ocn.pointer
+      Mode: read
+      Precision: double
+      Freq: 1
+      FreqUnits: OnStartup
+      UseStartEnd: true
+      StartTime: 0001-01-01_00:00:01
+      EndTime: 99999-12-31_00:00:00
+      Contents:
+      - Restart
+    RestartWrite:
+      UsePointerFile: true
+      PointerFilename: ocn.pointer
+      Filename: ocn.restart.$Y-$M-$D_$h.$m.$s
+      Mode: write
+      IfExists: replace
+      Precision: double
+      Freq: 6
+      FreqUnits: months
+      UseStartEnd: false
+      Contents:
+      - Restart
diff --git a/polaris/tasks/ocean/single_column/ideal_age/__init__.py b/polaris/tasks/ocean/single_column/ideal_age/__init__.py
index 614eff6da1..6c14f4dd2c 100644
--- a/polaris/tasks/ocean/single_column/ideal_age/__init__.py
+++ b/polaris/tasks/ocean/single_column/ideal_age/__init__.py
@@ -2,7 +2,6 @@
 
 from polaris import Task as Task
 from polaris.tasks.ocean.single_column.forward import Forward as Forward
-from polaris.tasks.ocean.single_column.init import Init as Init
 from polaris.tasks.ocean.single_column.viz import Viz as Viz
 
 
@@ -13,7 +12,7 @@ class IdealAge(Task):
     on 1 core.
     """
 
-    def __init__(self, component, ideal_age=True):
+    def __init__(self, component, config, init, indir, ideal_age=True):
         """
         Create the test case
         Parameters
@@ -22,15 +21,14 @@ def __init__(self, component, ideal_age=True):
             The ocean component that this task belongs to
         """
         name = 'ideal_age'
-        subdir = os.path.join('single_column', name)
+        subdir = os.path.join(indir, name)
         super().__init__(component=component, name=name, subdir=subdir)
+        config_filename = 'ideal_age.cfg'
+        self.set_shared_config(config, link=config_filename)
         self.config.add_from_package(
-            'polaris.tasks.ocean.single_column', 'single_column.cfg'
-        )
-
-        self.add_step(
-            Init(component=component, indir=self.subdir, ideal_age=ideal_age)
+            'polaris.tasks.ocean.single_column.ideal_age', config_filename
         )
+        self.add_step(init, symlink='init')
 
         validate_vars = ['temperature', 'salinity', 'iAge']
         step = Forward(
diff --git a/polaris/tasks/ocean/single_column/ideal_age/forward.yaml b/polaris/tasks/ocean/single_column/ideal_age/forward.yaml
index 487efcc5a1..62e8dbcdf5 100644
--- a/polaris/tasks/ocean/single_column/ideal_age/forward.yaml
+++ b/polaris/tasks/ocean/single_column/ideal_age/forward.yaml
@@ -1,6 +1,10 @@
 mpas-ocean:
   time_management:
     config_run_duration: 0010_00:00:00
+  tracer_forcing_activeTracers:
+    config_use_activeTracers_surface_bulk_forcing: true
+    config_use_activeTracers_surface_restoring: true
+    config_use_activeTracers_interior_restoring: true
   tracer_forcing_idealAgeTracers:
     config_use_idealAgeTracers: true
     config_use_idealAgeTracers_idealAge_forcing: true
diff --git a/polaris/tasks/ocean/single_column/ideal_age/ideal_age.cfg b/polaris/tasks/ocean/single_column/ideal_age/ideal_age.cfg
index 0252aced50..54db266069 100644
--- a/polaris/tasks/ocean/single_column/ideal_age/ideal_age.cfg
+++ b/polaris/tasks/ocean/single_column/ideal_age/ideal_age.cfg
@@ -1,33 +1,2 @@
-# config options for single column testcases
-[single_column]
-
-# Temperature gradient below the mixed layer
-temperature_gradient_interior = 0.01
-
-# Depth of the temperature mixed layer
-mixed_layer_depth_temperature =  25.0
-
-# The salinity below the mixed layer
-salinity_difference_across_mixed_layer = 1.0
-
-# config options for forcing single column testcases
-[single_column_forcing]
-
-# Net latent heat flux applied when bulk forcing is used. Positive values indicate a net
-# input of heat to ocean
-latent_heat_flux = -50.0
-
-# Net sensible heat flux applied when bulk forcing is used. Positive values indicate a
-# net input of heat to ocean
-sensible_heat_flux = -25.0
-
-# Net solar shortwave heat flux applied when bulk forcing is used. Positive values
-# indicate a net input of heat to ocean
-shortwave_heat_flux = 200.0
-
-# Net surface evaporation when bulk forcing is used. Positive values indicate a net
-# input of water to ocean
-evaporation_flux = 6.5E-4
-
-# Zonal surface wind stress over the domain
-wind_stress_zonal = 0.1
+# stable_stratification.cfg
+# wind.cfg
diff --git a/polaris/tasks/ocean/single_column/inertial/__init__.py b/polaris/tasks/ocean/single_column/inertial/__init__.py
index fb7bd85ac0..36c0f9f8c4 100644
--- a/polaris/tasks/ocean/single_column/inertial/__init__.py
+++ b/polaris/tasks/ocean/single_column/inertial/__init__.py
@@ -3,7 +3,6 @@
 from polaris import Task
 from polaris.tasks.ocean.single_column.forward import Forward
 from polaris.tasks.ocean.single_column.inertial.analysis import Analysis
-from polaris.tasks.ocean.single_column.init import Init
 from polaris.tasks.ocean.single_column.viz import Viz
 
 
@@ -13,7 +12,7 @@ class Inertial(Task):
     then performs a short forward run testing inertial oscillations on 1 core.
     """
 
-    def __init__(self, component):
+    def __init__(self, component, config, init, indir):
         """
         Create the test case
         Parameters
@@ -22,17 +21,14 @@ def __init__(self, component):
             The ocean component that this task belongs to
         """
         name = 'inertial'
-        subdir = os.path.join('single_column', name)
+        subdir = os.path.join(indir, name)
         super().__init__(component=component, name=name, subdir=subdir)
-
-        self.config.add_from_package(
-            'polaris.tasks.ocean.single_column', 'single_column.cfg'
-        )
+        config_filename = 'inertial.cfg'
+        self.set_shared_config(config, link=config_filename)
         self.config.add_from_package(
-            'polaris.tasks.ocean.single_column.inertial', 'inertial.cfg'
+            'polaris.tasks.ocean.single_column.inertial', config_filename
         )
-
-        self.add_step(Init(component=component, indir=self.subdir))
+        self.add_step(init, symlink='init')
 
         validate_vars = [
             'temperature',
diff --git a/polaris/tasks/ocean/single_column/inertial/forward.yaml b/polaris/tasks/ocean/single_column/inertial/forward.yaml
index ba694874f9..ca5fd87c5a 100644
--- a/polaris/tasks/ocean/single_column/inertial/forward.yaml
+++ b/polaris/tasks/ocean/single_column/inertial/forward.yaml
@@ -1,6 +1,8 @@
-mpas-ocean:
+ocean:
   time_management:
     config_run_duration: 0010_00:00:00
+
+mpas-ocean:
   tracer_forcing_activeTracers:
     config_use_activeTracers_surface_bulk_forcing: false
     config_use_activeTracers_surface_restoring: false
diff --git a/polaris/tasks/ocean/single_column/init.py b/polaris/tasks/ocean/single_column/init.py
index 4d59e7ac7f..be33671986 100644
--- a/polaris/tasks/ocean/single_column/init.py
+++ b/polaris/tasks/ocean/single_column/init.py
@@ -14,7 +14,7 @@ class Init(Step):
     test cases
     """
 
-    def __init__(self, component, indir, ideal_age=False):
+    def __init__(self, component, subdir):
         """
         Create the step
 
@@ -26,18 +26,13 @@ def __init__(self, component, indir, ideal_age=False):
         indir : str
             The subdirectory that the task belongs to, that this step will
             go into a subdirectory of
-
-        ideal_age : bool, optional
-            Whether the initial condition should include the ideal age tracer
         """
-        super().__init__(component=component, name='init', indir=indir)
-        self.ideal_age = ideal_age
+        super().__init__(component=component, name='init', subdir=subdir)
         for file in [
             'base_mesh.nc',
             'culled_mesh.nc',
             'culled_graph.info',
             'initial_state.nc',
-            'forcing.nc',
         ]:
             self.add_output_file(file)
 
@@ -48,7 +43,6 @@ def run(self):
         logger = self.logger
         config = self.config
         section = config['single_column']
-        ideal_age = self.ideal_age
         resolution = section.getfloat('resolution')
         nx = section.getint('nx')
         ny = section.getint('ny')
@@ -140,8 +134,8 @@ def run(self):
         normal_velocity = normal_velocity.transpose('nEdges', 'nVertLevels')
         normal_velocity = normal_velocity.expand_dims(dim='Time', axis=0)
 
-        if ideal_age:
-            ds['idealAgeTracers'] = xr.zeros_like(x_cell)
+        # We include this variable in initial conditions even when unused
+        ds['idealAgeTracers'] = xr.zeros_like(x_cell)
 
         ds['temperature'] = temperature
         ds['salinity'] = salinity
@@ -153,10 +147,9 @@ def run(self):
         ds.attrs['nx'] = nx
         ds.attrs['ny'] = ny
         ds.attrs['dc'] = dc
-        write_netcdf(ds, 'initial_state.nc')
 
         # create forcing stream
-        ds_forcing = xr.Dataset()
+        ds_forcing = ds.copy()
         forcing_array = xr.ones_like(temperature)
         forcing_array_surface = xr.ones_like(ds.bottomDepth)
         forcing_array_surface = forcing_array_surface.expand_dims(
@@ -239,4 +232,4 @@ def run(self):
         ds_forcing['icebergFreshWaterFlux'] = (
             iceberg_flux * forcing_array_surface
         )
-        write_netcdf(ds_forcing, 'forcing.nc')
+        write_netcdf(ds_forcing, 'initial_state.nc')
diff --git a/polaris/tasks/ocean/single_column/stable_stratification.cfg b/polaris/tasks/ocean/single_column/stable_stratification.cfg
new file mode 100644
index 0000000000..92de3c84cb
--- /dev/null
+++ b/polaris/tasks/ocean/single_column/stable_stratification.cfg
@@ -0,0 +1,18 @@
+
+# config options for single column testcases
+[single_column]
+
+# temperature gradient below the mixed layer
+temperature_gradient_interior = 0.05
+
+# Depth of the temperature mixed layer
+mixed_layer_depth_temperature =  0.0
+
+# The salinity below the mixed layer
+salinity_difference_across_mixed_layer = 0.0
+
+# Depth of the temperature mixed layer
+mixed_layer_depth_salinity =  0.0
+
+# Temperature gradient below the mixed layer
+salinity_gradient_interior = 0.0
diff --git a/polaris/tasks/ocean/single_column/unstable_stratification.cfg b/polaris/tasks/ocean/single_column/unstable_stratification.cfg
new file mode 100644
index 0000000000..d3831faafb
--- /dev/null
+++ b/polaris/tasks/ocean/single_column/unstable_stratification.cfg
@@ -0,0 +1,20 @@
+
+# config options for single column testcases
+[single_column]
+
+# Use this to test temperature and salinity
+#temperature_difference_across_mixed_layer = 1.0
+#Use this to test salinity only
+temperature_difference_across_mixed_layer = 0.0
+
+# Salinity gradient below the mixed layer
+salinity_gradient_interior = -0.01
+
+# Depth of the temperature mixed layer
+mixed_layer_depth_salinity =  25.0
+
+# Depth of the temperature mixed layer
+mixed_layer_depth_temperature =  25.0
+
+# The salinity below the mixed layer
+salinity_difference_across_mixed_layer = -1.0
diff --git a/polaris/tasks/ocean/single_column/viz.py b/polaris/tasks/ocean/single_column/viz.py
index 9d02337290..4a05a18393 100644
--- a/polaris/tasks/ocean/single_column/viz.py
+++ b/polaris/tasks/ocean/single_column/viz.py
@@ -1,17 +1,29 @@
+import os
+
 import matplotlib.pyplot as plt
 import numpy as np
-import xarray as xr
 
-from polaris import Step
+from polaris.ocean.model import OceanIOStep
+from polaris.ocean.time import get_days_since_start
+from polaris.ocean.vertical.diagnostics import depth_from_thickness
 from polaris.viz import use_mplstyle
 
+# TODO import rho_0 from constants
+
 
-class Viz(Step):
+class Viz(OceanIOStep):
     """
     A step for plotting the results of a single-column test
     """
 
-    def __init__(self, component, indir, ideal_age=False):
+    def __init__(
+        self,
+        component,
+        indir,
+        ideal_age=False,
+        comparisons=None,
+        variables=None,
+    ):
         """
         Create the step
 
@@ -26,69 +38,188 @@ def __init__(self, component, indir, ideal_age=False):
 
         ideal_age : bool, optional
             Whether the initial condition should include the ideal age tracer
+
+        comparisons : dict, optional
+            A dictionary of comparison datasets to use for validation
+
+        variables : dict, optional
+            A dictionary of variables to plot along with their units
         """
         super().__init__(component=component, name='viz', indir=indir)
-        self.ideal_age = ideal_age
-        self.add_input_file(
-            filename='initial_state.nc', target='../init/initial_state.nc'
+        self.comparisons = (
+            dict(comparisons)
+            if comparisons
+            else {'forward': '../forward/output.nc'}
         )
+        self.variables = (
+            dict(variables)
+            if variables
+            else {
+                'temperature': 'degC',
+                'salinity': 'PSU',
+                'velocity': 'm s$^{-1}$',
+            }
+        )
+        if ideal_age:
+            # Include age tracer
+            self.variables['iAge'] = 'seconds'
         self.add_input_file(
-            filename='output.nc', target='../forward/output.nc'
+            filename='initial_state.nc', target='../init/initial_state.nc'
         )
+        for comparison_name, comparison_path in self.comparisons.items():
+            self.add_input_file(
+                filename=f'{comparison_name}.nc',
+                target=f'{comparison_path}/output.nc',
+            )
+            self.add_input_file(
+                filename=f'{comparison_name}_diags.nc',
+                target=f'{comparison_path}/output/KPP_test.0001-01-01_00.00.00.nc',
+            )
 
     def run(self):
         """
         Run this step of the test case
         """
         use_mplstyle()
-        ideal_age = self.ideal_age
-        ds = xr.load_dataset('output.nc')
-        t_days = ds.daysSinceStartOfSim.values
-        t = t_days.astype('timedelta64[ns]')
-        t = t / np.timedelta64(1, 'D')
-        t_index = np.argmin(np.abs(t - 1.0))  # ds.sizes['Time'] - 1
-        t_days = t[t_index]
+
+        ds_init = self.open_model_dataset('initial_state.nc')
+        ds_init = ds_init.isel(Time=0)
+        if 'zMid' in ds_init:
+            z_mid_init = ds_init['zMid'].mean(dim='nCells')
+        else:
+            comparison_name = next(iter(self.comparisons))
+            ds_comp = self.open_model_dataset(
+                f'{comparison_name}.nc', decode_times=False
+            )
+            z_mid_init = ds_comp['zMid'].isel(Time=0).mean(dim='nCells')
+
+        section = self.config['single_column']
+        t_target = section.getfloat('run_duration')
+        # t_target = 0.
 
         # Plot temperature and salinity profiles
-        title = f'final time = {t_days} days'
-        fields = {'temperature': 'degC', 'salinity': 'PSU'}
-        if ideal_age:
-            # Include age tracer
-            fields['iAge'] = 'seconds'
-        z_mid = ds['zMid'].mean(dim='nCells')
-        z_mid_init = z_mid.isel(Time=0)
-        z_mid_final = z_mid.isel(Time=t_index)
-        for field_name, field_units in fields.items():
-            if field_name not in ds.keys():
-                raise ValueError(f'{field_name} not present in output.nc')
-            var = ds[field_name].mean(dim='nCells')
-            var_init = var.isel(Time=0)
-            var_final = var.isel(Time=t_index)
-            plt.figure(figsize=(3, 5))
-            ax = plt.subplot(111)
-            ax.plot(var_init, z_mid_init, '--k', label='initial')
-            ax.plot(var_final, z_mid_final, '-k', label='final')
-            ax.set_xlabel(f'{field_name} ({field_units})')
-            ax.set_ylabel('z (m)')
-            ax.legend()
-            plt.title(title)
-            plt.tight_layout(pad=0.5)
-            plt.savefig(f'{field_name}.png')
-            plt.close()
-
-        # Plot velocity profiles
-        u = ds['velocityZonal'].mean(dim='nCells')
-        v = ds['velocityMeridional'].mean(dim='nCells')
-        u_final = u.isel(Time=t_index)
-        v_final = v.isel(Time=t_index)
-        plt.figure(figsize=(3, 5))
-        ax = plt.subplot(111)
-        ax.plot(u_final, z_mid_final, '-k', label='u')
-        ax.plot(v_final, z_mid_final, '-b', label='v')
-        ax.set_xlabel('Velocity (m/s)')
-        ax.set_ylabel('z (m)')
-        ax.legend()
-        plt.title(title)
-        plt.tight_layout(pad=0.5)
-        plt.savefig('velocity.png')
-        plt.close()
+        for field_name, field_units in self.variables.items():
+            fig = plt.figure(figsize=(3, 5))
+            colors = ['b', 'r', 'darkgreen']
+            if field_name == 'velocity':
+                if (
+                    'velocityZonal' not in ds_comp.keys()
+                    and 'velocityMeridional' not in ds_comp.keys()
+                ):
+                    continue
+                ax = plt.subplot(111)
+                for comparison_name, color in zip(
+                    self.comparisons.keys(), colors, strict=False
+                ):
+                    if not os.path.exists(f'{comparison_name}.nc'):
+                        continue
+                    ds_comp = self.open_model_dataset(
+                        f'{comparison_name}.nc', decode_times=False
+                    )
+                    t_arr = get_days_since_start(ds_comp)
+                    t_index = np.argmin(np.abs(t_arr - t_target))
+                    t_days = float(t_arr[t_index])
+                    title = f'final time = {t_days:2.1g} days'
+                    self.logger.info(
+                        f'Plot {field_name} for {comparison_name} at {t_days} '
+                        'days'
+                    )
+                    ds_comp = ds_comp.isel(Time=t_index)
+                    z_mid_final = ds_comp['zMid'].mean(dim='nCells')
+                    u_final = ds_comp['velocityZonal'].mean(dim='nCells')
+                    v_final = ds_comp['velocityMeridional'].mean(dim='nCells')
+                    ax.plot(
+                        u_final,
+                        z_mid_final,
+                        '-',
+                        color=color,
+                        label=f'u {comparison_name}',
+                    )
+                    ax.plot(
+                        v_final,
+                        z_mid_final,
+                        '--',
+                        color=color,
+                        label=f'v {comparison_name}',
+                    )
+                ax.set_xlabel('Velocity (m/s)')
+                ax.set_ylabel('z (m)')
+                ax.legend()
+                plt.title(title)
+                plt.tight_layout(pad=0.5)
+                plt.savefig('velocity.png')
+                plt.close()
+            else:
+                # Plot initial state if available
+                if field_name in ds_init.keys():
+                    var_init = ds_init[field_name].mean(dim='nCells')
+                    plt.plot(var_init, z_mid_init, '--k', label='initial')
+
+                for comparison_name, color in zip(
+                    self.comparisons.keys(), colors, strict=False
+                ):
+                    if not os.path.exists(f'{comparison_name}.nc'):
+                        continue
+                    # Look for field_name in either output file
+                    ds_comp = self.open_model_dataset(
+                        f'{comparison_name}.nc', decode_times=False
+                    )
+                    ds_diags = self.open_model_dataset(
+                        f'{comparison_name}_diags.nc', decode_times=False
+                    )
+                    if field_name in ds_comp.keys():
+                        ds = ds_comp
+                    elif field_name in ds_diags.keys():
+                        ds = ds_diags
+                    else:
+                        self.logger.warn(
+                            f'{field_name} not present in {comparison_name}.nc'
+                        )
+                        continue
+                    t_arr = get_days_since_start(ds)
+                    t_index = np.argmin(np.abs(t_arr - t_target))
+                    t_days = float(t_arr[t_index])
+                    if abs(t_days - t_target) > (1 / 24):
+                        self.logger.warn(
+                            f'{comparison_name}: Time mismatch \n'
+                            f'expected {t_target}, got {t_days}'
+                        )
+                    ds_final = ds.isel(Time=t_index)
+                    var_comp = ds_final[field_name].mean(dim='nCells')
+                    if 'nVertLevelsP1' in var_comp.dims:
+                        var_comp = var_comp.isel(nVertLevelsP1=slice(0, -1))
+                    # TODO delete this line when MPAS-O bug is fixed
+                    if field_name == 'RiTopOfCell':
+                        var_comp[0] = np.nan
+                    # TODO use this line when Omega zMid is correct
+                    # z_mid_final = ds_comp['zMid'].mean(dim='nCells')
+                    if 'layerThickness' not in ds.keys():
+                        z_mid_final = z_mid_init
+                        self.logger.warn(
+                            'Using initial zMid values; may not represent '
+                            'plotted state'
+                        )
+                    else:
+                        z_mid_final = depth_from_thickness(ds_final).mean(
+                            dim='nCells'
+                        )
+                    plt.plot(
+                        var_comp,
+                        z_mid_final,
+                        '-',
+                        color=color,
+                        label=comparison_name,
+                    )
+                title = f'final time = {t_days:2.1g} days'
+                plt.ylim([-100, 0])
+                if field_name == 'temperature':
+                    plt.xlim([15, 20])
+                else:
+                    plt.xlim(auto=True)
+                plt.xlabel(f'{field_name} ({field_units})')
+                plt.ylabel('z (m)')
+                fig.legend(loc='center right')
+                plt.title(title)
+                plt.tight_layout(pad=0.5)
+                plt.savefig(f'{field_name}.png')
+                plt.close()
diff --git a/polaris/tasks/ocean/single_column/vmix/__init__.py b/polaris/tasks/ocean/single_column/vmix/__init__.py
new file mode 100644
index 0000000000..8aa1c4dfe9
--- /dev/null
+++ b/polaris/tasks/ocean/single_column/vmix/__init__.py
@@ -0,0 +1,103 @@
+from polaris import Task
+from polaris.tasks.ocean.single_column.forward import Forward
+from polaris.tasks.ocean.single_column.viz import Viz
+from polaris.tasks.ocean.single_column.vmix.analysis import Analysis
+
+
+class VMix(Task):
+    """
+    The VMix single-column test case creates the mesh and initial condition,
+    then performs a short forward run testing vertical mixing on 1 core.
+    """
+
+    def __init__(self, component, config, init, indir, name='vmix'):
+        """
+        Create the test case
+        Parameters
+        ----------
+        component : polaris.tasks.ocean.Ocean
+            The ocean component that this task belongs to
+        """
+        super().__init__(component=component, name=name, indir=indir)
+        config_filename = 'vmix.cfg'
+        self.set_shared_config(config, link=config_filename)
+        self.config.add_from_package(
+            'polaris.tasks.ocean.single_column.vmix', config_filename
+        )
+        self.add_step(init, symlink='init')
+
+        validate_vars = [
+            'temperature',
+            'salinity',
+            'layerThickness',
+            'normalVelocity',
+        ]
+        self.add_step(
+            Forward(
+                component=component,
+                indir=f'{indir}/{name}',
+                ntasks=1,
+                min_tasks=1,
+                openmp_threads=1,
+                validate_vars=validate_vars,
+                task_name='vmix',
+                enable_vadv=True,
+            ),
+            run_by_default=False,
+        )
+        self.add_step(
+            Forward(
+                component=component,
+                indir=f'{indir}/{name}',
+                ntasks=1,
+                min_tasks=1,
+                openmp_threads=1,
+                validate_vars=validate_vars,
+                task_name='vmix',
+                enable_vadv=False,
+            )
+        )
+        self.add_step(
+            Forward(
+                component=component,
+                indir=f'{indir}/{name}',
+                ntasks=1,
+                min_tasks=1,
+                openmp_threads=1,
+                validate_vars=validate_vars,
+                task_name='vmix',
+                constant_diff=True,
+                enable_vadv=False,
+            )
+        )
+
+        self.add_step(
+            Viz(
+                component=component,
+                indir=f'{indir}/{name}',
+                comparisons={
+                    'standard': '../forward',
+                    'no_vadv': '../forward_no_vadv',
+                    'constant': '../forward_no_vadv_constant',
+                },
+                variables={
+                    'temperature': 'degC',
+                    'salinity': 'PSU',
+                    'velocity': 'm s$^{-1}$',
+                    'RiTopOfCell': '',
+                    'BruntVaisalaFreqTop': '$s^{-2}$',
+                },
+            )
+        )
+
+        self.add_step(
+            Analysis(
+                component=component,
+                indir=f'{indir}/{name}',
+                comparisons={
+                    'standard': '../forward',
+                    'no_vadv': '../forward_no_vadv',
+                    'constant': '../forward_no_vadv_constant',
+                },
+            )
+        )
diff --git a/polaris/tasks/ocean/single_column/vmix/analysis.py b/polaris/tasks/ocean/single_column/vmix/analysis.py
new file mode 100644
index 0000000000..d9f4dc7289
--- /dev/null
+++ b/polaris/tasks/ocean/single_column/vmix/analysis.py
@@ -0,0 +1,72 @@
+import numpy as np
+
+from polaris.ocean.model import OceanIOStep
+from polaris.ocean.time import get_days_since_start
+
+# TODO import rho_0 from constants
+
+
+class Analysis(OceanIOStep):
+    """
+    A step for analyzing the results of a single-column wind-forced test
+    """
+
+    def __init__(
+        self,
+        component,
+        indir,
+        comparisons=None,
+    ):
+        super().__init__(component=component, name='analysis', indir=indir)
+        self.comparisons = (
+            dict(comparisons)
+            if comparisons
+            else {'forward': '../forward/output.nc'}
+        )
+        for comparison_name, comparison_path in self.comparisons.items():
+            self.add_input_file(
+                filename=f'{comparison_name}_diags.nc',
+                target=f'{comparison_path}/output/KPP_test.0001-01-01_00.00.00.nc',
+            )
+
+    def run(self):
+        """
+        Run this step of the test case
+        """
+        section = self.config['single_column_forcing']
+        wind_stress = np.sqrt(
+            section.getfloat('wind_stress_zonal') ** 2.0
+            + section.getfloat('wind_stress_meridional') ** 2.0
+        )
+        # u_star = 0.01
+        rho_0 = 1026.0
+        u_star = wind_stress / rho_0
+        # TODO why is this 0
+        # u_star = ds_diags_1day['surfaceFrictionVelocity'].mean(dim='nCells')
+        # TODO compute this based on config parameters
+        N_sq_init = 1.0e-4
+        for comparison_name in self.comparisons.keys():
+            ds_diags = self.open_model_dataset(
+                f'{comparison_name}_diags.nc', decode_times=False
+            )
+            t_target = 1.0  # empirical relationship hold for up to 30h
+            t_arr = get_days_since_start(ds_diags)
+            t_index = np.argmin(np.abs(t_arr - t_target))
+            t_days = float(t_arr[t_index])
+            if abs(t_days - t_target) > (1 / 24):
+                self.logger.warn(
+                    f'{comparison_name}: Time mismatch \n'
+                    f'expected {t_target}, got {t_days}'
+                )
+            bld_theory = -u_star * (15.0 * t_days * 86400.0 / N_sq_init) ** (
+                1 / 3
+            )
+            ds_diags_1day = ds_diags.isel(Time=t_index)
+            z_top_final = ds_diags_1day['zTop'].mean(dim='nCells')
+            N_sq = ds_diags_1day['BruntVaisalaFreqTop'].mean(dim='nCells')
+            index_bld = int(np.nanargmax(N_sq.values))
+            bld = z_top_final.isel(nVertLevels=index_bld)
+            self.logger.info(
+                f'{comparison_name}: boundary layer depth '
+                f'expected {bld_theory}, actual {bld.values}'
+            )
diff --git a/polaris/tasks/ocean/single_column/vmix/forward.yaml b/polaris/tasks/ocean/single_column/vmix/forward.yaml
new file mode 100644
index 0000000000..ee492cb073
--- /dev/null
+++ b/polaris/tasks/ocean/single_column/vmix/forward.yaml
@@ -0,0 +1,69 @@
+ocean:
+  time_management:
+    config_run_duration: 0010_00:00:00
+  hmix_del4:
+    config_use_mom_del4: false
+  eos:
+    config_eos_type: linear
+  eos_linear:
+    config_eos_linear_beta: 0.01
+
+mpas-ocean:
+  frazil_ice:
+    config_use_frazil_ice_formation: true
+  tracer_forcing_activeTracers:
+    config_use_activeTracers_surface_bulk_forcing: true
+    config_use_activeTracers_surface_restoring: false
+    config_use_activeTracers_interior_restoring: false
+  cvmix:
+    config_use_cvmix_kpp: false
+    config_cvmix_background_diffusion: 0.
+    config_cvmix_background_viscosity: 0.
+    config_cvmix_shear_mixing_scheme: PP
+    config_cvmix_num_ri_smooth_loops: 0
+  streams:
+    KPP_testing:
+      contents:
+      - mesh
+      - layerThickness
+      - tracers
+      - vertNonLocalFlux
+      - xtime
+      - zMid
+      - zTop
+      - velocityZonal
+      - velocityMeridional
+      - bulkRichardsonNumber
+      - bulkRichardsonNumberBuoy
+      - potentialDensity
+      - unresolvedShear
+      - boundaryLayerDepth
+      - surfaceFrictionVelocity
+      - surfaceBuoyancyForcing
+      - windStressZonal
+      - windStressMeridional
+      - transportVelocityZonal
+      - transportVelocityMeridional
+      - RiTopOfCell
+      - BruntVaisalaFreqTop
+      - vertViscTopOfCell
+      - vertDiffTopOfCell
+
+Omega:
+  Tendencies:
+    ThicknessFluxTendencyEnable: true
+    PVTendencyEnable: false
+    KETendencyEnable: false
+    SSHTendencyEnable: false
+    VelDiffTendencyEnable: false
+    VelHyperDiffTendencyEnable: false
+    VelVertMixTendencyEnable: true
+    PresForceTendencyEnable: false
+    PresGradForceTendencyEnable: false
+    GeoptGradTendencyEnable: false
+    WindForcingTendencyEnable: true
+    BottomDragTendencyEnable: false
+    TracerHorzAdvTendencyEnable: false
+    TracerDiffTendencyEnable: false
+    TracerHyperDiffTendencyEnable: false
+    TracerVertMixTendencyEnable: true
diff --git a/polaris/tasks/ocean/single_column/vmix/vmix.cfg b/polaris/tasks/ocean/single_column/vmix/vmix.cfg
new file mode 100644
index 0000000000..5f68cf311c
--- /dev/null
+++ b/polaris/tasks/ocean/single_column/vmix/vmix.cfg
@@ -0,0 +1,8 @@
+# stable_stratification.cfg
+# evap.cfg
+# wind.cfg
+
+[single_column]
+
+# Run duration in days
+run_duration = 1.
diff --git a/polaris/tasks/ocean/single_column/wind.cfg b/polaris/tasks/ocean/single_column/wind.cfg
new file mode 100644
index 0000000000..057c8c1258
--- /dev/null
+++ b/polaris/tasks/ocean/single_column/wind.cfg
@@ -0,0 +1,5 @@
+[single_column_forcing]
+
+# Zonal surface wind stress over the domain
+# corresponding to friction velocity of 0.01
+wind_stress_zonal = 10.26