diff --git a/examples/dalf3/dalf3.cxx b/examples/dalf3/dalf3.cxx
index 46a9a4286d..b27ae97113 100644
--- a/examples/dalf3/dalf3.cxx
+++ b/examples/dalf3/dalf3.cxx
@@ -71,8 +71,6 @@ class DALF3 : public PhysicsModel {
 
   BoutReal viscosity, hyper_viscosity;
 
-  bool smooth_separatrix;
-
   FieldGroup comms;
 
   std::unique_ptr<Laplacian> phiSolver{nullptr};  // Laplacian solver in X-Z
@@ -132,7 +130,6 @@ class DALF3 : public PhysicsModel {
     viscosity = options["viscosity"].withDefault(-1.0);
     hyper_viscosity = options["hyper_viscosity"].withDefault(-1.0);
     viscosity_par = options["viscosity_par"].withDefault(-1.0);
-    smooth_separatrix = options["smooth_separatrix"].withDefault(false);
 
     filter_z = options["filter_z"].withDefault(false);
 
@@ -466,11 +463,6 @@ class DALF3 : public PhysicsModel {
     ddt(Pe) = -bracket(phi, Pet, bm)
               + Pet * (Kappa(phi - Pe) + B0 * Grad_parP(jpar - Vpar) / B0);
 
-    if (smooth_separatrix) {
-      // Experimental smoothing across separatrix
-      ddt(Vort) += mesh->smoothSeparatrix(Vort);
-    }
-
     if (filter_z) {
       ddt(Pe) = filter(ddt(Pe), 1);
     }
diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx
index 2e27bed777..da8ac51127 100644
--- a/include/bout/mesh.hxx
+++ b/include/bout/mesh.hxx
@@ -505,9 +505,6 @@ public:
   virtual void addBoundaryPar(std::shared_ptr<BoundaryRegionPar> UNUSED(bndry),
                               BoundaryParType UNUSED(type)) {}
 
-  /// Branch-cut special handling (experimental)
-  virtual Field3D smoothSeparatrix(const Field3D& f) { return f; }
-
   virtual BoutReal GlobalX(int jx) const = 0;      ///< Continuous X index between 0 and 1
   virtual BoutReal GlobalY(int jy) const = 0;      ///< Continuous Y index (0 -> 1)
   virtual BoutReal GlobalX(BoutReal jx) const = 0; ///< Continuous X index between 0 and 1
diff --git a/src/mesh/impls/bout/boutmesh.cxx b/src/mesh/impls/bout/boutmesh.cxx
index 8917fa3753..11ea2f5c89 100644
--- a/src/mesh/impls/bout/boutmesh.cxx
+++ b/src/mesh/impls/bout/boutmesh.cxx
@@ -3115,47 +3115,6 @@ void BoutMesh::addBoundaryPar(std::shared_ptr<BoundaryRegionPar> bndry,
   par_boundary[static_cast<int>(BoundaryParType::all)].push_back(bndry);
 }
 
-Field3D BoutMesh::smoothSeparatrix(const Field3D& f) {
-  Field3D result{emptyFrom(f)};
-  if ((ixseps_inner > 0) && (ixseps_inner < nx - 1)) {
-    if (XPROC(ixseps_inner) == PE_XIND) {
-      int x = getLocalXIndex(ixseps_inner);
-      for (int y = 0; y < LocalNy; y++) {
-        for (int z = 0; z < LocalNz; z++) {
-          result(x, y, z) = 0.5 * (f(x, y, z) + f(x - 1, y, z));
-        }
-      }
-    }
-    if (XPROC(ixseps_inner - 1) == PE_XIND) {
-      int x = getLocalXIndex(ixseps_inner - 1);
-      for (int y = 0; y < LocalNy; y++) {
-        for (int z = 0; z < LocalNz; z++) {
-          result(x, y, z) = 0.5 * (f(x, y, z) + f(x + 1, y, z));
-        }
-      }
-    }
-  }
-  if ((ixseps_outer > 0) && (ixseps_outer < nx - 1) && (ixseps_outer != ixseps_inner)) {
-    if (XPROC(ixseps_outer) == PE_XIND) {
-      int x = getLocalXIndex(ixseps_outer);
-      for (int y = 0; y < LocalNy; y++) {
-        for (int z = 0; z < LocalNz; z++) {
-          result(x, y, z) = 0.5 * (f(x, y, z) + f(x - 1, y, z));
-        }
-      }
-    }
-    if (XPROC(ixseps_outer - 1) == PE_XIND) {
-      int x = getLocalXIndex(ixseps_outer - 1);
-      for (int y = 0; y < LocalNy; y++) {
-        for (int z = 0; z < LocalNz; z++) {
-          result(x, y, z) = 0.5 * (f(x, y, z) + f(x + 1, y, z));
-        }
-      }
-    }
-  }
-  return result;
-}
-
 BoutReal BoutMesh::GlobalX(int jx) const {
   if (symmetricGlobalX) {
     // With this definition the boundary sits dx/2 away form the first/last inner points
diff --git a/src/mesh/impls/bout/boutmesh.hxx b/src/mesh/impls/bout/boutmesh.hxx
index a22e0e4473..876edab1da 100644
--- a/src/mesh/impls/bout/boutmesh.hxx
+++ b/src/mesh/impls/bout/boutmesh.hxx
@@ -167,8 +167,6 @@ public:
                       BoundaryParType type) override;
   std::set<std::string> getPossibleBoundaries() const override;
 
-  Field3D smoothSeparatrix(const Field3D& f) override;
-
   int getNx() const { return nx; }
   int getNy() const { return ny; }