diff --git a/docs/make.jl b/docs/make.jl
index 85a1476d..91d58d9b 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -17,7 +17,7 @@ example_scripts = [
     "melting_in_spring.jl",
     "freezing_of_a_lake.jl",
     "ice_advected_by_anticyclone.jl",
-    # "ice_advected_on_coastline.jl",
+    "ice_advected_on_coastline.jl",
     "arctic_basin_seasonal_cycle.jl"
 ]
 
@@ -31,7 +31,7 @@ example_pages = [
     "Melting in Spring" => "literated/melting_in_spring.md",
     "Freezing of a Lake" => "literated/freezing_of_a_lake.md",
     "Ice advected by anticyclone" => "literated/ice_advected_by_anticyclone.md",
-    # "Ice advected on coastline" => "literated/ice_advected_on_coastline.md",
+    "Ice advected on coastline" => "literated/ice_advected_on_coastline.md",
     "Arctic basin seasonal cycle" => "literated/arctic_basin_seasonal_cycle.md"
 ]
 
diff --git a/examples/ice_advected_on_coastline.jl b/examples/ice_advected_on_coastline.jl
index 8e238d26..5e42870e 100644
--- a/examples/ice_advected_on_coastline.jl
+++ b/examples/ice_advected_on_coastline.jl
@@ -64,15 +64,31 @@ dynamics = SeaIceMomentumEquation(grid;
                                   rheology = ElastoViscoPlasticRheology(),
                                   solver = SplitExplicitSolver(substeps=150))
 
-u_bcs = FieldBoundaryConditions(top = nothing, bottom = nothing,
+@inline immersed_u_drag(i, j, k, grid, clock, fields, D) = @inbounds - D * fields.u[i, j, k]
+@inline immersed_v_drag(i, j, k, grid, clock, fields, D) = @inbounds - D * fields.v[i, j, k]
+
+immersed_u_bc = FluxBoundaryCondition(immersed_u_drag, discrete_form=true, parameters=3e-1)
+immersed_v_bc = FluxBoundaryCondition(immersed_v_drag, discrete_form=true, parameters=3e-1)
+
+immersed_u_bc = ImmersedBoundaryCondition(top=nothing, bottom=nothing, west=nothing, east=nothing, 
+                                          south=immersed_u_bc, north=immersed_u_bc)
+
+immersed_v_bc = ImmersedBoundaryCondition(top=nothing, bottom=nothing, south=nothing, north=nothing,
+                                          west=immersed_v_bc, east=immersed_v_bc)
+
+u_bcs = FieldBoundaryConditions(grid, (Face, Center, Nothing); 
                                 north = ValueBoundaryCondition(0),
-                                south = ValueBoundaryCondition(0))
+                                south = ValueBoundaryCondition(0),
+                                immersed = immersed_u_bc)
+
+v_bcs = FieldBoundaryConditions(grid, (Center, Face, Nothing); 
+                                immersed = immersed_v_bc)
 
 #Define the model! 
 model = SeaIceModel(grid; 
                     advection = WENO(order=7),
                     dynamics = dynamics,
-                    boundary_conditions = (; u=u_bcs),
+                    boundary_conditions = (; u=u_bcs, v=v_bcs),
                     ice_thermodynamics = nothing)
 
 # We start with a concentration of ℵ = 1 everywhere
@@ -159,4 +175,4 @@ CairoMakie.record(fig, "sea_ice_advected_on_coastline.mp4", 1:Nt, framerate = 8)
 end
 nothing #hide
 
-# ![](sea_ice_advected_on_coastline.mp4)
\ No newline at end of file
+# ![](sea_ice_advected_on_coastline.mp4)
diff --git a/src/Rheologies/ice_stress_divergence.jl b/src/Rheologies/ice_stress_divergence.jl
index 3515c1f0..be49b5b7 100644
--- a/src/Rheologies/ice_stress_divergence.jl
+++ b/src/Rheologies/ice_stress_divergence.jl
@@ -1,6 +1,7 @@
 using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, 
                                        ImmersedBoundaryCondition
 
+using Oceananigans.BoundaryConditions: FBC, getbc
 using Oceananigans.Advection:conditional_flux_ccc, conditional_flux_ffc
 using Oceananigans.Operators: index_left, index_right
 
@@ -52,44 +53,58 @@ end
 
 @inline function immersed_∂ⱼ_σ₁ⱼ(i, j, k, ibg::IBG, u_bc::IBC, rheology, clock, fields)
     # Fetch fluxes across immersed boundary
