diff --git a/amr-wind/equation_systems/icns/icns_advection.cpp b/amr-wind/equation_systems/icns/icns_advection.cpp
index 6bffd88035..8a715a33ea 100644
--- a/amr-wind/equation_systems/icns/icns_advection.cpp
+++ b/amr-wind/equation_systems/icns/icns_advection.cpp
@@ -55,6 +55,10 @@ MacProjOp::MacProjOp(
 {
     amrex::ParmParse pp("incflo");
     pp.query("density", m_rho_0);
+    amrex::ParmParse pp_ovst("Overset");
+    bool disable_ovst_mac = false;
+    pp_ovst.query("disable_coupled_mac_proj", disable_ovst_mac);
+    m_has_overset = m_has_overset && !disable_ovst_mac;
 }
 
 void MacProjOp::init_projector(const MacProjOp::FaceFabPtrVec& beta) noexcept
diff --git a/amr-wind/equation_systems/vof/volume_fractions.H b/amr-wind/equation_systems/vof/volume_fractions.H
index b67988a7aa..26ab361ba9 100644
--- a/amr-wind/equation_systems/vof/volume_fractions.H
+++ b/amr-wind/equation_systems/vof/volume_fractions.H
@@ -16,27 +16,27 @@ namespace amr_wind::multiphase {
  *  to zero.
  */
 
-AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void youngs_fd_normal(
-    int i,
-    int j,
-    int k,
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void youngs_finite_difference_normal(
+    const int i,
+    const int j,
+    const int k,
     amrex::Array4<amrex::Real const> const& volfrac,
     amrex::Real& mx,
     amrex::Real& my,
     amrex::Real& mz) noexcept
 {
-    amrex::Real mm1, mm2;
-
-    mm1 = volfrac(i - 1, j - 1, k - 1) + volfrac(i - 1, j - 1, k + 1) +
-          volfrac(i - 1, j + 1, k - 1) + volfrac(i - 1, j + 1, k + 1) +
-          2.0 * (volfrac(i - 1, j - 1, k) + volfrac(i - 1, j + 1, k) +
-                 volfrac(i - 1, j, k - 1) + volfrac(i - 1, j, k + 1)) +
-          4.0 * volfrac(i - 1, j, k);
-    mm2 = volfrac(i + 1, j - 1, k - 1) + volfrac(i + 1, j - 1, k + 1) +
-          volfrac(i + 1, j + 1, k - 1) + volfrac(i + 1, j + 1, k + 1) +
-          2.0 * (volfrac(i + 1, j - 1, k) + volfrac(i + 1, j + 1, k) +
-                 volfrac(i + 1, j, k - 1) + volfrac(i + 1, j, k + 1)) +
-          4.0 * volfrac(i + 1, j, k);
+    amrex::Real mm1 =
+        volfrac(i - 1, j - 1, k - 1) + volfrac(i - 1, j - 1, k + 1) +
+        volfrac(i - 1, j + 1, k - 1) + volfrac(i - 1, j + 1, k + 1) +
+        2.0 * (volfrac(i - 1, j - 1, k) + volfrac(i - 1, j + 1, k) +
+               volfrac(i - 1, j, k - 1) + volfrac(i - 1, j, k + 1)) +
+        4.0 * volfrac(i - 1, j, k);
+    amrex::Real mm2 =
+        volfrac(i + 1, j - 1, k - 1) + volfrac(i + 1, j - 1, k + 1) +
+        volfrac(i + 1, j + 1, k - 1) + volfrac(i + 1, j + 1, k + 1) +
+        2.0 * (volfrac(i + 1, j - 1, k) + volfrac(i + 1, j + 1, k) +
+               volfrac(i + 1, j, k - 1) + volfrac(i + 1, j, k + 1)) +
+        4.0 * volfrac(i + 1, j, k);
     mx = mm1 - mm2;
 
     mm1 = volfrac(i - 1, j - 1, k - 1) + volfrac(i - 1, j - 1, k + 1) +
@@ -64,25 +64,82 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void youngs_fd_normal(
     mz = mm1 - mm2;
 }
 
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
+youngs_finite_difference_normal_neumann(
+    const int i,
+    const int j,
+    const int k,
+    const int ibdy,
+    const int jbdy,
+    const int kbdy,
+    amrex::Array4<amrex::Real const> const& volfrac,
+    amrex::Real& mx,
+    amrex::Real& my,
+    amrex::Real& mz) noexcept
+{
+    // Do neumann condition via indices
+    const int im1 = ibdy == -1 ? i : i - 1;
+    const int jm1 = jbdy == -1 ? j : j - 1;
+    const int km1 = kbdy == -1 ? k : k - 1;
+    const int ip1 = ibdy == +1 ? i : i + 1;
+    const int jp1 = jbdy == +1 ? j : j + 1;
+    const int kp1 = kbdy == +1 ? k : k + 1;
+
+    amrex::Real mm1 = volfrac(im1, jm1, km1) + volfrac(im1, jm1, kp1) +
+                      volfrac(im1, jp1, km1) + volfrac(im1, jp1, kp1) +
+                      2.0 * (volfrac(im1, jm1, k) + volfrac(im1, jp1, k) +
+                             volfrac(im1, j, km1) + volfrac(im1, j, kp1)) +
+                      4.0 * volfrac(im1, j, k);
+    amrex::Real mm2 = volfrac(ip1, jm1, km1) + volfrac(ip1, jm1, kp1) +
+                      volfrac(ip1, jp1, km1) + volfrac(ip1, jp1, kp1) +
+                      2.0 * (volfrac(ip1, jm1, k) + volfrac(ip1, jp1, k) +
+                             volfrac(ip1, j, km1) + volfrac(ip1, j, kp1)) +
+                      4.0 * volfrac(ip1, j, k);
+    mx = mm1 - mm2;
+
+    mm1 = volfrac(im1, jm1, km1) + volfrac(im1, jm1, kp1) +
+          volfrac(ip1, jm1, km1) + volfrac(ip1, jm1, kp1) +
+          2.0 * (volfrac(im1, jm1, k) + volfrac(ip1, jm1, k) +
+                 volfrac(i, jm1, km1) + volfrac(i, jm1, kp1)) +
+          4.0 * volfrac(i, jm1, k);
+    mm2 = volfrac(im1, jp1, km1) + volfrac(im1, jp1, kp1) +
+          volfrac(ip1, jp1, km1) + volfrac(ip1, jp1, kp1) +
+          2.0 * (volfrac(im1, jp1, k) + volfrac(ip1, jp1, k) +
+                 volfrac(i, jp1, km1) + volfrac(i, jp1, kp1)) +
+          4.0 * volfrac(i, jp1, k);
+    my = mm1 - mm2;
+
+    mm1 = volfrac(im1, jm1, km1) + volfrac(im1, jp1, km1) +
+          volfrac(ip1, jm1, km1) + volfrac(ip1, jp1, km1) +
+          2.0 * (volfrac(im1, j, km1) + volfrac(ip1, j, km1) +
+                 volfrac(i, jm1, km1) + volfrac(i, jp1, km1)) +
+          4.0 * volfrac(i, j, km1);
+    mm2 = volfrac(im1, jm1, kp1) + volfrac(im1, jp1, kp1) +
+          volfrac(ip1, jm1, kp1) + volfrac(ip1, jp1, kp1) +
+          2.0 * (volfrac(im1, j, kp1) + volfrac(ip1, j, kp1) +
+                 volfrac(i, jm1, kp1) + volfrac(i, jp1, kp1)) +
+          4.0 * volfrac(i, j, kp1);
+    mz = mm1 - mm2;
+}
+
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mixed_youngs_central_normal(
-    int i,
-    int j,
-    int k,
+    const int i,
+    const int j,
+    const int k,
     amrex::Array4<amrex::Real const> const& volfrac,
     amrex::Real& mx,
     amrex::Real& my,
     amrex::Real& mz) noexcept
 {
     amrex::Array2D<amrex::Real, 0, AMREX_SPACEDIM + 1, 0, AMREX_SPACEDIM> m;
-    amrex::Real m1, m2;
     // write the plane as: sgn(mx) X =  my Y +  mz Z + alpha
     //                           m00 X = m01 Y + m02 Z + alpha
-    m1 = volfrac(i - 1, j, k - 1) + volfrac(i - 1, j, k + 1) +
-         volfrac(i - 1, j - 1, k) + volfrac(i - 1, j + 1, k) +
-         volfrac(i - 1, j, k);
-    m2 = volfrac(i + 1, j, k - 1) + volfrac(i + 1, j, k + 1) +
-         volfrac(i + 1, j - 1, k) + volfrac(i + 1, j + 1, k) +
-         volfrac(i + 1, j, k);
+    amrex::Real m1 = volfrac(i - 1, j, k - 1) + volfrac(i - 1, j, k + 1) +
+                     volfrac(i - 1, j - 1, k) + volfrac(i - 1, j + 1, k) +
+                     volfrac(i - 1, j, k);
+    amrex::Real m2 = volfrac(i + 1, j, k - 1) + volfrac(i + 1, j, k + 1) +
+                     volfrac(i + 1, j - 1, k) + volfrac(i + 1, j + 1, k) +
+                     volfrac(i + 1, j, k);
 
     m(0, 0) = (m1 > m2) ? 1.0 : -1.0;
 
@@ -180,7 +237,8 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mixed_youngs_central_normal(
     }
 
     // Youngs-CIAM scheme */
-    youngs_fd_normal(i, j, k, volfrac, m(3, 0), m(3, 1), m(3, 2));
+    youngs_finite_difference_normal(
+        i, j, k, volfrac, m(3, 0), m(3, 1), m(3, 2));
     // normalize the set (mx,my,mz): |mx|+|my|+|mz| = 1
     constexpr amrex::Real tiny = 1e-20;
     t0 = std::abs(m(3, 0)) + std::abs(m(3, 1)) + std::abs(m(3, 2)) + tiny;
@@ -214,7 +272,10 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mixed_youngs_central_normal(
  * given that m1+m2+m3=1 (m1,m2,m3>0) and the volumetric fraction volF
  */
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real volume_intercept(
-    amrex::Real b1, amrex::Real b2, amrex::Real b3, amrex::Real volF) noexcept
+    const amrex::Real b1,
+    const amrex::Real b2,
+    const amrex::Real b3,
+    const amrex::Real volF) noexcept
 {
     using namespace amrex;
 
@@ -296,15 +357,15 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real volume_intercept(
  * (5) calculate volume (NOTE: it is assumed:s0=t0=0; ds0=dt0=1.)
  */
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real cut_volume(
-    amrex::Real m1,
-    amrex::Real m2,
-    amrex::Real m3,
-    amrex::Real alpha,
-    amrex::Real r0,
-    amrex::Real dr0) noexcept
+    const amrex::Real m1,
+    const amrex::Real m2,
+    const amrex::Real m3,
+    const amrex::Real alpha,
+    const amrex::Real r0,
+    const amrex::Real dr0) noexcept
 {
 
-    amrex::Real const_tiny = std::numeric_limits<amrex::Real>::epsilon();
+    const amrex::Real const_tiny = std::numeric_limits<amrex::Real>::epsilon();
     amrex::Real al;
 
     // move origin to x0
@@ -428,11 +489,11 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void fit_plane(
 }
 
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool interface_band(
-    int i,
-    int j,
-    int k,
+    const int i,
+    const int j,
+    const int k,
     amrex::Array4<amrex::Real const> const& volfrac,
-    int n_band = 1) noexcept
+    const int n_band = 1) noexcept
 {
     // n_band must be <= number of vof ghost cells (3)
     constexpr amrex::Real tiny = 1e-12;
@@ -456,14 +517,14 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool interface_band(
 }
 
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real levelset_to_vof(
-    int i,
-    int j,
-    int k,
-    amrex::Real eps,
+    const int i,
+    const int j,
+    const int k,
+    const amrex::Real eps,
     amrex::Array4<amrex::Real const> const& phi) noexcept
 {
     amrex::Real mx, my, mz;
-    youngs_fd_normal(i, j, k, phi, mx, my, mz);
+    youngs_finite_difference_normal(i, j, k, phi, mx, my, mz);
     mx = std::abs(mx / 32.);
     my = std::abs(my / 32.);
     mz = std::abs(mz / 32.);
diff --git a/amr-wind/incflo.cpp b/amr-wind/incflo.cpp
index 750cfe2e0e..4cc56b030a 100644
--- a/amr-wind/incflo.cpp
+++ b/amr-wind/incflo.cpp
@@ -319,6 +319,9 @@ void incflo::do_advance()
     } else {
         advance();
     }
+    if (m_sim.has_overset()) {
+        m_ovst_ops.post_advance_work();
+    }
 }
 
 // Make a new level from scratch using provided BoxArray and
diff --git a/amr-wind/overset/OversetOps.H b/amr-wind/overset/OversetOps.H
index 1220e99570..0658019d4d 100644
--- a/amr-wind/overset/OversetOps.H
+++ b/amr-wind/overset/OversetOps.H
@@ -13,10 +13,52 @@ public:
 
     void initialize(CFDSim& sim);
     void pre_advance_work();
+    void post_advance_work();
 
     void update_gradp();
 
 private:
+    // Functions called within public functions
+    void parameter_output() const;
+    void sharpen_nalu_data();
+    void set_hydrostatic_gradp();
+    void replace_masked_gradp();
+
+    // Check for multiphase sim
+    bool m_vof_exists{false};
+    // Check for perturbational pressure
+    // (will be removed soon)
+    bool m_perturb_p{false};
+
+    // This is the only option for now
+    const bool m_use_hydrostatic_gradp{true};
+
+    // Coupling options
+    bool m_disable_nodal_proj{false};
+    bool m_disable_mac_proj{false};
+    bool m_replace_gradp_postsolve{false};
+
+    // Verbosity
+    int m_verbose{0};
+
+    // Reinitialization parameters
+    int m_n_iterations{10};
+    int m_calc_convg_interval{1}; // calctolniter, cconv
+    amrex::Real m_convg_tol = 1e-12;
+    amrex::Real m_relative_length_scale = 1.5;
+    amrex::Real m_upw_margin = 0.1;
+    amrex::Real m_target_cutoff = 0.0; // proc_tgvof_tol
+
+    // Tolerance for VOF-related bound checks
+    const amrex::Real m_vof_tol = 1e-12;
+    // Small number for approximate signed distance function
+    const amrex::Real m_asdf_tiny = 1e-12;
+
+    // Pointer for pressure gradient copy field
+    amr_wind::Field* m_gp_copy{nullptr};
+    // Pointer for MultiPhase physics
+    amr_wind::MultiPhase* m_mphase{nullptr};
+
     CFDSim* m_sim_ptr;
 };
 
diff --git a/amr-wind/overset/OversetOps.cpp b/amr-wind/overset/OversetOps.cpp
index 31f2f1e94e..3c03ec7498 100644
--- a/amr-wind/overset/OversetOps.cpp
+++ b/amr-wind/overset/OversetOps.cpp
@@ -1,4 +1,5 @@
 #include "amr-wind/overset/OversetOps.H"
+#include "amr-wind/overset/overset_ops_routines.H"
 #include "amr-wind/core/field_ops.H"
 #include "AMReX_ParmParse.H"
 #include "AMReX_MultiFabUtil.H"
@@ -8,12 +9,67 @@
 
 namespace amr_wind {
 
-void OversetOps::initialize(CFDSim& sim) { m_sim_ptr = &sim; }
+void OversetOps::initialize(CFDSim& sim)
+{
+    m_sim_ptr = &sim;
+    // Queries for reinitialization options
+    amrex::ParmParse pp("Overset");
+    pp.query("reinit_iterations", m_n_iterations);
+    pp.query("reinit_convg_interval", m_calc_convg_interval);
+    pp.query("reinit_convg_tolerance", m_convg_tol);
+    pp.query("reinit_rlscale", m_relative_length_scale);
+    pp.query("reinit_upw_margin", m_upw_margin);
+    pp.query("reinit_target_cutoff", m_target_cutoff);
+
+    // Queries for coupling options
+    pp.query("replace_gradp_postsolve", m_replace_gradp_postsolve);
+    // OversetOps does not control these coupling options, merely reports them
+    pp.query("disable_coupled_nodal_proj", m_disable_nodal_proj);
+    pp.query("disable_coupled_mac_proj", m_disable_mac_proj);
+
+    pp.query("verbose", m_verbose);
+
+    amrex::ParmParse pp_icns("ICNS");
+    pp_icns.query("use_perturb_pressure", m_perturb_p);
+
+    m_vof_exists = (*m_sim_ptr).repo().field_exists("vof");
+    if (m_vof_exists) {
+        m_mphase = &(*m_sim_ptr).physics_manager().get<MultiPhase>();
+    }
+    if (m_replace_gradp_postsolve) {
+        m_gp_copy = &(*m_sim_ptr).repo().declare_field("gp_copy", 3);
+    }
+
+    parameter_output();
+}
 
 void OversetOps::pre_advance_work()
 {
-    // Update pressure gradient using updated overset pressure field
-    update_gradp();
+    // Pressure gradient not updated for current multiphase approach
+    if (!(m_vof_exists && m_use_hydrostatic_gradp)) {
+        // Update pressure gradient using updated overset pressure field
+        update_gradp();
+    }
+
+    if (m_vof_exists) {
+        // Reinitialize fields
+        sharpen_nalu_data();
+        if (m_use_hydrostatic_gradp) {
+            // Use hydrostatic pressure gradient
+            set_hydrostatic_gradp();
+        }
+    }
+
+    // If pressure gradient will be replaced, store current pressure gradient
+    if (m_replace_gradp_postsolve) {
+        auto& gp = (*m_sim_ptr).repo().get_field("gp");
+        for (int lev = 0; lev < (*m_sim_ptr).repo().num_active_levels();
+             ++lev) {
+            amrex::MultiFab::Copy(
+                (*m_gp_copy)(lev), gp(lev), 0, 0, gp(lev).nComp(),
+                (m_gp_copy)->num_grow());
+        }
+    }
 }
 
 void OversetOps::update_gradp()
@@ -73,4 +129,269 @@ void OversetOps::update_gradp()
     // Averaging down here would be unnecessary; it is built into calcGradPhi
 }
 
+void OversetOps::post_advance_work()
+{
+    // Replace and reapply pressure gradient if requested
+    if (m_replace_gradp_postsolve) {
+        replace_masked_gradp();
+    }
+}
+
+/* ----------------------------------------------- */
+/* PUBLIC FUNCTIONS ABOVE, PRIVATE FUNCTIONS BELOW */
+/* ----------------------------------------------- */
+
+void OversetOps::parameter_output() const
+{
+    // Print the details
+    if (m_verbose > 0) {
+        // Important parameters
+        amrex::Print() << "\nOverset Coupling Parameters: \n"
+                       << "---- Coupled nodal projection : "
+                       << !m_disable_nodal_proj << std::endl
+                       << "---- Coupled MAC projection   : "
+                       << !m_disable_mac_proj << std::endl
+                       << "---- Replace overset pres grad: "
+                       << m_replace_gradp_postsolve << std::endl;
+        if (m_vof_exists) {
+            amrex::Print() << "Overset Reinitialization Parameters:\n"
+                           << "---- Maximum iterations   : " << m_n_iterations
+                           << std::endl
+                           << "---- Convergence tolerance: " << m_convg_tol
+                           << std::endl
+                           << "---- Relative length scale: "
+                           << m_relative_length_scale << std::endl
+                           << "---- Upwinding VOF margin : " << m_upw_margin
+                           << std::endl;
+            if (m_verbose > 1) {
+                // Less important or less used parameters
+                amrex::Print()
+                    << "---- Calc. conv. interval : " << m_calc_convg_interval
+                    << std::endl
+                    << "---- Target field cutoff  : " << m_target_cutoff
+                    << std::endl;
+            }
+        }
+        amrex::Print() << std::endl;
+    }
+}
+
+void OversetOps::sharpen_nalu_data()
+{
+    const auto& repo = (*m_sim_ptr).repo();
+    const auto nlevels = repo.num_active_levels();
+    const auto geom = (*m_sim_ptr).mesh().Geom();
+
+    // Get blanking for cells
+    const auto& iblank_cell = repo.get_int_field("iblank_cell");
+
+    // Get fields that will be modified
+    auto& vof = repo.get_field("vof");
+    auto& levelset = repo.get_field("levelset");
+    auto& rho = repo.get_field("density");
+    auto& velocity = repo.get_field("velocity");
+
+    // Create scratch fields for fluxes
+    // 5 components are vof, density, and 3 of velocity
+    auto flux_x = repo.create_scratch_field(5, 0, amr_wind::FieldLoc::XFACE);
+    auto flux_y = repo.create_scratch_field(5, 0, amr_wind::FieldLoc::YFACE);
+    auto flux_z = repo.create_scratch_field(5, 0, amr_wind::FieldLoc::ZFACE);
+    // Create scratch field for approximate signed distance function and grad
+    // (components 0-2 are gradient, 3 is asdf)
+    auto normal_vec = repo.create_scratch_field(3, vof.num_grow()[0] - 1);
+
+    auto target_vof = repo.create_scratch_field(1, vof.num_grow()[0]);
+
+    // Create levelset field
+    for (int lev = 0; lev < nlevels; ++lev) {
+        // Thickness used here is user parameter, whatever works best
+        auto dx = (geom[lev]).CellSizeArray();
+        const amrex::Real i_th =
+            m_relative_length_scale * std::cbrt(dx[0] * dx[1] * dx[2]);
+
+        // Populate approximate signed distance function
+        overset_ops::populate_psi(levelset(lev), vof(lev), i_th, m_asdf_tiny);
+    }
+    amrex::Gpu::synchronize();
+
+    // Convert levelset to vof to get target_vof
+    m_mphase->levelset2vof(iblank_cell, *target_vof);
+
+    // Process target vof for tiny margins from single-phase
+    for (int lev = 0; lev < nlevels; ++lev) {
+        // A tolerance of 0 should do nothing
+        overset_ops::process_vof((*target_vof)(lev), m_target_cutoff);
+    }
+    amrex::Gpu::synchronize();
+
+    // Replace vof with original values in amr domain
+    for (int lev = 0; lev < nlevels; ++lev) {
+        overset_ops::harmonize_vof(
+            (*target_vof)(lev), vof(lev), iblank_cell(lev));
+    }
+    amrex::Gpu::synchronize();
+
+    // Put fluxes in vector for averaging down during iterations
+    amrex::Vector<amrex::Array<amrex::MultiFab*, AMREX_SPACEDIM>> fluxes(
+        repo.num_active_levels());
+    for (int lev = 0; lev < nlevels; ++lev) {
+        fluxes[lev][0] = &(*flux_x)(lev);
+        fluxes[lev][1] = &(*flux_y)(lev);
+        fluxes[lev][2] = &(*flux_z)(lev);
+    }
+
+    // Pseudo-time loop
+    amrex::Real err = 100.0 * m_convg_tol;
+    int n = 0;
+    while (n < m_n_iterations && err > m_convg_tol) {
+        ++n;
+        bool calc_convg = n % m_calc_convg_interval == 0;
+        // Zero error if being calculated this step
+        err = calc_convg ? 0.0 : err;
+
+        // Maximum possible value of pseudo time factor (dtau)
+        amrex::Real ptfac = 1.0;
+        // Maximum pseudoCFL, 0.5 seems to work well
+        const amrex::Real pCFL = 0.5;
+
+        for (int lev = 0; lev < nlevels; ++lev) {
+            // Populate normal vector
+            overset_ops::populate_normal_vector(
+                (*normal_vec)(lev), vof(lev), iblank_cell(lev));
+
+            // Sharpening fluxes for vof, density, and momentum
+            overset_ops::populate_sharpen_fluxes(
+                (*flux_x)(lev), (*flux_y)(lev), (*flux_z)(lev), vof(lev),
+                (*target_vof)(lev), (*normal_vec)(lev), velocity(lev),
+                m_upw_margin, m_mphase->rho1(), m_mphase->rho2());
+
+            // Process fluxes
+            overset_ops::process_fluxes(
+                (*flux_x)(lev), (*flux_y)(lev), (*flux_z)(lev),
+                iblank_cell(lev));
+
+            // Measure convergence to determine if loop can stop
+            if (calc_convg) {
+                // Update error at specified interval of steps
+                const amrex::Real err_lev = overset_ops::measure_convergence(
+                    (*flux_x)(lev), (*flux_y)(lev), (*flux_z)(lev));
+                err = amrex::max(err, err_lev);
+            }
+        }
+        amrex::Gpu::synchronize();
+
+        // Average down fluxes across levels for consistency
+        for (int lev = nlevels - 1; lev > 0; --lev) {
+            amrex::IntVect rr =
+                geom[lev].Domain().size() / geom[lev - 1].Domain().size();
+            amrex::average_down_faces(
+                GetArrOfConstPtrs(fluxes[lev]), fluxes[lev - 1], rr,
+                geom[lev - 1]);
+        }
+
+        // Get pseudo dt (dtau)
+        for (int lev = 0; lev < nlevels; ++lev) {
+            // Compare vof fluxes to vof in source cells
+            // Convergence tolerance determines what size of fluxes matter
+            const amrex::Real ptfac_lev = overset_ops::calculate_pseudo_dt_flux(
+                (*flux_x)(lev), (*flux_y)(lev), (*flux_z)(lev), vof(lev),
+                m_convg_tol);
+            ptfac = amrex::min(ptfac, ptfac_lev);
+        }
+        amrex::ParallelDescriptor::ReduceRealMin(ptfac);
+
+        // Conform pseudo dt (dtau) to pseudo CFL
+        ptfac = pCFL * ptfac;
+
+        // Apply fluxes
+        for (int lev = 0; lev < nlevels; ++lev) {
+            overset_ops::apply_fluxes(
+                (*flux_x)(lev), (*flux_y)(lev), (*flux_z)(lev), vof(lev),
+                rho(lev), velocity(lev), ptfac, m_vof_tol);
+        }
+        amrex::Gpu::synchronize();
+
+        // Fillpatch for ghost cells
+        vof.fillpatch((*m_sim_ptr).time().current_time());
+        velocity.fillpatch((*m_sim_ptr).time().current_time());
+
+        // Update density (fillpatch built in)
+        m_mphase->set_density_via_vof();
+
+        // Ensure that err is same across processors
+        if (calc_convg) {
+            amrex::ParallelDescriptor::ReduceRealMax(err);
+        }
+
+        if (m_verbose > 0) {
+            amrex::Print() << "OversetOps: sharpen step " << n << "  conv. err "
+                           << err << "  tol " << m_convg_tol << std::endl;
+        }
+    }
+
+    // Purely for debugging via visualization, should be removed later
+    // Currently set up to overwrite the levelset field (not used as time
+    // evolves) with the post-sharpening velocity magnitude
+    for (int lev = 0; lev < nlevels; ++lev) {
+        overset_ops::equate_field(levelset(lev), velocity(lev));
+    }
+    amrex::Gpu::synchronize();
+}
+
+void OversetOps::set_hydrostatic_gradp()
+{
+    const auto& repo = (*m_sim_ptr).repo();
+    const auto nlevels = repo.num_active_levels();
+    const auto geom = (*m_sim_ptr).mesh().Geom();
+
+    // Get blanking for cells
+    const auto& iblank_cell = repo.get_int_field("iblank_cell");
+
+    // Get fields that will be modified or used
+    Field* rho0{nullptr};
+    auto& rho = repo.get_field("density");
+    auto& gp = repo.get_field("gp");
+    if (m_perturb_p) {
+        rho0 = &((*m_sim_ptr).repo().get_field("reference_density"));
+    } else {
+        // Point to existing field, won't be used
+        rho0 = &rho;
+    }
+
+    // Replace initial gp with best guess (hydrostatic)
+    for (int lev = 0; lev < nlevels; ++lev) {
+        overset_ops::replace_gradp_hydrostatic(
+            gp(lev), rho(lev), (*rho0)(lev), iblank_cell(lev),
+            m_mphase->gravity()[2], m_perturb_p);
+    }
+}
+
+void OversetOps::replace_masked_gradp()
+{
+    const auto& repo = (*m_sim_ptr).repo();
+    const auto nlevels = repo.num_active_levels();
+
+    // Get timestep
+    const amrex::Real dt = (*m_sim_ptr).time().deltaT();
+
+    // Get blanking for cells
+    const auto& iblank_cell = repo.get_int_field("iblank_cell");
+
+    // Get fields that will be modified or used
+    auto& vel = repo.get_field("velocity");
+    auto& rho = repo.get_field("density");
+    auto& gp = repo.get_field("gp");
+
+    // For iblanked cells, replace gp with original gp, to get original vel
+    for (int lev = 0; lev < nlevels; ++lev) {
+        // Remove pressure gradient term
+        overset_ops::apply_pressure_gradient(vel(lev), rho(lev), gp(lev), -dt);
+        // Modify pressure gradient
+        overset_ops::replace_gradp(
+            gp(lev), (*m_gp_copy)(lev), iblank_cell(lev));
+        // Reapply pressure gradient term
+        overset_ops::apply_pressure_gradient(vel(lev), rho(lev), gp(lev), dt);
+    }
+}
+
 } // namespace amr_wind
\ No newline at end of file
diff --git a/amr-wind/overset/TiogaInterface.cpp b/amr-wind/overset/TiogaInterface.cpp
index 59e42d414d..ec1a75280f 100644
--- a/amr-wind/overset/TiogaInterface.cpp
+++ b/amr-wind/overset/TiogaInterface.cpp
@@ -127,6 +127,9 @@ void TiogaInterface::post_overset_conn_work()
     for (int lev = 0; lev < nlevels; ++lev) {
         htod_memcpy(m_iblank_cell(lev), (*m_iblank_cell_host)(lev), 0, 0, 1);
         htod_memcpy(m_iblank_node(lev), (*m_iblank_node_host)(lev), 0, 0, 1);
+
+        m_iblank_cell(lev).FillBoundary(m_sim.mesh().Geom()[lev].periodicity());
+        m_iblank_node(lev).FillBoundary(m_sim.mesh().Geom()[lev].periodicity());
     }
 
     iblank_to_mask(m_iblank_cell, m_mask_cell);
diff --git a/amr-wind/overset/overset_ops_K.H b/amr-wind/overset/overset_ops_K.H
new file mode 100644
index 0000000000..6b03c6a63e
--- /dev/null
+++ b/amr-wind/overset/overset_ops_K.H
@@ -0,0 +1,108 @@
+#ifndef OVERSET_OPS_K_H_
+#define OVERSET_OPS_K_H_
+
+#include <AMReX_FArrayBox.H>
+
+namespace amr_wind::overset_ops {
+
+// Approximate signed distance function
+amrex::Real AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+asdf(const amrex::Real a_vof, const amrex::Real i_th, const amrex::Real tiny)
+{
+    // function of local vof value and interface thickness
+    return (i_th * log((a_vof + tiny) / (1. - a_vof + tiny)));
+}
+
+amrex::Real AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE alpha_flux(
+    const int i,
+    const int j,
+    const int k,
+    const int dir,
+    const amrex::Real margin,
+    amrex::Array4<amrex::Real const> const& vof,
+    amrex::Array4<amrex::Real const> const& tg_vof,
+    amrex::Array4<amrex::Real const> const& normal)
+{
+    // Set up neighbor indices
+    const amrex::IntVect iv{i, j, k};
+    const amrex::IntVect dv{(int)(dir == 0), (int)(dir == 1), (int)(dir == 2)};
+    const amrex::IntVect ivm = iv - dv;
+
+    // Gradient of phi normal to interface
+    const amrex::Real gphi = (vof(iv) - vof(ivm));
+    // Normal vector in each cell (already normalized)
+    const amrex::Real norm_ = normal(iv, dir);
+    const amrex::Real norm_nb = normal(ivm, dir);
+
+    // Determine which delta_phi (and multiply by normal)
+    // The sign depends on side of flux face (like upwinding)
+    const amrex::Real dphi_ = (tg_vof(iv) - vof(iv)) * (-norm_);
+    const amrex::Real dphi_nb = (tg_vof(ivm) - vof(ivm)) * norm_nb;
+    // Average value used across the interface
+    amrex::Real dphi_eval = 0.5 * (dphi_ + dphi_nb);
+    // Upwinding when on the gas side, downwinding on the liquid
+    // Across the interface defined as crossing 0.5 or within margin of 0.5
+    if ((std::abs(vof(iv) - 0.5) > margin ||
+         std::abs(vof(ivm) - 0.5) > margin)) {
+        if (gphi > 0.0) {
+            dphi_eval = (vof(ivm) < 0.5 && vof(iv) <= 0.5 + margin) ? dphi_nb
+                                                                    : dphi_eval;
+            dphi_eval =
+                (vof(ivm) >= 0.5 - margin && vof(iv) > 0.5) ? dphi_ : dphi_eval;
+        }
+        if (gphi < 0.0) {
+            dphi_eval =
+                (vof(iv) < 0.5 && vof(ivm) <= 0.5 + margin) ? dphi_ : dphi_eval;
+            dphi_eval = (vof(iv) >= 0.5 - margin && vof(ivm) > 0.5) ? dphi_nb
+                                                                    : dphi_eval;
+        }
+    }
+    return dphi_eval;
+}
+
+void AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE velocity_face(
+    const int i,
+    const int j,
+    const int k,
+    const int dir,
+    amrex::Array4<amrex::Real const> const& vof,
+    amrex::Array4<amrex::Real const> const& velocity,
+    amrex::Real& uface,
+    amrex::Real& vface,
+    amrex::Real& wface)
+{
+    // Set up neighbor indices
+    const amrex::IntVect iv{i, j, k};
+    const amrex::IntVect dv{(int)(dir == 0), (int)(dir == 1), (int)(dir == 2)};
+    const amrex::IntVect ivm = iv - dv;
+
+    // Gradient of phi normal to interface
+    const amrex::Real gphi = (vof(iv) - vof(ivm));
+
+    // Get velocities on both sides
+    const amrex::Real u_ = velocity(iv, 0);
+    const amrex::Real v_ = velocity(iv, 1);
+    const amrex::Real w_ = velocity(iv, 2);
+    const amrex::Real u_nb = velocity(ivm, 0);
+    const amrex::Real v_nb = velocity(ivm, 1);
+    const amrex::Real w_nb = velocity(ivm, 2);
+    // Average value when gphi = 0
+    uface = 0.5 * (u_ + u_nb);
+    vface = 0.5 * (v_ + v_nb);
+    wface = 0.5 * (w_ + w_nb);
+    // Use simple upwinding
+    if (gphi > 0.0) {
+        uface = u_nb;
+        vface = v_nb;
+        wface = w_nb;
+    }
+    if (gphi < 0.0) {
+        uface = u_;
+        vface = v_;
+        wface = w_;
+    }
+}
+
+} // namespace amr_wind::overset_ops
+
+#endif
\ No newline at end of file
diff --git a/amr-wind/overset/overset_ops_routines.H b/amr-wind/overset/overset_ops_routines.H
new file mode 100644
index 0000000000..9e289f8ac2
--- /dev/null
+++ b/amr-wind/overset/overset_ops_routines.H
@@ -0,0 +1,453 @@
+#ifndef OVERSET_OPS_ROUTINES_H_
+#define OVERSET_OPS_ROUTINES_H_
+
+#include "amr-wind/equation_systems/vof/volume_fractions.H"
+#include "amr-wind/overset/overset_ops_K.H"
+
+namespace amr_wind::overset_ops {
+
+// Populate approximate signed distance function using vof field
+void populate_psi(
+    amrex::MultiFab& mf_psi,
+    const amrex::MultiFab& mf_vof,
+    const amrex::Real i_th,
+    const amrex::Real asdf_tiny)
+{
+    const auto& psi = mf_psi.arrays();
+    const auto& vof = mf_vof.const_arrays();
+    amrex::ParallelFor(
+        mf_psi, mf_psi.n_grow,
+        [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept {
+            psi[nbx](i, j, k) = asdf(vof[nbx](i, j, k), i_th, asdf_tiny);
+        });
+}
+
+// Modify a vof field to not have values that barely differ from 0 or 1
+void process_vof(amrex::MultiFab& mf_vof, const amrex::Real vof_tol)
+{
+    const auto& vof = mf_vof.arrays();
+    amrex::ParallelFor(
+        mf_vof, mf_vof.n_grow,
+        [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept {
+            vof[nbx](i, j, k) =
+                vof[nbx](i, j, k) < vof_tol
+                    ? 0.0
+                    : (vof[nbx](i, j, k) > 1. - vof_tol ? 1.
+                                                        : vof[nbx](i, j, k));
+        });
+}
+
+// Combine overset target vof field with current non-overset vof field
+void harmonize_vof(
+    amrex::MultiFab& mf_vof_target,
+    const amrex::MultiFab& mf_vof_original,
+    const amrex::iMultiFab& mf_iblank)
+{
+    const auto& tg_vof = mf_vof_target.arrays();
+    const auto& og_vof = mf_vof_original.const_arrays();
+    const auto& iblank = mf_iblank.const_arrays();
+    amrex::ParallelFor(
+        mf_vof_target,
+        [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept {
+            // Replace amr-wind vof values with originals
+            if (iblank[nbx](i, j, k) != -1) {
+                tg_vof[nbx](i, j, k) = og_vof[nbx](i, j, k);
+            }
+        });
+}
+
+// Populate normal vector with special treatment of overset boundary
+void populate_normal_vector(
+    amrex::MultiFab& mf_normvec,
+    const amrex::MultiFab& mf_vof,
+    const amrex::iMultiFab& mf_iblank)
+{
+    const auto& normvec = mf_normvec.arrays();
+    const auto& vof = mf_vof.const_arrays();
+    const auto& iblank = mf_iblank.const_arrays();
+    // Calculate gradients in each direction with centered diff
+    amrex::ParallelFor(
+        mf_normvec, mf_normvec.n_grow - amrex::IntVect(1),
+        [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept {
+            // Neumann condition across nalu bdy
+            int ibdy =
+                (iblank[nbx](i, j, k) != iblank[nbx](i - 1, j, k)) ? -1 : 0;
+            int jbdy =
+                (iblank[nbx](i, j, k) != iblank[nbx](i, j - 1, k)) ? -1 : 0;
+            int kbdy =
+                (iblank[nbx](i, j, k) != iblank[nbx](i, j, k - 1)) ? -1 : 0;
+            // no cell should be isolated such that -1 and 1 are needed
+            ibdy =
+                (iblank[nbx](i, j, k) != iblank[nbx](i + 1, j, k)) ? +1 : ibdy;
+            jbdy =
+                (iblank[nbx](i, j, k) != iblank[nbx](i, j + 1, k)) ? +1 : jbdy;
+            kbdy =
+                (iblank[nbx](i, j, k) != iblank[nbx](i, j, k + 1)) ? +1 : kbdy;
+            // Calculate normal
+            amrex::Real mx, my, mz, mmag;
+            multiphase::youngs_finite_difference_normal_neumann(
+                i, j, k, ibdy, jbdy, kbdy, vof[nbx], mx, my, mz);
+            // Normalize normal
+            mmag = std::sqrt(mx * mx + my * my + mz * mz + 1e-20);
+            // Save normal
+            normvec[nbx](i, j, k, 0) = mx / mmag;
+            normvec[nbx](i, j, k, 1) = my / mmag;
+            normvec[nbx](i, j, k, 2) = mz / mmag;
+        });
+}
+
+// Calculate fluxes for reinitialization over entire domain without concern for
+// overset bdy
+void populate_sharpen_fluxes(
+    amrex::MultiFab& mf_fx,
+    amrex::MultiFab& mf_fy,
+    amrex::MultiFab& mf_fz,
+    const amrex::MultiFab& mf_vof,
+    const amrex::MultiFab& mf_target_vof,
+    const amrex::MultiFab& mf_norm,
+    const amrex::MultiFab& mf_velocity,
+    const amrex::Real margin,
+    const amrex::Real rho1,
+    const amrex::Real rho2)
+{
+    const auto& fx = mf_fx.arrays();
+    const auto& fy = mf_fy.arrays();
+    const auto& fz = mf_fz.arrays();
+    const auto& vof = mf_vof.const_arrays();
+    const auto& tg_vof = mf_target_vof.const_arrays();
+    const auto& norm = mf_norm.const_arrays();
+    const auto& vel = mf_velocity.const_arrays();
+    amrex::ParallelFor(
+        mf_fx, mf_fx.n_grow,
+        [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept {
+            // vof flux
+            amrex::Real flux = alpha_flux(
+                i, j, k, 0, margin, vof[nbx], tg_vof[nbx], norm[nbx]);
+            fx[nbx](i, j, k, 0) = flux;
+            // density flux
+            flux *= (rho1 - rho2);
+            fx[nbx](i, j, k, 1) = flux;
+            // momentum fluxes (dens flux * face vel)
+            amrex::Real uf, vf, wf;
+            velocity_face(i, j, k, 0, vof[nbx], vel[nbx], uf, vf, wf);
+            fx[nbx](i, j, k, 2) = flux * uf;
+            fx[nbx](i, j, k, 3) = flux * vf;
+            fx[nbx](i, j, k, 4) = flux * wf;
+        });
+    amrex::ParallelFor(
+        mf_fy, mf_fy.n_grow,
+        [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept {
+            amrex::Real flux = alpha_flux(
+                i, j, k, 1, margin, vof[nbx], tg_vof[nbx], norm[nbx]);
+            fy[nbx](i, j, k, 0) = flux;
+            flux *= (rho1 - rho2);
+            fy[nbx](i, j, k, 1) = flux;
+            amrex::Real uf, vf, wf;
+            velocity_face(i, j, k, 1, vof[nbx], vel[nbx], uf, vf, wf);
+            fy[nbx](i, j, k, 2) = flux * uf;
+            fy[nbx](i, j, k, 3) = flux * vf;
+            fy[nbx](i, j, k, 4) = flux * wf;
+        });
+    amrex::ParallelFor(
+        mf_fz, mf_fz.n_grow,
+        [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept {
+            amrex::Real flux = alpha_flux(
+                i, j, k, 2, margin, vof[nbx], tg_vof[nbx], norm[nbx]);
+            fz[nbx](i, j, k, 0) = flux;
+            flux *= (rho1 - rho2);
+            fz[nbx](i, j, k, 1) = flux;
+            amrex::Real uf, vf, wf;
+            velocity_face(i, j, k, 2, vof[nbx], vel[nbx], uf, vf, wf);
+            fz[nbx](i, j, k, 2) = flux * uf;
+            fz[nbx](i, j, k, 3) = flux * vf;
+            fz[nbx](i, j, k, 4) = flux * wf;
+        });
+}
+
+// Process reinitialization fluxes - zero non-internal to overset region
+void process_fluxes(
+    amrex::MultiFab& mf_fx,
+    amrex::MultiFab& mf_fy,
+    amrex::MultiFab& mf_fz,
+    const amrex::iMultiFab& mf_iblank)
+{
+    const auto& fx = mf_fx.arrays();
+    const auto& fy = mf_fy.arrays();
+    const auto& fz = mf_fz.arrays();
+    const auto& iblank = mf_iblank.const_arrays();
+    amrex::ParallelFor(
+        mf_fx, mf_fx.n_grow, mf_fx.n_comp,
+        [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k, int n) noexcept {
+            bool zero_all =
+                (iblank[nbx](i - 1, j, k) + iblank[nbx](i, j, k) > -2);
+            fx[nbx](i, j, k, n) *= zero_all ? 0. : 1.;
+        });
+    amrex::ParallelFor(
+        mf_fy, mf_fy.n_grow, mf_fy.n_comp,
+        [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k, int n) noexcept {
+            bool zero_all =
+                (iblank[nbx](i, j - 1, k) + iblank[nbx](i, j, k) > -2);
+            fy[nbx](i, j, k, n) *= zero_all ? 0. : 1.;
+        });
+    amrex::ParallelFor(
+        mf_fz, mf_fz.n_grow, mf_fz.n_comp,
+        [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k, int n) noexcept {
+            bool zero_all =
+                (iblank[nbx](i, j, k - 1) + iblank[nbx](i, j, k) > -2);
+            fz[nbx](i, j, k, n) *= zero_all ? 0. : 1.;
+        });
+}
+
+// Calculate a type of CFL by measuring how much % VOF is being removed per cell
+amrex::Real calculate_pseudo_dt_flux(
+    amrex::MultiFab& mf_fx,
+    amrex::MultiFab& mf_fy,
+    amrex::MultiFab& mf_fz,
+    amrex::MultiFab& mf_vof,
+    amrex::Real tol)
+{
+    // Get the maximum flux magnitude, but just for vof fluxes
+    const amrex::Real pdt_fx = amrex::ReduceMin(
+        mf_fx, mf_vof, 0,
+        [=] AMREX_GPU_HOST_DEVICE(
+            amrex::Box const& bx, amrex::Array4<amrex::Real const> const& fx,
+            amrex::Array4<amrex::Real const> const& vof) -> amrex::Real {
+            amrex::Real pdt_fab = 1.0;
+            amrex::Loop(bx, [=, &pdt_fab](int i, int j, int k) noexcept {
+                amrex::Real pdt_lim = 1.0;
+                if (fx(i, j, k, 0) > tol && vof(i, j, k) > tol) {
+                    // VOF is removed from cell i
+                    pdt_lim = vof(i, j, k) / fx(i, j, k, 0);
+                } else if (fx(i, j, k, 0) < -tol && vof(i - 1, j, k) > tol) {
+                    // VOF is removed from cell i-1
+                    pdt_lim = vof(i - 1, j, k) / -fx(i, j, k, 0);
+                }
+                pdt_fab = amrex::min(pdt_fab, pdt_lim);
+            });
+            return pdt_fab;
+        });
+    const amrex::Real pdt_fy = amrex::ReduceMin(
+        mf_fy, mf_vof, 0,
+        [=] AMREX_GPU_HOST_DEVICE(
+            amrex::Box const& bx, amrex::Array4<amrex::Real const> const& fy,
+            amrex::Array4<amrex::Real const> const& vof) -> amrex::Real {
+            amrex::Real pdt_fab = 1.0;
+            amrex::Loop(bx, [=, &pdt_fab](int i, int j, int k) noexcept {
+                amrex::Real pdt_lim = 1.0;
+                if (fy(i, j, k, 0) > tol && vof(i, j, k) > tol) {
+                    // VOF is removed from cell j
+                    pdt_lim = vof(i, j, k) / fy(i, j, k, 0);
+                } else if (fy(i, j, k, 0) < -tol && vof(i, j - 1, k) > tol) {
+                    // VOF is removed from cell j-1
+                    pdt_lim = vof(i, j - 1, k) / -fy(i, j, k, 0);
+                }
+                pdt_fab = amrex::min(pdt_fab, pdt_lim);
+            });
+            return pdt_fab;
+        });
+    const amrex::Real pdt_fz = amrex::ReduceMin(
+        mf_fz, mf_vof, 0,
+        [=] AMREX_GPU_HOST_DEVICE(
+            amrex::Box const& bx, amrex::Array4<amrex::Real const> const& fz,
+            amrex::Array4<amrex::Real const> const& vof) -> amrex::Real {
+            amrex::Real pdt_fab = 1.0;
+            amrex::Loop(bx, [=, &pdt_fab](int i, int j, int k) noexcept {
+                amrex::Real pdt_lim = 1.0;
+                if (fz(i, j, k, 0) > tol && vof(i, j, k) > tol) {
+                    // VOF is removed from cell k
+                    pdt_lim = vof(i, j, k) / fz(i, j, k, 0);
+                } else if (fz(i, j, k, 0) < -tol && vof(i, j, k - 1) > tol) {
+                    // VOF is removed from cell k-1
+                    pdt_lim = vof(i, j, k - 1) / -fz(i, j, k, 0);
+                }
+                pdt_fab = amrex::min(pdt_fab, pdt_lim);
+            });
+            return pdt_fab;
+        });
+    const amrex::Real pdt = amrex::min(pdt_fx, amrex::min(pdt_fy, pdt_fz));
+    return pdt;
+}
+
+// Apply reinitialization fluxes to modify fields
+void apply_fluxes(
+    amrex::MultiFab& mf_fx,
+    amrex::MultiFab& mf_fy,
+    amrex::MultiFab& mf_fz,
+    amrex::MultiFab& mf_vof,
+    amrex::MultiFab& mf_dens,
+    amrex::MultiFab& mf_vel,
+    amrex::Real ptfac,
+    amrex::Real vof_tol)
+{
+    const auto& fx = mf_fx.const_arrays();
+    const auto& fy = mf_fy.const_arrays();
+    const auto& fz = mf_fz.const_arrays();
+    const auto& vof = mf_vof.arrays();
+    const auto& dens = mf_dens.arrays();
+    const auto& vel = mf_vel.arrays();
+
+    amrex::ParallelFor(
+        mf_vof, [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept {
+            const amrex::Real olddens = dens[nbx](i, j, k);
+            vof[nbx](i, j, k) +=
+                ptfac * ((fx[nbx](i + 1, j, k, 0) - fx[nbx](i, j, k, 0)) +
+                         (fy[nbx](i, j + 1, k, 0) - fy[nbx](i, j, k, 0)) +
+                         (fz[nbx](i, j, k + 1, 0) - fz[nbx](i, j, k, 0)));
+            dens[nbx](i, j, k) +=
+                ptfac * ((fx[nbx](i + 1, j, k, 1) - fx[nbx](i, j, k, 1)) +
+                         (fy[nbx](i, j + 1, k, 1) - fy[nbx](i, j, k, 1)) +
+                         (fz[nbx](i, j, k + 1, 1) - fz[nbx](i, j, k, 1)));
+            vel[nbx](i, j, k, 0) =
+                1.0 / dens[nbx](i, j, k) *
+                (olddens * vel[nbx](i, j, k, 0) +
+                 ptfac * ((fx[nbx](i + 1, j, k, 2) - fx[nbx](i, j, k, 2)) +
+                          (fy[nbx](i, j + 1, k, 2) - fy[nbx](i, j, k, 2)) +
+                          (fz[nbx](i, j, k + 1, 2) - fz[nbx](i, j, k, 2))));
+            vel[nbx](i, j, k, 1) =
+                1.0 / dens[nbx](i, j, k) *
+                (olddens * vel[nbx](i, j, k, 1) +
+                 ptfac * ((fx[nbx](i + 1, j, k, 3) - fx[nbx](i, j, k, 3)) +
+                          (fy[nbx](i, j + 1, k, 3) - fy[nbx](i, j, k, 3)) +
+                          (fz[nbx](i, j, k + 1, 3) - fz[nbx](i, j, k, 3))));
+            vel[nbx](i, j, k, 2) =
+                1.0 / dens[nbx](i, j, k) *
+                (olddens * vel[nbx](i, j, k, 2) +
+                 ptfac * ((fx[nbx](i + 1, j, k, 4) - fx[nbx](i, j, k, 4)) +
+                          (fy[nbx](i, j + 1, k, 4) - fy[nbx](i, j, k, 4)) +
+                          (fz[nbx](i, j, k + 1, 4) - fz[nbx](i, j, k, 4))));
+
+            // Ensure vof is bounded
+            vof[nbx](i, j, k) =
+                vof[nbx](i, j, k) < vof_tol
+                    ? 0.0
+                    : (vof[nbx](i, j, k) > 1. - vof_tol ? 1.
+                                                        : vof[nbx](i, j, k));
+            // Density bounds are enforced elsewhere
+        });
+}
+
+// Get the size of the smallest VOF flux to quantify convergence
+amrex::Real measure_convergence(
+    amrex::MultiFab& mf_fx, amrex::MultiFab& mf_fy, amrex::MultiFab& mf_fz)
+{
+    // Get the maximum flux magnitude, but just for vof fluxes
+    const amrex::Real err_fx = amrex::ReduceMax(
+        mf_fx, 0,
+        [=] AMREX_GPU_HOST_DEVICE(
+            amrex::Box const& bx,
+            amrex::Array4<amrex::Real const> const& fx) -> amrex::Real {
+            amrex::Real err_fab = -1.0;
+            amrex::Loop(bx, [=, &err_fab](int i, int j, int k) noexcept {
+                err_fab = amrex::max(err_fab, std::abs(fx(i, j, k, 0)));
+            });
+            return err_fab;
+        });
+    const amrex::Real err_fy = amrex::ReduceMax(
+        mf_fy, 0,
+        [=] AMREX_GPU_HOST_DEVICE(
+            amrex::Box const& bx,
+            amrex::Array4<amrex::Real const> const& fy) -> amrex::Real {
+            amrex::Real err_fab = -1.0;
+            amrex::Loop(bx, [=, &err_fab](int i, int j, int k) noexcept {
+                err_fab = amrex::max(err_fab, std::abs(fy(i, j, k, 0)));
+            });
+            return err_fab;
+        });
+    const amrex::Real err_fz = amrex::ReduceMax(
+        mf_fz, 0,
+        [=] AMREX_GPU_HOST_DEVICE(
+            amrex::Box const& bx,
+            amrex::Array4<amrex::Real const> const& fz) -> amrex::Real {
+            amrex::Real err_fab = -1.0;
+            amrex::Loop(bx, [=, &err_fab](int i, int j, int k) noexcept {
+                err_fab = amrex::max(err_fab, std::abs(fz(i, j, k, 0)));
+            });
+            return err_fab;
+        });
+    const amrex::Real err = amrex::max(err_fx, amrex::max(err_fy, err_fz));
+    return err;
+}
+
+// Set levelset field to another quantity to view in plotfile for debugging
+void equate_field(amrex::MultiFab& mf_dest, const amrex::MultiFab& mf_src)
+{
+    const auto& dest = mf_dest.arrays();
+    const auto& src = mf_src.const_arrays();
+    amrex::ParallelFor(
+        mf_dest, [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept {
+            dest[nbx](i, j, k) = std::sqrt(
+                src[nbx](i, j, k, 0) * src[nbx](i, j, k, 0) +
+                src[nbx](i, j, k, 1) * src[nbx](i, j, k, 1) +
+                src[nbx](i, j, k, 2) * src[nbx](i, j, k, 2));
+        });
+}
+
+// Replace pressure gradient with hydrostatic field in overset regions
+void replace_gradp_hydrostatic(
+    amrex::MultiFab& mf_gp,
+    const amrex::MultiFab& mf_density,
+    const amrex::MultiFab& mf_refdens,
+    const amrex::iMultiFab& mf_iblank,
+    const amrex::Real grav_z,
+    const bool is_pptb)
+{
+    const auto& gp = mf_gp.arrays();
+    const auto& rho = mf_density.const_arrays();
+    const auto& rho0 = mf_refdens.const_arrays();
+    const auto& iblank = mf_iblank.const_arrays();
+    amrex::ParallelFor(
+        mf_gp, mf_gp.n_grow,
+        [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept {
+            if (iblank[nbx](i, j, k) == -1) {
+                const amrex::Real dfac =
+                    is_pptb ? rho[nbx](i, j, k) - rho0[nbx](i, j, k)
+                            : rho[nbx](i, j, k);
+                gp[nbx](i, j, k, 0) = 0.;
+                gp[nbx](i, j, k, 1) = 0.;
+                gp[nbx](i, j, k, 2) = dfac * grav_z;
+            }
+        });
+}
+
+// Swap pressure gradient values in overset region
+void replace_gradp(
+    amrex::MultiFab& mf_gp,
+    const amrex::MultiFab& mf_gp0,
+    const amrex::iMultiFab& mf_iblank)
+{
+    const auto gp = mf_gp.arrays();
+    const auto gp0 = mf_gp0.const_arrays();
+    const auto iblank = mf_iblank.const_arrays();
+    amrex::ParallelFor(
+        mf_gp, mf_gp.n_grow,
+        [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept {
+            if (iblank[nbx](i, j, k) == -1) {
+                gp[nbx](i, j, k, 0) = gp0[nbx](i, j, k, 0);
+                gp[nbx](i, j, k, 1) = gp0[nbx](i, j, k, 1);
+                gp[nbx](i, j, k, 2) = gp0[nbx](i, j, k, 2);
+            }
+        });
+}
+
+// Apply pressure gradient to velocity field
+void apply_pressure_gradient(
+    amrex::MultiFab& mf_vel,
+    const amrex::MultiFab& mf_density,
+    const amrex::MultiFab& mf_gp,
+    const amrex::Real scaling_factor)
+{
+    const auto& vel = mf_vel.arrays();
+    const auto& rho = mf_density.const_arrays();
+    const auto& gp = mf_gp.const_arrays();
+    amrex::ParallelFor(
+        mf_vel, [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept {
+            const amrex::Real soverrho = scaling_factor / rho[nbx](i, j, k);
+            vel[nbx](i, j, k, 0) -= gp[nbx](i, j, k, 0) * soverrho;
+            vel[nbx](i, j, k, 1) -= gp[nbx](i, j, k, 1) * soverrho;
+            vel[nbx](i, j, k, 2) -= gp[nbx](i, j, k, 2) * soverrho;
+        });
+}
+
+} // namespace amr_wind::overset_ops
+
+#endif
\ No newline at end of file
diff --git a/amr-wind/physics/multiphase/MultiPhase.H b/amr-wind/physics/multiphase/MultiPhase.H
index a900acda88..8b43eef505 100644
--- a/amr-wind/physics/multiphase/MultiPhase.H
+++ b/amr-wind/physics/multiphase/MultiPhase.H
@@ -3,6 +3,8 @@
 
 #include "amr-wind/core/Physics.H"
 #include "amr-wind/core/Field.H"
+#include "amr-wind/core/IntField.H"
+#include "amr-wind/core/ScratchField.H"
 
 /** Multiphase physics
  *
@@ -51,6 +53,8 @@ public:
 
     void levelset2vof();
 
+    void levelset2vof(const IntField& iblank_cell, ScratchField& vof_scr);
+
     void favre_filtering();
 
     amrex::Real volume_fraction_sum();
diff --git a/amr-wind/physics/multiphase/MultiPhase.cpp b/amr-wind/physics/multiphase/MultiPhase.cpp
index 33bb7140b0..efc390ddfc 100644
--- a/amr-wind/physics/multiphase/MultiPhase.cpp
+++ b/amr-wind/physics/multiphase/MultiPhase.cpp
@@ -488,7 +488,8 @@ void MultiPhase::levelset2vof()
             amrex::ParallelFor(
                 vbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
                     amrex::Real mx, my, mz;
-                    multiphase::youngs_fd_normal(i, j, k, phi, mx, my, mz);
+                    multiphase::youngs_finite_difference_normal(
+                        i, j, k, phi, mx, my, mz);
                     mx = std::abs(mx / 32.);
                     my = std::abs(my / 32.);
                     mz = std::abs(mz / 32.);
@@ -498,13 +499,73 @@ void MultiPhase::levelset2vof()
                     mz = mz / normL1;
                     // Make sure that alpha is negative far away from the
                     // interface
-                    amrex::Real alpha;
-                    if (phi(i, j, k) < -eps) {
-                        alpha = -1.0;
+                    const amrex::Real alpha = (phi(i, j, k) < -eps)
+                                                  ? -1.0
+                                                  : phi(i, j, k) / normL1 + 0.5;
+                    if (alpha >= 1.0) {
+                        volfrac(i, j, k) = 1.0;
+                    } else if (alpha <= 0.0) {
+                        volfrac(i, j, k) = 0.0;
                     } else {
-                        alpha = phi(i, j, k) / normL1;
-                        alpha = alpha + 0.5;
+                        volfrac(i, j, k) =
+                            multiphase::cut_volume(mx, my, mz, alpha, 0.0, 1.0);
                     }
+                });
+        }
+    }
+    // Fill ghost and boundary cells before simulation begins
+    (*m_vof).fillpatch(0.0);
+}
+
+// Do levelset2vof with iblank neumann and into supplied scratch field
+void MultiPhase::levelset2vof(
+    const IntField& iblank_cell, ScratchField& vof_scr)
+{
+    const int nlevels = m_sim.repo().num_active_levels();
+    (*m_levelset).fillpatch(m_sim.time().current_time());
+    const auto& geom = m_sim.mesh().Geom();
+
+    for (int lev = 0; lev < nlevels; ++lev) {
+        auto& levelset = (*m_levelset)(lev);
+        auto& vof = vof_scr(lev);
+        const auto& dx = geom[lev].CellSizeArray();
+
+        for (amrex::MFIter mfi(levelset); mfi.isValid(); ++mfi) {
+            const auto& vbx = mfi.validbox();
+            const amrex::Array4<amrex::Real>& phi = levelset.array(mfi);
+            const amrex::Array4<amrex::Real>& volfrac = vof.array(mfi);
+            const amrex::Array4<const int>& iblank =
+                iblank_cell(lev).const_array(mfi);
+            const amrex::Real eps = 2. * std::cbrt(dx[0] * dx[1] * dx[2]);
+            amrex::ParallelFor(
+                vbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+                    // Neumann of levelset across iblank boundaries
+                    int ibdy =
+                        (iblank(i, j, k) != iblank(i - 1, j, k)) ? -1 : 0;
+                    int jbdy =
+                        (iblank(i, j, k) != iblank(i, j - 1, k)) ? -1 : 0;
+                    int kbdy =
+                        (iblank(i, j, k) != iblank(i, j, k - 1)) ? -1 : 0;
+                    // no cell should be isolated such that -1 and 1 are
+                    // needed
+                    ibdy = (iblank(i, j, k) != iblank(i + 1, j, k)) ? +1 : ibdy;
+                    jbdy = (iblank(i, j, k) != iblank(i, j + 1, k)) ? +1 : jbdy;
+                    kbdy = (iblank(i, j, k) != iblank(i, j, k + 1)) ? +1 : kbdy;
+                    amrex::Real mx, my, mz;
+                    multiphase::youngs_finite_difference_normal_neumann(
+                        i, j, k, ibdy, jbdy, kbdy, phi, mx, my, mz);
+                    mx = std::abs(mx / 32.);
+                    my = std::abs(my / 32.);
+                    mz = std::abs(mz / 32.);
+                    amrex::Real normL1 = (mx + my + mz);
+                    mx = mx / normL1;
+                    my = my / normL1;
+                    mz = mz / normL1;
+                    // Make sure that alpha is negative far away from the
+                    // interface
+                    const amrex::Real alpha = (phi(i, j, k) < -eps)
+                                                  ? -1.0
+                                                  : phi(i, j, k) / normL1 + 0.5;
                     if (alpha >= 1.0) {
                         volfrac(i, j, k) = 1.0;
                     } else if (alpha <= 0.0) {
@@ -517,7 +578,7 @@ void MultiPhase::levelset2vof()
         }
     }
     // Fill ghost and boundary cells before simulation begins
-    (*m_vof).fillpatch(0.0);
+    vof_scr.fillpatch(0.0);
 }
 
 } // namespace amr_wind
diff --git a/amr-wind/physics/multiphase/SloshingTank.H b/amr-wind/physics/multiphase/SloshingTank.H
index bf7b24b763..4d467ac9bc 100644
--- a/amr-wind/physics/multiphase/SloshingTank.H
+++ b/amr-wind/physics/multiphase/SloshingTank.H
@@ -35,6 +35,7 @@ public:
 private:
     Field& m_velocity;
     Field& m_levelset;
+    Field& m_pressure;
 
     //! Initial free surface amplitude magnitude
     amrex::Real m_amplitude{0.1};
@@ -44,6 +45,14 @@ private:
 
     //! Initial zero-level free-surface water depth
     amrex::Real m_waterlevel{0.0};
+
+    //! Stuff to get from MultiPhase physics
+    amrex::Real m_rho1{std::numeric_limits<amrex::Real>::quiet_NaN()};
+    amrex::Real m_rho2{std::numeric_limits<amrex::Real>::quiet_NaN()};
+    amrex::Vector<amrex::Real> m_gravity{
+        std::numeric_limits<amrex::Real>::quiet_NaN(),
+        std::numeric_limits<amrex::Real>::quiet_NaN(),
+        std::numeric_limits<amrex::Real>::quiet_NaN()};
 };
 
 } // namespace amr_wind
diff --git a/amr-wind/physics/multiphase/SloshingTank.cpp b/amr-wind/physics/multiphase/SloshingTank.cpp
index 7558141c7c..167911812f 100644
--- a/amr-wind/physics/multiphase/SloshingTank.cpp
+++ b/amr-wind/physics/multiphase/SloshingTank.cpp
@@ -1,4 +1,5 @@
 #include "amr-wind/physics/multiphase/SloshingTank.H"
+#include "amr-wind/physics/multiphase/MultiPhase.H"
 #include "amr-wind/CFDSim.H"
 #include "AMReX_ParmParse.H"
 #include "amr-wind/fvm/gradient.H"
@@ -9,16 +10,22 @@ namespace amr_wind {
 SloshingTank::SloshingTank(CFDSim& sim)
     : m_velocity(sim.repo().get_field("velocity"))
     , m_levelset(sim.repo().get_field("levelset"))
+    , m_pressure(sim.repo().get_field("p"))
 {
     amrex::ParmParse pp(identifier());
     pp.query("amplitude", m_amplitude);
     pp.query("peak_enhance", m_kappa);
     pp.query("water_level", m_waterlevel);
+
+    const auto& mphase = sim.physics_manager().get<MultiPhase>();
+    m_rho1 = mphase.rho1();
+    m_rho2 = mphase.rho2();
+    m_gravity = mphase.gravity();
 }
 
 /** Initialize the velocity and levelset fields at the beginning of the
- *  simulation.
- *
+ *  simulation. Initializing pressure has little effect for pure AMR-Wind
+ *  simulations, but it improves initialization in the overset solver.
  */
 void SloshingTank::initialize_fields(int level, const amrex::Geometry& geom)
 {
@@ -26,6 +33,7 @@ void SloshingTank::initialize_fields(int level, const amrex::Geometry& geom)
     velocity.setVal(0.0);
 
     auto& levelset = m_levelset(level);
+    auto& pressure = m_pressure(level);
     const auto& dx = geom.CellSizeArray();
     const auto& problo = geom.ProbLoArray();
     const auto& probhi = geom.ProbHiArray();
@@ -34,10 +42,14 @@ void SloshingTank::initialize_fields(int level, const amrex::Geometry& geom)
     const amrex::Real water_level = m_waterlevel;
     const amrex::Real Lx = probhi[0] - problo[0];
     const amrex::Real Ly = probhi[1] - problo[1];
+    const amrex::Real rho1 = m_rho1;
+    const amrex::Real rho2 = m_rho2;
+    const amrex::Real grav_z = m_gravity[2];
 
     for (amrex::MFIter mfi(levelset); mfi.isValid(); ++mfi) {
         const auto& vbx = mfi.validbox();
         auto phi = levelset.array(mfi);
+        auto p = pressure.array(mfi);
 
         amrex::ParallelFor(
             vbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
@@ -51,6 +63,29 @@ void SloshingTank::initialize_fields(int level, const amrex::Geometry& geom)
                                         std::pow(y - problo[1] - 0.5 * Ly, 2)));
                 phi(i, j, k) = z0 - z;
             });
+        amrex::Box const& nbx = mfi.grownnodaltilebox();
+        amrex::ParallelFor(
+            nbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+                // For pressure nodes, no offset
+                const amrex::Real x = problo[0] + i * dx[0];
+                const amrex::Real y = problo[1] + j * dx[1];
+                const amrex::Real z = problo[2] + k * dx[2];
+                const amrex::Real z0 =
+                    water_level +
+                    Amp * std::exp(
+                              -kappa * (std::pow(x - problo[0] - 0.5 * Lx, 2) +
+                                        std::pow(y - problo[1] - 0.5 * Ly, 2)));
+                // Integrated (top-down in z) phase heights to pressure node
+                const amrex::Real ih_g =
+                    amrex::max(0.0, amrex::min(probhi[2] - z0, probhi[2] - z));
+                const amrex::Real ih_l =
+                    amrex::max(0.0, amrex::min(z0 - z, z0 - problo[2]));
+                // Integrated rho at pressure node
+                const amrex::Real irho = rho1 * ih_l + rho2 * ih_g;
+
+                // Add term to reference pressure
+                p(i, j, k) = -irho * grav_z;
+            });
     }
 }
 
diff --git a/amr-wind/physics/multiphase/hydrostatic_ops.H b/amr-wind/physics/multiphase/hydrostatic_ops.H
index 4f4eeb1942..0cea6e1af8 100644
--- a/amr-wind/physics/multiphase/hydrostatic_ops.H
+++ b/amr-wind/physics/multiphase/hydrostatic_ops.H
@@ -46,15 +46,12 @@ static void define_p0(
             auto p0_arr = p0(lev).array(mfi);
             amrex::ParallelFor(
                 nbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
-                    // Height of pressure node
-                    const amrex::Real hnode = k * dx[2];
-                    // Liquid height
-                    const amrex::Real hliq = wlev - problo[2];
+                    const amrex::Real z = problo[2] + k * dx[2];
                     // Integrated (top-down in z) phase heights to pressure node
                     amrex::Real ih_g = amrex::max(
-                        0.0, amrex::min(probhi[2] - hliq, probhi[2] - hnode));
-                    amrex::Real ih_l = amrex::max(
-                        0.0, amrex::min(hliq - hnode, hliq - problo[2]));
+                        0.0, amrex::min(probhi[2] - wlev, probhi[2] - z));
+                    amrex::Real ih_l =
+                        amrex::max(0.0, amrex::min(wlev - z, wlev - problo[2]));
                     // Integrated rho at pressure node
                     const amrex::Real irho = rho1 * ih_l + rho2 * ih_g;
 
diff --git a/amr-wind/projection/incflo_apply_nodal_projection.cpp b/amr-wind/projection/incflo_apply_nodal_projection.cpp
index fe1cac1bfc..4592caaf64 100644
--- a/amr-wind/projection/incflo_apply_nodal_projection.cpp
+++ b/amr-wind/projection/incflo_apply_nodal_projection.cpp
@@ -358,8 +358,15 @@ void incflo::ApplyProjection(
         nodal_projector->setCustomRHS(div_vel_rhs->vec_const_ptrs());
     }
 
-    // Setup masking for overset simulations
+    // Determine if nodal projection should be coupled for overset
+    bool disable_ovst_nodal = false;
     if (sim().has_overset()) {
+        amrex::ParmParse pp("Overset");
+        pp.query("disable_coupled_nodal_proj", disable_ovst_nodal);
+    }
+
+    // Setup masking for overset simulations
+    if (sim().has_overset() && !disable_ovst_nodal) {
         auto& linop = nodal_projector->getLinOp();
         const auto& imask_node = repo().get_int_field("mask_node");
         for (int lev = 0; lev <= finest_level; ++lev) {
@@ -367,7 +374,7 @@ void incflo::ApplyProjection(
         }
     }
 
-    if (m_sim.has_overset()) {
+    if (m_sim.has_overset() && !disable_ovst_nodal) {
         auto phif = m_repo.create_scratch_field(1, 1, amr_wind::FieldLoc::NODE);
         if (incremental) {
             for (int lev = 0; lev <= finestLevel(); ++lev) {
diff --git a/unit_tests/aw_test_utils/iter_tools.H b/unit_tests/aw_test_utils/iter_tools.H
index f8a67ea13d..908937d1c2 100644
--- a/unit_tests/aw_test_utils/iter_tools.H
+++ b/unit_tests/aw_test_utils/iter_tools.H
@@ -20,6 +20,25 @@ namespace amr_wind_tests {
  *  @param field MultiFabs for each of the `nlevels` levels
  *  @param func Functor to execute over the MultiFabs
  */
+
+template <typename Functor>
+void run_algorithm_single_lev(
+    const int lev, amrex::MultiFab& lfab, const Functor& func)
+{
+    for (amrex::MFIter mfi(lfab); mfi.isValid(); ++mfi) {
+        func(lev, mfi);
+    }
+}
+
+template <typename Functor>
+void run_algorithm_single_lev(
+    const int lev, amrex::iMultiFab& lfab, const Functor& func)
+{
+    for (amrex::MFIter mfi(lfab); mfi.isValid(); ++mfi) {
+        func(lev, mfi);
+    }
+}
+
 template <typename Functor>
 void run_algorithm(
     const int nlevels,
@@ -28,10 +47,7 @@ void run_algorithm(
 {
     for (int lev = 0; lev < nlevels; ++lev) {
         auto& lfab = *field[lev];
-
-        for (amrex::MFIter mfi(lfab); mfi.isValid(); ++mfi) {
-            func(lev, mfi);
-        }
+        run_algorithm_single_lev(lev, lfab, func);
     }
 }
 
@@ -41,10 +57,29 @@ void run_algorithm(amr_wind::Field& field, const Functor& func)
     const int nlevels = field.repo().num_active_levels();
     for (int lev = 0; lev < nlevels; ++lev) {
         auto& lfab = field(lev);
+        run_algorithm_single_lev(lev, lfab, func);
+    }
+}
 
-        for (amrex::MFIter mfi(lfab); mfi.isValid(); ++mfi) {
-            func(lev, mfi);
-        }
+template <typename Functor>
+void run_algorithm(
+    const int nlevels,
+    amrex::Vector<amrex::iMultiFab*>& intfield,
+    const Functor& func)
+{
+    for (int lev = 0; lev < nlevels; ++lev) {
+        auto& lfab = *intfield[lev];
+        run_algorithm_single_lev(lev, lfab, func);
+    }
+}
+
+template <typename Functor>
+void run_algorithm(amr_wind::IntField& intfield, const Functor& func)
+{
+    const int nlevels = intfield.repo().num_active_levels();
+    for (int lev = 0; lev < nlevels; ++lev) {
+        auto& lfab = intfield(lev);
+        run_algorithm_single_lev(lev, lfab, func);
     }
 }
 
diff --git a/unit_tests/multiphase/CMakeLists.txt b/unit_tests/multiphase/CMakeLists.txt
index 7bc2bc6de7..8386a857d8 100644
--- a/unit_tests/multiphase/CMakeLists.txt
+++ b/unit_tests/multiphase/CMakeLists.txt
@@ -7,4 +7,5 @@ target_sources(
   test_vof_BCs.cpp
   test_mflux_schemes.cpp
   test_reference_fields.cpp
+  test_vof_overset_ops.cpp
   )
diff --git a/unit_tests/multiphase/test_reference_fields.cpp b/unit_tests/multiphase/test_reference_fields.cpp
index 257c024a6e..df80ff4976 100644
--- a/unit_tests/multiphase/test_reference_fields.cpp
+++ b/unit_tests/multiphase/test_reference_fields.cpp
@@ -108,7 +108,7 @@ class MultiPhaseHydroStatic : public MeshTest
         }
         {
             amrex::ParmParse pp("geometry");
-            amrex::Vector<amrex::Real> problo{{0.0, 0.0, 0.0}};
+            amrex::Vector<amrex::Real> problo{{-1.0, -1.0, -1.0}};
             amrex::Vector<amrex::Real> probhi{{1.0, 1.0, 1.0}};
 
             pp.addarr("prob_lo", problo);
@@ -120,7 +120,7 @@ class MultiPhaseHydroStatic : public MeshTest
     const amrex::Real m_rho2 = 1.0;
     const amrex::Real m_wlev = 0.5;
     const amrex::Real m_gz = -9.81;
-    const int m_nx = 3;
+    const int m_nx = 16;
 };
 
 TEST_F(MultiPhaseHydroStatic, reference_density)
diff --git a/unit_tests/multiphase/test_vof_overset_ops.cpp b/unit_tests/multiphase/test_vof_overset_ops.cpp
new file mode 100644
index 0000000000..82b80083fe
--- /dev/null
+++ b/unit_tests/multiphase/test_vof_overset_ops.cpp
@@ -0,0 +1,406 @@
+#include "aw_test_utils/MeshTest.H"
+#include "aw_test_utils/iter_tools.H"
+#include "aw_test_utils/test_utils.H"
+#include "amr-wind/overset/overset_ops_K.H"
+
+namespace amr_wind_tests {
+
+class VOFOversetOps : public MeshTest
+{
+protected:
+    void populate_parameters() override
+    {
+        MeshTest::populate_parameters();
+
+        {
+            amrex::ParmParse pp("amr");
+            amrex::Vector<int> ncell{{8, 8, 8}};
+            pp.add("max_level", 0);
+            pp.add("max_grid_size", 4);
+            pp.addarr("n_cell", ncell);
+        }
+        {
+            amrex::ParmParse pp("geometry");
+            amrex::Vector<amrex::Real> problo{{0.0, 0.0, 0.0}};
+            amrex::Vector<amrex::Real> probhi{{1.0, 1.0, 1.0}};
+
+            pp.addarr("prob_lo", problo);
+            pp.addarr("prob_hi", probhi);
+        }
+    }
+};
+
+namespace {
+void init_vof_etc(
+    amr_wind::Field& vof,
+    amr_wind::Field& tg_vof,
+    amr_wind::Field& norm,
+    const int& dir)
+{
+    run_algorithm(vof, [&](const int lev, const amrex::MFIter& mfi) {
+        auto vof_arr = vof(lev).array(mfi);
+        // tgvof and norm values are set randomly to make sure they are involved
+        auto tgvof_arr = tg_vof(lev).array(mfi);
+        auto norm_arr = norm(lev).array(mfi);
+        const auto& bx = mfi.validbox();
+        amrex::ParallelFor(
+            grow(bx, 2), [=] AMREX_GPU_DEVICE(int i, int j, int k) {
+                vof_arr(i, j, k) = 0.0;
+                tgvof_arr(i, j, k) = 0.;
+                norm_arr(i, j, k, dir) = 1.0;
+                const int idx = (dir == 0 ? i : (dir == 1 ? j : k));
+                if (idx == -1) {
+                    // Within margin
+                    vof_arr(i, j, k) = 0.41;
+                    tgvof_arr(i, j, k) = 0.5;
+                    norm_arr(i, j, k, dir) = -1.0;
+                } else if (idx == 0) {
+                    // Within margin
+                    vof_arr(i, j, k) = 0.55;
+                    tgvof_arr(i, j, k) = 0.5;
+                    norm_arr(i, j, k, dir) = 1.0;
+                } else if (idx == 1) {
+                    // Outside margin, low
+                    vof_arr(i, j, k) = 0.1;
+                    tgvof_arr(i, j, k) = 0.2;
+                    norm_arr(i, j, k, dir) = -1.0;
+                } else if (idx == 2) {
+                    // Also low
+                    vof_arr(i, j, k) = 0.2;
+                    tgvof_arr(i, j, k) = 0.22;
+                    norm_arr(i, j, k, dir) = 1.0;
+                } else if (idx == 3) {
+                    // Above half
+                    vof_arr(i, j, k) = 0.7;
+                    tgvof_arr(i, j, k) = 0.8;
+                    norm_arr(i, j, k, dir) = -1.0;
+                } else if (idx == 4) {
+                    // Also high
+                    vof_arr(i, j, k) = 0.65;
+                    tgvof_arr(i, j, k) = 0.91;
+                    norm_arr(i, j, k, dir) = 1.0;
+                } else if (idx == 5) {
+                    // Also high, positive gradient
+                    vof_arr(i, j, k) = 0.8;
+                    tgvof_arr(i, j, k) = 0.75;
+                    norm_arr(i, j, k, dir) = -1.0;
+                } else if (idx == 6) {
+                    // Within margin
+                    vof_arr(i, j, k) = 0.45;
+                    tgvof_arr(i, j, k) = 0.5;
+                    norm_arr(i, j, k, dir) = 1.0;
+                } else if (idx == 7) {
+                    // Low, negative gradient
+                    vof_arr(i, j, k) = 0.3;
+                    tgvof_arr(i, j, k) = 0.2;
+                    norm_arr(i, j, k, dir) = -1.0;
+                } else if (idx == 8) {
+                    // High
+                    vof_arr(i, j, k) = 0.7;
+                    tgvof_arr(i, j, k) = 0.6;
+                    norm_arr(i, j, k, dir) = 1.0;
+                }
+            });
+    });
+}
+
+void init_velocity_etc(
+    amr_wind::Field& velocity, amr_wind::Field& vof, const int& dir)
+{
+    run_algorithm(velocity, [&](const int lev, const amrex::MFIter& mfi) {
+        auto vel_arr = velocity(lev).array(mfi);
+        auto vof_arr = vof(lev).array(mfi);
+        const auto& bx = mfi.validbox();
+        amrex::ParallelFor(
+            grow(bx, 2), [=] AMREX_GPU_DEVICE(int i, int j, int k) {
+                vel_arr(i, j, k) = 0.0;
+                vof_arr(i, j, k) = 0.0;
+                const int idx = (dir == 0 ? i : (dir == 1 ? j : k));
+                if (idx == -1) {
+                    vof_arr(i, j, k) = 0.41;
+                    vel_arr(i, j, k, 0) = 1.0;
+                    vel_arr(i, j, k, 1) = 2.0;
+                    vel_arr(i, j, k, 2) = 3.0;
+                } else if (idx == 0) {
+                    vof_arr(i, j, k) = 0.55;
+                    vel_arr(i, j, k, 0) = 1.5;
+                    vel_arr(i, j, k, 1) = 2.5;
+                    vel_arr(i, j, k, 2) = 3.5;
+                } else if (idx == 1) {
+                    vof_arr(i, j, k) = 0.2;
+                    vel_arr(i, j, k, 0) = 2.0;
+                    vel_arr(i, j, k, 1) = 3.0;
+                    vel_arr(i, j, k, 2) = 4.0;
+                } else if (idx == 2) {
+                    vof_arr(i, j, k) = 0.2;
+                    vel_arr(i, j, k, 0) = 2.5;
+                    vel_arr(i, j, k, 1) = 3.5;
+                    vel_arr(i, j, k, 2) = 4.5;
+                }
+            });
+    });
+}
+
+void calc_alpha_flux(
+    amr_wind::Field& flux,
+    amr_wind::Field& vof,
+    amr_wind::Field& tg_vof,
+    amr_wind::Field& norm,
+    const int& dir,
+    const amrex::Real& margin)
+{
+    run_algorithm(flux, [&](const int lev, const amrex::MFIter& mfi) {
+        auto f_arr = flux(lev).array(mfi);
+        auto vof_arr = vof(lev).array(mfi);
+        auto tgvof_arr = tg_vof(lev).array(mfi);
+        auto norm_arr = norm(lev).array(mfi);
+        const auto& bx = mfi.validbox();
+        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
+            f_arr(i, j, k) = amr_wind::overset_ops::alpha_flux(
+                i, j, k, dir, margin, vof_arr, tgvof_arr, norm_arr);
+        });
+    });
+}
+
+void calc_velocity_face(
+    amr_wind::Field& flux,
+    amr_wind::Field& velocity,
+    amr_wind::Field& vof,
+    const int& dir)
+{
+    run_algorithm(flux, [&](const int lev, const amrex::MFIter& mfi) {
+        auto f_arr = flux(lev).array(mfi);
+        auto vof_arr = vof(lev).array(mfi);
+        auto vel_arr = velocity(lev).array(mfi);
+        const auto& bx = mfi.validbox();
+        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
+            amrex::Real u_f, v_f, w_f;
+            amr_wind::overset_ops::velocity_face(
+                i, j, k, dir, vof_arr, vel_arr, u_f, v_f, w_f);
+            f_arr(i, j, k, 0) = u_f;
+            f_arr(i, j, k, 1) = v_f;
+            f_arr(i, j, k, 2) = w_f;
+        });
+    });
+}
+
+amrex::Real check_alpha_flux_impl(amr_wind::Field& flux, const int& dir)
+{
+    amrex::Real error_total = 0;
+
+    for (int lev = 0; lev < flux.repo().num_active_levels(); ++lev) {
+
+        error_total += amrex::ReduceSum(
+            flux(lev), 0,
+            [=] AMREX_GPU_HOST_DEVICE(
+                amrex::Box const& bx,
+                amrex::Array4<amrex::Real const> const& f_arr) -> amrex::Real {
+                amrex::Real error = 0;
+
+                amrex::Loop(bx, [=, &error](int i, int j, int k) noexcept {
+                    const int idx = (dir == 0 ? i : (dir == 1 ? j : k));
+                    // Difference between actual and expected
+                    if (idx == 0) {
+                        // Both within margin, will average
+                        const amrex::Real flux_answer =
+                            0.5 * ((0.5 - 0.41) * -1.0 + (0.5 - 0.55) * -1.0);
+                        error += std::abs(f_arr(i, j, k) - flux_answer);
+                    } else if (idx == 1) {
+                        // Current is low, neighbor within margin
+                        const amrex::Real flux_answer = (0.2 - 0.1) * 1.0;
+                        error += std::abs(f_arr(i, j, k) - flux_answer);
+                    } else if (idx == 2) {
+                        // Low on both sides, gphi > 0
+                        const amrex::Real flux_answer = (0.2 - 0.1) * -1.0;
+                        error += std::abs(f_arr(i, j, k) - flux_answer);
+                    } else if (idx == 3) {
+                        // Opposite sides of margin, no conditional fits
+                        const amrex::Real flux_answer =
+                            0.5 * ((0.8 - 0.7) * 1.0 + (0.22 - 0.2) * 1.0);
+                        error += std::abs(f_arr(i, j, k) - flux_answer);
+                    } else if (idx == 4) {
+                        // High on both sides, gphi < 0
+                        const amrex::Real flux_answer = (0.8 - 0.7) * -1.0;
+                        error += std::abs(f_arr(i, j, k) - flux_answer);
+                    } else if (idx == 5) {
+                        // High on both sides, gphi > 0
+                        const amrex::Real flux_answer = (0.75 - 0.8) * 1.0;
+                        error += std::abs(f_arr(i, j, k) - flux_answer);
+                    } else if (idx == 6) {
+                        // Current is within margin, neighbor is high
+                        const amrex::Real flux_answer = (0.75 - 0.8) * -1.0;
+                        error += std::abs(f_arr(i, j, k) - flux_answer);
+                    } else if (idx == 7) {
+                        // Current is low, neighbor is within margin
+                        const amrex::Real flux_answer = (0.2 - 0.3) * 1.0;
+                        error += std::abs(f_arr(i, j, k) - flux_answer);
+                    } else if (idx == 8) {
+                        // Opposite sides of margin, no conditional fits
+                        const amrex::Real flux_answer =
+                            0.5 * ((0.6 - 0.7) * -1.0 + (0.2 - 0.3) * -1.0);
+                        error += std::abs(f_arr(i, j, k) - flux_answer);
+                    }
+                });
+
+                return error;
+            });
+    }
+    return error_total;
+}
+
+amrex::Real check_velocity_face_impl(amr_wind::Field& flux, const int& dir)
+{
+    amrex::Real error_total = 0;
+
+    for (int lev = 0; lev < flux.repo().num_active_levels(); ++lev) {
+
+        error_total += amrex::ReduceSum(
+            flux(lev), 0,
+            [=] AMREX_GPU_HOST_DEVICE(
+                amrex::Box const& bx,
+                amrex::Array4<amrex::Real const> const& f_arr) -> amrex::Real {
+                amrex::Real error = 0;
+
+                amrex::Loop(bx, [=, &error](int i, int j, int k) noexcept {
+                    const int idx = (dir == 0 ? i : (dir == 1 ? j : k));
+                    if (idx == 0) {
+                        // gphi > 0, uwpind from the "left"
+                        amrex::Real flux_answer = 1.0;
+                        error += std::abs(f_arr(i, j, k, 0) - flux_answer);
+                        flux_answer = 2.0;
+                        error += std::abs(f_arr(i, j, k, 1) - flux_answer);
+                        flux_answer = 3.0;
+                        error += std::abs(f_arr(i, j, k, 2) - flux_answer);
+                    } else if (idx == 1) {
+                        // gphi < 0, upwind from the "right"
+                        amrex::Real flux_answer = 2.0;
+                        error += std::abs(f_arr(i, j, k, 0) - flux_answer);
+                        flux_answer = 3.0;
+                        error += std::abs(f_arr(i, j, k, 1) - flux_answer);
+                        flux_answer = 4.0;
+                        error += std::abs(f_arr(i, j, k, 2) - flux_answer);
+                    } else if (idx == 2) {
+                        // gphi = 0, average both sides
+                        amrex::Real flux_answer = 0.5 * (2.5 + 2.0);
+                        error += std::abs(f_arr(i, j, k, 0) - flux_answer);
+                        flux_answer = 0.5 * (3.5 + 3.0);
+                        error += std::abs(f_arr(i, j, k, 1) - flux_answer);
+                        flux_answer = 0.5 * (4.5 + 4.0);
+                        error += std::abs(f_arr(i, j, k, 2) - flux_answer);
+                    }
+                });
+
+                return error;
+            });
+    }
+    return error_total;
+}
+} // namespace
+
+TEST_F(VOFOversetOps, alpha_flux)
+{
+    populate_parameters();
+    initialize_mesh();
+
+    auto& repo = sim().repo();
+    const int nghost = 3;
+    auto& vof = repo.declare_field("vof", 1, nghost);
+    auto& tg_vof = repo.declare_field("target_vof", 1, nghost);
+    auto& norm = repo.declare_field("int_normal", 3, nghost);
+    const amrex::Real margin = 0.1;
+
+    // Create flux fields
+    auto& flux_x =
+        repo.declare_field("flux_x", 1, 0, 1, amr_wind::FieldLoc::XFACE);
+    auto& flux_y =
+        repo.declare_field("flux_y", 1, 0, 1, amr_wind::FieldLoc::YFACE);
+    auto& flux_z =
+        repo.declare_field("flux_z", 1, 0, 1, amr_wind::FieldLoc::ZFACE);
+
+    // -- Variations in x direction -- //
+    int dir = 0;
+    // Initialize vof and other fields
+    init_vof_etc(vof, tg_vof, norm, dir);
+    // Populate flux field
+    calc_alpha_flux(flux_x, vof, tg_vof, norm, dir, margin);
+    // Check results
+    amrex::Real error_total = check_alpha_flux_impl(flux_x, dir);
+    amrex::ParallelDescriptor::ReduceRealSum(error_total);
+    EXPECT_NEAR(error_total, 0.0, 1e-15);
+
+    // -- Variations in y direction -- //
+    dir = 1;
+    // Initialize vof and other fields
+    init_vof_etc(vof, tg_vof, norm, dir);
+    // Populate flux field
+    calc_alpha_flux(flux_y, vof, tg_vof, norm, dir, margin);
+    // Check results
+    error_total = check_alpha_flux_impl(flux_y, dir);
+    amrex::ParallelDescriptor::ReduceRealSum(error_total);
+    EXPECT_NEAR(error_total, 0.0, 1e-15);
+
+    // -- Variations in z direction -- //
+    dir = 2;
+    // Initialize vof and other fields
+    init_vof_etc(vof, tg_vof, norm, dir);
+    // Populate flux field
+    calc_alpha_flux(flux_z, vof, tg_vof, norm, dir, margin);
+    // Check results
+    error_total = check_alpha_flux_impl(flux_z, dir);
+    amrex::ParallelDescriptor::ReduceRealSum(error_total);
+    EXPECT_NEAR(error_total, 0.0, 1e-15);
+}
+
+TEST_F(VOFOversetOps, velocity_face)
+{
+    populate_parameters();
+    initialize_mesh();
+
+    auto& repo = sim().repo();
+    const int nghost = 3;
+    auto& velocity = repo.declare_field("velocity", 3, nghost);
+    auto& vof = repo.declare_field("vof", 1, nghost);
+
+    // Create flux fields
+    auto& flux_x =
+        repo.declare_field("flux_x", 3, 0, 1, amr_wind::FieldLoc::XFACE);
+    auto& flux_y =
+        repo.declare_field("flux_y", 3, 0, 1, amr_wind::FieldLoc::YFACE);
+    auto& flux_z =
+        repo.declare_field("flux_z", 3, 0, 1, amr_wind::FieldLoc::ZFACE);
+
+    // -- Variations in x direction -- //
+    int dir = 0;
+    // Initialize velocity and vof
+    init_velocity_etc(velocity, vof, dir);
+    // Populate flux field
+    calc_velocity_face(flux_x, velocity, vof, dir);
+    // Check results
+    amrex::Real error_total = check_velocity_face_impl(flux_x, dir);
+    amrex::ParallelDescriptor::ReduceRealSum(error_total);
+    EXPECT_NEAR(error_total, 0.0, 1e-15);
+
+    // -- Variations in y direction -- //
+    dir = 1;
+    // Initialize velocity and vof
+    init_velocity_etc(velocity, vof, dir);
+    // Populate flux field
+    calc_velocity_face(flux_y, velocity, vof, dir);
+    // Check results
+    error_total = check_velocity_face_impl(flux_y, dir);
+    amrex::ParallelDescriptor::ReduceRealSum(error_total);
+    EXPECT_NEAR(error_total, 0.0, 1e-15);
+
+    // -- Variations in z direction -- //
+    dir = 2;
+    // Initialize velocity and vof
+    init_velocity_etc(velocity, vof, dir);
+    // Populate flux field
+    calc_velocity_face(flux_z, velocity, vof, dir);
+    // Check results
+    error_total = check_velocity_face_impl(flux_z, dir);
+    amrex::ParallelDescriptor::ReduceRealSum(error_total);
+    EXPECT_NEAR(error_total, 0.0, 1e-15);
+}
+
+} // namespace amr_wind_tests
\ No newline at end of file
diff --git a/unit_tests/multiphase/test_vof_plic.cpp b/unit_tests/multiphase/test_vof_plic.cpp
index 71cefbbbcd..8d12e54632 100644
--- a/unit_tests/multiphase/test_vof_plic.cpp
+++ b/unit_tests/multiphase/test_vof_plic.cpp
@@ -39,11 +39,10 @@ void initialize_volume_fractions(
 {
     // grow the box by 1 so that x,y,z go out of bounds and min(max()) corrects
     // it and it fills the ghosts with wall values
-    const int d = dir;
     amrex::ParallelFor(grow(bx, 1), [=] AMREX_GPU_DEVICE(int i, int j, int k) {
-        int ii = (d != 0 ? i : 0);
-        int jj = (d != 1 ? j : 0);
-        int kk = (d != 2 ? k : 0);
+        int ii = (dir != 0 ? i : 0);
+        int jj = (dir != 1 ? j : 0);
+        int kk = (dir != 2 ? k : 0);
         if (ii + jj + kk > 3) {
             vof_arr(i, j, k) = 0.0;
         }
@@ -73,10 +72,9 @@ void initialize_volume_fractions_horizontal(
 {
     // grow the box by 1 so that x,y,z go out of bounds and min(max()) corrects
     // it and it fills the ghosts with wall values
-    const int d = dir;
     const amrex::Real vv = vof_val;
     amrex::ParallelFor(grow(bx, 1), [=] AMREX_GPU_DEVICE(int i, int j, int k) {
-        int ii = (d == 0 ? i : (d == 1 ? j : k));
+        int ii = (dir == 0 ? i : (dir == 1 ? j : k));
         if (ii > 1) {
             vof_arr(i, j, k) = 0.0;
         }
@@ -98,10 +96,105 @@ void init_vof_h(amr_wind::Field& vof, const amrex::Real vof_val, const int dir)
     });
 }
 
+void initialize_volume_fractions(
+    const amrex::Box& bx, const amrex::Array4<amrex::Real>& vof_arr)
+{
+    amrex::ParallelFor(grow(bx, 1), [=] AMREX_GPU_DEVICE(int i, int j, int k) {
+        vof_arr(i, j, k) = 0.13 * ((amrex::Real)i - 1.5) +
+                           0.04 * std::pow((amrex::Real)j - 1., 2) +
+                           0.01 * std::pow((amrex::Real)k - 2., 3) + 0.5;
+    });
+}
+
+void init_vof(amr_wind::Field& vof)
+{
+    run_algorithm(vof, [&](const int lev, const amrex::MFIter& mfi) {
+        auto vof_arr = vof(lev).array(mfi);
+        const auto& bx = mfi.validbox();
+        initialize_volume_fractions(bx, vof_arr);
+    });
+}
+
+void initialize_iblank_distribution(
+    const int dir, const amrex::Box& bx, const amrex::Array4<int>& iblk_arr)
+{
+    amrex::ParallelFor(grow(bx, 1), [=] AMREX_GPU_DEVICE(int i, int j, int k) {
+        int ii = (dir != 0 ? i : 0);
+        int jj = (dir != 1 ? j : 0);
+        int kk = (dir != 2 ? k : 0);
+        if (ii + jj + kk > 4 || ii + jj + kk < 2) {
+            iblk_arr(i, j, k) = -1;
+        } else {
+            iblk_arr(i, j, k) = 1;
+        }
+    });
+}
+
+void init_iblank(amr_wind::IntField& iblank, const int dir)
+{
+    run_algorithm(iblank, [&](const int lev, const amrex::MFIter& mfi) {
+        auto iblk_arr = iblank(lev).array(mfi);
+        const auto& bx = mfi.validbox();
+        initialize_iblank_distribution(dir, bx, iblk_arr);
+    });
+}
+
+void initialize_volume_fractions_iblank(
+    const int dir,
+    const amrex::Box& bx,
+    const amrex::Array4<amrex::Real>& vof_arr,
+    const amrex::Array4<int>& iblk_arr)
+{
+    // Does a horizontal interface that is different outside of iblank region
+    amrex::ParallelFor(grow(bx, 1), [=] AMREX_GPU_DEVICE(int i, int j, int k) {
+        // Turns on index tangent to interface
+        int iit = (dir != 0 ? i : 1);
+        int jjt = (dir != 1 ? j : 1);
+        int kkt = (dir != 2 ? k : 1);
+        // Turns on index normal to interface
+        int iin = (dir == 0 ? i : 0);
+        int jjn = (dir == 1 ? j : 0);
+        int kkn = (dir == 2 ? k : 0);
+        // Ordinary vof distribution
+        if (iin + jjn + kkn > 2) {
+            vof_arr(i, j, k) = 0.0;
+        }
+        if (iin + jjn + kkn == 2) {
+            vof_arr(i, j, k) = 0.5;
+        }
+        if (iin + jjn + kkn < 2) {
+            vof_arr(i, j, k) = 1.0;
+        }
+        // iblank distribution
+        iblk_arr(i, j, k) = 1;
+        if (iit > 2 || iit < 1 || jjt > 2 || jjt < 1 || kkt > 2 || kkt < 1) {
+            // Restrict iblank = 1 to central block
+            iblk_arr(i, j, k) = -1;
+            // Change vof where iblank = -1
+            if (iin + jjn + kkn == 3) {
+                vof_arr(i, j, k) = 0.25;
+            }
+            if (iin + jjn + kkn == 1) {
+                vof_arr(i, j, k) = 0.75;
+            }
+        }
+    });
+}
+
+void init_vof_iblank(
+    amr_wind::Field& vof, amr_wind::IntField& iblank, const int dir)
+{
+    run_algorithm(vof, [&](const int lev, const amrex::MFIter& mfi) {
+        auto vof_arr = vof(lev).array(mfi);
+        auto iblk_arr = iblank(lev).array(mfi);
+        const auto& bx = mfi.validbox();
+        initialize_volume_fractions_iblank(dir, bx, vof_arr, iblk_arr);
+    });
+}
+
 amrex::Real normal_vector_test_impl(amr_wind::Field& vof, const int dir)
 {
     amrex::Real error_total = 0.0;
-    const int d = dir;
 
     for (int lev = 0; lev < vof.repo().num_active_levels(); ++lev) {
 
@@ -118,15 +211,141 @@ amrex::Real normal_vector_test_impl(amr_wind::Field& vof, const int dir)
                     amr_wind::multiphase::mixed_youngs_central_normal(
                         i, j, k, vof_arr, mx, my, mz);
 
-                    int ii = (d != 0 ? i : 0);
-                    int jj = (d != 1 ? j : 0);
-                    int kk = (d != 2 ? k : 0);
+                    int ii = (dir != 0 ? i : 0);
+                    int jj = (dir != 1 ? j : 0);
+                    int kk = (dir != 2 ? k : 0);
 
                     // Use L1 norm, check cells where slope is known
                     if (ii + jj + kk == 3) {
-                        error += std::abs(mx - (d != 0 ? 0.5 : 0.0));
-                        error += std::abs(my - (d != 1 ? 0.5 : 0.0));
-                        error += std::abs(mz - (d != 2 ? 0.5 : 0.0));
+                        error += std::abs(mx - (dir != 0 ? 0.5 : 0.0));
+                        error += std::abs(my - (dir != 1 ? 0.5 : 0.0));
+                        error += std::abs(mz - (dir != 2 ? 0.5 : 0.0));
+                    }
+                });
+
+                return error;
+            });
+    }
+    return error_total;
+}
+
+amrex::Real normal_vector_neumann_test_impl(
+    amr_wind::Field& vof, amr_wind::IntField& iblk_fld)
+{
+    amrex::Real error_total = 0.0;
+
+    for (int lev = 0; lev < vof.repo().num_active_levels(); ++lev) {
+
+        error_total += amrex::ReduceSum(
+            vof(lev), iblk_fld(lev), 0,
+            [=] AMREX_GPU_HOST_DEVICE(
+                amrex::Box const& bx,
+                amrex::Array4<amrex::Real const> const& vof_arr,
+                amrex::Array4<int const> const& iblank) -> amrex::Real {
+                amrex::Real error = 0.0;
+
+                amrex::Loop(bx, [=, &error](int i, int j, int k) noexcept {
+                    int ibdy =
+                        (iblank(i, j, k) != iblank(i - 1, j, k)) ? -1 : 0;
+                    int jbdy =
+                        (iblank(i, j, k) != iblank(i, j - 1, k)) ? -1 : 0;
+                    int kbdy =
+                        (iblank(i, j, k) != iblank(i, j, k - 1)) ? -1 : 0;
+                    ibdy = (iblank(i, j, k) != iblank(i + 1, j, k)) ? +1 : ibdy;
+                    jbdy = (iblank(i, j, k) != iblank(i, j + 1, k)) ? +1 : jbdy;
+                    kbdy = (iblank(i, j, k) != iblank(i, j, k + 1)) ? +1 : kbdy;
+                    amrex::Real mxn, myn, mzn;
+                    amr_wind::multiphase::
+                        youngs_finite_difference_normal_neumann(
+                            i, j, k, ibdy, jbdy, kbdy, vof_arr, mxn, myn, mzn);
+                    amrex::Real mx, my, mz;
+                    amr_wind::multiphase::youngs_finite_difference_normal(
+                        i, j, k, vof_arr, mx, my, mz);
+
+                    // Use L1 norm, check against non-neumann implementation
+                    // Slope across overset boundary should be different
+                    constexpr amrex::Real slp_tol = 1e-8;
+                    if (ibdy != 0) {
+                        // x slope should be different
+                        error += std::abs(mx - mxn) > slp_tol ? 0. : 1.0;
+                    }
+                    if (jbdy != 0) {
+                        // y slope should be different
+                        error += std::abs(my - myn) > slp_tol ? 0. : 1.0;
+                    }
+                    if (kbdy != 0) {
+                        // z slope should be different
+                        error += std::abs(mz - mzn) > slp_tol ? 0. : 1.0;
+                    }
+                    // Slope should otherwise be the same
+                    if (ibdy == 0 && jbdy == 0 && kbdy == 0) {
+                        error += std::abs(mx - mxn);
+                        error += std::abs(my - myn);
+                        error += std::abs(mz - mzn);
+                    }
+                });
+
+                return error;
+            });
+    }
+    return error_total;
+}
+
+amrex::Real normal_vector_neumann_test_impl(
+    amr_wind::Field& vof, amr_wind::IntField& iblk_fld, const int& dir)
+{
+    const amrex::Real ref_m = 16.0 * (1.0 - 0.0);
+    amrex::Real error_total = 0.0;
+
+    for (int lev = 0; lev < vof.repo().num_active_levels(); ++lev) {
+
+        error_total += amrex::ReduceSum(
+            vof(lev), iblk_fld(lev), 0,
+            [=] AMREX_GPU_HOST_DEVICE(
+                amrex::Box const& bx,
+                amrex::Array4<amrex::Real const> const& vof_arr,
+                amrex::Array4<int const> const& iblank) -> amrex::Real {
+                amrex::Real error = 0.0;
+
+                amrex::Loop(bx, [=, &error](int i, int j, int k) noexcept {
+                    int ibdy =
+                        (iblank(i, j, k) != iblank(i - 1, j, k)) ? -1 : 0;
+                    int jbdy =
+                        (iblank(i, j, k) != iblank(i, j - 1, k)) ? -1 : 0;
+                    int kbdy =
+                        (iblank(i, j, k) != iblank(i, j, k - 1)) ? -1 : 0;
+                    ibdy = (iblank(i, j, k) != iblank(i + 1, j, k)) ? +1 : ibdy;
+                    jbdy = (iblank(i, j, k) != iblank(i, j + 1, k)) ? +1 : jbdy;
+                    kbdy = (iblank(i, j, k) != iblank(i, j, k + 1)) ? +1 : kbdy;
+                    amrex::Real mxn, myn, mzn;
+                    amr_wind::multiphase::
+                        youngs_finite_difference_normal_neumann(
+                            i, j, k, ibdy, jbdy, kbdy, vof_arr, mxn, myn, mzn);
+
+                    // Use L1 norm, check for 0
+                    if (iblank(i, j, k) == 1) {
+                        // Only interested in normals from field cells
+
+                        // Slope in non-normal directions should be 0
+                        if (dir != 0) {
+                            error += std::abs(mxn - 0.);
+                        }
+                        if (dir != 1) {
+                            error += std::abs(myn - 0.);
+                        }
+                        if (dir != 2) {
+                            error += std::abs(mzn - 0.);
+                        }
+                        // Slope in normal direction, at center, should be same
+                        if (dir == 0 && i == 2) {
+                            error += std::abs(mxn - ref_m);
+                        }
+                        if (dir == 1 && j == 2) {
+                            error += std::abs(myn - ref_m);
+                        }
+                        if (dir == 2 && k == 2) {
+                            error += std::abs(mzn - ref_m);
+                        }
                     }
                 });
 
@@ -139,7 +358,6 @@ amrex::Real normal_vector_test_impl(amr_wind::Field& vof, const int dir)
 amrex::Real fit_plane_test_impl(amr_wind::Field& vof, const int dir)
 {
     amrex::Real error_total = 0.0;
-    const int d = dir;
 
     for (int lev = 0; lev < vof.repo().num_active_levels(); ++lev) {
 
@@ -152,9 +370,9 @@ amrex::Real fit_plane_test_impl(amr_wind::Field& vof, const int dir)
                 amrex::Real error = 0.0;
 
                 amrex::Loop(bx, [=, &error](int i, int j, int k) noexcept {
-                    int ii = (d != 0 ? i : 0);
-                    int jj = (d != 1 ? j : 0);
-                    int kk = (d != 2 ? k : 0);
+                    int ii = (dir != 0 ? i : 0);
+                    int jj = (dir != 1 ? j : 0);
+                    int kk = (dir != 2 ? k : 0);
                     // Check multiphase cells
                     if (ii + jj + kk == 3) {
                         amrex::Real mx, my, mz, alpha;
@@ -162,9 +380,9 @@ amrex::Real fit_plane_test_impl(amr_wind::Field& vof, const int dir)
                             i, j, k, vof_arr, mx, my, mz, alpha);
 
                         // Check slope
-                        error += std::abs(mx - (d != 0 ? 0.5 : 0.0));
-                        error += std::abs(my - (d != 1 ? 0.5 : 0.0));
-                        error += std::abs(mz - (d != 2 ? 0.5 : 0.0));
+                        error += std::abs(mx - (dir != 0 ? 0.5 : 0.0));
+                        error += std::abs(my - (dir != 1 ? 0.5 : 0.0));
+                        error += std::abs(mz - (dir != 2 ? 0.5 : 0.0));
                         // Check intercept
                         error += std::abs(alpha - 0.5);
                     }
@@ -180,7 +398,6 @@ amrex::Real fit_plane_test_impl_h(
     amr_wind::Field& vof, const amrex::Real vof_val, const int dir)
 {
     amrex::Real error_total = 0.0;
-    const int d = dir;
     const amrex::Real vv = vof_val;
 
     for (int lev = 0; lev < vof.repo().num_active_levels(); ++lev) {
@@ -194,7 +411,7 @@ amrex::Real fit_plane_test_impl_h(
                 amrex::Real error = 0.0;
 
                 amrex::Loop(bx, [=, &error](int i, int j, int k) noexcept {
-                    int ii = (d == 0 ? i : (d == 1 ? j : k));
+                    int ii = (dir == 0 ? i : (dir == 1 ? j : k));
                     // Check multiphase cells
                     if (ii == 1) {
                         amrex::Real mx, my, mz, alpha;
@@ -202,9 +419,9 @@ amrex::Real fit_plane_test_impl_h(
                             i, j, k, vof_arr, mx, my, mz, alpha);
 
                         // Check slope
-                        error += std::abs(mx - (d == 0 ? 1.0 : 0.0));
-                        error += std::abs(my - (d == 1 ? 1.0 : 0.0));
-                        error += std::abs(mz - (d == 2 ? 1.0 : 0.0));
+                        error += std::abs(mx - (dir == 0 ? 1.0 : 0.0));
+                        error += std::abs(my - (dir == 1 ? 1.0 : 0.0));
+                        error += std::abs(mz - (dir == 2 ? 1.0 : 0.0));
                         // Check intercept
                         error += std::abs(alpha - vv);
                     }
@@ -367,4 +584,63 @@ TEST_F(VOFOpTest, interface_plane)
     }
 }
 
+TEST_F(VOFOpTest, interface_normal_neumann)
+{
+
+    constexpr double tol = 1.0e-11;
+
+    populate_parameters();
+    {
+        amrex::ParmParse pp("geometry");
+        amrex::Vector<int> periodic{{0, 0, 0}};
+        pp.addarr("is_periodic", periodic);
+    }
+
+    initialize_mesh();
+
+    auto& repo = sim().repo();
+    const int ncomp = 1;
+    const int nghost = 3;
+    auto& vof = repo.declare_field("vof", ncomp, nghost);
+    auto& iblank = repo.declare_int_field("iblank", ncomp, 1);
+
+    // Check agreement / disagreement with ordinary calculation
+    amrex::Real error_total = 0.0;
+    // iblank constant in x, vof varies in every direction
+    init_vof(vof);
+    init_iblank(iblank, 0);
+    error_total = normal_vector_neumann_test_impl(vof, iblank);
+    amrex::ParallelDescriptor::ReduceRealSum(error_total);
+    EXPECT_NEAR(error_total, 0.0, tol);
+    // iblank constant in y, vof varies in every direction
+    init_vof(vof);
+    init_iblank(iblank, 1);
+    error_total = normal_vector_neumann_test_impl(vof, iblank);
+    amrex::ParallelDescriptor::ReduceRealSum(error_total);
+    EXPECT_NEAR(error_total, 0.0, tol);
+    // iblank constant in z, vof varies in every direction
+    init_vof(vof);
+    init_iblank(iblank, 2);
+    error_total = normal_vector_neumann_test_impl(vof, iblank);
+    amrex::ParallelDescriptor::ReduceRealSum(error_total);
+    EXPECT_NEAR(error_total, 0.0, tol);
+
+    // Confirm that neumann gets 0 when it should
+    // iblank varies in y and z, vof varies in x
+    init_vof_iblank(vof, iblank, 0);
+    error_total = normal_vector_neumann_test_impl(vof, iblank, 0);
+    amrex::ParallelDescriptor::ReduceRealSum(error_total);
+    EXPECT_NEAR(error_total, 0.0, tol);
+    // iblank varies in x and z, vof varies in y
+    init_vof_iblank(vof, iblank, 1);
+    error_total = normal_vector_neumann_test_impl(vof, iblank, 1);
+    amrex::ParallelDescriptor::ReduceRealSum(error_total);
+    EXPECT_NEAR(error_total, 0.0, tol);
+    // iblank varies in x and y, vof varies in z
+    init_vof_iblank(vof, iblank, 2);
+    error_total = normal_vector_neumann_test_impl(vof, iblank, 2);
+    amrex::ParallelDescriptor::ReduceRealSum(error_total);
+    EXPECT_NEAR(error_total, 0.0, tol);
+}
+
 } // namespace amr_wind_tests