diff --git a/docs/advanced/benchmarks/_ti_vep_bdf2_isolation.py b/docs/advanced/benchmarks/_ti_vep_bdf2_isolation.py
index eacffedea..4c5157d65 100644
--- a/docs/advanced/benchmarks/_ti_vep_bdf2_isolation.py
+++ b/docs/advanced/benchmarks/_ti_vep_bdf2_isolation.py
@@ -66,6 +66,7 @@ def build(label, *, spatial_yield, yield_mode):
     cm.Parameters.shear_viscosity_1 = ETA_1
     cm.Parameters.shear_modulus = MU
     cm.Parameters.yield_stress = ty
+    cm.Parameters.yield_stress_min = -sympy.oo  # the fault ty = 1/weakness reads a field but stays >= TAU_Y > 0
     cm.Parameters.director = sympy.Matrix([0.0, 1.0])  # θ=0 throughout
     cm.Parameters.strainrate_inv_II_min = 1.0e-6
     cm._yield_mode = yield_mode
diff --git a/docs/advanced/benchmarks/bench_ti_vep_harmonic.py b/docs/advanced/benchmarks/bench_ti_vep_harmonic.py
index cd03c7323..d6484eb01 100644
--- a/docs/advanced/benchmarks/bench_ti_vep_harmonic.py
+++ b/docs/advanced/benchmarks/bench_ti_vep_harmonic.py
@@ -130,6 +130,7 @@ def build_ti_stokes(label, theta_deg, tau_y, bdf_order):
     cm.Parameters.shear_viscosity_1 = ETA_1
     cm.Parameters.shear_modulus = MU
     cm.Parameters.yield_stress = tau_y_field
+    cm.Parameters.yield_stress_min = -sympy.oo  # reads the distance field, but 1/weakness >= tau_y > 0 by construction
     cm.Parameters.director = director
     cm.Parameters.shear_viscosity_min = ETA_0 * 1.0e-3
     cm.Parameters.strainrate_inv_II_min = 1.0e-6
diff --git a/docs/advanced/fault-networks.md b/docs/advanced/fault-networks.md
index f87cd1774..1aa4b3004 100644
--- a/docs/advanced/fault-networks.md
+++ b/docs/advanced/fault-networks.md
@@ -91,6 +91,7 @@ yield. Read the next paragraph before using it as glue.
 ```python
 stokes.constitutive_model = uw.constitutive_models.ViscoPlasticFlowModel
 stokes.constitutive_model.Parameters.yield_stress = net.band_yield(4.0)
+stokes.constitutive_model.Parameters.yield_stress_min = -sympy.oo  # a field, but tau_y or tau_far: positive by construction
 stokes.consistent_jacobian = True
 ```
 
diff --git a/docs/advanced/listric-fault-extension.md b/docs/advanced/listric-fault-extension.md
index 790887390..781de2b1f 100644
--- a/docs/advanced/listric-fault-extension.md
+++ b/docs/advanced/listric-fault-extension.md
@@ -69,6 +69,7 @@ net.prepare(h=CELL, ligament=1.0)
 net.mesh = base_mesh().add_fault(net.prepared)   # same base as the ribbons
 ...
 stokes.constitutive_model.Parameters.yield_stress = net.damage_yield(v, dial=0.05)
+stokes.constitutive_model.Parameters.yield_stress_min = -sympy.oo  # dial*(1 + 2 edot_II) > 0 by construction
 net.apply_contact(stokes, conds=ETA_F)
 ```
 
diff --git a/docs/developer/guides/nonlinear-solver.md b/docs/developer/guides/nonlinear-solver.md
index 73b83a7ff..31a760c64 100644
--- a/docs/developer/guides/nonlinear-solver.md
+++ b/docs/developer/guides/nonlinear-solver.md
@@ -84,6 +84,7 @@ case felt like whack-a-mole. Check them first.
 | **A pressure-dependent yield stress makes the Newton Schur complement INDEFINITE far from the solution**: the consistent tangent differentiates η by p, so J_vp is no longer −Bᵀ (Spiegelman notch at η_bg 5e24: 41% off, 40 of 1968 Schur eigenvalues change sign). In the basin of convergence the full tangent is benign again (no wrong-sign eigenvalues, no negative pivots). | fieldsplit Schur diverges at ~0.99 per iteration even with the velocity block solved exactly; `DIVERGED_LINEAR_SOLVE` on the first Newton steps of a cold Drucker-Prager solve | `newton_pressure_coupling=False` (the Newton velocity block with the frozen pressure column, a partial tangent with the standard saddle-point structure) on the approach, or `consistent_jacobian="continuation"`; the default keeps the full tangent, which is right once inside the basin |
 | **Perfect plasticity's consistent tangent is SINGULAR along the flow**: on the hard-`Min` plastic branch η = τ_y/2ε̇_II, so 2η + 2η′ε̇_II = 0 — the velocity block is symmetric but semi-definite in every yielded cell. (An earlier version of this row blamed *asymmetry*; that is wrong for any η(ε̇_II) law — the rank-one term η′ ε̇⊗ε̇/ε̇_II is symmetric. Pressure-dependent yield adds a non-symmetric v–p coupling, not a non-symmetric velocity block. Corrected 2026-08-26, maintainer review.) | benign while yielded cells are few (the viscous neighbours regularise); with a large yielded fraction the velocity sub-solve caps out and Newton stalls at ~1e-3, no failure reason | give the plastic branch a positive tangent: the declared `Parameters.plastic_rate_strengthening` (η_reg adds 2η_reg along the flow in every yielded cell; its scale is the model's), a rounded viscosity floor (clips only the fastest cells). The default soft-min (sqrt, `"onset"`, δ = 0.1) does NOT do it — its yielded tangent is ≥ −1e-4 of the viscous one, still effectively singular; the sqrt law under `yield_anchor="yield"` is negative (−0.029). Picard converges regardless (full 2η stiffness) but is linear-rate. The FMG bundle's `gmres`+`sor` smoother is Newton-safe either way |
 | `preconditioner="fmg"` (vs explicit `pc_type=mg` + manual mg opts) | outer KSP "converges" in **1 iteration** → no real Newton correction → stall → `DIVERGED_LINE_SEARCH` | use explicit `pc_type=mg` with the smoother opts above; bound the outer KSP (`ksp_max_it`~80) so a hostile step fails fast |
+| A yield stress that reads a field (pressure, damage) left at the zero `yield_stress_min` | `ValueError: ... the yield stress reads ... residual strength` at setup (the cut-off would make η_pl = 0, a singular operator, #759) | state the residual strength (`yield_stress_min > 0`, stress units), or `-sympy.oo` if the yield stress is positive by construction |
 | Cold plastic start `v=0`, or any rigid/unyielded point | `DIVERGED_FNORM_NAN` at iteration 0 | **Not** a div/0: `ε̇=0` gives `η_pl=+inf`, which `Min` and the sqrt soft-min carry correctly to the viscous branch. Only a soft-min form that computes `η_ve·η_pl/(η_ve+η_pl)` breaks (`inf/inf`). Fixed in the power-mean; if you hand-roll a blend, write the harmonic mean as `η_ve/(1+η_ve/η_pl)`. **Do not reach for a strain-rate floor** — it hides this rather than fixing it |
 | LU velocity block with all-Dirichlet-ish BC | pressure nullspace singular | attach the Stokes nullspace / avoid a bare LU there |
 | Hand-rolled `snes_monitor` to "see what's happening" | you read residuals but miss the tell | use `solve_with_diagnostics` / `get_snes_diagnostics` instead (below) |
diff --git a/docs/developer/guides/plasticity-solvers.md b/docs/developer/guides/plasticity-solvers.md
index b32328b7d..2dc138be9 100644
--- a/docs/developer/guides/plasticity-solvers.md
+++ b/docs/developer/guides/plasticity-solvers.md
@@ -17,6 +17,7 @@ and records what was retired.
 ```python
 stokes.constitutive_model = cm           # ViscoPlastic / ViscoElasticPlastic / TI-VEP
 cm.Parameters.yield_stress = tau_y       # finite -> plasticity active
