diff --git a/docs/src/library/public.md b/docs/src/library/public.md
index 53379acf..1973d4a1 100644
--- a/docs/src/library/public.md
+++ b/docs/src/library/public.md
@@ -39,9 +39,9 @@ Modules = [ClimaSeaIce.Rheologies]
 Private = false
 ```
 
-## ClimaSeaIce.SeaIceMomentumEquations
+## ClimaSeaIce.SeaIceDynamics
 
 ```@autodocs
-Modules = [ClimaSeaIce.SeaIceMomentumEquations]
+Modules = [ClimaSeaIce.SeaIceDynamics]
 Private = false
 ```
diff --git a/examples/ice_advected_by_anticyclone.jl b/examples/ice_advected_by_anticyclone.jl
index e71984d5..9298ac21 100644
--- a/examples/ice_advected_by_anticyclone.jl
+++ b/examples/ice_advected_by_anticyclone.jl
@@ -4,7 +4,7 @@
 #
 #
 using ClimaSeaIce
-using ClimaSeaIce.SeaIceMomentumEquations
+using ClimaSeaIce.SeaIceDynamics
 using ClimaSeaIce.Rheologies
 using Oceananigans
 using Oceananigans.Units
@@ -86,18 +86,13 @@ fill_halo_regions!((Uₐ, Vₐ))
 # We use an elasto-visco-plastic rheology and WENO seventh order 
 # for advection of h and ℵ
 
-momentum_equations = SeaIceMomentumEquation(grid; 
-                                            top_momentum_stress = (u=τₐu, v=τₐv),
-                                            bottom_momentum_stress = τₒ,
-                                            coriolis = FPlane(f=1e-4),
-                                            ocean_velocities = (u = Uₒ, v = Vₒ),
-                                            rheology = ElastoViscoPlasticRheology(),
-                                            solver   = SplitExplicitSolver(substeps=150))
-
-# Define the model!
+dynamics = SeaIceMomentumEquation(grid; 
+                                  top_momentum_stress = (u=τₐu, v=τₐv),
+                                  bottom_momentum_stress = τₒ,
+                                  coriolis     = FPlane(f=1e-4))
 
 model = SeaIceModel(grid; 
-                    dynamics = momentum_equations,
+                    dynamics,
                     ice_thermodynamics = nothing, # No ice_thermodynamics here
                     advection = WENO(order=7),
                     boundary_conditions = (u=u_bcs, v=v_bcs))
diff --git a/examples/ice_advected_on_coastline.jl b/examples/ice_advected_on_coastline.jl
index 734d466e..8e238d26 100644
--- a/examples/ice_advected_on_coastline.jl
+++ b/examples/ice_advected_on_coastline.jl
@@ -1,5 +1,5 @@
 using ClimaSeaIce
-using ClimaSeaIce.SeaIceMomentumEquations
+using ClimaSeaIce.SeaIceDynamics
 using ClimaSeaIce.Rheologies
 using Oceananigans
 using Oceananigans.Units
diff --git a/src/ClimaSeaIce.jl b/src/ClimaSeaIce.jl
index 4ca3488d..abb6acb3 100644
--- a/src/ClimaSeaIce.jl
+++ b/src/ClimaSeaIce.jl
@@ -39,7 +39,7 @@ include("forward_euler_timestepper.jl")
 
 include("SeaIceThermodynamics/SeaIceThermodynamics.jl")
 include("Rheologies/Rheologies.jl")
-include("SeaIceMomentumEquations/SeaIceMomentumEquations.jl")
+include("SeaIceDynamics/SeaIceDynamics.jl")
 include("sea_ice_model.jl")
 include("sea_ice_advection.jl")
 include("tracer_tendency_kernel_functions.jl")
@@ -47,7 +47,7 @@ include("sea_ice_time_stepping.jl")
 include("EnthalpyMethodSeaIceModel.jl")
 
 using .SeaIceThermodynamics
-using .SeaIceMomentumEquations
+using .SeaIceDynamics
 using .Rheologies
 
 # Advection timescale for a `SeaIceModel`. Sea ice dynamics are two-dimensional so 
diff --git a/src/Rheologies/Rheologies.jl b/src/Rheologies/Rheologies.jl
index 276ff632..dc82048b 100644
--- a/src/Rheologies/Rheologies.jl
+++ b/src/Rheologies/Rheologies.jl
@@ -20,6 +20,11 @@ compute_stresses!(model, dynamics, rheology, Δt) = nothing
 @inline sum_of_forcing_u(i, j, k, grid, rheology, u_forcing, fields, Δt) = u_forcing(i, j, k, grid, fields)
 @inline sum_of_forcing_v(i, j, k, grid, rheology, v_forcing, fields, Δt) = v_forcing(i, j, k, grid, fields)
 
+@inline ice_stress_ux(i, j, k, grid, ::Nothing, args...) = zero(grid)
+@inline ice_stress_uy(i, j, k, grid, ::Nothing, args...) = zero(grid)
+@inline ice_stress_vx(i, j, k, grid, ::Nothing, args...) = zero(grid)
+@inline ice_stress_vy(i, j, k, grid, ::Nothing, args...) = zero(grid)
+
 include("ice_stress_divergence.jl")
 include("viscous_rheology.jl")
 include("elasto_visco_plastic_rheology.jl")
diff --git a/src/SeaIceMomentumEquations/SeaIceMomentumEquations.jl b/src/SeaIceDynamics/SeaIceDynamics.jl
similarity index 94%
rename from src/SeaIceMomentumEquations/SeaIceMomentumEquations.jl
rename to src/SeaIceDynamics/SeaIceDynamics.jl
index 1f4b600b..4782e835 100644
--- a/src/SeaIceMomentumEquations/SeaIceMomentumEquations.jl
+++ b/src/SeaIceDynamics/SeaIceDynamics.jl
@@ -1,8 +1,8 @@
-module SeaIceMomentumEquations
+module SeaIceDynamics
 
 # The only functions provided by the module
 export compute_momentum_tendencies!, time_step_momentum!
-export SeaIceMomentumEquation, ExplicitSolver, SplitExplicitSolver, SemiImplicitStress
+export SeaIceMomentumEquation, ExplicitSolver, SplitExplicitSolver, SemiImplicitStress, StressBalanceFreeDrift
 
 using ClimaSeaIce
 
@@ -46,8 +46,9 @@ import Oceananigans: fields
 time_step_momentum!(model, dynamics, Δt) = nothing
 compute_momentum_tendencies!(model, dynamics, Δt) = nothing
 
-include("sea_ice_momentum_equations.jl")
 include("sea_ice_external_stress.jl")
+include("stress_balance_free_drift.jl")
+include("sea_ice_momentum_equations.jl")
 include("momentum_tendencies_kernel_functions.jl")
 include("explicit_momentum_equations.jl")
 include("split_explicit_momentum_equations.jl")
diff --git a/src/SeaIceMomentumEquations/explicit_momentum_equations.jl b/src/SeaIceDynamics/explicit_momentum_equations.jl
similarity index 92%
rename from src/SeaIceMomentumEquations/explicit_momentum_equations.jl
rename to src/SeaIceDynamics/explicit_momentum_equations.jl
index 00dab8d9..d83ab0fe 100644
--- a/src/SeaIceMomentumEquations/explicit_momentum_equations.jl
+++ b/src/SeaIceDynamics/explicit_momentum_equations.jl
@@ -18,7 +18,7 @@ function time_step_momentum!(model, ::ExplicitMomentumEquation, Δt)
                          ρ = model.ice_density))
 
 
-    ocean_velocities = dynamics.ocean_velocities
+    free_drift = dynamics.free_drift
     clock = model.clock
     minimum_mass = dynamics.minimum_mass
     minimum_concentration = dynamics.minimum_concentration
@@ -27,7 +27,7 @@ function time_step_momentum!(model, ::ExplicitMomentumEquation, Δt)
     bottom_stress = dynamics.external_momentum_stresses.bottom
 
     launch!(arch, grid, :xy, _step_velocities!, u, v, grid, Gⁿ, Δt, 
-            top_stress, bottom_stress, ocean_velocities, 
+            top_stress, bottom_stress, free_drift, 
             minimum_mass, minimum_concentration, clock, model_fields)
 
     return nothing
@@ -35,7 +35,7 @@ end
 
 @kernel function _step_velocities!(u, v, grid, Gⁿ, Δt, 
                                    top_stress, bottom_stress, 
-                                   ocean_velocities, minimum_mass, minimum_concentration, clock, fields)
+                                   free_drift, minimum_mass, minimum_concentration, clock, fields)
 
     i, j = @index(Global, NTuple)
     kᴺ   = size(grid, 3)
@@ -56,8 +56,8 @@ end
         uᴰ = (u[i, j, 1] + Δt * Gⁿ.u[i, j, 1]) / (1 + Δt * τuᵢ)
         vᴰ = (v[i, j, 1] + Δt * Gⁿ.v[i, j, 1]) / (1 + Δt * τvᵢ)
 
-        uᶠ = free_drift_u(i, j, kᴺ, grid, ocean_velocities)
-        vᶠ = free_drift_v(i, j, kᴺ, grid, ocean_velocities)
+        uᶠ = free_drift_u(i, j, kᴺ, grid, free_drift, clock, fields)
+        vᶠ = free_drift_v(i, j, kᴺ, grid, free_drift, clock, fields)
 
         sea_ice = (mᶠᶜ ≥ minimum_mass) & (ℵᶠᶜ ≥ minimum_concentration)
         u[i, j, 1] = ifelse(sea_ice, uᴰ, uᶠ)
diff --git a/src/SeaIceMomentumEquations/momentum_tendencies_kernel_functions.jl b/src/SeaIceDynamics/momentum_tendencies_kernel_functions.jl
similarity index 100%
rename from src/SeaIceMomentumEquations/momentum_tendencies_kernel_functions.jl
rename to src/SeaIceDynamics/momentum_tendencies_kernel_functions.jl
diff --git a/src/SeaIceMomentumEquations/sea_ice_external_stress.jl b/src/SeaIceDynamics/sea_ice_external_stress.jl
similarity index 100%
rename from src/SeaIceMomentumEquations/sea_ice_external_stress.jl
rename to src/SeaIceDynamics/sea_ice_external_stress.jl
diff --git a/src/SeaIceMomentumEquations/sea_ice_momentum_equations.jl b/src/SeaIceDynamics/sea_ice_momentum_equations.jl
similarity index 65%
rename from src/SeaIceMomentumEquations/sea_ice_momentum_equations.jl
rename to src/SeaIceDynamics/sea_ice_momentum_equations.jl
index 842c3627..c46b4925 100644
--- a/src/SeaIceMomentumEquations/sea_ice_momentum_equations.jl
+++ b/src/SeaIceDynamics/sea_ice_momentum_equations.jl
@@ -1,12 +1,12 @@
 using ClimaSeaIce.Rheologies
 using Adapt
 
-struct SeaIceMomentumEquation{S, C, R, V, A, ES, FT}
+struct SeaIceMomentumEquation{S, C, R, F, A, ES, FT}
     coriolis :: C
     rheology :: R
     auxiliary_fields :: A
     solver :: S
-    ocean_velocities :: V
+    free_drift :: F
     external_momentum_stresses :: ES
     minimum_concentration :: FT
     minimum_mass :: FT
@@ -17,16 +17,19 @@ struct ExplicitSolver end
 
 """
     SeaIceMomentumEquation(grid; 
-                           coriolis=nothing,
-                           rheology=ElastoViscoPlasticRheology(eltype(grid)),
-                           auxiliary_fields=NamedTuple(),
-                           ocean_velocities=nothing,
-                           solver=ExplicitSolver(),
-                           minimum_concentration=1e-3,
-                           minimum_mass=1.0)
+                           coriolis = nothing,
+                           rheology = ElastoViscoPlasticRheology(eltype(grid)),
+                           auxiliary_fields = NamedTuple(),
+                           top_momentum_stress    = (u=nothing, v=nothing),
+                           bottom_momentum_stress = (u=nothing, v=nothing),
+                           free_drift = StressBalanceFreeDrift(top_momentum_stress, bottom_momentum_stress),
+                           solver = SplitExplicitSolver(150),
+                           minimum_concentration = 1e-3,
+                           minimum_mass = 1.0)
 
 Constructs a `SeaIceMomentumEquation` object that controls the dynamical evolution of sea-ice momentum.
 The sea-ice momentum obey the following evolution equation:
+
 ```math
     ∂u                   τₒ    τₐ
     -- + f x u = ∇ ⋅ σ + --  + -- 
