diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx
index 8cbe26c54..fa08f6e62 100644
--- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx
+++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx
@@ -4374,7 +4374,7 @@ class SNES_Vector(SolverBaseClass):
 
 
     def add_nitsche_bc(self, conds=None, boundary=None, direction=None,
-                       normal=None, gamma=10.0, theta=1, mask=None,
+                       normal=None, gamma=12.5, theta=1, mask=None,
                        local_h=True, g=None):
         r"""Add Nitsche weak enforcement of a velocity constraint along a direction.
 
@@ -4396,8 +4396,31 @@ class SNES_Vector(SolverBaseClass):
             terms — the same geometric-normal override as on the Stokes
             variant. Default ``None`` uses the per-boundary,
             deformation-tracking ``mesh.boundary_normal(boundary)``.
-        gamma : float, default=10.0
-            Dimensionless stabilisation parameter.
+        gamma : float, default=12.5
+            Dimensionless stabilisation parameter. The penalty is
+            ``gamma*mu/h``, so this is calibrated against the definition of
+            ``h``. It was 10.0 while ``h`` came from a kd-tree of neighbouring
+            centroids; ``mesh.cell_size()`` is now PETSc's ``volume**(1/dim)``
+            (#694), and 12.5 is calibrated against THAT definition on the Zhong
+            spherical shell — the benchmark whose 0.2% response drifted to
+            2.4-5.7% when ``h`` was last redefined without recalibrating
+            (#734).
+
+            The shift in ``h`` is not one number: it changes sign with the
+            dimension. Measured as new/old per cell,
+
+              2-D simplex box (unstructured)   +12.2%
+              2-D simplex box (regular)         +6.1%   (closed form: 6.07%)
+              2-D annulus                      +10.0%
+              3-D simplex box                  -29.8%
+              3-D spherical shell              -33.3%
+
+            so ``h`` grows by about a tenth on 2-D triangles and SHRINKS by
+            about a third on tetrahedra. Since the penalty is ``gamma*mu/h``, no
+            single gamma can reproduce the old enforcement in both. 12.5 is the
+            3-D number; **a 2-D sweep against an independent benchmark has not
+            been done**, and if one is wanted it belongs with #734 rather than
+            in this docstring.
         theta : {-1, 0, 1}, default=1
             Symmetry parameter (1=symmetric, -1=skew-symmetric).
         mask : sympy expression, optional
@@ -6591,7 +6614,7 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
                      remove_mean=remove_mean)
 
     def add_nitsche_bc(self, conds=None, boundary=None, direction=None, normal=None,
-                       gamma=10.0, theta=1, mask=None, local_h=True, g=None):
+                       gamma=12.5, theta=1, mask=None, local_h=True, g=None):
         r"""Add Nitsche weak enforcement of a velocity constraint along a direction.
 
         Nitsche's method provides a variationally consistent alternative to
@@ -6628,9 +6651,13 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
             Boundary unit normal used in the Nitsche consistency, symmetry,
             and pressure-coupling terms. Default ``None`` uses the per-boundary,
             deformation-tracking ``mesh.boundary_normal(boundary)``.
-        gamma : float, default=10.0
+        gamma : float, default=12.5
             Dimensionless stabilisation parameter. Typical values 5--20
-            for P2 elements.
+            for P2 elements. The penalty is ``gamma*mu/h``, so this is
+            calibrated against the definition of ``h``: it was 10.0 while
+            ``h`` came from a kd-tree of neighbouring centroids, and moved
+            with ``mesh.cell_size()`` becoming PETSc's ``volume**(1/dim)``
+            (#694).
         theta : {-1, 0, 1}, default=1
             Symmetry parameter:
              1: symmetric (default — optimal convergence and solver efficiency)
@@ -6650,13 +6677,15 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
             on every facet. Set ``False`` to restore the legacy global-h
             behaviour exactly.
 
-            The two coincide on **tensor** cells only. On a uniform **simplex**
-            mesh they differ by exactly :math:`\sqrt{2}` — for congruent
-            right-isosceles cells of legs :math:`h`, :meth:`Mesh.cell_size` is
-            :math:`2h/3` while :meth:`Mesh.get_min_radius` is
-            :math:`\sqrt{2}h/3` — so the penalty :math:`\gamma\mu/h` differs
-            between the two settings on the simplex meshes the free-slip and
-            fault models use. See ``tests/test_0010_cell_size_geometry.py``.
+            Since #694 both read the same quantity, PETSc's
+            :math:`\mathrm{volume}^{1/d}`, so this flag is now a choice between
+            the **local** cell and the **global minimum** and nothing else: on a
+            uniform mesh the two coincide exactly, for simplices as well as
+            tensor cells. They did not before — ``cell_size`` was a vertex-RMS
+            about the centroid and differed from ``get_min_radius`` by
+            :math:`\sqrt{2}` on simplices — so the flag silently rescaled the
+            penalty :math:`\gamma\mu/h` by cell type, which is how #734
+            happened. See ``tests/test_0010_cell_size_geometry.py``.
         g : sympy expression or float, optional
             Deprecated keyword alias for ``conds`` (one DeprecationWarning).
 
@@ -6755,12 +6784,13 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
         # local_h=False to restore the legacy single global-minimum scalar
         # (mesh.get_min_radius()).
         #
-        # The two coincide on TENSOR cells only. On a uniform SIMPLEX mesh --
-        # which is what the free-slip and fault models are built on -- they
-        # differ by exactly sqrt(2): on congruent right-isosceles cells of legs
-        # h, cell_size is 2h/3 and get_min_radius is sqrt(2)h/3. The penalty
-        # gamma*mu/h moves with that, so the two settings are NOT interchangeable
-        # there (see #734 and tests/test_0010_cell_size_geometry.py).
+        # Since #694 both read PETSc's volume**(1/dim), so this is a choice
+        # between the LOCAL cell and the GLOBAL minimum and nothing else -- on a
+        # uniform mesh they coincide exactly, simplices included. Before #694
+        # cell_size was a vertex-RMS about the centroid and differed from
+        # get_min_radius by sqrt(2) on simplices, so the flag silently rescaled
+        # gamma*mu/h by cell type (see #734 and
+        # tests/test_0010_cell_size_geometry.py).
         if local_h:
             h_sym = mesh.cell_size()
         else:
diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py
index 7f045997e..229e930e6 100644
--- a/src/underworld3/discretisation/discretisation_mesh.py
+++ b/src/underworld3/discretisation/discretisation_mesh.py
@@ -2840,12 +2840,10 @@ def nuke_coords_and_rebuild(
                 flush=True,
             )
 
-        (
-            self._min_size,
-            self._radii,
-            self._centroids,
-            self._search_lengths,
-        ) = self._get_mesh_sizes()
+        # `_min_size` and `_search_lengths` used to be unpacked here and were
+        # never read anywhere in src/ or tests/ -- the kd-tree loop computed
+        # three distance statistics per cell and two were discarded.
+        self._cell_radii, self._centroids = self._get_cell_radii()
 
         # Skip self-copy when hierarchy is trivial (issue #96 investigation)
         if self.dm is not self.dm_hierarchy[-1]:
@@ -3234,8 +3232,8 @@ def cell_size(self):
 
         Returns the ``.sym`` of a cell-constant (degree-0, discontinuous)
         scalar MeshVariable holding each cell's characteristic length (the
-        RMS distance of its vertices from their own centroid). This is a
-        purely cell-local quantity, independent of the MPI partition. Unlike
+        ``volume**(1/dim)`` equivalent radius, i.e. ``self._cell_radii``,
+        which comes from PETSc and is independent of the MPI partition). Unlike
         the single *global* scalar from :meth:`get_min_radius` (the smallest
         cell anywhere), this varies cell to cell, so a stabilisation that
         scales as :math:`1/h` — e.g. the Nitsche free-slip penalty