+cm.Parameters.yield_stress_min = tau_r   # residual strength: REQUIRED if tau_y reads a field (-sympy.oo if positive by construction)
 stokes.consistent_jacobian = True        # Newton tangent (non-elastic DP; see table)
 stokes.solve()                           # cold: ONE Picard step automatically, then Newton
 stokes.solve(picard=3)                   # the recipe: a few Picard steps, then Newton
@@ -235,6 +236,11 @@ could.
   plastic branch wherever `m * eta_reg >= shear_viscosity_min`; the floor remains a
   safeguard for the creep branch. `yield_stress_min` (the POSITIVE residual strength in
   tension) is a statement of the yield law and stays.
+- **`yield_stress_min` must be stated when the yield stress reads a field** (a
+  Drucker-Prager `C + sin(phi) p`, a cohesion weakening with damage). At the zero
+  default a yield stress cut off in tension gives `eta_pl = 0`, a singular operator
+  (`DIVERGED_LINEAR_SOLVE`, #759); the model refuses it. State the residual strength
+  in stress units, or `-sympy.oo` for a yield stress positive by construction.
 
 - **`shear_viscosity_min`** (default `-oo` = off) is applied through
   `uw.maths.smooth_max`, but the default rounding scale is zero under
@@ -258,6 +264,7 @@ could.
 |---------|-------|-----|
 | `DIVERGED_LINEAR_SOLVE`, 0 iters, VEP | consistent Newton over the elastic block → indefinite | Picard (`consistent_jacobian=False`) |
 | `DIVERGED_LINEAR_SOLVE` at nl=0 with a viscosity floor set | δ→0 leaves the floor's `Max` corner exact | set `viscosity_min_rounding` |
+| `ValueError: the yield stress depends on a field ... residual strength` | a field-dependent yield stress with the zero default floor (#759) | set `yield_stress_min` to the residual strength, or `-sympy.oo` if the yield stress cannot reach zero |
 | Newton stalls with no divergence reason | admissible uselessness — steps accepted, residual flat | revert to best iterate, Picard block (`solve(picard=N)` or `"continuation"`); consider grid sequencing |
 | Converges but σ sits **below** τ_y | a fixed δ>0 soft-min under the default `"onset"` anchor is a WEAKER law (4.5 % at the yield point at δ = 0.1) | that is the modelling choice you made — a smaller δ, or `yield_mode="min"` for the exact surface (`yield_anchor="yield"` reaches τ_y but stiffens every unyielded cell and softens above yield) |
 | Linear (~20-iter) convergence | Picard tangent when you wanted Newton | `consistent_jacobian=True` on a non-elastic model (see "Confirm" above) |
diff --git a/docs/examples/WIP/VEP_Fault_Investigation.py b/docs/examples/WIP/VEP_Fault_Investigation.py
index 3ea9dd169..f47309cd6 100644
--- a/docs/examples/WIP/VEP_Fault_Investigation.py
+++ b/docs/examples/WIP/VEP_Fault_Investigation.py
@@ -87,6 +87,7 @@
 cm.Parameters.shear_viscosity_0 = ETA
 cm.Parameters.shear_modulus = MU
 cm.Parameters.yield_stress = tau_y_field
+cm.Parameters.yield_stress_min = -sympy.oo  # reads the distance field, but 1/weakness >= TAU_Y_FAULT > 0 by construction
 # yield_mode="softmin" is the default (smooth approximation to Min)
 cm.Parameters.shear_viscosity_min = ETA * 1.0e-2
 cm.Parameters.strainrate_inv_II_min = 1.0e-5
diff --git a/docs/examples/solid_mechanics/advanced/Ex_Shear_Band_Plasticity_SS.py b/docs/examples/solid_mechanics/advanced/Ex_Shear_Band_Plasticity_SS.py
index c29fce27f..11e1fd664 100644
--- a/docs/examples/solid_mechanics/advanced/Ex_Shear_Band_Plasticity_SS.py
+++ b/docs/examples/solid_mechanics/advanced/Ex_Shear_Band_Plasticity_SS.py
@@ -182,7 +182,7 @@ class boundaries(Enum):
 stokes.constitutive_model = uw.constitutive_models.ViscoPlasticFlowModel
 stokes.constitutive_model.Parameters.shear_viscosity_0 = 1
 stokes.constitutive_model.Parameters.yield_stress = mu * p_soln.sym[0] + C
-stokes.constitutive_model.Parameters.yield_min = 1000
+stokes.constitutive_model.Parameters.yield_stress_min = 1000  # residual strength (= C): the Drucker-Prager yield reads p
 stokes.penalty = 1
 
 stokes.petsc_options["ksp_monitor"] = None
diff --git a/docs/examples/solid_mechanics/advanced/Ex_Sheared_Layer_Test.py b/docs/examples/solid_mechanics/advanced/Ex_Sheared_Layer_Test.py
index b4fc78fe7..931a7bfff 100644
--- a/docs/examples/solid_mechanics/advanced/Ex_Sheared_Layer_Test.py
+++ b/docs/examples/solid_mechanics/advanced/Ex_Sheared_Layer_Test.py
@@ -427,6 +427,7 @@ class boundaries(Enum):
 C = 2 * (y * 2) ** 2 + (C0 - Cinf) * (1 - sympy.tanh((strain.sym[0] / eps_ref - 1) * scale)) / 2 + Cinf
 
 stokes.constitutive_model.Parameters.yield_stress = C + mu * p_soln.sym[0]
+stokes.constitutive_model.Parameters.yield_stress_min = 0.1 * Cinf  # residual strength: 10% of the softened cohesion (yield reads p and strain)
 stokes.constitutive_model.Parameters.edot_II_fn = stokes.Unknowns.Einv2
 stokes.constitutive_model.Parameters.min_viscosity = 0.1
 stokes.saddle_preconditioner = 1 / stokes.constitutive_model.viscosity
diff --git a/docs/examples/solid_mechanics/intermediate/Ex_TI_VEP_Angled_Fault.py b/docs/examples/solid_mechanics/intermediate/Ex_TI_VEP_Angled_Fault.py
index 4c045ba0f..90d148663 100644
--- a/docs/examples/solid_mechanics/intermediate/Ex_TI_VEP_Angled_Fault.py
+++ b/docs/examples/solid_mechanics/intermediate/Ex_TI_VEP_Angled_Fault.py
@@ -165,6 +165,7 @@
 cm.Parameters.shear_viscosity_1 = ETA_1
 cm.Parameters.shear_modulus = MU
 cm.Parameters.yield_stress = tau_y_field
+cm.Parameters.yield_stress_min = -sympy.oo  # reads the distance field, but 1/weakness >= TAU_Y > 0 by construction
 cm.Parameters.director = director
 cm.Parameters.shear_viscosity_min = ETA_0 * 1.0e-3
 cm.Parameters.strainrate_inv_II_min = 1.0e-6
diff --git a/src/underworld3/constitutive_models.py b/src/underworld3/constitutive_models.py
index 0b44198e2..eac21f590 100644
--- a/src/underworld3/constitutive_models.py
+++ b/src/underworld3/constitutive_models.py
@@ -75,6 +75,72 @@ def _parameter_is_zero(par):
         return False
 
 
+def _fields_read_by(par):
+    """The names of the mesh or swarm variables a parameter reads (the pressure, a
+    damage field, ...): a value not known until the solve, which can change sign."""
+    from underworld3.function.expressions import unwrap_expression
+    try:
+        value = unwrap_expression(par, mode="nondimensional")
+    except (TypeError, ValueError, AttributeError):
+        # Sanctioned: a parameter that cannot be unwrapped is not inspected, and the
+        # residual-strength check then trusts it (fails open, as before #759)
+        return []
+    names = set()
+    for f in sympy.sympify(value).atoms(sympy.Function):
+        if hasattr(f, "meshvar"):
+            variable = f.meshvar() if callable(f.meshvar) else f.meshvar
+            names.add(getattr(variable, "name", None) or str(f.func))
+    return sorted(names)
+
+
+def _numeric_value(par):
+    """A parameter's non-dimensional value as a float, or None when it is not a
+    number (an expression of a field)."""
+    from underworld3.function.expressions import unwrap_expression
+    try:
+        value = unwrap_expression(par, mode="nondimensional")
+        return float(getattr(value, "magnitude", value))
+    except (TypeError, ValueError, AttributeError):
+        # Sanctioned: an expression has no single value; None says so
+        return None
+
+
+def _floor_is_off(floor):
+    """True when the yield-stress floor is disabled: -oo, bare or wrapped in an
+    expression (a wrapped -oo applied as a floor evaluates to NaN)."""
+    if floor.sym == -sympy.oo:
+        return True
+    value = _numeric_value(floor)
+    return value is not None and value == -float("inf")
+
+
+def _require_residual_strength(constitutive_model, parameters):
+    """Refuse a yield stress that reads a field (Drucker-Prager through the pressure,
+    a cohesion that weakens with damage) with a floor that is not positive: where it
+    reaches the floor the plastic viscosity tau_y / (2 edot_II) is zero or negative
+    and the operator singular (#759). The residual strength must be stated, > 0 in
+    stress units; an explicit -oo removes the floor for a yield stress known to stay
+    positive. A floor that is itself a field expression is trusted as stated."""
+    floor = parameters.yield_stress_min
+    if _floor_is_off(floor):
+        return
+    value = _numeric_value(floor)
+    if value is not None and value > 0.0:
+        return
+    if value is None and not _parameter_is_zero(floor):
+        return                                     # a floor that is itself a field
+    fields = _fields_read_by(parameters.yield_stress)
+    if fields:
+        material = getattr(constitutive_model, "_material_name", None)
+        where = type(constitutive_model).__name__ + (f" ({material})" if material else "")
+        raise ValueError(
+            f"{where}: the yield stress reads {', '.join(fields)} and can fall to the "
+            f"floor yield_stress_min = {floor.sym}, where the plastic viscosity is zero or "
+            "negative and the operator singular. Set Parameters.yield_stress_min to the "
+            "residual strength (> 0, in stress units), or to -sympy.oo if the yield "
+            "stress is known to stay positive.")
+
+
 def _plastic_rate_strengthening_parameter():
     """The declared plastic rate strengthening, shared by every yielding model's
     Parameters (ViscoPlastic, ViscoElasticPlastic, TransverseIsotropicVEP)."""
@@ -1379,9 +1445,11 @@ def yield_smoother(self):
         Which SIDE of ``Min`` a family sits on is set by ``yield_anchor``, not by the
         family — see that property.
 
-        Selecting ``"powermean"`` bumps a zero ``yield_softness`` to ``1.0``
-        (``s≈1``, the parameter-free harmonic mean) since ``δ=0`` is the sharp
-        hard-``Min`` limit it only *approaches*.
+        Selecting a family never changes ``yield_softness``: the power mean runs
+        at the δ the model states, ``s = 1/(δ + 0.001)`` -- ``δ = 0`` is
+        ``s = 1000``, close to the exact ``Min``; ``δ = 1`` is the harmonic mean,
+        a distinct, strongly smoothed law (#782: selecting it used to bump a zero
+        δ to 1 silently).
         """
         return getattr(self, "_yield_smoother", "sqrt")
 
@@ -1392,10 +1460,6 @@ def yield_smoother(self, value):
                 f"yield_smoother must be 'sqrt' or 'powermean', got '{value}'"
             )
         self._yield_smoother = value
-        if value == "powermean" and getattr(self, "_yield_softness", 0.0) == 0.0:
-            # δ=0 ⇒ s=1/δ=∞ is the singular Min limit; default to the
-            # parameter-free harmonic mean (s=1) instead.
-            self.yield_softness = 1.0
         self._reset()
 
     # --- Smooth lower bounds (viscosity / yield floors) -----------------------------
@@ -1706,7 +1770,7 @@ class _Parameters(_ParameterBase, _ViscousParameterAlias):
         yield_stress_min = api_tools.Parameter(
             R"{\tau_{y, \mathrm{min}}}",
             lambda inner_self: 0,
-            "Yield stress (DP) minimum cutoff",
+            "Residual strength: the floor of the yield stress (stated, > 0, when the yield stress reads a field)",
             units="Pa",
         )
 
@@ -1751,8 +1815,10 @@ def viscosity(self):
         # than propagated into tau_y/(2 edot_II). (The ±oo defaults elsewhere in
         # Parameters exist so sympy can cancel an unused term away; that trick is
         # wrong here, so this one defaults to 0 — maintainer ruling 2026-07-26.)
-        # An explicit -oo still disables the floor.
-        if inner_self.yield_stress_min.sym != -sympy.oo:
+        # An explicit -oo still disables the floor. A yield stress that reads a field
+        # needs a stated, positive floor (_require_residual_strength).
+        _require_residual_strength(self, inner_self)
+        if not _floor_is_off(inner_self.yield_stress_min):
             yield_stress = self._apply_floor(
                 inner_self.yield_stress, inner_self.yield_stress_min
             )
@@ -1812,7 +1878,7 @@ def _rate_strengthening_probe_terms(self):
         # the viscosity above, so the ladder reads the law the solve runs.
         inner_self = self.Parameters
         tau = inner_self.yield_stress
-        if inner_self.yield_stress_min.sym != -sympy.oo:
+        if not _floor_is_off(inner_self.yield_stress_min):
             tau = self._apply_floor(tau, inner_self.yield_stress_min)
         edot = self._strainrate_inv_II
         if not _parameter_is_zero(inner_self.strainrate_inv_II_min):
@@ -2189,7 +2255,7 @@ def dt_elastic(inner_self, value):
         yield_stress_min = api_tools.Parameter(
             R"{\tau_{y, \mathrm{min}}}",
             lambda inner_self: 0,
-            "Yield stress (DP) minimum cutoff",
+            "Residual strength: the floor of the yield stress (stated, > 0, when the yield stress reads a field)",
             units="Pa",
         )
 
@@ -2835,8 +2901,10 @@ def _plastic_effective_viscosity(self):
         # Guard on the DISABLING sentinel, not on zero: zero is the default and a
         # physically meaningful floor (the yield stress is compared against the second
         # invariant of the stress, so a negative tau_y is meaningless). Only an
-        # explicit -oo turns the floor off.
-        if parameters.yield_stress_min.sym != -sympy.oo:
+        # explicit -oo turns the floor off. A yield stress that reads a field needs a
+        # stated, positive floor (_require_residual_strength).
+        _require_residual_strength(self, parameters)
+        if not _floor_is_off(parameters.yield_stress_min):
             # Literal 0 for the default floor, not the parameter atom: sympy cannot
             # fuzzy-compare an opaque UWexpression and Max canonicalisation recurses.
             # smooth_max keeps both operands symbolic (the JIT routes the floor
@@ -3910,7 +3978,7 @@ def dt_elastic(inner_self, value):
         yield_stress_min = api_tools.Parameter(
             R"{\tau_{y, \mathrm{min}}}",
             lambda inner_self: 0,
-            "Yield stress minimum cutoff", units="Pa",
+            "Residual strength: the floor of the yield stress (stated, > 0, when the yield stress reads a field)", units="Pa",
         )
         strainrate_inv_II_min = api_tools.Parameter(
             R"{\dot\varepsilon_{II,\mathrm{min}}}",
@@ -4270,7 +4338,9 @@ def _plastic_effective_viscosity(self):
         tau_y = parameters.yield_stress
         # Guard on the DISABLING sentinel, not on zero — zero is the default and a
         # real floor (a negative yield stress is meaningless against an invariant).
-        if parameters.yield_stress_min.sym != -sympy.oo:
+        # A yield stress that reads a field needs a stated, positive floor.
+        _require_residual_strength(self, parameters)
+        if not _floor_is_off(parameters.yield_stress_min):
             # Literal 0 for the default floor (sympy fuzzy-compare recursion on atoms).
             # smooth_max: keeps the floor symbolic, no ordering test (see the note
             # in Constitutive_Model._apply_floor). eps = 0 is exactly Max.
@@ -4794,6 +4864,8 @@ def _eta_par_eff_lagged(self):
             except (TypeError, ValueError):
                 edot_min_val = 1.0e-6
 
+        # TODO(BUG): the lagged estimate uses the unfloored yield stress, so it ignores
+        # yield_stress_min and the residual-strength rule (#759); see _require_residual_strength.
         tau_y_sym = params.yield_stress.sym
         vp_eff_lag = tau_y_sym / (
             2 * (gamma_dot_abs_lag + sympy.Float(edot_min_val))
diff --git a/src/underworld3/function/expressions.py b/src/underworld3/function/expressions.py
index b119f98c2..15172ba87 100644
--- a/src/underworld3/function/expressions.py
+++ b/src/underworld3/function/expressions.py
@@ -706,6 +706,17 @@ def live_expressions():
     return [UWexpression._expr_names[k] for k in sorted(UWexpression._expr_names)]
 
 
+def _quantity_as_expression(quantity):
+    """Wrap a bare UWQuantity as its own UWexpression atom, so that a sum or
+    difference with a symbolic operand keeps it in the tree and unwrap
+    non-dimensionalises it by its own units rather than by its magnitude."""
+    if quantity._pint_qty is not None:
+        latex_name = f"{quantity._pint_qty:~L}"
+    else:
+        latex_name = str(quantity.value)
+    return UWexpression(latex_name, quantity, _unique_name_generation=True)
+
+
 class UWexpression(MathematicalMixin, uw_object, Symbol):
     """
     A SymPy Symbol that wraps a value for lazy evaluation.
@@ -1483,33 +1494,12 @@ def __add__(self, other):
                     result_qty,  # Pass full UWQuantity - Transparent Container stores it
                     _unique_name_generation=True,
                 )
-            elif isinstance(self._sym, UWQuantity):
-                # self has units, other doesn't - convert other to UWQuantity
-                other_qty = UWQuantity(other.value, self.units) if self.units else other.value
-                if isinstance(other_qty, UWQuantity):
-                    result_qty = self._sym + other_qty
-                    return UWexpression(
-                        f"({self.name}+{other.name})",
-                        result_qty,
-                        _unique_name_generation=True,
-                    )
-                else:
-                    result_value = self.value + other.value
-                    return UWexpression(
-                        f"({self.name}+{other.name})",
-                        result_value,
-                        _unique_name_generation=True,
-                        units=self.units
-                    )
-            else:
-                # Neither has UWQuantity in _sym - use simple value arithmetic
-                result_value = self.value + other.value
-                return UWexpression(
-                    f"({self.name}+{other.name})",
-                    result_value,
-                    _unique_name_generation=True,
-                    units=self.units or other.units
-                )
+            # Otherwise keep both symbols in the tree, as __sub__ does, so that
+            # unwrap non-dimensionalises each operand by its own units. Adding the
+            # operands' .value eagerly mixed systems: a quantity's .value is its
+            # dimensional magnitude (0.35 for 0.35 MPa), which then entered a
+            # non-dimensional sum unscaled (#758).
+            return Symbol.__add__(self, other)
 
         # Handle UWexpression + UWQuantity
         if isinstance(other, UWQuantity):
@@ -1522,25 +1512,8 @@ def __add__(self, other):
                     result_qty,
                     _unique_name_generation=True,
                 )