-    q̃ᵂ = ib_ice_stress_ux(i,   j,   k, ibg, u_bc.west,  rheology, clock, fields)
-    q̃ᴱ = ib_ice_stress_ux(i+1, j,   k, ibg, u_bc.east,  rheology, clock, fields)
-    q̃ˢ = ib_ice_stress_uy(i,   j-1, k, ibg, u_bc.south, rheology, clock, fields)
-    q̃ᴺ = ib_ice_stress_uy(i,   j,   k, ibg, u_bc.north, rheology, clock, fields)
+    q̃ᵂ =  ib_ice_stress_ux_west(i,   j,   k, ibg, u_bc.west,  rheology, clock, fields)
+    q̃ᴱ =  ib_ice_stress_ux_east(i+1, j,   k, ibg, u_bc.east,  rheology, clock, fields)
+    q̃ˢ = ib_ice_stress_uy_south(i,   j-1, k, ibg, u_bc.south, rheology, clock, fields)
+    q̃ᴺ = ib_ice_stress_uy_north(i,   j,   k, ibg, u_bc.north, rheology, clock, fields)
 
     iᵂ, jˢ, _ = map(index_left,  (i, j, k), (c, c, c)) # Broadcast instead of map causes inference failure
     iᴱ, jᴺ, _ = map(index_right, (i, j, k), (f, f, c))
 
-    # Impose i) immersed fluxes if we're on an immersed boundary or ii) zero otherwise.
-    qᵂ = conditional_flux_ccc(iᵂ, j, k, ibg, q̃ᵂ, zero(ibg)) * Axᶜᶜᶜ(iᵂ, j, k, grid)
-    qᴱ = conditional_flux_ccc(iᴱ, j, k, ibg, q̃ᴱ, zero(ibg)) * Axᶜᶜᶜ(iᴱ, j, k, grid)
-    qˢ = conditional_flux_ffc(i, jˢ, k, ibg, q̃ˢ, zero(ibg)) * Ayᶠᶠᶜ(i, jˢ, k, grid)
-    qᴺ = conditional_flux_ffc(i, jᴺ, k, ibg, q̃ᴺ, zero(ibg)) * Ayᶠᶠᶜ(i, jᴺ, k, grid)
+    # Impose i) immersed fluxes if we're on an immersed boundary or ii) zero oibgwise.
+    qᵂ = conditional_flux_ccc(iᵂ, j, k, ibg, q̃ᵂ, zero(ibg)) * Axᶜᶜᶜ(iᵂ, j, k, ibg)
+    qᴱ = conditional_flux_ccc(iᴱ, j, k, ibg, q̃ᴱ, zero(ibg)) * Axᶜᶜᶜ(iᴱ, j, k, ibg)
+    qˢ = conditional_flux_ffc(i, jˢ, k, ibg, q̃ˢ, zero(ibg)) * Ayᶠᶠᶜ(i, jˢ, k, ibg)
+    qᴺ = conditional_flux_ffc(i, jᴺ, k, ibg, q̃ᴺ, zero(ibg)) * Ayᶠᶠᶜ(i, jᴺ, k, ibg)
 
-    return (qᴱ - qᵂ + qᴺ - qˢ) / Vᶠᶜᶜ(i, j, k, grid)
+    return (qᴱ - qᵂ + qᴺ - qˢ) / Vᶠᶜᶜ(i, j, k, ibg)
 end
 
 @inline function immersed_∂ⱼ_σ₂ⱼ(i, j, k, ibg::IBG, v_bc::IBC, rheology, clock, fields)
     # Fetch fluxes across immersed boundary
-    q̃ᵂ = ib_ice_stress_vx(i-1, j,   k, ibg, v_bc.west,  rheology, clock, fields)
-    q̃ᴱ = ib_ice_stress_vx(i,   j,   k, ibg, v_bc.east,  rheology, clock, fields)
-    q̃ˢ = ib_ice_stress_vy(i,   j,   k, ibg, v_bc.south, rheology, clock, fields)
-    q̃ᴺ = ib_ice_stress_vy(i,   j+1, k, ibg, v_bc.north, rheology, clock, fields)
+    q̃ᵂ =  ib_ice_stress_vx_west(i-1, j,   k, ibg, v_bc.west,  rheology, clock, fields)
+    q̃ᴱ =  ib_ice_stress_vx_east(i,   j,   k, ibg, v_bc.east,  rheology, clock, fields)
+    q̃ˢ = ib_ice_stress_vy_south(i,   j-1, k, ibg, v_bc.south, rheology, clock, fields)
+    q̃ᴺ = ib_ice_stress_vy_north(i,   j,   k, ibg, v_bc.north, rheology, clock, fields)
 
     iᵂ, jˢ, _ = map(index_left,  (i, j, k), (f, f, c)) # Broadcast instead of map causes inference failure
     iᴱ, jᴺ, _ = map(index_right, (i, j, k), (c, c, c))
 
     # Impose i) immersed fluxes if we're on an immersed boundary or ii) zero otherwise.