@@ -46,8 +49,8 @@ Keyword Arguments
 - `coriolis`: Parameters for the background rotation rate of the model.
 - `rheology`: The sea ice rheology model, default is `ElastoViscoPlasticRheology(eltype(grid))`.
 - `auxiliary_fields`: A named tuple of auxiliary fields, default is an empty `NamedTuple()`.
-- `ocean_velocities`: The ocean surface velocities used to limit the sea ice momentum when the mass or the concentration are
-                      below a certain threshold. default is `nothing` (indicating that the free drift velocities are zero).
+- `free_drift`: The free drift velocities used to limit sea ice momentum when the mass or the concentration are
+                below a certain threshold. default is `nothing` (indicating that the free drift velocities are zero).
 - `solver`: The momentum solver to be used.
 - `minimum_concentration`: The minimum sea ice concentration above which the sea ice velocity is dynamically calculated, default is `1e-3`.
 - `minimum_mass`: The minimum sea ice mass per area above which the sea ice velocity is dynamically calculated, default is `1.0 kg/m²`.
@@ -56,10 +59,10 @@ function SeaIceMomentumEquation(grid;
                                 coriolis = nothing,
                                 rheology = ElastoViscoPlasticRheology(eltype(grid)),
                                 auxiliary_fields = NamedTuple(),
-                                ocean_velocities = nothing,
-                                top_momentum_stress    = (u = nothing, v = nothing),
-                                bottom_momentum_stress = (u = nothing, v = nothing),
-                                solver = ExplicitSolver(),
+                                top_momentum_stress    = nothing,
+                                bottom_momentum_stress = nothing,
+                                free_drift = StressBalanceFreeDrift(top_momentum_stress, bottom_momentum_stress),
+                                solver = SplitExplicitSolver(150),
                                 minimum_concentration = 1e-3,
                                 minimum_mass = 1.0)
 