-            else:
-                # self.sym is not UWQuantity - convert and add
-                if self.units is not None:
-                    self_qty = UWQuantity(self.value, self.units)
-                    result_qty = self_qty + other
-                    return UWexpression(
-                        f"({self.name}+qty)",
-                        result_qty,
-                        _unique_name_generation=True,
-                    )
-                else:
-                    # No units on self - just add values
-                    result_value = self.value + other.value
-                    return UWexpression(
-                        f"({self.name}+qty)",
-                        result_value,
-                        _unique_name_generation=True,
-                        units=other.units
-                    )
+            # self is symbolic: keep the quantity as its own atom (#758), as above
+            return Symbol.__add__(self, _quantity_as_expression(other))
         return Symbol.__add__(self, other)
 
     def __radd__(self, other):
@@ -1575,18 +1548,7 @@ def __sub__(self, other):
 
         # Handle UWexpression - UWQuantity: Wrap in UWexpression first (LAZY)
         if isinstance(other, UWQuantity):
-            # Wrap the UWQuantity in a UWexpression to preserve it as a symbol
-            # This allows unwrap_for_evaluate to find and nondimensionalize it
-            if other._pint_qty is not None:
-                latex_name = f"{other._pint_qty:~L}"
-            else:
-                latex_name = str(other.value)
-            wrapped_other = UWexpression(
-                latex_name,
-                other,  # Store the full UWQuantity - Transparent Container
-                _unique_name_generation=True
-            )
-            return Symbol.__sub__(self, wrapped_other)
+            return Symbol.__sub__(self, _quantity_as_expression(other))
 
         return Symbol.__sub__(self, other)
 
