diff --git a/compass/ocean/tests/hurricane/__init__.py b/compass/ocean/tests/hurricane/__init__.py
index 15b5536299..6d59f39c3e 100644
--- a/compass/ocean/tests/hurricane/__init__.py
+++ b/compass/ocean/tests/hurricane/__init__.py
@@ -19,7 +19,7 @@ def __init__(self, mpas_core):
         storm = 'sandy'
         mesh_name = 'DEQU120at30cr10rr2'
 
-        for use_lts in [False, True]:
+        for use_lts in [False, 'LTS', 'FB_LTS']:
 
             mesh = Mesh(test_group=self,
                         mesh_name=mesh_name,
diff --git a/compass/ocean/tests/hurricane/fblts/__init_.py b/compass/ocean/tests/hurricane/fblts/__init_.py
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/compass/ocean/tests/hurricane/fblts/forward/__init__.py b/compass/ocean/tests/hurricane/fblts/forward/__init__.py
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/compass/ocean/tests/hurricane/fblts/forward/namelist.ocean b/compass/ocean/tests/hurricane/fblts/forward/namelist.ocean
new file mode 100644
index 0000000000..3850c97f8c
--- /dev/null
+++ b/compass/ocean/tests/hurricane/fblts/forward/namelist.ocean
@@ -0,0 +1,39 @@
+config_start_time = '2012-10-10_00:00:00'
+config_stop_time = '2012-11-03_00:00:00'
+config_run_duration = 'none'
+config_dt = '72'
+config_time_integrator = 'FB_LTS'
+config_dt_scaling_LTS = 4
+config_number_of_time_levels = 5
+config_hmix_scaleWithMesh = .true.
+config_use_mom_del4 = .true.
+config_mom_del4 = 4.0e8
+config_use_cvmix = .false.
+config_use_bulk_wind_stress = .true.
+config_use_time_varying_atmospheric_forcing = .true.
+config_time_varying_atmospheric_forcing_start_time = '2012-10-10_00:00:00'
+config_time_varying_atmospheric_forcing_reference_time = '2012-10-10_00:00:00'
+config_time_varying_atmospheric_forcing_cycle_start = 'none'
+config_time_varying_atmospheric_forcing_cycle_duration = '2-00-00_00:00:00'
+config_use_tidal_potential_forcing = .true.
+config_tidal_potential_reference_time = '2012-10-10_00:00:00'
+config_bottom_drag_mode = 'explicit'
+config_use_topographic_wave_drag = .true.
+config_topographic_wave_drag_coeff = 5.0e-3
+config_pressure_gradient_type = 'ssh_gradient'
+config_vert_coord_movement = 'impermeable_interfaces'
+config_ALE_thickness_proportionality = 'weights_only'
+config_disable_thick_vadv = .true.
+config_disable_thick_sflux = .true.
+config_disable_vel_hmix = .true.
+config_disable_vel_vmix = .true.
+config_disable_vel_vadv = .true.
+config_disable_tr_adv = .true.
+config_disable_tr_hmix = .true.
+config_disable_tr_vmix = .true.
+config_disable_tr_sflux = .true.
+config_disable_tr_nonlocalflux = .true.
+config_disable_redi_k33 = .true.
+config_open_ocean_freezing_temperature_coeff_0 = -1.8
+config_disable_tr_all_tend = .true.
+config_AM_pointwiseStats_enable = .true.
diff --git a/compass/ocean/tests/hurricane/fblts/forward/streams.ocean b/compass/ocean/tests/hurricane/fblts/forward/streams.ocean
new file mode 100644
index 0000000000..deb78f154b
--- /dev/null
+++ b/compass/ocean/tests/hurricane/fblts/forward/streams.ocean
@@ -0,0 +1,108 @@
+<streams>
+
+<immutable_stream name="mesh"
+                  type="input"
+                  filename_template="input.nc"
+                  input_interval="initial_only"/>
+
+<immutable_stream name="input"
+                  type="input"
+                  filename_template="input.nc"
+                  input_interval="initial_only"/>
+
+<immutable_stream name="restart"
+                  type="input;output"
+                  filename_template="restarts/restart.$Y-$M-$D_$h.$m.$s.nc"
+                  filename_interval="output_interval"
+                  reference_time="0001-01-01_00:00:00"
+                  clobber_mode="truncate"
+                  input_interval="initial_only"
+                  output_interval="0005_00:00:00"/>
+
+<immutable_stream name="atmospheric_forcing"
+                  type="input"
+                  filename_template="atmospheric_forcing.nc"
+                  filename_interval="none"
+                  packages="timeVaryingAtmosphericForcingPKG"
+                  input_interval="none"/>
+
+<stream name="pointLocationsInput"
+        filename_template="points.nc"
+        type="input"
+        input_interval="initial_only"
+        runtime_format="single_file"
+        mode="forward;analysis">
+
+    <var name="pointCellGlobalID"/>
+</stream>
+
+<stream name="topographic_wave_drag"
+                  type="input"
+                  filename_template="topographic_wave_drag.nc"
+                  filename_interval="none"
+                  packages="topographicWaveDragPKG"
+                  input_interval="initial_only">
+
+    <var name="topographic_wave_drag"/>
+</stream>
+
+<stream name="ltsregion"
+                  type="input"
+                  filename_template="input.nc"
+                  filename_interval="none"
+                  input_interval="initial_only">
+    <var name="LTSRegion"/>
+</stream>
+
+<stream name="output"
+        type="output"
+        filename_template="output.nc"
+        output_interval="12:00:00"
+        precision="double"
+        clobber_mode="truncate">
+
+    <var name="xtime"/>
+    <var name="ssh"/>
+    <var name="kineticEnergyCell"/>
+    <var name="layerThickness"/>
+    <var name="normalVelocity"/>
+    <var name="daysSinceStartOfSim"/>
+    <var name="LTSRegion"/>
+</stream>
+
+<stream name="pointwiseStatsOutput"
+        filename_template="pointwiseStats.nc"
+        type="output"
+        mode="forward;analysis"
+        output_interval="00:30:00"
+        packages="pointwiseStatsAMPKG"
+        clobber_mode="truncate"
+        io_type="netcdf"
+        runtime_format="single_file">
+
+    <var name="pointCellGlobalID"/>
+    <var name="xtime"/>
+    <var name="latCell"/>
+    <var name="lonCell"/>
+    <var name="ssh"/>
+    <var name="layerThickness"/>
+</stream>
+
+<stream name="globalStatsOutput"
+        type="none"
+        filename_template="globalStats.nc"
+        filename_interval="01-00-00_00:00:00"
+        clobber_mode="truncate"
+        packages="globalStatsAMPKG"
+        output_interval="none" >
+
+    <var_array name="minGlobalStats"/>
+    <var_array name="maxGlobalStats"/>
+    <var_array name="avgGlobalStats"/>
+    <var name="xtime"/>
+    <var name="daysSinceStartOfSim"/>
+    <var name="volumeCellGlobal"/>
+    <var name="CFLNumberGlobal"/>
+</stream>
+
+</streams>
diff --git a/compass/ocean/tests/hurricane/forward/__init__.py b/compass/ocean/tests/hurricane/forward/__init__.py
index 289b9292e7..6654e467a4 100644
--- a/compass/ocean/tests/hurricane/forward/__init__.py
+++ b/compass/ocean/tests/hurricane/forward/__init__.py
@@ -41,8 +41,10 @@ def __init__(self, test_group, mesh, storm, init, use_lts):
         """
         mesh_name = mesh.mesh_name
 
-        if use_lts:
+        if use_lts == 'LTS':
             name = f'{storm}_lts'
+        elif use_lts == 'FB_LTS':
+            name = f'{storm}_fblts'
         else:
             name = storm
         subdir = os.path.join(mesh_name, name)
diff --git a/compass/ocean/tests/hurricane/forward/forward.py b/compass/ocean/tests/hurricane/forward/forward.py
index 1d2d392d33..8ce12580e9 100644
--- a/compass/ocean/tests/hurricane/forward/forward.py
+++ b/compass/ocean/tests/hurricane/forward/forward.py
@@ -49,13 +49,22 @@ def __init__(self, test_case, mesh, init, use_lts,
 
         super().__init__(test_case=test_case, name=name)
 
-        if use_lts:
+        if use_lts == 'LTS':
 
             self.add_namelist_file(
                 'compass.ocean.tests.hurricane.lts.forward', 'namelist.ocean')
             self.add_streams_file(
                 'compass.ocean.tests.hurricane.lts.forward', 'streams.ocean')
 
+        elif use_lts == 'FB_LTS':
+
+            self.add_namelist_file(
+                'compass.ocean.tests.hurricane.fblts.forward',
+                'namelist.ocean')
+            self.add_streams_file(
+                'compass.ocean.tests.hurricane.fblts.forward',
+                'streams.ocean')
+
         else:
 
             self.add_namelist_file(
diff --git a/compass/ocean/tests/hurricane/init/__init__.py b/compass/ocean/tests/hurricane/init/__init__.py
index cc24063ae1..1c14fc59fa 100644
--- a/compass/ocean/tests/hurricane/init/__init__.py
+++ b/compass/ocean/tests/hurricane/init/__init__.py
@@ -45,8 +45,10 @@ def __init__(self, test_group, mesh, storm, use_lts):
         self.mesh = mesh
         self.use_lts = use_lts
 
-        if use_lts:
+        if use_lts == 'LTS':
             name = 'init_lts'
+        elif use_lts == 'FB_LTS':
+            name = 'init_fblts'
         else:
             name = 'init'
         mesh_name = mesh.mesh_name
diff --git a/compass/ocean/tests/hurricane/lts/mesh/lts_regions.py b/compass/ocean/tests/hurricane/lts/mesh/lts_regions.py
index 87e7025946..fbb95b0cb6 100644
--- a/compass/ocean/tests/hurricane/lts/mesh/lts_regions.py
+++ b/compass/ocean/tests/hurricane/lts/mesh/lts_regions.py
@@ -22,6 +22,7 @@ class LTSRegionsStep(Step):
         The culled mesh step containing input files to this step
     """
     def __init__(self, test_case, cull_mesh_step,
+                 num_interface_adjacent,
                  name='lts_regions', subdir='lts_regions'):
         """
         Create a new step
@@ -46,6 +47,7 @@ def __init__(self, test_case, cull_mesh_step,
             self.add_output_file(filename=file)
 
         self.cull_mesh_step = cull_mesh_step
+        self.num_interface_adjacent = num_interface_adjacent
 
     def setup(self):
         """
@@ -90,11 +92,12 @@ def run(self):
         use_progress_bar = self.log_filename is None
         label_mesh(fine_region, mesh='culled_mesh.nc',
                    graph_info='culled_graph.info', num_interface=2,
+                   num_interface_adjacent=self.num_interface_adjacent,
                    logger=self.logger, use_progress_bar=use_progress_bar)
 
 
 def label_mesh(fine_region, mesh, graph_info, num_interface,  # noqa: C901
-               logger, use_progress_bar):
+               num_interface_adjacent, logger, use_progress_bar):
 
     # read in mesh data
     ds = xr.open_dataset(mesh)
