diff --git a/src/underworld3/utilities/boundary_flux.py b/src/underworld3/utilities/boundary_flux.py
index c0b8e5b2..9ac32c70 100644
--- a/src/underworld3/utilities/boundary_flux.py
+++ b/src/underworld3/utilities/boundary_flux.py
@@ -633,10 +633,38 @@ def vcoord(q): return cvec[csec.getOffset(q) // dim]
     return np.array([sig[gi[_key(x, dim)]] for x in xs])
 
 
+def _refuse_a_constrained_boundary(solver, boundary):
+    """CBF recovery reads the VELOCITY residual. On a boundary held by
+    ``add_constraint_bc`` the traction has been moved into the multiplier block, so
+    what remains in the velocity rows is a CONSTANT traction — the reaction comes back
+    exactly proportional to the nodal boundary mass, and de-smearing hands that constant
+    straight back (#614).
+
+    That is not an approximation, it is the wrong quantity: measured on SolCx the
+    recovered flux had a standard deviation of 2e-15 across the boundary while the
+    equivalent ``Stokes`` solve varied correctly and correlated -0.997 with the exact
+    topography. It returns a plausible-looking number, which is why it went unnoticed.
+
+    Refuse rather than mislead. The traction on such a boundary is the multiplier
+    (plus its augmented-Lagrangian share — #607)."""
+    constraints = getattr(solver, "_block_constraint_bcs", None)
+    if not constraints:
+        return
+    if any(getattr(c, "boundary", None) == boundary for c in constraints):
+        raise NotImplementedError(
+            f"boundary_flux cannot recover the traction on {boundary!r}: it is held "
+            "by add_constraint_bc, so the traction is carried by the multiplier and "
+            "the velocity residual this reads holds only a constant (#614). Use "
+            f"solver.multiplier({boundary!r}) or solver.topography({boundary!r}) "
+            "instead."
+        )
+
+
 def boundary_flux(solver, boundary, mass="auto", remove_mean=False, normal=None):
     """See ``SolverBaseClass.boundary_flux``. Returns ``(xs, flux)`` for this rank's
     boundary nodes; scalar solver → normal flux, vector solver → traction (or its normal
     component if ``normal`` is given)."""
+    _refuse_a_constrained_boundary(solver, boundary)
     dm = solver.dm; dim = solver.mesh.dim
     ra = np.asarray(solver._assemble_volume_reaction()).ravel()
     nodes, lsec, csec, cvec, v0, v1, edge_nodes = _boundary_field_nodes(
diff --git a/tests/test_1019_boundary_flux.py b/tests/test_1019_boundary_flux.py
index 86e67faa..805f8675 100644
--- a/tests/test_1019_boundary_flux.py
+++ b/tests/test_1019_boundary_flux.py
@@ -11,6 +11,7 @@
 import sympy
 import pytest
 import underworld3 as uw
+from underworld3 import analytic as A
 
 pytestmark = [pytest.mark.level_1, pytest.mark.tier_a]
 
@@ -316,3 +317,45 @@ def test_volume_residual_fields_insert_essential_values():
             f"{name}: compute_volume_residual_fields diverges from "
             f"_assemble_volume_reaction on g != 0 walls — essential values not inserted"
         )
+
+
+@pytest.mark.level_2
+def test_cbf_refuses_a_multiplier_constrained_boundary():
+    """CBF must not silently return a constant on a constrained boundary (#614).
+
+    The back-calculation reads the VELOCITY residual. Where `add_constraint_bc`
+    holds the boundary, the traction has been moved into the multiplier block and
+    the velocity rows retain only a constant one — the reaction comes back exactly
+    proportional to the nodal boundary mass (ratio 4:2:1 at midpoint / interior
+    vertex / end vertex, the P2 line lumping), so de-smearing returns that constant.
+    Measured std across the boundary was 2e-15 where the equivalent `Stokes` solve
+    varied correctly. The danger is that the value LOOKS reasonable.
+
+    The unconstrained boundary of the same solver must keep working, or this is a
+    guard on the solver rather than on the boundary.
+    """
+    mesh = uw.meshing.StructuredQuadBox(
+        elementRes=(16, 16), minCoords=(0, 0), maxCoords=(1, 1), qdegree=3)
+    sol = A.SolCx(mesh, eta_A=1.0, eta_B=1.0e3, x_c=0.5, n=1)
+    solver = uw.systems.Stokes_Constrained(mesh)
+    solver.constitutive_model = uw.constitutive_models.ViscousFlowModel
+    solver.constitutive_model.Parameters.shear_viscosity_0 = sol.fn_viscosity
+    solver.bodyforce = sol.fn_bodyforce
+    solver.tolerance = 1.0e-9
+    solver.petsc_use_pressure_nullspace = True
+    solver.add_constraint_bc(0.0, "Top")
+    solver.add_dirichlet_bc((None, 0.0), "Bottom")
+    solver.add_dirichlet_bc((0.0, None), "Left")
+    solver.add_dirichlet_bc((0.0, None), "Right")
+    solver.solve()
+
+    with pytest.raises(NotImplementedError, match="add_constraint_bc"):
+        solver.boundary_flux("Top", normal=[[0.0, 1.0]])
+
+    # negative control: the Dirichlet boundary of the SAME solver still recovers,
+    # and recovers something with spatial content rather than a constant.
+    xs, flux = solver.boundary_flux("Bottom", normal=[[0.0, -1.0]])
+    assert len(xs) > 0
+    assert np.asarray(flux).std() > 1.0e-3, (
+        "the unconstrained boundary came back constant too — the guard is on the "
+        "solver, not on the boundary")