@@ -73,18 +76,10 @@ function SeaIceMomentumEquation(grid;
                                   rheology, 
                                   auxiliary_fields, 
                                   solver,
-                                  ocean_velocities,
+                                  free_drift,
                                   external_momentum_stresses,
                                   convert(FT, minimum_concentration),
                                   convert(FT, minimum_mass))
 end
 
 fields(mom::SeaIceMomentumEquation) = mom.auxiliary_fields
-
-# Just passing ocean velocities without mitigation
-@inline free_drift_u(i, j, k, grid, f) = @inbounds f.u[i, j, k] 
-@inline free_drift_v(i, j, k, grid, f) = @inbounds f.v[i, j, k] 
-
-# Passing no velocities
-@inline free_drift_u(i, j, k, grid, ::Nothing) = zero(grid)
-@inline free_drift_v(i, j, k, grid, ::Nothing) = zero(grid)
diff --git a/src/SeaIceMomentumEquations/split_explicit_momentum_equations.jl b/src/SeaIceDynamics/split_explicit_momentum_equations.jl
similarity index 90%
rename from src/SeaIceMomentumEquations/split_explicit_momentum_equations.jl
rename to src/SeaIceDynamics/split_explicit_momentum_equations.jl
index 9d30c0fe..39884aeb 100644
--- a/src/SeaIceMomentumEquations/split_explicit_momentum_equations.jl
+++ b/src/SeaIceDynamics/split_explicit_momentum_equations.jl
@@ -35,7 +35,7 @@ function time_step_momentum!(model, dynamics::SplitExplicitMomentumEquation, Δt
 
     u, v = model.velocities
   
-    ocean_velocities = dynamics.ocean_velocities
+    free_drift = dynamics.free_drift
     clock = model.clock
     coriolis = dynamics.coriolis
 
@@ -55,8 +55,10 @@ function time_step_momentum!(model, dynamics::SplitExplicitMomentumEquation, Δt
                          ℵ = model.ice_concentration, 
                          ρ = model.ice_density))
 
