diff --git a/examples/performance/bracket/bracket.cxx b/examples/performance/bracket/bracket.cxx
index 50012510ff..0ef9e241fb 100644
--- a/examples/performance/bracket/bracket.cxx
+++ b/examples/performance/bracket/bracket.cxx
@@ -68,9 +68,6 @@ int main(int argc, char** argv) {
     ITERATOR_TEST_BLOCK("Bracket [2D,3D] ARAKAWA",
                         result = bracket(a, c, BRACKET_ARAKAWA););
 
-    ITERATOR_TEST_BLOCK("Bracket [2D,3D] ARAKAWA_OLD",
-                        result = bracket(a, c, BRACKET_ARAKAWA_OLD););
-
     ITERATOR_TEST_BLOCK("Bracket [2D,3D] SIMPLE",
                         result = bracket(a, c, BRACKET_SIMPLE););
 
@@ -81,21 +78,12 @@ int main(int argc, char** argv) {
     ITERATOR_TEST_BLOCK("Bracket [3D,3D] ARAKAWA",
                         result = bracket(a, b, BRACKET_ARAKAWA););
 
-    ITERATOR_TEST_BLOCK("Bracket [3D,3D] ARAKAWA_OLD",
-                        result = bracket(a, b, BRACKET_ARAKAWA_OLD););
-
     ITERATOR_TEST_BLOCK("Bracket [3D,3D] SIMPLE",
                         result = bracket(a, b, BRACKET_SIMPLE););
 
     ITERATOR_TEST_BLOCK("Bracket [3D,3D] DEFAULT", result = bracket(a, b, BRACKET_STD););
   }
 
-  // Uncomment below for a "correctness" check
-  // Field3D resNew = bracket(a, b, BRACKET_ARAKAWA); mesh->communicate(resNew);
-  // Field3D resOld = bracket(a, b, BRACKET_ARAKAWA_OLD); mesh->communicate(resOld);
-  // time_output << "Max abs diff is
-  // "<<max(abs(resNew-resOld),true)/max(abs(resOld),true)<<std::endl;
-
   if (profileMode) {
     int nthreads = 0;
 #if BOUT_USE_OPENMP
diff --git a/include/bout/difops.hxx b/include/bout/difops.hxx
index 71053d454a..8b502002f2 100644
--- a/include/bout/difops.hxx
+++ b/include/bout/difops.hxx
@@ -271,18 +271,16 @@ Field3D b0xGrad_dot_Grad(const Field3D& phi, const Field3D& A,
  * Poisson bracket methods
  */
 enum class BRACKET_METHOD {
-  standard,   ///< Use b0xGrad_dot_Grad
-  simple,     ///< Keep only terms in X-Z
-  arakawa,    ///< Arakawa method in X-Z (optimised)
-  ctu,        ///< Corner Transport Upwind (CTU) method. Explicit method only, needs the
-              ///  timestep from the solver
-  arakawa_old ///< Older version, for regression testing of optimised version.
+  standard, ///< Use b0xGrad_dot_Grad
+  simple,   ///< Keep only terms in X-Z
+  arakawa,  ///< Arakawa method in X-Z
+  ctu,      ///< Corner Transport Upwind (CTU) method. Explicit method only, needs the
+            ///  timestep from the solver
 };
 constexpr BRACKET_METHOD BRACKET_STD = BRACKET_METHOD::standard;
 constexpr BRACKET_METHOD BRACKET_SIMPLE = BRACKET_METHOD::simple;
 constexpr BRACKET_METHOD BRACKET_ARAKAWA = BRACKET_METHOD::arakawa;
 constexpr BRACKET_METHOD BRACKET_CTU = BRACKET_METHOD::ctu;
-constexpr BRACKET_METHOD BRACKET_ARAKAWA_OLD = BRACKET_METHOD::arakawa_old;
 
 /*!
  * Compute advection operator terms, which can be cast as
diff --git a/manual/sphinx/user_docs/python_boutpp.rst b/manual/sphinx/user_docs/python_boutpp.rst
index 939a817ffe..e79facbf55 100644
--- a/manual/sphinx/user_docs/python_boutpp.rst
+++ b/manual/sphinx/user_docs/python_boutpp.rst
@@ -160,7 +160,7 @@ A real example - check derivative contributions:
               phi[-1,:,z]=phi_arr
       with open(path+"/equilibrium/phi_eq.dat","rb") as inf:
           phi_arr=np.fromfile(inf,dtype=np.double)
-          bm="BRACKET_ARAKAWA_OLD"
+          bm="BRACKET_ARAKAWA"
 
           for tind in range(len(times)):
               vort     = Field3D.fromCollect("vort"     ,path=path,tind=tind,info=False)
diff --git a/src/mesh/difops.cxx b/src/mesh/difops.cxx
index 42fa4d6ca5..53ba6f5ad3 100644
--- a/src/mesh/difops.cxx
+++ b/src/mesh/difops.cxx
@@ -771,46 +771,6 @@ Field3D bracket(const Field3D& f, const Field2D& g, BRACKET_METHOD method,
 
     break;
   }
-  case BRACKET_ARAKAWA_OLD: {
-#if not(BOUT_USE_METRIC_3D)
-    const int ncz = mesh->LocalNz;
-    BOUT_OMP_PERF(parallel for)
-    for (int jx = mesh->xstart; jx <= mesh->xend; jx++) {
-      for (int jy = mesh->ystart; jy <= mesh->yend; jy++) {
-        const BoutReal partialFactor = 1.0 / (12 * metric->dz(jx, jy));
-        const BoutReal spacingFactor = partialFactor / metric->dx(jx, jy);
-        for (int jz = 0; jz < mesh->LocalNz; jz++) {
-          const int jzp = jz + 1 < ncz ? jz + 1 : 0;
-          // Above is alternative to const int jzp = (jz + 1) % ncz;
-          const int jzm = jz - 1 >= 0 ? jz - 1 : ncz - 1;
-          // Above is alternative to const int jzmTmp = (jz - 1 + ncz) % ncz;
-
-          // J++ = DDZ(f)*DDX(g) - DDX(f)*DDZ(g)
-          BoutReal Jpp =
-              ((f(jx, jy, jzp) - f(jx, jy, jzm)) * (g(jx + 1, jy) - g(jx - 1, jy))
-               - (f(jx + 1, jy, jz) - f(jx - 1, jy, jz)) * (g(jx, jy) - g(jx, jy)));
-
-          // J+x
-          BoutReal Jpx = (g(jx + 1, jy) * (f(jx + 1, jy, jzp) - f(jx + 1, jy, jzm))
-                          - g(jx - 1, jy) * (f(jx - 1, jy, jzp) - f(jx - 1, jy, jzm))
-                          - g(jx, jy) * (f(jx + 1, jy, jzp) - f(jx - 1, jy, jzp))
-                          + g(jx, jy) * (f(jx + 1, jy, jzm) - f(jx - 1, jy, jzm)));
-
-          // Jx+
-          BoutReal Jxp = (g(jx + 1, jy) * (f(jx, jy, jzp) - f(jx + 1, jy, jz))
-                          - g(jx - 1, jy) * (f(jx - 1, jy, jz) - f(jx, jy, jzm))
-                          - g(jx - 1, jy) * (f(jx, jy, jzp) - f(jx - 1, jy, jz))
-                          + g(jx + 1, jy) * (f(jx + 1, jy, jz) - f(jx, jy, jzm)));
-
-          result(jx, jy, jz) = (Jpp + Jpx + Jxp) * spacingFactor;
-        }
-      }
-    }
-#else
-    throw BoutException("BRACKET_ARAKAWA_OLD not valid with 3D metrics yet.");
-#endif
-    break;
-  }
   case BRACKET_SIMPLE: {
     // Use a subset of terms for comparison to BOUT-06
     result = VDDX(DDZ(f, outloc), g, outloc);
@@ -1090,63 +1050,6 @@ Field3D bracket(const Field3D& f, const Field3D& g, BRACKET_METHOD method,
     }
     break;
   }
-  case BRACKET_ARAKAWA_OLD: {
-    // Arakawa scheme for perpendicular flow
-
-    const int ncz = mesh->LocalNz;
-
-    // We need to discard const qualifier in order to manipulate
-    // storage array directly
-    Field3D f_temp = f;
-    Field3D g_temp = g;
-
-    BOUT_OMP_PERF(parallel for)
-    for (int jx = mesh->xstart; jx <= mesh->xend; jx++) {
-      for (int jy = mesh->ystart; jy <= mesh->yend; jy++) {
-#if not(BOUT_USE_METRIC_3D)
-        const BoutReal spacingFactor =
-            1.0 / (12 * metric->dz(jx, jy) * metric->dx(jx, jy));
-#endif
-        const BoutReal* Fxm = f_temp(jx - 1, jy);
-        const BoutReal* Fx = f_temp(jx, jy);
-        const BoutReal* Fxp = f_temp(jx + 1, jy);
-        const BoutReal* Gxm = g_temp(jx - 1, jy);
-        const BoutReal* Gx = g_temp(jx, jy);
-        const BoutReal* Gxp = g_temp(jx + 1, jy);
-        for (int jz = 0; jz < mesh->LocalNz; jz++) {
-#if BOUT_USE_METRIC_3D
-          const BoutReal spacingFactor =
-              1.0 / (12 * metric->dz(jx, jy, jz) * metric->dx(jx, jy, jz));
-#endif
-          const int jzp = jz + 1 < ncz ? jz + 1 : 0;
-          // Above is alternative to const int jzp = (jz + 1) % ncz;
-          const int jzm = jz - 1 >= 0 ? jz - 1 : ncz - 1;
-          // Above is alternative to const int jzm = (jz - 1 + ncz) % ncz;
-
-          // J++ = DDZ(f)*DDX(g) - DDX(f)*DDZ(g)
-          // NOLINTNEXTLINE
-          BoutReal Jpp = ((Fx[jzp] - Fx[jzm]) * (Gxp[jz] - Gxm[jz])
-                          - (Fxp[jz] - Fxm[jz]) * (Gx[jzp] - Gx[jzm]));
-
-          // J+x
-          // NOLINTNEXTLINE
-          BoutReal Jpx =
-              (Gxp[jz] * (Fxp[jzp] - Fxp[jzm]) - Gxm[jz] * (Fxm[jzp] - Fxm[jzm])
-               - Gx[jzp] * (Fxp[jzp] - Fxm[jzp]) + Gx[jzm] * (Fxp[jzm] - Fxm[jzm]));
-
-          // Jx+
-          // NOLINTNEXTLINE
-          BoutReal Jxp =
-              (Gxp[jzp] * (Fx[jzp] - Fxp[jz]) - Gxm[jzm] * (Fxm[jz] - Fx[jzm])
-               - Gxm[jzp] * (Fx[jzp] - Fxm[jz]) + Gxp[jzm] * (Fxp[jz] - Fx[jzm]));
-
-          result(jx, jy, jz) = (Jpp + Jpx + Jxp) * spacingFactor;
-        }
-      }
-    }
-
-    break;
-  }
   case BRACKET_SIMPLE: {
     // Use a subset of terms for comparison to BOUT-06
     result = VDDX(DDZ(f, outloc), g, outloc) + VDDZ(-DDX(f, outloc), g, outloc);
diff --git a/tests/MMS/spatial/advection/runtest b/tests/MMS/spatial/advection/runtest
index 9114664b47..fa06c478ee 100755
--- a/tests/MMS/spatial/advection/runtest
+++ b/tests/MMS/spatial/advection/runtest
@@ -43,7 +43,6 @@ build_and_log("MMS steady-state advection test")
 # List of options to be passed for each test
 options = [
     ("method=2", "Arakawa", "-^", 2),
-    ("method=4", "Arakawa-old", "-.", 2),
     (
         "method=1 mesh:ddx:upwind=u1 mesh:ddz:upwind=u1",
         "SIMPLE: 1st order upwind",
diff --git a/tools/pylib/_boutpp_build/other_enums.hxx b/tools/pylib/_boutpp_build/other_enums.hxx
index 30480a1007..a1da5ded7e 100644
--- a/tools/pylib/_boutpp_build/other_enums.hxx
+++ b/tools/pylib/_boutpp_build/other_enums.hxx
@@ -1,4 +1,4 @@
 #ifdef YOU_SHOULDNT_READ_THIS
 #error YOU_SHOULDNT_BE_HERE
-enum class BRACKET_METHOD { standard, simple, arakawa, ctu, arakawa_old };
+enum class BRACKET_METHOD { standard, simple, arakawa, ctu };
 #endif
diff --git a/tools/pylib/_boutpp_build/scan_enums.py b/tools/pylib/_boutpp_build/scan_enums.py
index 837380c828..a403a92625 100644
--- a/tools/pylib/_boutpp_build/scan_enums.py
+++ b/tools/pylib/_boutpp_build/scan_enums.py
@@ -78,8 +78,6 @@ def __str__(self):
     "ARAKAWA": "arakawa",
     "BRACKET_CTU": "ctu",
     "CTU": "ctu",
-    "BRACKET_ARAKAWA_OLD": "arakawa_old",
-    "ARAKAWA_OLD": "arakawa_old",
 }
 enums["CELL_LOC"].extra = {
     "CELL_DEFAULT": "deflt",