@@ -1594,32 +1556,16 @@ def __rsub__(self, other):
         """Right subtraction - handle UWQuantity specially."""
         from .quantities import UWQuantity
         if isinstance(other, UWQuantity):
-            # UWQuantity - UWexpression → UWexpression
-            # Convert self to other's units (other is the "base" unit here)
-            self_units = self.units
-            other_units = other.units
-
-            if self_units is not None and other_units is not None:
-                # Convert self's value to other's units for correct subtraction
-                # other - self: result is in other's units
-                try:
-                    from ..scaling import units as ureg
-                    self_in_other_units = (self.value * self_units).to(other_units).magnitude
-                except Exception:
-                    self_in_other_units = self.value
-                result_value = other.value - self_in_other_units
-                result_units = other_units
-            else:
-                # Use self.value (not self.sym) for arithmetic
-                result_value = other.value - self.value
-                result_units = self_units or other_units
-
-            return UWexpression(
-                f"(qty-{self.name})",
-                result_value,
-                _unique_name_generation=True,
-                units=result_units
-            )
+            # Both are quantities: Pint handles the unit conversion
+            if isinstance(self._sym, UWQuantity):
+                return UWexpression(
+                    f"(qty-{self.name})",
+                    other - self._sym,
+                    _unique_name_generation=True,
+                )
+            # self is symbolic: keep the quantity as its own atom (#758), as in
+            # __add__ — converting self.value by Pint scaled a symbolic magnitude
+            return Symbol.__sub__(_quantity_as_expression(other), self)
         return Symbol.__rsub__(self, other)
 
     def __pow__(self, other):