-    qᵂ = conditional_flux_ffc(iᵂ, j, k, ibg, q̃ᵂ, zero(ibg)) * Axᶠᶠᶜ(iᵂ, j, k, grid)
-    qᴱ = conditional_flux_ffc(iᴱ, j, k, ibg, q̃ᴱ, zero(ibg)) * Axᶠᶠᶜ(iᴱ, j, k, grid)
-    qˢ = conditional_flux_ccc(i, jˢ, k, ibg, q̃ˢ, zero(ibg)) * Ayᶜᶜᶜ(i, jˢ, k, grid)
-    qᴺ = conditional_flux_ccc(i, jᴺ, k, ibg, q̃ᴺ, zero(ibg)) * Ayᶜᶜᶜ(i, jᴺ, k, grid)
+    qᵂ = conditional_flux_ffc(iᵂ, j, k, ibg, q̃ᵂ, zero(ibg)) * Axᶠᶠᶜ(iᵂ, j, k, ibg)
+    qᴱ = conditional_flux_ffc(iᴱ, j, k, ibg, q̃ᴱ, zero(ibg)) * Axᶠᶠᶜ(iᴱ, j, k, ibg)
+    qˢ = conditional_flux_ccc(i, jˢ, k, ibg, q̃ˢ, zero(ibg)) * Ayᶜᶜᶜ(i, jˢ, k, ibg)
+    qᴺ = conditional_flux_ccc(i, jᴺ, k, ibg, q̃ᴺ, zero(ibg)) * Ayᶜᶜᶜ(i, jᴺ, k, ibg)
 
-    return (qᴱ - qᵂ + qᴺ - qˢ) / Vᶜᶠᶜ(i, j, k, grid)
+    return (qᴱ - qᵂ + qᴺ - qˢ) / Vᶜᶠᶜ(i, j, k, ibg)
 end
 
 # TODO: Implement immersed fluxes (0 for the moment)
-@inline ib_ice_stress_ux(i, j, k, grid, args...) = zero(grid)
-@inline ib_ice_stress_vx(i, j, k, grid, args...) = zero(grid)
-@inline ib_ice_stress_uy(i, j, k, grid, args...) = zero(grid)
-@inline ib_ice_stress_vy(i, j, k, grid, args...) = zero(grid)
\ No newline at end of file
+@inline ib_ice_stress_ux_west(i, j, k, ibg, args...) = zero(ibg)
+@inline ib_ice_stress_ux_east(i, j, k, ibg, args...) = zero(ibg)
+@inline ib_ice_stress_vx_west(i, j, k, ibg, args...) = zero(ibg)
+@inline ib_ice_stress_vx_east(i, j, k, ibg, args...) = zero(ibg)
+@inline ib_ice_stress_uy_south(i, j, k, ibg, args...) = zero(ibg)
+@inline ib_ice_stress_uy_north(i, j, k, ibg, args...) = zero(ibg)
+@inline ib_ice_stress_vy_south(i, j, k, ibg, args...) = zero(ibg)
+@inline ib_ice_stress_vy_north(i, j, k, ibg, args...) = zero(ibg)
+
+# Only supporting FluxBoundaryConditions for now
+@inline ib_ice_stress_ux_west(i, j, k, ibg, bc::FBC, rheology, args...) = + getbc(bc, i, j, k, ibg, args...)
+@inline ib_ice_stress_ux_east(i, j, k, ibg, bc::FBC, rheology, args...) = - getbc(bc, i, j, k, ibg, args...)
+@inline ib_ice_stress_vx_west(i, j, k, ibg, bc::FBC, rheology, args...) = + getbc(bc, i, j, k, ibg, args...)
+@inline ib_ice_stress_vx_east(i, j, k, ibg, bc::FBC, rheology, args...) = - getbc(bc, i, j, k, ibg, args...)
+@inline ib_ice_stress_uy_south(i, j, k, ibg, bc::FBC, rheology, args...) = + getbc(bc, i, j, k, ibg, args...)
+@inline ib_ice_stress_uy_north(i, j, k, ibg, bc::FBC, rheology, args...) = - getbc(bc, i, j, k, ibg, args...)
+@inline ib_ice_stress_vy_south(i, j, k, ibg, bc::FBC, rheology, args...) = + getbc(bc, i, j, k, ibg, args...)
+@inline ib_ice_stress_vy_north(i, j, k, ibg, bc::FBC, rheology, args...) = - getbc(bc, i, j, k, ibg, args...)