-    u_velocity_kernel!, _ = configure_kernel(arch, grid, :xy, _u_velocity_step!)
-    v_velocity_kernel!, _ = configure_kernel(arch, grid, :xy, _v_velocity_step!)
+    active_cells_map = Oceananigans.Grids.get_active_column_map(grid)
+
+    u_velocity_kernel!, _ = configure_kernel(arch, grid, :xy, _u_velocity_step!; active_cells_map)
+    v_velocity_kernel!, _ = configure_kernel(arch, grid, :xy, _v_velocity_step!; active_cells_map)
 
     substeps = dynamics.solver.substeps
     
@@ -73,24 +75,24 @@ function time_step_momentum!(model, dynamics::SplitExplicitMomentumEquation, Δt
         # In odd substeps we switch and calculate vⁿ⁺¹ = f(uⁿ) and uⁿ⁺¹ = f(vⁿ⁺¹).
         if iseven(substep) 
             u_velocity_kernel!(u, grid, Δt, substeps, rheology, model_fields, 
-                               ocean_velocities, clock, coriolis,
+                               free_drift, clock, coriolis,
                                minimum_mass, minimum_concentration, 
                                u_immersed_bc, top_stress, bottom_stress, u_forcing)
 
             v_velocity_kernel!(v, grid, Δt, substeps, rheology, model_fields, 
-                               ocean_velocities, clock, coriolis, 
+                               free_drift, clock, coriolis, 
                                minimum_mass, minimum_concentration,
                                v_immersed_bc, top_stress, bottom_stress, v_forcing)
 
         else
             v_velocity_kernel!(v, grid, Δt, substeps, rheology, model_fields, 
-                               ocean_velocities, clock, coriolis, 
+                               free_drift, clock, coriolis, 
                                minimum_mass, minimum_concentration,
                                v_immersed_bc, top_stress, bottom_stress, v_forcing)
             
 
             u_velocity_kernel!(u, grid, Δt, substeps, rheology, model_fields, 
-                               ocean_velocities, clock, coriolis,
+                               free_drift, clock, coriolis,
                                minimum_mass, minimum_concentration, 
                                u_immersed_bc, top_stress, bottom_stress, u_forcing)
         end
@@ -107,7 +109,7 @@ end
 
 @kernel function _u_velocity_step!(u, grid, Δt, 
                                    substeps, rheology, 
-                                   model_fields, ocean_velocities, 
+                                   model_fields, free_drift, 
                                    clock, coriolis, 
                                    minimum_mass, minimum_concentration,
                                    u_immersed_bc, u_top_stress, u_bottom_stress, u_forcing)
@@ -127,7 +129,7 @@ end
 
     τuᵢ = ifelse(mᵢ ≤ 0, zero(grid), τuᵢ)
     uᴰ  = @inbounds (u[i, j, 1] + Δτ * Gu) / (1 + Δτ * τuᵢ) # dynamical velocity 
-    uᶠ  = free_drift_u(i, j, kᴺ, grid, ocean_velocities) # free drift velocity
+    uᶠ  = free_drift_u(i, j, kᴺ, grid, free_drift, clock, model_fields) # free drift velocity
 
     # If the ice mass or the ice concentration are below a certain threshold, 
     # the sea ice velocity is set to the free drift velocity
@@ -138,7 +140,7 @@ end
 
 @kernel function _v_velocity_step!(v, grid, Δt, 
                                    substeps, rheology, 
-                                   model_fields, ocean_velocities, 
+                                   model_fields, free_drift, 
                                    clock, coriolis, 
                                    minimum_mass, minimum_concentration,
                                    v_immersed_bc, v_top_stress, v_bottom_stress, v_forcing)
@@ -159,7 +161,7 @@ end
     τvᵢ = ifelse(mᵢ ≤ 0, zero(grid), τvᵢ)
 
     vᴰ = @inbounds (v[i, j, 1] + Δτ * Gv) / (1 + Δτ * τvᵢ)# dynamical velocity 
-    vᶠ = free_drift_v(i, j, kᴺ, grid, ocean_velocities) # free drift velocity
+    vᶠ = free_drift_v(i, j, kᴺ, grid, free_drift, clock, model_fields)  # free drift velocity
 
     # If the ice mass or the ice concentration are below a certain threshold, 
     # the sea ice velocity is set to the free drift velocity
diff --git a/src/SeaIceDynamics/stress_balance_free_drift.jl b/src/SeaIceDynamics/stress_balance_free_drift.jl
new file mode 100644
index 00000000..f56d79d7
--- /dev/null
+++ b/src/SeaIceDynamics/stress_balance_free_drift.jl
@@ -0,0 +1,119 @@
+abstract type AbstractFreeDriftDynamics end
+
+"""
+    StressBalanceFreeDrift{T, B}
+
+A free drift parameterization that computes the free drift velocities as a balance between top and bottom stresses ``τa ≈ τo``.
+
+In case the only one of the stresses is a `SemiImplicitStress`, the model will compute the free drift velocity exactly 
+assuming that the other stress does not depend on the sea ice velocity. Stresses other than `SemiImplicitStress` 
+are assumed to be ice-velocity independent.
+
+Can be used to limit the sea ice velocity when the mass or the concentration are below a certain threshold, or
+as a `dynamics` model itself that substitutes the sea ice momentum equation calculation everywhere.
+"""
+struct StressBalanceFreeDrift{T, B} <: AbstractFreeDriftDynamics
+    top_momentum_stress :: T
+    bottom_momentum_stress :: B
+end
+
+Adapt.adapt_structure(to, s::StressBalanceFreeDrift) = 
+    StressBalanceFreeDrift(Adapt.adapt(to, s.top_momentum_stress),
+                           Adapt.adapt(to, s.bottom_momentum_stress))
+
+fields(::StressBalanceFreeDrift) = NamedTuple()
+
+# Stress balance when either the top or the bottom stresses do not depend on ice velocity
+# In this case we have a simplified form of the free drift velocity
+# Otherwise, to avoid a nonlinear solve, we assume the stress is only lineary dependent on the velocity at time-step
+# n+1 and use the ice velocities at time-step n to compute the nonlinear term. 
+# Note that this is the same formulation we use to solve for stresses in the `SeaIceMomentumEquation` dynamics.
+const TISB = StressBalanceFreeDrift{<:Any, <:SemiImplicitStress}
+const BISB = StressBalanceFreeDrift{<:SemiImplicitStress, <:Any}
+
+# Stress balance when only the bottom stress is ice-velocity dependent:
+# Then: 𝒰ᵢ = 𝒰ᴮ - τᵀ / sqrt(Cᴮ * ||τᵀ||)
+@inline function free_drift_u(i, j, k, grid, f::TISB, clock, fields) 
+    τxᵀ = x_momentum_stress(i, j, k, grid, f.top_momentum_stress, clock, fields)
+    τyᵀ = ℑxyᶠᶜᵃ(i, j, k, grid, y_momentum_stress, f.top_momentum_stress, clock, fields)
+    τᵀ  = sqrt(τxᵀ^2 + τyᵀ^2)
+
+    τᴮ = f.bottom_momentum_stress
+    uᴮ = @inbounds τᴮ.uₑ[i, j, k]
+    Cᴮ = τᴮ.ρₑ * τᴮ.Cᴰ
+
+    return uᴮ - ifelse(τᵀ == 0, τᵀ, τxᵀ / sqrt(Cᴮ * τᵀ))
+end
+
+@inline function free_drift_v(i, j, k, grid, f::TISB, clock, fields) 
+    τxᵀ = ℑxyᶜᶠᵃ(i, j, k, grid, x_momentum_stress, f.top_momentum_stress, clock, fields)
+    τyᵀ = y_momentum_stress(i, j, k, grid, f.top_momentum_stress, clock, fields)
+    τᵀ  = sqrt(τxᵀ^2 + τyᵀ^2)
+
+    τᴮ = f.bottom_momentum_stress
+    vᴮ = @inbounds τᴮ.vₑ[i, j, k]
+    Cᴮ = τᴮ.ρₑ * τᴮ.Cᴰ
+
+    return vᴮ - ifelse(τᵀ == 0, τᵀ, τyᵀ / sqrt(Cᴮ * τᵀ))
+end
+
+# Stress balance when only the top stress is ice-velocity dependent:
+# Then: 𝒰ᵢ = 𝒰ᵀ - τᴮ / sqrt(Cᵀ * ||τᴮ||)
+@inline function free_drift_u(i, j, k, grid, f::BISB, clock, fields) 
+    τxᴮ = x_momentum_stress(i, j, k, grid, f.bottom_momentum_stress, clock, fields)
+    τyᴮ = ℑxyᶠᶜᵃ(i, j, k, grid, y_momentum_stress, f.bottom_momentum_stress, clock, fields)
+    τᴮ  = sqrt(τxᴮ^2 + τyᴮ^2)
+
+    τᵀ = f.top_momentum_stess
+    uᵀ = @inbounds τᵀ.uₑ[i, j, k]
+    Cᵀ = τᵀ.ρₑ * τᵀ.Cᴰ
+
+    return uᵀ - ifelse(τᴮ == 0, τᴮ, τxᴮ / sqrt(Cᵀ * τᴮ))
+end
+
+@inline function free_drift_v(i, j, k, grid, f::BISB, clock, fields) 
+    τxᴮ = ℑxyᶜᶠᵃ(i, j, k, grid, x_momentum_stress, f.bottom_momentum_stress, clock, fields)
+    τyᴮ = y_momentum_stress(i, j, k, grid, f.bottom_momentum_stress, clock, fields)
+    τᴮ  = sqrt(τxᴮ^2 + τyᴮ^2)
+
+    τᵀ = f.top_momentum_stess
+    vᵀ = @inbounds τᵀ.vₑ[i, j, k]
+    Cᵀ = τᵀ.ρₑ * τᵀ.Cᴰ
+
+    return vᵀ - ifelse(τᴮ == 0, τᴮ, τyᴮ / sqrt(Cᵀ * τᴮ))
+end
+
+const NoFreeDrift = StressBalanceFreeDrift{<:Nothing, <:Nothing}
+
+@inline free_drift_u(i, j, k, grid, ::NoFreeDrift, clock, fields) = zero(grid)
+@inline free_drift_v(i, j, k, grid, ::NoFreeDrift, clock, fields) = zero(grid)
+
+# Fallbacks for a given velocity field.
+@inline free_drift_u(i, j, k, grid, f::NamedTuple, clock, fields)  = @inbounds f.u[i, j, k] 
+@inline free_drift_v(i, j, k, grid, f::NamedTuple, clock, fields)  = @inbounds f.v[i, j, k] 
+
+# Passing no velocities
+@inline free_drift_u(i, j, k, grid, ::Nothing, clock, fields) = zero(grid)
+@inline free_drift_v(i, j, k, grid, ::Nothing, clock, fields) = zero(grid)
+
+# What if we want to use _only_ the free drift velocities? (not advised)
+function time_step_momentum!(model, dynamics::AbstractFreeDriftDynamics, args...)
+    
+    model_fields = fields(model)
+    clock = model.clock
+    grid  = model.grid
+    arch  = architecture(grid)
+    u, v  = model.velocities
+
+    launch!(arch, grid, :xy, _free_drift_velocity_step!, u, v, grid, dynamics, clock, model_fields)
+
+    return nothing
+end
+
+@kernel function _free_drift_velocity_step!(u, v, grid, dynamics, clock, fields)
+    i, j = @index(Global, NTuple)
+    kᴺ   = size(grid, 3)
+
+    @inbounds u[i, j, 1] = free_drift_u(i, j, kᴺ, grid, dynamics, clock, fields)
+    @inbounds v[i, j, 1] = free_drift_v(i, j, kᴺ, grid, dynamics, clock, fields)
+end
diff --git a/src/sea_ice_time_stepping.jl b/src/sea_ice_time_stepping.jl
index df4ef8ba..d041aa26 100644
--- a/src/sea_ice_time_stepping.jl
+++ b/src/sea_ice_time_stepping.jl
@@ -3,7 +3,7 @@ using Oceananigans.Fields: flattened_unique_values
 using Oceananigans.OutputReaders: extract_field_time_series, update_field_time_series!
 using Oceananigans.ImmersedBoundaries: mask_immersed_field_xy!
 
-using ClimaSeaIce.SeaIceMomentumEquations: time_step_momentum!
+using ClimaSeaIce.SeaIceDynamics: time_step_momentum!
 using ClimaSeaIce.SeaIceThermodynamics: thermodynamic_time_step!
 
 import Oceananigans.Models: update_model_field_time_series!
diff --git a/src/tracer_tendency_kernel_functions.jl b/src/tracer_tendency_kernel_functions.jl
index 12272436..fb3fb196 100644
--- a/src/tracer_tendency_kernel_functions.jl
+++ b/src/tracer_tendency_kernel_functions.jl
@@ -1,5 +1,5 @@
 using Oceananigans.Advection
-using ClimaSeaIce.SeaIceMomentumEquations: compute_momentum_tendencies!
+using ClimaSeaIce.SeaIceDynamics: compute_momentum_tendencies!
 
 function compute_tendencies!(model::SIM, Δt)
     compute_tracer_tendencies!(model)
diff --git a/test/test_sea_ice_advection.jl b/test/test_sea_ice_advection.jl
index 6f129106..112d284a 100644
--- a/test/test_sea_ice_advection.jl
+++ b/test/test_sea_ice_advection.jl
@@ -1,5 +1,5 @@
 using ClimaSeaIce
-using ClimaSeaIce.SeaIceMomentumEquations
+using ClimaSeaIce.SeaIceDynamics
 using ClimaSeaIce.Rheologies
 using Oceananigans