@@ -118,6 +121,7 @@ def label_mesh(fine_region, mesh, graph_info, num_interface,  # noqa: C901
 
     # first layer of cells with label 5
     changed_cells = [[], []]
+    border_cells = []
     for iedge in range(0, n_edges):
         cell1 = cells_on_edge[iedge, 0] - 1
         cell2 = cells_on_edge[iedge, 1] - 1
@@ -125,20 +129,20 @@ def label_mesh(fine_region, mesh, graph_info, num_interface,  # noqa: C901
         if (cell1 != -1 and cell2 != -1):
             if (lts_rgn[cell1] == 1 and lts_rgn[cell2] == 2):
 
-                lts_rgn[cell2] = 5
-                changed_cells[0].append(cell2)
+                lts_rgn[cell1] = 5
+                changed_cells[0].append(cell1)
 
             elif (lts_rgn[cell1] == 2 and lts_rgn[cell2] == 1):
 
-                lts_rgn[cell1] = 5
-                changed_cells[0].append(cell1)
+                lts_rgn[cell2] = 5
+                changed_cells[0].append(cell2)
+
+    border_cells = changed_cells[0]
 
-    # second and third layer of cells with label 5
+    # num_interface_adjacent - 1 layers of cells with label 5
     # only looping over cells changed during loop for previous layer
-    # at the end of this loop, changed_cells[0] will have the list of cells
-    # sharing edegs with the coarse cells
     logger.info('Labeling interface-adjacent fine cells...')
-    for i in range(0, 2):  # this loop creates 2 layers
+    for i in range(0, num_interface_adjacent - 1):
         changed_cells[(i + 1) % 2] = []
 
         for icell in changed_cells[i % 2]:
@@ -149,17 +153,19 @@ def label_mesh(fine_region, mesh, graph_info, num_interface,  # noqa: C901
                     cell2 = cells_on_edge[iedge - 1, 1] - 1
 
                     if (cell1 != -1 and cell2 != -1):
-                        if (lts_rgn[cell1] == 5 and lts_rgn[cell2] == 2):
+                        if (lts_rgn[cell1] == 5 and lts_rgn[cell2] == 1):
 
                             lts_rgn[cell2] = 5
                             changed_cells[(i + 1) % 2].append(cell2)
 
-                        elif (lts_rgn[cell1] == 2 and lts_rgn[cell2] == 5):
+                        elif (lts_rgn[cell1] == 1 and lts_rgn[cell2] == 5):
 
                             lts_rgn[cell1] = 5
                             changed_cells[(i + 1) % 2].append(cell1)
 
-    # n layers of interface region with label 4
+    changed_cells[0] = border_cells
+
+    # num_interface layers of interface region with label 4
     logger.info('Labeling interface cells...')
     for i in range(0, num_interface):
         changed_cells[(i + 1) % 2] = []
@@ -200,7 +206,7 @@ def label_mesh(fine_region, mesh, graph_info, num_interface,  # noqa: C901
 
     changed_cells[0] = changed_cells[num_interface % 2]
 
-    # n layers of interface region with label 3
+    # num_interface layers of interface region with label 3
     for i in range(0, num_interface):
         changed_cells[(i + 1) % 2] = []
 
@@ -213,9 +219,9 @@ def label_mesh(fine_region, mesh, graph_info, num_interface,  # noqa: C901
 
                     if (cell1 != -1 and cell2 != -1):
                         # for the first layer, need to check neighbors are
-                        # 5 and 2
-                        # for further layers, need to check neighbors are
                         # 3 and 2
+                        # for further layers, need to check neighbors are
+                        # 4 and 2
                         if (i == 0):
                             if (lts_rgn[cell1] == 3 and lts_rgn[cell2] == 2):
 
diff --git a/compass/ocean/tests/hurricane/mesh/__init__.py b/compass/ocean/tests/hurricane/mesh/__init__.py
index 23f7964604..94f8188d53 100644
--- a/compass/ocean/tests/hurricane/mesh/__init__.py
+++ b/compass/ocean/tests/hurricane/mesh/__init__.py
@@ -29,8 +29,10 @@ def __init__(self, test_group, mesh_name, use_lts):
         self.mesh_name = mesh_name
         self.use_lts = use_lts
 
-        if use_lts:
+        if use_lts == 'LTS':
             name = 'mesh_lts'
+        elif use_lts == 'FB_LTS':
+            name = 'mesh_fblts'
         else:
             name = 'mesh'
         subdir = '{}/{}'.format(mesh_name, name)
@@ -60,10 +62,14 @@ def __init__(self, test_group, mesh_name, use_lts):
 
         self.add_step(cull_mesh_step)
 
-        if use_lts:
-
+        if use_lts == 'LTS':
+            self.add_step(LTSRegionsStep(
+                          test_case=self, cull_mesh_step=cull_mesh_step,
+                          num_interface_adjacent=2))
+        elif use_lts == 'FB_LTS':
             self.add_step(LTSRegionsStep(
-                          test_case=self, cull_mesh_step=cull_mesh_step))
+                          test_case=self, cull_mesh_step=cull_mesh_step,
+                          num_interface_adjacent=10))
 
     def configure(self):
         """
diff --git a/docs/users_guide/ocean/test_groups/hurricane.rst b/docs/users_guide/ocean/test_groups/hurricane.rst
index be50d7e36e..b9976022d7 100644
--- a/docs/users_guide/ocean/test_groups/hurricane.rst
+++ b/docs/users_guide/ocean/test_groups/hurricane.rst
@@ -19,13 +19,16 @@ which can be used to simulate coastal inundation using MPAS-Ocean's
 wetting and drying scheme.
 
 The time stepping options to run the simulations include the fourth
-order Runge-Kutta scheme (RK4), which is the default, and a local
-time-stepping (LTS) scheme, see 
+order Runge-Kutta scheme (RK4), and two local time-stepping schemes.
+The first LTS scheme is based on a strong stability preserving Runge-Kutta
+scheme of order three and is called LTS3, see 
 `Lilly et al. (2023) <https://doi.org/10.1029/2022MS003327>`_
-for details. Each test case in the ``ocean/hurricane`` test group has an 
-LTS counterpart which is identified by the ``_lts`` appendix in the test case name.
-Note that the executable to be used with LTS has to be compiled with
-the following options: ``USE_LTS=true OPENMP=false``.
+for details.
+The second LTS scheme is based on a forward-backward Runge-Kutta scheme
+of order two and is called FB-LTS.
+Each test case in the ``ocean/hurricane`` test group has a counterpart 
+for each LTS scheme which is identified by appending the test case name
+with ``_lts`` for LTS3 and ``_fblts`` for FB-LTS.
 
 Shared config options
 ---------------------
@@ -96,7 +99,7 @@ since this data is not carried over from the cell culling process.
 
 .. _hurricane_mesh_lts:
 
-If the LTS option is selected for the mesh test case, an additional step
+If either LTS option is selected for the mesh test case, an additional step
 is carried out after the mesh culling. This step appropriately flags 
 the cells of the mesh according to a user defined criterion in order to
 use time-steps of different sizes on different regions of the mesh.
@@ -136,7 +139,7 @@ pointWiseStats analysis member for the forward run.
 
 compute topographic wave drag step
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-This step is carried out only if the LTS option is selected for the init test case.
+This step is carried out only if either LTS option is selected for the init test case.
 
 The reciprocal of the e-folding time, ``r_inv``, from the HyCOM model,
 is computed in this step. See 
@@ -165,6 +168,6 @@ is used for the validation.
 
 .. _hurricane_sandy_lts:
 
-If the LTS option is selected for the sandy test case, the LTS scheme
+If either LTS option is selected for the sandy test case, the LTS scheme
 is used to advance the solution in time rather than the default RK4 scheme.