diff --git a/src/underworld3/function/quantities.py b/src/underworld3/function/quantities.py
index bf29a1f78..0fdd595d9 100644
--- a/src/underworld3/function/quantities.py
+++ b/src/underworld3/function/quantities.py
@@ -385,31 +385,11 @@ def __sub__(self, other: Union['UWQuantity', float, int]) -> 'UWQuantity':
             else:
                 return UWQuantity(self._value - other._value, self._pint_unit or other._pint_unit)
 
-        # Handle UWexpression: UWQuantity - UWexpression → UWexpression
+        # Handle UWexpression: delegate to UWexpression.__rsub__, as __add__ does,
+        # which keeps a symbolic operand lazy so each side is scaled by its own
+        # units (#758)
         if isinstance(other, UWexpression):
-            from ..units import get_units
-            other_units = other.units
-
-            # For subtraction, units must be compatible - convert other to self's units
-            if self._pint_unit is not None and other_units is not None:
-                try:
-                    from ..scaling import units as ureg
-                    # Convert other's value to self's units
-                    other_converted = (other.value * other_units).to(self._pint_unit).magnitude
-                except Exception:
-                    # Units incompatible - use raw value (will be wrong but won't crash)
-                    other_converted = other.value
-            else:
-                other_converted = other.value
-
-            result_sym = self._value - other_converted
-            combined_units = self._pint_unit  # Result in self's units
-            return UWexpression(
-                f"(qty-{other.name})",
-                result_sym,
-                _unique_name_generation=True,
-                units=combined_units
-            )
+            return NotImplemented
 
         # Handle SymPy expressions - use LAZY EVALUATION like __add__
         if isinstance(other, sympy.Basic):
diff --git a/src/underworld3/meshing/fault_network.py b/src/underworld3/meshing/fault_network.py
index 2d27e6a8e..57137750f 100644
--- a/src/underworld3/meshing/fault_network.py
+++ b/src/underworld3/meshing/fault_network.py
@@ -55,6 +55,8 @@
 >>> stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0
 >>> stokes.constitutive_model.Parameters.yield_stress = \\
 ...     net.damage_yield(v, dial=0.05)
+>>> # positive by construction (dial + rate term): no floor needed
+>>> stokes.constitutive_model.Parameters.yield_stress_min = -sympy.oo
 >>> stokes.consistent_jacobian = True
 >>> net.apply(stokes)
 >>> # ... boundary conditions ...
@@ -1171,6 +1173,8 @@ def band_yield(self, tau_y, tau_far=1.0e8, tag=""):
 
             stokes.constitutive_model.Parameters.yield_stress = \
                 net.band_yield(tau_y=4.0)
+            # tau_y or tau_far everywhere: positive, no floor needed
+            stokes.constitutive_model.Parameters.yield_stress_min = -sympy.oo
 
         This is the damage the SPLIT realisation wants. A released fault
         flank sits far below ``tau_y`` and is untouched; the places that
diff --git a/tests/test_0106_a_sum_of_expressions_stays_live.py b/tests/test_0106_a_sum_of_expressions_stays_live.py
new file mode 100644
index 000000000..c2d770a07
--- /dev/null
+++ b/tests/test_0106_a_sum_of_expressions_stays_live.py
@@ -0,0 +1,37 @@
+"""A sum of two expressions keeps both as live atoms.
+
+The sum used to be formed eagerly from the operands' current values, so the
+compiled kernel froze the number: changing one operand afterwards changed nothing
+(and in a units model a quantity entered the sum unscaled, #758). Checked on a
+Poisson solve with diffusivity a + b: the solution is inversely proportional to
+the diffusivity, so changing a must rescale it exactly.
+"""
+import numpy as np
+import pytest
+
+import underworld3 as uw
+
+pytestmark = [pytest.mark.level_1, pytest.mark.tier_b]
+
+
+def test_changing_an_operand_of_a_sum_reaches_the_solve():
+    uw.reset_default_model()
+    mesh = uw.meshing.UnstructuredSimplexBox(minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.2)
+    u = uw.discretisation.MeshVariable("u_sum", mesh, 1, degree=2)
+    poisson = uw.systems.Poisson(mesh, u_Field=u)
+    a = uw.expression(r"a_{\mathrm{sum}}", 2.0, "first operand")
+    b = uw.expression(r"b_{\mathrm{sum}}", 1.0, "second operand")
+    poisson.constitutive_model = uw.constitutive_models.DiffusionModel
+    poisson.constitutive_model.Parameters.diffusivity = a + b
+    poisson.f = 1.0
+    poisson.tolerance = 1.0e-12
+    for boundary in ("Left", "Right", "Top", "Bottom"):
+        poisson.add_dirichlet_bc(0.0, boundary)
+    poisson.solve()
+    u3 = np.asarray(u.array).copy()
+    a.sym = 7.0                                     # diffusivity 3 -> 8
+    poisson.solve()
+    u8 = np.asarray(u.array).copy()
+    # a frozen sum leaves u unchanged (an error of order one); the live sum rescales
+    # it to the solver's tolerance
+    assert np.abs(u8 * 8.0 - u3 * 3.0).max() <= 1e-8 * np.abs(u3 * 3.0).max()
diff --git a/tests/test_0850_fault_network_toolkit.py b/tests/test_0850_fault_network_toolkit.py
index d3f6eb589..6277ea308 100644
--- a/tests/test_0850_fault_network_toolkit.py
+++ b/tests/test_0850_fault_network_toolkit.py
@@ -60,6 +60,7 @@ def test_network_end_to_end_transmits_slip():
     cm.yield_mode = "min"
     cm.Parameters.shear_viscosity_0 = 1.0
     cm.Parameters.yield_stress = net.damage_yield(v, dial=0.05)
