diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py
index 25c85b16..56058ce1 100644
--- a/src/underworld3/systems/solvers.py
+++ b/src/underworld3/systems/solvers.py
@@ -2737,20 +2737,22 @@ def __init__(
         ### sets up DuDt and DFDt
         ## ._setup_history_terms()
 
-        # If DuDt is not provided, then we can build a SLCN version
+        # If DuDt is not provided, build it here. We prefer the Eulerian_DDt
+        # formulation over the SemiLagrangian alternative to keep the solver
+        # initialisation simpler and easier to maintain.
         if self.Unknowns.DuDt is None:
-            self.Unknowns.DuDt = uw.systems.ddt.SemiLagrangian(
+            self.Unknowns.DuDt = uw.systems.Eulerian_DDt(
                 self.mesh,
-                self.u.sym,  # Symbolic expression - SemiLagrangian evaluates this at each update
-                self.u.sym,
-                vtype=uw.VarType.VECTOR,
+                self.u,
+                vtype=self.u.vtype,
                 degree=self.u.degree,
                 continuous=self.u.continuous,
+                V_fn=self.u.sym,
                 varsymbol=self.u.symbol,
-                verbose=self.verbose,
                 bcs=self.essential_bcs,
                 order=self._order,
-                smoothing=0.0001,
+                smoothing=0.0,
+                verbose=self.verbose,
             )
 
         # F (at least for N-S) is a nodal point variable so there is no benefit
@@ -2760,19 +2762,33 @@ def __init__(
         # Maybe u.degree-1. The scalar equivalent seems to show
         # little benefit from specific choices here other than
         # discontinuous flux variables amplifying instabilities.
+        if self.Unknowns.DFDt is None:
+            # self.Unknowns.DFDt = uw.systems.ddt.SemiLagrangian(
+            #     self.mesh,
+            #     sympy.Matrix.zeros(self.mesh.dim, self.mesh.dim),
+            #     self.u.sym,
+            #     vtype=uw.VarType.SYM_TENSOR,
+            #     degree=self.u.degree,
+            #     continuous=self.u.continuous,
+            #     varsymbol=rf"{{ F[ {self.u.symbol} ] }}",
+            #     verbose=self.verbose,
+            #     bcs=None,
+            #     order=self._order,
+            # )
+            self.Unknowns.DFDt = uw.systems.ddt.Eulerian(
+                self.mesh,
+                sympy.Matrix.zeros(self.mesh.dim, self.mesh.dim),
+                vtype=uw.VarType.SYM_TENSOR,
+                degree=self.u.degree,
+                continuous=self.u.continuous,
+                varsymbol=rf"{{ F[ {self.u.symbol} ] }}",
+                bcs=None,
+                order=self._order,
+                smoothing=0.0,
+                verbose=self.verbose,
+            )
+
 
-        self.Unknowns.DFDt = uw.systems.ddt.SemiLagrangian(
-            self.mesh,
-            sympy.Matrix.zeros(self.mesh.dim, self.mesh.dim),
-            self.u.sym,
-            vtype=uw.VarType.SYM_TENSOR,
-            degree=self.u.degree,
-            continuous=self.u.continuous,
-            varsymbol=rf"{{ F[ {self.u.symbol} ] }}",
-            verbose=self.verbose,
-            bcs=None,
-            order=self._order,
-        )
 
         ## Add in the history terms provided ...