@@ -3301,27 +3299,45 @@ def _refresh():
     def _assemble_cell_size(self, var):
         """Fill ``var`` (degree-0 scalar) with each cell's characteristic size.
 
-        Uses the cell-geometry characteristic lengths ``self._cell_radii`` computed by
-        :meth:`_get_mesh_sizes` on the *current* geometry. A degree-0
+        Uses the per-cell characteristic lengths ``self._cell_radii`` computed
+        by :meth:`_get_cell_radii` on the *current* geometry. A degree-0
         discontinuous variable's local DOFs and ``self._cell_radii`` are BOTH
         indexed by this rank's cell-stratum order, so a direct assignment is
         correct on every rank.
 
-        This is deliberately a purely RANK-LOCAL operation (no ``var.coords``
-        access, no collective): mixing a rank-local fast path with a
-        collective fallback would diverge across ranks and deadlock, because
-        ``var.coords`` triggers the collective ``_get_coords_for_basis``."""
-        # Own-cell radii fix #687 without changing the legacy kd-tree radii
-        # used by global timestep estimates, adaptivity, and mesh relaxation.
+        This routine has TWO collectives in it and therefore no early return,
+        which is the opposite of what its previous docstring claimed (#698):
+
+        * the first ``var.data`` access lazily reaches ``MeshVariable._set_vec``,
+          which calls ``dm.createSubDM`` and ``createGlobalVector``;
+        * assigning into ``var.data`` fires the array's write-back callback,
+          ``pack_raw_data_to_petsc``.
+
+        A rank owning no cells used to short-circuit on ``radii.size == 0`` and
+        return before either, while its populated peers made both. Measured on a
+        region submesh at np=8 with cells per rank ``[12, 11, 0, 19, 0, 0, 0, 0]``:
+        the populated ranks sat in ``pack_raw_data_to_petsc`` and the job never
+        finished. A starved rank now walks the same path writing a zero-length
+        slice.
+
+        ``var.coords`` is still deliberately not read: it triggers the collective
+        ``_get_coords_for_basis``, and reading it only on some ranks would put a
+        third conditional collective back in."""
+        # `_cell_radii` is PETSc's volume**(1/dim), a property of each cell, so
+        # the values here do not depend on the partition -- and neither does the
+        # Nitsche penalty gamma*mu/h that consumes them under the default
+        # local_h=True. It was a kd-tree distance to the nearest centroid among
+        # THIS RANK's centroids, which near a seam could simply be absent (#694).
+        # There is NO early return here, and that is the point. Both steps below
+        # are collective, so a rank owning no cells has to walk through them
+        # writing nothing rather than skipping them (#698).
+        data = var.data                                   # allocates: createSubDM + createGlobalVector
         radii = numpy.asarray(self._cell_radii).reshape(-1)
-        # Empty partition (no local cells): nothing to fill on this rank.
-        if radii.size == 0 or var.data.shape[0] == 0:
-            return
-        # Assign over the common length. In practice these match exactly (same
-        # local cell set / ordering); the slice only guards a stray off-by-ghost
-        # mismatch without ever taking a collective path on a subset of ranks.
-        n = min(var.data.shape[0], radii.shape[0])
-        var.data[:n, 0] = radii[:n]
+
+        # `n` is 0 on a starved rank. The assignment still fires the array's
+        # write-back callback, which is the second collective.
+        n = min(data.shape[0], radii.shape[0])
+        data[:n, 0] = radii[:n]
 
     @property
     def Gamma_P1(self):
@@ -6934,56 +6950,44 @@ def _eval_use_robust_location(self) -> bool:
         """
         return (uw.mpi.size > 1) and (self._location_capability() != "none")
 
-    def _get_mesh_sizes(self, verbose=False):
-        """
-        Cache own-cell radii for cell_size and return legacy kd-tree radii.
-
-        Own-cell sizes use current DM vertices, so neither partition-local
-        neighbours nor stale coordinate views affect stabilization (#687).
-        Legacy radii remain unchanged for their other consumers.
-        """
-
+    def _get_cell_radii(self):
+        """Each cell's characteristic length, and the cell centroids.
+
+        The length is PETSc's ``volume**(1/dim)`` from
+        ``DMPlexComputeGeometryFVM``. A cell's volume is a property of that
+        cell, so this cannot depend on how the mesh was partitioned — which is
+        the point.
+
+        It replaces a kd-tree of THIS RANK's centroids queried with each cell's
+        vertices. Near a partition boundary the true nearest centroid can belong
+        to a cell owned by another rank and be absent from the tree, so the
+        answer moved with the rank count: per-cell by 3.3e-03 at np=2 and
+        4.1e-03 at np=4, `get_max_radius()` by 4.9% at np=8, `get_mean_radius()`
+        at every rank count, and `mesh.cell_size()` with them -- which scales
+        the Nitsche penalty under the DEFAULT ``local_h=True`` (#569, #687,
+        #694).
+
+        The FVM routine had been abandoned with a note that it "does not
+        compute all cells". That does not reproduce: measured on 2-D simplex,
+        2-D quad, 3-D tetrahedra, 3-D hexahedra and a deformed mesh, it returns
+        one finite positive value per local cell and is bit-identical across
+        rank counts in every case. (The note also named ``DMPlexGetMinRadius``,
+        which is a different call and is not used here.)
+        """
+        from underworld3.cython import petsc_discretisation
+
+        radii, _fvm_centroids = petsc_discretisation.petsc_fvm_get_local_cell_sizes(self)
+
+        # The FVM centroids are discarded: `_get_coords_for_basis(0, False)` is
+        # the degree-0 coordinate array the rest of the mesh indexes by cell,
+        # and mixing the two orderings would misalign every per-cell lookup.
         centroids = self._get_coords_for_basis(0, False)
-        centroids_kd_tree = uw.kdtree.KDTree(centroids)
 
-        import numpy as np
-
-        cStart, cEnd = self.dm.getHeightStratum(0)
-        pStart, pEnd = self.dm.getDepthStratum(0)
-        cell_length = np.empty(centroids.shape[0])
-        cell_min_r = np.empty(centroids.shape[0])
-        cell_r = np.empty(centroids.shape[0])
-        cell_radii = np.empty(centroids.shape[0])
-        coordinate_section = self.dm.getCoordinateDM().getLocalSection()
-        vertex_coordinates = self.dm.getCoordinatesLocal().array
-
-        for cell in range(cEnd - cStart):
-            cell_num_points = self.dm.getConeSize(cell)
-            cell_points = self.dm.getTransitiveClosure(cell)[0][-cell_num_points:]
-            # Use raw internal array for internal mesh operations (avoid unit-aware wrapping)
-            cell_coords = self._coords[cell_points - pStart]
-
-            distsq, _ = centroids_kd_tree.query(cell_coords, k=1, sqr_dists=True)
-
-            cell_length[cell] = np.sqrt(distsq.max())
-            cell_r[cell] = np.sqrt(distsq.mean())
-            cell_min_r[cell] = np.sqrt(distsq.min())
-
-            # A hex has six faces but eight vertices: select the vertex
-            # stratum, not a cone-sized suffix of its transitive closure.
-            closure = self.dm.getTransitiveClosure(cStart + cell)[0]
-            vertices = closure[(closure >= pStart) & (closure < pEnd)]
-            offsets = np.array([coordinate_section.getOffset(int(v)) for v in vertices])
-            own_coords = vertex_coordinates[offsets[:, None] + np.arange(self.cdim)]
-            delta = own_coords - own_coords.mean(axis=0)
-            cell_radii[cell] = np.sqrt(np.mean(np.sum(delta ** 2, axis=1)))
-
-        self._cell_radii = cell_radii
-        return cell_min_r, cell_r, centroids, cell_length
+        return radii, centroids
 
     # ==========
 
-    # Deprecated in favour of _get_mesh_sizes (above)
+    # Deprecated in favour of _get_cell_radii (above)
     def _get_mesh_centroids(self):
         """
         Obtain and cache the (local) mesh centroids using underworld swarm technology.
@@ -7073,7 +7077,7 @@ def get_min_radius(self) -> float:
         import numpy as np
         from mpi4py import MPI
 
-        radii = np.asarray(self._radii).reshape(-1)
+        radii = np.asarray(self._cell_radii).reshape(-1)
         local_min = float(radii.min()) if radii.size else float("inf")
         if uw.mpi.size > 1:
             local_min = uw.mpi.comm.allreduce(local_min, op=MPI.MIN)
@@ -7095,7 +7099,7 @@ def get_max_radius(self) -> float:
         import numpy as np
         from mpi4py import MPI
 
-        radii = np.asarray(self._radii).reshape(-1)
+        radii = np.asarray(self._cell_radii).reshape(-1)
         local_max = float(radii.max()) if radii.size else float("-inf")
         if uw.mpi.size > 1:
             local_max = uw.mpi.comm.allreduce(local_max, op=MPI.MAX)
@@ -7114,15 +7118,22 @@ def get_mean_radius(self) -> float:
         this is the canonical "mesh length" API. Use this anywhere you
         need a representative h0 (smoothing-length defaults, diffusion-
         stability heuristics, problem-scale normalisation) rather than
-        reaching for the rank-local ``self._radii`` array, which gives
-        different answers on different MPI ranks and leaks downstream
-        (e.g. into JIT C source via per-rank pointwise-function inputs).
+        reducing a per-rank array by hand, which gives different answers on
+        different MPI ranks and leaks downstream (e.g. into JIT C source via
+        per-rank pointwise-function inputs).
+
+        The value is the same at every RANK COUNT as well as on every rank:
+        ``self._cell_radii`` is PETSc's ``volume**(1/dim)``, a property of each
+        cell rather than of the partition. That was not true while these
+        reduced over a kd-tree of this rank's centroids -- an allreduce made
+        the answer agree across ranks without making it agree across rank
+        counts, and ``get_max_radius()`` moved 4.9% at np=8 (#694).
         """
 
         import numpy as np
         from mpi4py import MPI
 
-        radii = np.asarray(self._radii)
+        radii = np.asarray(self._cell_radii)
         local_sum = float(radii.sum())
         local_n = int(radii.size)
         if uw.mpi.size > 1:
diff --git a/src/underworld3/meshing/smoothing/api.py b/src/underworld3/meshing/smoothing/api.py
index 145e7977d..5b4160085 100644
--- a/src/underworld3/meshing/smoothing/api.py
+++ b/src/underworld3/meshing/smoothing/api.py
@@ -713,7 +713,7 @@ def follow_metric(
             mesh, T,
             refinement=2.0, coarsening=2.0,
             metric="gradient-uniform",
-            gradient_smoothing_length=2.0 * mesh._radii.mean(),
+            gradient_smoothing_length=2.0 * mesh.get_mean_radius(),
         )
 
     See Also
diff --git a/src/underworld3/swarm.py b/src/underworld3/swarm.py
index 7fc16b643..60590fe37 100644
--- a/src/underworld3/swarm.py
+++ b/src/underworld3/swarm.py
@@ -5908,7 +5908,7 @@ def estimate_dt(self, V_fn):
         vel = uw.function.evaluate(V_fn, self._particle_coordinates.data, evalf=True)
 
         # If vel is unit-aware (UnitAwareArray), nondimensionalise it to get
-        # consistent nondimensional values that match mesh._radii
+        # consistent nondimensional values that match mesh._cell_radii
         # Note: .magnitude returns physical units, which would be wrong here
         if hasattr(vel, "units") and vel.units is not None:
             vel = uw.non_dimensionalise(vel)
diff --git a/src/underworld3/systems/free_surface.py b/src/underworld3/systems/free_surface.py
index 6fd131772..dc6f803c9 100644
--- a/src/underworld3/systems/free_surface.py
+++ b/src/underworld3/systems/free_surface.py
@@ -121,7 +121,7 @@ def __init__(
         driving_buoyancy=None,
         smooth_length=0.0,
         mass="lumped",
-        max_surface_cfl=0.5,
+        max_surface_cfl=0.437,
         tangent_advect=None,
         tangent_spectral_modes=0,
         surface_mask=None,
diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py
index a9f59bbd7..c7fb18cde 100644
--- a/src/underworld3/systems/solvers.py
+++ b/src/underworld3/systems/solvers.py
@@ -256,7 +256,7 @@ def _global_max_diffusivity(constitutive_K, mesh):
         diffusivity = K
 
     # If unit-aware (UnitAwareArray), nondimensionalise so the value is
-    # consistent with mesh._radii. Note: .magnitude alone would keep the
+    # consistent with mesh._cell_radii. Note: .magnitude alone would keep the
     # physical-units number, which would be wrong here.
     if hasattr(diffusivity, "units") and diffusivity.units is not None:
         diffusivity = uw.non_dimensionalise(diffusivity)
@@ -276,7 +276,7 @@ def _centroid_velocities_nd(V_fn, mesh, basis=None, ensure_2d=True):
 
     Shared by the ``estimate_dt`` implementations: the advective CFL limit
     needs per-element centroid velocities in the same (nondimensional)
-    scale as ``mesh._radii``.
+    scale as ``mesh._cell_radii``.
 
     Parameters
     ----------
@@ -303,7 +303,7 @@ def _centroid_velocities_nd(V_fn, mesh, basis=None, ensure_2d=True):
         vel = uw.function.evaluate(V_fn, mesh._centroids)
 
     # If unit-aware (UnitAwareArray), nondimensionalise so the values are
-    # consistent with mesh._radii. Note: .magnitude alone would keep the
+    # consistent with mesh._cell_radii. Note: .magnitude alone would keep the
     # physical-units numbers, which would be wrong here.
     if hasattr(vel, "units") and vel.units is not None:
         vel = uw.non_dimensionalise(vel)
@@ -361,7 +361,7 @@ def _advective_diffusive_dt(constitutive_K, V_fn, mesh, direction_aware=False,
     diffusivity_glob = _global_max_diffusivity(constitutive_K, mesh)
     vel = _centroid_velocities_nd(V_fn, mesh)
     vel_magnitudes = np.linalg.norm(vel, axis=1)
-    element_radii = mesh._radii
+    element_radii = mesh._cell_radii
 
     def _reduce_dt(per_elem):
         fin = per_elem[np.isfinite(per_elem)] if len(per_elem) else per_elem
@@ -2422,7 +2422,7 @@ def estimate_dt(self):
         vel_magnitudes = np.linalg.norm(vel, axis=1)
 
         # Get per-element radii (characteristic element size)
-        element_radii = self.mesh._radii
+        element_radii = self.mesh._cell_radii
 
         # Compute per-element advective timestep: dt_i = h_i / |v_i|
         # Avoid division by zero for elements with zero velocity
@@ -4459,7 +4459,7 @@ def estimate_dt(self, direction_aware: bool = False, percentile: float = 0.0):
             centroid) · v̂` over the cell vertices. This is the
             distance material actually traverses through the cell
             per unit ``|v|``, and is **always ≥ the isotropic
-            mesh._radii estimate**, by 1.5–3× for equant cells
+            mesh._cell_radii estimate**, by 1.5–3× for equant cells
             (geometric factor) and up to ~10× for cells that the
             mover has stretched along the flow direction. On
             adapted meshes the gain is substantial; on uniform
diff --git a/tests/parallel/test_0774_empty_rank_reductions_mpi.py b/tests/parallel/test_0774_empty_rank_reductions_mpi.py
index e763ddb8d..f609ca752 100644
--- a/tests/parallel/test_0774_empty_rank_reductions_mpi.py
+++ b/tests/parallel/test_0774_empty_rank_reductions_mpi.py
@@ -1,7 +1,7 @@
 """Parallel regression tests for issue #405 — reductions on a zero-cell rank.
 
 A rank that owns NO CELLS used to raise a rank-local ``ValueError`` from an
-unguarded local reduction (``self._radii.min()`` and friends) while its
+unguarded local reduction (``self._cell_radii.min()`` and friends) while its
 populated peers sat in the matching collective. The job then hung or aborted
 asymmetrically. Every global quantity computed from rank-local data must
 instead reduce across ranks, with the starved rank contributing the identity
@@ -34,10 +34,11 @@
     pytest.mark.tier_a,
 ]
 
-# Cell is 1.0 wide x 0.5 tall; the characteristic length UW3 reports is the
-# centroid-to-corner half-diagonal. Analytic, so this is a real oracle rather
-# than a recorded number.
-SERIAL_RADIUS = math.sqrt(0.5 ** 2 + 0.25 ** 2)
+# Cell is 1.0 wide x 0.5 tall, so its area is 0.5 and the characteristic
+# length UW3 reports is area**(1/2). Analytic, so this is a real oracle rather
+# than a recorded number. (It was the centroid-to-corner half-diagonal before
+# the cell size came from PETSc's DMPlexComputeGeometryFVM, issue #694.)
+SERIAL_RADIUS = (1.0 * 0.5) ** (1.0 / 2.0)
 
 
 def _starved_box():
@@ -89,7 +90,7 @@ def test_negative_control_rank_local_minimum_would_be_caught():
     """
     mesh = _starved_box()
 
-    radii = np.asarray(mesh._radii).reshape(-1)
+    radii = np.asarray(mesh._cell_radii).reshape(-1)
     rank_local_min = float(radii.min()) if radii.size else float("inf")
     gathered = uw.mpi.comm.allgather(rank_local_min)
 
diff --git a/tests/parallel/test_0779_submesh_cell_size_starved_rank.py b/tests/parallel/test_0779_submesh_cell_size_starved_rank.py
new file mode 100644
index 000000000..80358dcb5
--- /dev/null
+++ b/tests/parallel/test_0779_submesh_cell_size_starved_rank.py
@@ -0,0 +1,69 @@
+"""``cell_size()`` must survive a rank that owns none of the submesh.
+
+A region submesh keeps only the cells carrying its label, so a partition whose
+share of the parent lies outside that region legitimately owns nothing. That is
+routine for submeshes and is the case this guards -- not an over-decomposed
+full mesh, which is a different (and degenerate) situation.
+
+The defect (#698): ``_assemble_cell_size`` returned early on ``radii.size == 0``,
+in front of TWO collectives --- the first ``var.data`` access, which allocates
+through ``dm.createSubDM`` and ``createGlobalVector``, and the assignment into
+it, which fires ``pack_raw_data_to_petsc``. Starved ranks skipped both while
+their populated peers made them, and the job never finished. Measured at np=8
+with cells per rank ``[12, 11, 0, 19, 0, 0, 0, 0]``.
+"""
+
+import pytest
+
+import underworld3 as uw
+
+pytestmark = [pytest.mark.mpi(min_size=2), pytest.mark.level_1, pytest.mark.tier_a]
+
+
+def _thin_slab_parent():
+    """A box split near the top, so the outer region is a few cells deep."""
+    return uw.meshing.BoxInternalBoundary(
+        elementRes=(12, 12), zelementRes=(10, 2), zintCoord=0.85,
+        minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), simplex=True, qdegree=2)
+
+
+@pytest.mark.mpi(min_size=2)
+def test_cell_size_completes_when_ranks_own_none_of_the_submesh():
+    """The whole test is that this returns. A regression hangs rather than fails.
+
+    The starved-rank count is reported, not asserted: how the parent is split is
+    the partitioner's business, and at np=2 the slab may well reach both ranks.
+    What must hold at every rank count is that every rank comes back.
+    """
+    submesh = _thin_slab_parent().extract_region("Outer")
+    start, end = submesh.dm.getHeightStratum(0)
+    counts = uw.mpi.comm.allgather(end - start)
+
+    submesh.cell_size()
+    uw.mpi.comm.barrier()
+
+    uw.mpi.pprint(f"SUBMESH_CELL_SIZE ranks={uw.mpi.size} cells={counts} "
+                  f"starved={counts.count(0)}")
+    assert sum(counts) > 0, "the region submesh is empty everywhere; test is vacuous"
+
+
+@pytest.mark.mpi(min_size=2)
+def test_the_size_field_is_filled_where_there_are_cells():
+    """A starved rank writing nothing must not stop the others writing.
+
+    Without this, the fix above could be satisfied by never filling the field
+    at all.
+    """
+    import numpy as np
+
+    submesh = _thin_slab_parent().extract_region("Outer")
+    start, end = submesh.dm.getHeightStratum(0)
+    submesh.cell_size()
+
+    local = np.asarray(submesh._cell_size_variable.array[:, 0, 0])
+    assert local.shape[0] == end - start
+    if local.size:
+        assert np.all(local > 0.0), "populated rank has non-positive cell sizes"
+    from mpi4py import MPI
+    filled = uw.mpi.comm.allreduce(int(local.size), op=MPI.SUM)
+    assert filled > 0, "no rank filled the field at all"
diff --git a/tests/parallel/test_1069_boundary_normal_parallel.py b/tests/parallel/test_1069_boundary_normal_parallel.py
index 036f39a45..152b25101 100644
--- a/tests/parallel/test_1069_boundary_normal_parallel.py
+++ b/tests/parallel/test_1069_boundary_normal_parallel.py
@@ -333,12 +333,14 @@ def _nitsche_annulus_diagnostics():
     leakage. Both are stable from tolerance 1e-9 to 1e-12, so neither is the linear
     solve.
 
-    This test now leaves ``local_h`` at its default ``True``. Before #569/#687,
-    doing so mixed the boundary-normal regression with a second partition-dependent
-    input from ``mesh.cell_size()``; this test therefore had to disable the public
-    default. The cell-local geometric size is now partition independent, so retaining
-    the default jointly guards the normal assembly and the Nitsche penalty path users
-    actually run.
+    This runs with the DEFAULT ``local_h=True``. It used to pass ``local_h=False``
+    because ``mesh.cell_size()`` was itself partition-dependent -- a kd-tree over
+    THIS RANK's cell centroids, so on this mesh the field's sum was 26.0822 at
+    np=1, 26.1211 at np=2 and 26.1386 at np=4 -- and leaving the default on would
+    have made this test measure two defects at once. ``cell_size()`` is now
+    PETSc's ``volume**(1/dim)`` and is partition independent (#569, #687, #694),
+    so keeping the default guards the normal assembly and the Nitsche penalty
+    path users actually run, together.
     """
     RI, RO = 0.5, 1.0
     mesh = uw.meshing.Annulus(radiusInner=RI, radiusOuter=RO, cellSize=0.12, qdegree=3)
diff --git a/tests/parallel/test_1077_cell_size_partition_independence.py b/tests/parallel/test_1079_cell_size_partition_independence.py
similarity index 100%
rename from tests/parallel/test_1077_cell_size_partition_independence.py
rename to tests/parallel/test_1079_cell_size_partition_independence.py
diff --git a/tests/parallel/test_1080_radius_accessors_partition_independence.py b/tests/parallel/test_1080_radius_accessors_partition_independence.py
new file mode 100644
index 000000000..fdcaf5f6b
--- /dev/null
+++ b/tests/parallel/test_1080_radius_accessors_partition_independence.py
@@ -0,0 +1,119 @@
+"""The radius ACCESSORS must not depend on how the mesh was partitioned.
+
+``test_1079`` covers the per-cell field cell by cell. This file covers
+``get_min_radius()``, ``get_max_radius()`` and ``get_mean_radius()``, which
+reduce that field and advertise a global mesh length -- and which were the half
+left partition-dependent when only ``cell_size()`` was fixed: ``get_max_radius()``
+moved 4.9% at np=8 and ``get_mean_radius()`` at every rank count.
+
+The defect (#569, #687, #694): the per-cell length came from a kd-tree over
+THIS RANK's centroids, queried with each cell's vertices. Near a partition
+boundary the true nearest centroid can belong to a cell owned by another rank
+and so be absent from the tree, and the answer moved with the rank count --
+per-cell by 3.3e-03 at np=2, ``get_max_radius()`` by 4.9% at np=8,
+``get_mean_radius()`` at every rank count.
+
+It reached users through ``mesh.cell_size()``, which scales the Nitsche penalty
+under the default ``local_h=True``, and through the three radius accessors,
+whose docstrings advertise a global mesh length.
+
+Two things this file is careful about, both learned the hard way:
+
+* **It compares two rank counts.** Partition independence is a statement about
+  two runs agreeing, so nothing measured at a single rank count establishes it.
+  The tests that shipped with the first attempt at this fix asserted a
+  within-rank oracle -- each cell matching its own vertices -- which is true of
+  a partition-dependent field as well.
+* **The reference is computed here, not recorded.** ``serial_reference`` runs
+  this module's own ``__main__`` at np=1 in this environment and asserts the
+  mesh fingerprints match, so a host that triangulates differently is reported
+  as that rather than as partition dependence.
+"""
+
+import numpy as np
+import pytest
+
+import underworld3 as uw
+from mpi4py import MPI
+
+from serial_reference import compare, emit, mesh_fingerprint, serial_reference
+
+pytestmark = [pytest.mark.mpi(min_size=2), pytest.mark.level_1, pytest.mark.tier_a]
+
+LABELS = ("min radius", "max radius", "mean radius", "sum of cell radii")
+
+# min and max are exact reductions of identical per-cell values, so they must
+# agree to the bit. The mean is a distributed sum and reduces in partition
+# order, so it is allowed the last couple of bits -- that ordering difference
+# is not the defect under test.
+RTOLS = (0.0, 0.0, 1.0e-12, 1.0e-12)
+
+
+def _box():
+    mesh = uw.meshing.UnstructuredSimplexBox(
+        minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.12, qdegree=2)
+    # mesh_fingerprint integrates over the mesh, and integration needs at
+    # least one variable to exist on it.
+    uw.discretisation.MeshVariable("cell_size_probe", mesh, 1, degree=1)
+    return mesh
+
+
+def _cell_size_diagnostics():
+    """The three accessors, plus a global sum that no single one of them sees.
+
+    The sum is the sharp one: a per-cell change that leaves the extremes alone
+    still moves it, and the old field's per-cell values differed while its
+    minimum happened not to.
+    """
+    mesh = _box()
+    radii = np.asarray(mesh._cell_radii).reshape(-1)
+    local_sum = float(radii.sum())
+    global_sum = uw.mpi.comm.allreduce(local_sum, op=MPI.SUM) \
+        if uw.mpi.size > 1 else local_sum
+
+    values = (
+        mesh.get_min_radius(),
+        mesh.get_max_radius(),
+        mesh.get_mean_radius(),
+        global_sum,
+    )
+    return values, mesh_fingerprint(mesh)
+
+
+@pytest.mark.mpi(min_size=2)
+def test_cell_size_is_the_same_at_every_rank_count():
+    """The accessors and the per-cell sum reproduce their own np=1 answer."""
+    values, fingerprint = _cell_size_diagnostics()
+    compare(values, serial_reference(__file__, "cell_size"),
+            rtols=RTOLS, labels=LABELS, fingerprint=fingerprint,
+            what="cell size / radius accessors")
+
+
+@pytest.mark.mpi(min_size=2)
+def test_every_rank_agrees_on_the_accessors():
+    """The weaker property, kept because it is the one the allreduce gives.
+
+    An allreduce makes an answer identical on every rank; it does not make it
+    identical at every rank count if the values being reduced are themselves
+    partition-dependent. Failing this while passing the test above would mean
+    the reduction is broken rather than its input, so keeping both separates
+    the two.
+    """
+    mesh = _box()
+    for name in ("min", "max", "mean"):
+        got = getattr(mesh, f"get_{name}_radius")()
+        everyone = uw.mpi.comm.allgather(got)
+        assert len(set(everyone)) == 1, (
+            f"get_{name}_radius() differs between ranks: {everyone}"
+        )
+
+
+if __name__ == "__main__":
+    import sys
+
+    _kind = sys.argv[1] if len(sys.argv) > 1 else "cell_size"
+    if _kind == "cell_size":
+        _values, _fingerprint = _cell_size_diagnostics()
+        emit(_values, _fingerprint)
+    else:
+        raise SystemExit(f"unknown kind {_kind!r}")
diff --git a/tests/test_0010_cell_size_geometry.py b/tests/test_0010_cell_size_geometry.py
index cd2de0d08..2a685c8b9 100644
--- a/tests/test_0010_cell_size_geometry.py
+++ b/tests/test_0010_cell_size_geometry.py
@@ -1,7 +1,14 @@
-"""Issue #687: cell_size is an own-cell geometric quantity, including after deform.
+"""``mesh.cell_size()`` is each cell's ``volume**(1/dim)``, and it tracks deformation.
 
-The independent oracle reads vertex coordinates through the coordinate section;
-it does not use the mesh's cached radii or centroid kd-tree. Run serial and MPI.
+The oracle here is computed from the vertex coordinates, not from PETSc, so this
+is a check and not a restatement of the implementation. For simplices that is
+the determinant volume of the cell; for the structured boxes it is the analytic
+cell volume, scaled by the determinant of the affine map when the mesh is
+deformed.
+
+Partition independence is NOT tested here -- it cannot be, at one rank count.
+``tests/parallel/test_1079`` compares the field cell by cell against its own
+serial answer, and ``test_1080`` does the same for the three radius accessors.
 """
 
 import numpy as np
@@ -9,81 +16,113 @@
 
 import underworld3 as uw
 
-pytestmark = [pytest.mark.level_1, pytest.mark.tier_b]
+pytestmark = [pytest.mark.level_1, pytest.mark.tier_a]
+
+# The deformation applied below, and its Jacobian determinant. Volumes scale by
+# |det|, so lengths scale by |det|**(1/dim).
+AFFINE = np.array([[1.7, 0.2], [0.0, 1.0]])
+AFFINE_DET = 1.7
+
+
+def _simplex_volume_from_vertices(mesh):
+    """Each simplex's volume from its own vertex coordinates.
+
+    Triangle: ``|det[v1-v0, v2-v0]| / 2``. Tetrahedron: ``|det[...]| / 6``.
+    """
+    dim = mesh.dim
+    cell_start, cell_end = mesh.dm.getHeightStratum(0)
+    point_start, _point_end = mesh.dm.getDepthStratum(0)
+
+    volumes = np.empty(cell_end - cell_start)
+    factorial = 2.0 if dim == 2 else 6.0
+    for cell in range(cell_end - cell_start):
+        corners = mesh.dm.getTransitiveClosure(cell)[0][-(dim + 1):]
+        coords = mesh._coords[corners - point_start]
+        edges = coords[1:] - coords[0]
+        volumes[cell] = abs(np.linalg.det(edges)) / factorial
+    return volumes
+
 
+@pytest.mark.parametrize("dim", [2, 3])
+def test_simplex_cell_size_is_the_cube_root_of_its_own_volume(dim):
+    """Against a determinant volume computed from the cell's vertices."""
+    mesh = uw.meshing.UnstructuredSimplexBox(
+        minCoords=(0.0,) * dim, maxCoords=(1.0,) * dim, qdegree=3,
+        cellSize=0.25, regular=False)
+    mesh.cell_size()
+    actual = np.asarray(mesh._cell_size_variable.array[:, 0, 0])
+    expected = _simplex_volume_from_vertices(mesh) ** (1.0 / dim)
 
-def _vertex_rms(mesh):
-    dm = mesh.dm
-    section = dm.getCoordinateDM().getLocalSection()
-    coordinates = dm.getCoordinatesLocal().array
-    start, end = dm.getHeightStratum(0)
-    first_vertex, last_vertex = dm.getDepthStratum(0)
-    radii = []
-    for cell in range(start, end):
-        vertices = [int(point) for point in dm.getTransitiveClosure(cell)[0]
-                    if first_vertex <= point < last_vertex]
-        points = np.array([coordinates[section.getOffset(v):section.getOffset(v) + mesh.cdim]
-                           for v in vertices])
-        radii.append(np.sqrt(np.mean(np.sum((points - points.mean(axis=0)) ** 2, axis=1))))
-    return np.asarray(radii)
+    error = float(np.abs(actual - expected).max(initial=0.0))
+    uw.pprint(f"CELL_SIZE_GEOMETRY simplex dim={dim} max_error={error:.12g}")
+    assert error < 1.0e-12, error
 
 
 @pytest.mark.parametrize("dim", [2, 3])
-@pytest.mark.parametrize("simplex", [True, False], ids=["simplex", "tensor"])
-def test_cell_size_matches_own_vertices_and_tracks_deform(dim, simplex):
-    geometry = dict(minCoords=(0.0,) * dim, maxCoords=(1.0,) * dim, qdegree=3)
-    mesh = (uw.meshing.UnstructuredSimplexBox(**geometry, cellSize=0.25, regular=False)
-            if simplex else uw.meshing.StructuredQuadBox(**geometry, elementRes=(4,) * dim))
+def test_structured_cell_size_is_the_analytic_value_and_tracks_deformation(dim):
+    """A regular box has an exact answer, and an affine deform scales it.
+
+    ``elementRes=4`` on the unit box gives cells of side 0.25, so
+    ``volume**(1/dim)`` is 0.25 whatever the dimension. Applying a linear map
+    multiplies every cell volume by ``|det|``, hence every length by
+    ``|det|**(1/dim)`` -- which the field must follow after ``deform``.
+    """
+    mesh = uw.meshing.StructuredQuadBox(
+        minCoords=(0.0,) * dim, maxCoords=(1.0,) * dim, qdegree=3,
+        elementRes=(4,) * dim)
     mesh.cell_size()
     field = mesh._cell_size_variable
-    errors = []
-    for phase in ("initial", "deformed"):
-        if phase == "deformed":
-            coordinates = np.array(mesh.X.coords)
-            coordinates[:, 0] = 1.7 * coordinates[:, 0] + 0.2 * coordinates[:, 1]
-            mesh.deform(coordinates)
-        expected = _vertex_rms(mesh)
-        actual = np.asarray(field.array[:, 0, 0])
-        shapes_match = actual.shape == expected.shape
-        assert all(uw.mpi.comm.allgather(shapes_match)), (actual.shape, expected.shape)
-        local_error = float(np.abs(actual - expected).max(initial=0.0))
-        error = max(uw.mpi.comm.allgather(local_error))
-        errors.append(error)
-        uw.pprint(f"CELL_SIZE_GEOMETRY dim={dim} simplex={simplex} phase={phase} "
-                  f"ranks={uw.mpi.size} max_error={error:.12g}")
-    assert max(errors) < 1e-12, errors
-
-
-def test_regular_square_cell_size_keeps_global_radius():
+
+    actual = np.asarray(field.array[:, 0, 0])
+    error = float(np.abs(actual - 0.25).max(initial=0.0))
+    assert max(uw.mpi.comm.allgather(error)) < 1.0e-12, error
+
+    coordinates = np.array(mesh.X.coords)
+    coordinates[:, 0] = AFFINE[0, 0] * coordinates[:, 0] + AFFINE[0, 1] * coordinates[:, 1]
+    mesh.deform(coordinates)
+
+    expected = 0.25 * AFFINE_DET ** (1.0 / dim)
+    actual = np.asarray(field.array[:, 0, 0])
+    error = float(np.abs(actual - expected).max(initial=0.0))
+    uw.pprint(f"CELL_SIZE_GEOMETRY tensor dim={dim} deformed max_error={error:.12g}")
+    assert max(uw.mpi.comm.allgather(error)) < 1.0e-12, (error, expected)
+
+
+def test_the_global_minimum_agrees_with_the_field():
+    """``get_min_radius()`` reduces the same field ``cell_size()`` exposes.
+
+    On a regular box every cell is the same size, so the global minimum is that
+    size -- which also pins the reduction to the analytic value rather than to
+    whatever the field happens to hold.
+    """
     mesh = uw.meshing.StructuredQuadBox(elementRes=(4, 4), qdegree=2)
-    legacy = np.array(mesh._radii)
-    global_radius = mesh.get_min_radius()
     mesh.cell_size()
-    expected = np.sqrt(2.0) / 8.0
-    error = float(np.abs(np.asarray(mesh._cell_size_variable.array) - expected).max(initial=0.0))
-    assert max(uw.mpi.comm.allgather(error)) < 1e-12
-    assert global_radius == pytest.approx(expected, rel=1e-12)
-    assert all(uw.mpi.comm.allgather(np.array_equal(mesh._radii, legacy)))
+    assert mesh.get_min_radius() == pytest.approx(0.25, rel=1.0e-12)
+    assert float(np.asarray(mesh._cell_size_variable.array).max()) == pytest.approx(
+        0.25, rel=1.0e-12)
 
 
 @pytest.mark.parametrize("h", [0.5, 0.25])
 def test_regular_simplex_cell_size_is_the_closed_form(h):
-    """cell_size on congruent right-isosceles cells is 2h/3, exactly.
+    """cell_size on congruent right-isosceles cells is h/sqrt(2), exactly.
 
     `test_cell_size_matches_own_vertices_and_tracks_deform` above checks the
     implementation against an independent reading of the same DEFINITION, so it
     stays true if the definition itself is changed on both sides. This pins the
     VALUE against geometry instead: `regular=True` tiles the box with congruent
-    right-isosceles triangles of legs h, whose vertices sit at (0,0), (h,0),
-    (0,h) up to rigid motion, so the RMS distance to the centroid is
+    right-isosceles triangles of legs h, whose area is h^2/2, and #694 defines
+    the characteristic length as PETSc's volume**(1/dim):
 
-        sqrt( ( 2(h/3)^2 + 2[(2h/3)^2 + (h/3)^2] ) / 3 ) = 2h/3
+        cell_size = (h^2 / 2)^(1/2) = h / sqrt(2)
 
     A redefinition of `cell_size` fails here immediately, at the quantity that
     changed. #692 changed it for simplices and nothing failed, so the Nitsche
     penalty (gamma*mu/h) moved 43% unnoticed and surfaced months later as a
-    5.67% miss on a spherical-shell benchmark (#734) — on one platform's
+    5.67% miss on a spherical-shell benchmark (#734) - on one platform's
     triangulation only, which is the hardest kind of failure to read backwards.
+    That is why this is written as geometry rather than as whatever the code
+    returns: the previous vertex-RMS definition gave 2h/3 here, and the number
+    moving from 2h/3 to h/sqrt(2) is exactly the event this test announces.
     """
     mesh = uw.meshing.UnstructuredSimplexBox(
         minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0),
@@ -92,27 +131,29 @@ def test_regular_simplex_cell_size_is_the_closed_form(h):
     mesh.cell_size()
     radii = np.asarray(mesh._cell_size_variable.array[:, 0, 0])
 
-    error = float(np.abs(radii - 2.0 * h / 3.0).max(initial=0.0))
+    expected = h / np.sqrt(2.0)
+    error = float(np.abs(radii - expected).max(initial=0.0))
     assert max(uw.mpi.comm.allgather(error)) < 1e-12, (
         f"cell_size on congruent legs-{h} right-isosceles cells is "
-        f"{radii.min():.12g}..{radii.max():.12g}, expected exactly {2.0 * h / 3.0:.12g}"
+        f"{radii.min():.12g}..{radii.max():.12g}, expected exactly {expected:.12g}"
     )
 
 
-def test_cell_size_and_min_radius_agree_on_tensor_cells_but_not_simplices():
-    """The two mesh-size measures coincide on TENSOR cells only, by sqrt(2).
+def test_cell_size_and_min_radius_are_now_one_measure():
+    """``cell_size()`` and ``get_min_radius()`` read the SAME quantity (#694).
 
-    `add_nitsche_bc` offers `local_h=False` to fall back from `mesh.cell_size()`
-    to the global `mesh.get_min_radius()`, and its docstring used to say the two
-    coincide "on a uniform mesh". They do on a regular quad box — which is
-    presumably where that was checked — and they do NOT on a uniform SIMPLEX
-    mesh, which is what every free-slip and fault model is built on.
+    ``add_nitsche_bc`` offers ``local_h=False`` to fall back from
+    ``mesh.cell_size()`` to the global ``mesh.get_min_radius()``. Before #694
+    those were different DEFINITIONS - vertex-RMS about the centroid versus
+    PETSc's face-distance radius - so the switch silently rescaled the penalty
+    ``gamma*mu/h`` by a factor that depended on the cell type: 1 on tensor cells
+    and sqrt(2) on simplices. That is how #734 happened.
 
-    Both measures have closed forms on congruent right-isosceles cells of legs h:
-    cell_size is 2h/3 (vertex RMS about the centroid) and get_min_radius is
-    sqrt(2)h/3, so the ratio is exactly sqrt(2). Pinning it means neither measure
-    can be redefined without this saying so, and says which way the Nitsche
-    penalty moves when it is.
+    #694 makes both read ``volume**(1/dim)``, so ``local_h`` is now a choice
+    between the LOCAL cell and the GLOBAL minimum and nothing else. On a uniform
+    mesh, where those two coincide, the measures must therefore agree EXACTLY -
+    on simplices as well as on tensor cells, which is the half that used not to
+    hold. Pinning both means neither can drift back apart silently.
     """
     quad = uw.meshing.StructuredQuadBox(elementRes=(4, 4), qdegree=2)
     quad.cell_size()
@@ -129,9 +170,11 @@ def test_cell_size_and_min_radius_agree_on_tensor_cells_but_not_simplices():
     simplex.cell_size()
     simplex_local = np.asarray(simplex._cell_size_variable.array[:, 0, 0])
 
-    assert simplex.get_min_radius() == pytest.approx(np.sqrt(2.0) * h / 3.0, rel=1e-12)
+    # both sides against the closed form, so a common drift cannot cancel
+    assert simplex.get_min_radius() == pytest.approx(h / np.sqrt(2.0), rel=1e-12)
     ratio = float(simplex_local.min()) / simplex.get_min_radius()
-    assert ratio == pytest.approx(np.sqrt(2.0), rel=1e-12), (
+    assert ratio == pytest.approx(1.0, rel=1e-12), (
         f"cell_size/get_min_radius on uniform simplices is {ratio:.12g}, "
-        "expected sqrt(2) — the Nitsche penalty gamma*mu/h scales with this"
+        "expected 1 - before #694 this was sqrt(2), and the Nitsche penalty "
+        "gamma*mu/h scaled with it"
     )
diff --git a/tests/test_1060_nitsche_freeslip.py b/tests/test_1060_nitsche_freeslip.py
index 2f849d635..641208293 100644
--- a/tests/test_1060_nitsche_freeslip.py
+++ b/tests/test_1060_nitsche_freeslip.py
@@ -67,8 +67,8 @@ def _solve_freeslip_box(method, res=8):
         stokes.add_natural_bc(1e4 * Gamma.dot(v.sym) * Gamma, "Top")
         stokes.add_natural_bc(1e4 * Gamma.dot(v.sym) * Gamma, "Bottom")
     elif method == "nitsche":
-        stokes.add_nitsche_bc(0.0, "Top", gamma=10.0)
-        stokes.add_nitsche_bc(0.0, "Bottom", gamma=10.0)
+        stokes.add_nitsche_bc(0.0, "Top")
+        stokes.add_nitsche_bc(0.0, "Bottom")
     else:
         raise ValueError(f"Unknown method: {method}")
 
diff --git a/tests/test_1065_nitsche_local_h.py b/tests/test_1065_nitsche_local_h.py
index e6c225611..b643b46b4 100644
--- a/tests/test_1065_nitsche_local_h.py
+++ b/tests/test_1065_nitsche_local_h.py
@@ -152,7 +152,7 @@ def _box_wobble(X0, amp):
 # --------------------------------------------------------------------------
 def test_cell_size_is_local_per_cell():
     """``mesh.cell_size()`` is a per-cell field equal to each cell's
-    own-vertex RMS size (``mesh._cell_radii``), not the single global minimum."""
+    per-cell size (``mesh._cell_radii``, PETSc's ``volume**(1/dim)``), not the single global minimum."""
     mesh = _graded_box()
     h = mesh.cell_size()  # sympy symbol -> backed by a P0 field
     field = np.asarray(mesh._cell_size_variable.array[:, 0, 0]).reshape(-1)
@@ -225,7 +225,7 @@ def test_cell_size_tracks_deformation():
 # --------------------------------------------------------------------------
 # 3. local-h still solves free-slip correctly (back-compat / correctness)
 # --------------------------------------------------------------------------
-def _solve_freeslip(mesh, method, gamma=10.0):
+def _solve_freeslip(mesh, method, gamma=12.5):
     v = uw.discretisation.MeshVariable(
         "U", mesh, mesh.dim, degree=2, vtype=uw.VarType.VECTOR)
     p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1)