+    cm.Parameters.yield_stress_min = -sympy.oo   # dial + rate term: positive by construction (#759)
     stokes.consistent_jacobian = True
     stokes.bodyforce = [0.0, 0.0]
     trend = float(np.degrees(np.arctan2(MAIN[1][1] - MAIN[0][1],
diff --git a/tests/test_0851_fault_network_3d.py b/tests/test_0851_fault_network_3d.py
index 7da2265d7..81a30f916 100644
--- a/tests/test_0851_fault_network_3d.py
+++ b/tests/test_0851_fault_network_3d.py
@@ -7,6 +7,7 @@
 """
 import numpy as np
 import pytest
+import sympy
 
 import underworld3 as uw
 from underworld3.meshing.fault_network_3d import (crossing_segment,
@@ -144,6 +145,7 @@ def test_network_3d_end_to_end():
     cm.yield_mode = "min"
     cm.Parameters.shear_viscosity_0 = 1.0
     cm.Parameters.yield_stress = net.damage_yield(v, dial=0.05)
+    cm.Parameters.yield_stress_min = -sympy.oo   # dial + rate term: positive by construction (#759)
     stokes.consistent_jacobian = True
     stokes.bodyforce = [0.0, 0.0, 0.0]
     for wall in ("Bottom", "Top", "Left", "Right", "Front", "Back"):
diff --git a/tests/test_1055_yield_smoother.py b/tests/test_1055_yield_smoother.py
index b17290cb8..48d114d63 100644
--- a/tests/test_1055_yield_smoother.py
+++ b/tests/test_1055_yield_smoother.py
@@ -93,10 +93,12 @@ def test_powermean_smoother_undershoots_min():
 
     cm.yield_mode = "softmin"              # dev: soft-min family lives under "softmin"
     cm.yield_anchor = "onset"              # the law under test; not left to the default
-    cm.yield_softness = 0.0                # start from δ=0 so the powermean select bumps it
-    cm.yield_smoother = "powermean"        # bumps δ 0 → 1 (s=1, harmonic mean)
+    cm.yield_softness = 0.0
+    cm.yield_smoother = "powermean"
     assert cm.yield_smoother == "powermean"
-    assert cm.yield_softness == 1.0
+    # selecting the family leaves the stated δ alone (#782: it used to bump 0 -> 1,
+    # the harmonic mean, without a word)
+    assert cm.yield_softness == 0.0
 
     # Includes a geodynamic-range pair (1e25, 1e21) that overflows the naive
     # η^(−s) form above s≈40 but is finite in the harmonic-normalised form.
@@ -200,10 +202,12 @@ def test_dp_model_smoother_optin():
     cm.yield_smoother = "powermean"
     assert cm.yield_softness == 0.1
     cm.yield_smoother = "sqrt"
-    # only the singular δ = 0 is bumped, to the power mean's entry value
+    # a family switch never changes δ, not even the singular δ = 0 (#782); the
+    # power mean's entry value, the harmonic mean, is stated explicitly
     cm.yield_softness = 0.0
-    cm.yield_smoother = "powermean"           # bumps δ→1
-    assert cm.yield_softness == 1.0
+    cm.yield_smoother = "powermean"
+    assert cm.yield_softness == 0.0
+    cm.yield_softness = 1.0
     for eta_ve, eta_pl in [(1.0, 2.0), (5.0, 0.3), (1e25, 1e21)]:
         comb = cm._combine_yield(sympy.Float(eta_ve), sympy.Float(eta_pl))
         val = float(unwrap_expression(comb, mode="nondimensional"))
@@ -239,6 +243,7 @@ def test_dp_pressure_yield_consistent_jacobian_builds():
     sphi = float(np.sin(np.deg2rad(30.0)))
     cm.Parameters.shear_viscosity_0 = 1.0e3
     cm.Parameters.yield_stress = sympy.Max(1.0 + sphi * p.sym[0], sympy.sympify(0))
+    cm.Parameters.yield_stress_min = 0.05   # a tension cut-off needs a stated residual strength (#759)
     cm.Parameters.strainrate_inv_II_min = 1.0e-8
     s.consistent_jacobian = True
     # Build the pointwise residual + Jacobian; on the bug this recurses to RecursionError.
diff --git a/tests/test_1057_yield_homotopy_solve.py b/tests/test_1057_yield_homotopy_solve.py
index a1da222db..c02f0729c 100644
--- a/tests/test_1057_yield_homotopy_solve.py
+++ b/tests/test_1057_yield_homotopy_solve.py
@@ -178,7 +178,7 @@ def test_yield_stress_is_floored_at_zero_by_default():
                          [("min", None), ("softmin", "powermean")])
 def test_pressure_dependent_yield_survives_tension(mode, smoother):
     """A Drucker-Prager yield that goes NEGATIVE in tension must not poison the
-    viscosity in any yield mode — the zero floor must catch it first."""
+    viscosity in any yield mode — the residual-strength floor must catch it first."""
     import numpy as np
 
     mesh, stokes = _viscoplastic_stokes(cellSize=0.5)
@@ -187,6 +187,8 @@ def test_pressure_dependent_yield_survives_tension(mode, smoother):
     # C + sin(phi)*p with the pressure driven strongly negative (tension) over part
     # of the domain, so the raw yield stress changes sign inside the mesh.
     cm.Parameters.yield_stress = 1.0 + 0.5 * stokes.p.sym[0]
+    # the residual strength left in tension, stated as a tension cut-off needs (#759)
+    cm.Parameters.yield_stress_min = 0.05
     cm.yield_mode = mode
     if smoother is not None:
         cm.yield_smoother = smoother
@@ -354,3 +356,47 @@ def test_homotopy_refuses_a_stress_history_solver():
 
     with pytest.raises(NotImplementedError, match="stress history"):
         stokes.solve(homotopy=True, timestep=0.1)
+
+
+def test_a_field_dependent_yield_needs_a_residual_strength():
+    """#759: a yield stress that reads a field (here Drucker-Prager through the
+    pressure) floored at zero gives eta_pl = tau_y/(2 edot_II) = 0 wherever it is cut
+    off in tension, a singular operator. It is refused until the residual strength is
+    stated; a stated strength bounds the plastic viscosity below by
+    tau_res/(2 edot_II) > 0; an explicit -oo is the opt-out for a yield stress known
+    to stay positive; a constant yield stress is unaffected."""
+    import numpy as np
+
+    mesh, stokes = _viscoplastic_stokes(cellSize=0.5)
+    cm = stokes.constitutive_model
+    cm.Parameters.yield_stress = 1.0 + 0.5 * stokes.p.sym[0]
+    with pytest.raises(ValueError, match="residual strength"):
+        cm.viscosity
+
+    cm.Parameters.yield_stress_min = -0.1         # a negative floor is no residual strength
+    with pytest.raises(ValueError, match="residual strength"):
+        cm.viscosity
+
+    # In tension everywhere (p = -4, so 1 + 0.5 p = -1) under a uniform shear
+    # (u = y: edot_II = 1/2), the stated strength sets eta = tau_res / (2 edot_II).
+    tau_res = 0.05
+    cm.Parameters.yield_stress_min = tau_res
+    cm.yield_mode = "min"
+    stokes.p.array[...] = -4.0
+    X = np.asarray(stokes.u.coords)
+    stokes.u.array[:, 0, 0] = X[:, 1]
+    stokes.u.array[:, 0, 1] = 0.0
+    eta = np.asarray(uw.function.evaluate(cm.viscosity.sym, mesh.X.coords)).ravel()
+    assert np.abs(eta - tau_res / (2 * 0.5)).max() < 1e-10, eta
+
+    cm.Parameters.yield_stress_min = -sympy.oo
+    cm.viscosity                                   # the opt-out is honoured
+    # the opt-out wrapped in an expression is the same opt-out, not a NaN floor
+    cm.Parameters.yield_stress_min = uw.expression(r"\tau_{\mathrm{off}}", -sympy.oo, "no floor")
+    cm.Parameters.yield_stress = 2.0 + 0.0 * stokes.p.sym[0]
+    eta = np.asarray(uw.function.evaluate(cm.viscosity.sym, mesh.X.coords)).ravel()
+    assert np.all(np.isfinite(eta)), eta
+
+    cm.Parameters.yield_stress_min = 0
+    cm.Parameters.yield_stress = 5.0
+    cm.viscosity                                   # a constant cannot fall to the floor
diff --git a/tests/test_1059_yield_stress_min_units.py b/tests/test_1059_yield_stress_min_units.py
new file mode 100644
index 000000000..40433d2d1
--- /dev/null
+++ b/tests/test_1059_yield_stress_min_units.py
@@ -0,0 +1,146 @@
+"""The yield-stress floor is scaled like the yield stress in a units model (#758).
+
+``ViscoPlasticFlowModel`` floors the yield stress, :math:`\\tau_y \\ge \\tau_{y,min}`,
+with ``uw.maths.smooth_max(tau_y, tau_y_min, 0)``, i.e.
+:math:`\\tfrac12(a + b + |a - b|)`. The floor can be written as a plain
+non-dimensional number, as a ``uw.quantity`` in any stress unit, or as a product of a
+dimensional expression and a number. All are the same physical stress, so all must
+give the same effective viscosity: the exact :math:`\\max(\\tau_y, \\tau_{y,min})`.
+
+Before the fix, ``UWexpression + UWexpression`` added the operands' dimensional
+magnitudes eagerly (``self.value + other.value``) whenever the left operand was not
+a bare quantity, so in :math:`a + b` the floor arrived as 0.35 (MPa) or 3.5e5 (Pa)
+in a non-dimensional sum while :math:`|a - b|` scaled it correctly. A floor 83 times
+below the yield stress raised it (quantity) or stopped yielding entirely (product).
+
+Baseline: pure shear with a linear velocity field, which the degree-2 velocity
+represents exactly, so :math:`\\dot\\varepsilon_{II} = R` and the yielded viscosity is
+:math:`\\eta = \\max(\\tau_y, \\tau_{y,min}) / (2R)` in closed form.
+"""
+
+import numpy as np
+import pytest
+import sympy
+
+import underworld3 as uw
+from underworld3.function.expressions import unwrap_expression
+
+pytestmark = [pytest.mark.level_2]
+
+H = 10.4e3                                  # m
+V = 44.496e-3 / (365.25 * 86400)            # m/s
+ETA = 5.8333e19                             # Pa s
+STRESS_SCALE = ETA * V / H                  # Pa
+
+FLOOR_PA = 0.35e6
+
+# every spelling of the same 0.35 MPa floor; C0 is the cohesion expression
+FLOOR_SPELLINGS = {
+    "nondimensional_float": lambda C0, tau_pa: FLOOR_PA / STRESS_SCALE,
+    "quantity_Pa": lambda C0, tau_pa: uw.quantity(FLOOR_PA, "Pa"),
+    "quantity_MPa": lambda C0, tau_pa: uw.quantity(FLOOR_PA / 1.0e6, "MPa"),
+    "product_with_expression": lambda C0, tau_pa: (FLOOR_PA / tau_pa) * C0,
+}
+
+
+def _pure_shear_viscosity(tau_pa, floor_spelling):
+    """Effective viscosity, in model units (eta / eta_0), of a ViscoPlasticFlowModel
+    in pure shear at strain rate R = 29.1 MPa / eta_0, the yield stress being the
+    cohesion expression times a unit field. Returns the evaluate() value, the
+    compiled (projected) nodal values, and R."""
+    uw.reset_default_model()
+    uw.get_default_model().set_reference_quantities(
+        length=uw.quantity(H, "m"),
+        velocity=uw.quantity(V, "m/s"),
+        viscosity=uw.quantity(ETA, "Pa*s"),
+    )
+    mesh = uw.meshing.StructuredQuadBox(
+        elementRes=(4, 4), minCoords=(0.0, 0.0),
+        maxCoords=(uw.quantity(H, "m"), uw.quantity(H, "m")),
+    )
+    v = uw.discretisation.MeshVariable("U", mesh, 2, degree=2, units="m/s")
+    p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1, units="Pa")
+    seed = uw.discretisation.MeshVariable("seed", mesh, 1, degree=1)
+    seed.array[...] = 1.0
+
+    rate = 29.1e6 / ETA                     # 1/s
+    x = np.asarray(v.coords.to("m"))
+    v.array[...] = uw.utilities.UnitAwareArray(
+        (rate * np.stack([x[:, 0], -x[:, 1]], axis=-1)).reshape(-1, 1, 2), units="m/s")
+
+    stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p)
+    stokes.constitutive_model = uw.constitutive_models.ViscoPlasticFlowModel
+    constitutive_model = stokes.constitutive_model
+    constitutive_model.Parameters.shear_viscosity_0 = uw.quantity(ETA, "Pa*s")
+    C0 = uw.expression(r"C_0", uw.quantity(tau_pa, "Pa"), "cohesion")
+    constitutive_model.Parameters.yield_stress = C0 * seed.sym[0]
+    constitutive_model.yield_mode = "min"
+    if floor_spelling is None:
+        # no floor at all: the yield stress reads a field, so the floor's absence
+        # is stated explicitly (#759)
+        constitutive_model.Parameters.yield_stress_min = -sympy.oo
+    else:
+        constitutive_model.Parameters.yield_stress_min = (
+            FLOOR_SPELLINGS[floor_spelling](C0, tau_pa))
+
+    viscosity = constitutive_model.viscosity
+    centre = np.array([[0.5, 0.5]])
+    evaluated = float(np.asarray(uw.function.evaluate(viscosity, centre)).ravel()[0])
+
+    eta_field = uw.discretisation.MeshVariable("eta", mesh, 1, degree=1)
+    projection = uw.systems.Projection(mesh, eta_field)
+    projection.uw_function = viscosity
+    projection.smoothing = 0.0
+    projection.tolerance = 1.0e-12
+    projection.solve()
+    compiled = np.asarray(eta_field.array).ravel()
+    return evaluated, compiled, rate
+
+
+@pytest.mark.tier_b
+@pytest.mark.parametrize("floor_spelling", [None, *FLOOR_SPELLINGS])
+def test_floor_far_below_the_yield_stress_is_inert(floor_spelling):
+    """tau_y = 29.1 MPa, floor 0.35 MPa: eta = tau_y / (2R) = eta_0 / 2 exactly."""
+    tau_pa = 29.1e6
+    evaluated, compiled, rate = _pure_shear_viscosity(tau_pa, floor_spelling)
+    expected = tau_pa / (2.0 * rate) / ETA
+    assert evaluated == pytest.approx(expected, rel=1.0e-10)
+    assert compiled == pytest.approx(np.full_like(compiled, expected), rel=1.0e-8)
+
+
+@pytest.mark.tier_b
+@pytest.mark.parametrize("floor_spelling", list(FLOOR_SPELLINGS))
+def test_floor_above_the_yield_stress_is_honoured(floor_spelling):
+    """tau_y = 0.2 MPa, floor 0.35 MPa: eta = floor / (2R), for every spelling."""
+    tau_pa = 0.2e6
+    evaluated, compiled, rate = _pure_shear_viscosity(tau_pa, floor_spelling)
+    expected = FLOOR_PA / (2.0 * rate) / ETA
+    assert evaluated == pytest.approx(expected, rel=1.0e-10)
+    assert compiled == pytest.approx(np.full_like(compiled, expected), rel=1.0e-8)
+
+
+@pytest.mark.tier_b
+def test_sum_of_field_expression_and_quantity_expression_scales_each_operand():
+    """The root cause, directly: a sum of a field expression and a dimensional
+    constant (expression or bare quantity) non-dimensionalises each term."""
+    uw.reset_default_model()
+    uw.get_default_model().set_reference_quantities(
+        length=uw.quantity(H, "m"),
+        velocity=uw.quantity(V, "m/s"),
+        viscosity=uw.quantity(ETA, "Pa*s"),
+    )
+    mesh = uw.meshing.StructuredQuadBox(
+        elementRes=(2, 2), minCoords=(0.0, 0.0),
+        maxCoords=(uw.quantity(H, "m"), uw.quantity(H, "m")),
+    )
+    seed = uw.discretisation.MeshVariable("seed_sum", mesh, 1, degree=1)
+    C0 = uw.expression(r"C_0", uw.quantity(29.1, "MPa"), "cohesion")
+    tau = uw.expression(r"\tau", C0 * seed.sym[0], "yield stress")
+    floor = uw.expression(r"f", uw.quantity(0.35, "MPa"), "floor")
+    for total in (tau + floor, floor + tau, tau + uw.quantity(0.35, "MPa"),
+                  tau - (-1 * uw.quantity(0.35, "MPa")),
+                  -(uw.quantity(-0.35, "MPa") - tau)):
+        nondimensional = unwrap_expression(total, mode="nondimensional")
+        at_unit_seed = float(nondimensional.subs(seed.sym[0], 1))
+        assert at_unit_seed == pytest.approx(
+            (29.1e6 + 0.35e6) / STRESS_SCALE, rel=1.0e-12)
diff --git a/tests/test_1068_picard_warmup_is_a_frozen_tangent_step.py b/tests/test_1068_picard_warmup_is_a_frozen_tangent_step.py
index 940329ea7..fbf321eb9 100644
--- a/tests/test_1068_picard_warmup_is_a_frozen_tangent_step.py
+++ b/tests/test_1068_picard_warmup_is_a_frozen_tangent_step.py
@@ -319,6 +319,7 @@ def _pressure_yield_box(tag, pressure_coupling):
     s.constitutive_model.Parameters.shear_viscosity_0 = 1.0
     # positive for every p, with a non-zero derivative in p (a smoothed |p|)
     s.constitutive_model.Parameters.yield_stress = 0.3 + 0.3 * sympy.sqrt(p.sym[0] ** 2 + 1.0e-2)
+    s.constitutive_model.Parameters.yield_stress_min = -sympy.oo   # positive by construction: no floor needed (#759)
     s.add_essential_bc((1.0, 0.0), "Top")
     s.add_essential_bc((0.0, 0.0), "Bottom")
     s.petsc_options.delValue("ksp_monitor")
diff --git a/tests/test_1069_pressure_gauge_is_recorded.py b/tests/test_1069_pressure_gauge_is_recorded.py
index 989da001d..61189546b 100644
--- a/tests/test_1069_pressure_gauge_is_recorded.py
+++ b/tests/test_1069_pressure_gauge_is_recorded.py
@@ -30,6 +30,7 @@ def _box(tag, closed, pressure_dependent, nullspace):
     s.constitutive_model.Parameters.shear_viscosity_0 = 1.0
     if pressure_dependent:
         s.constitutive_model.Parameters.yield_stress = 0.3 + 0.3 * sympy.sqrt(p.sym[0] ** 2 + 1.0e-2)
+        s.constitutive_model.Parameters.yield_stress_min = -sympy.oo   # positive by construction: no floor needed (#759)
         # a stated regularisation, so the only warning these tests can see is the gauge one
         s.constitutive_model.Parameters.plastic_rate_strengthening = 0.05
     else:
@@ -176,6 +177,7 @@ def test_switching_on_a_pressure_dependent_yield_is_seen_by_the_record():
     assert s.solve_report.config["gauge"]["pressure_in_rheology"] is False
     cm = s.constitutive_model
     cm.Parameters.yield_stress = 0.3 + 0.3 * sympy.sqrt(s.p.sym[0] ** 2 + 1.0e-2)
+    cm.Parameters.yield_stress_min = -sympy.oo   # positive by construction: no floor needed (#759)
     cm.Parameters.plastic_rate_strengthening = 0.05
     with pytest.warns(RuntimeWarning, match="GAUGE-fixed pressure"):
         s.solve(zero_init_guess=False)
diff --git a/tests/test_1071_plastic_rate_strengthening.py b/tests/test_1071_plastic_rate_strengthening.py
index 2d02eb79e..470b89a7e 100644
--- a/tests/test_1071_plastic_rate_strengthening.py
+++ b/tests/test_1071_plastic_rate_strengthening.py
@@ -29,6 +29,7 @@ def _box(tag, eta_reg):
     cm = s.constitutive_model
     cm.Parameters.shear_viscosity_0 = 1.0
     cm.Parameters.yield_stress = 0.3 + 0.3 * sympy.sqrt(p.sym[0] ** 2 + 1.0e-2)
+    cm.Parameters.yield_stress_min = -sympy.oo   # positive by construction: no floor needed (#759)
     if eta_reg is not None:          # None = never stated; 0.0 = stated as zero
         cm.Parameters.plastic_rate_strengthening = eta_reg
     s.add_essential_bc((1.0, 0.0), "Top")