diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py
index 37524079..9bbe79cd 100644
--- a/src/underworld3/systems/solvers.py
+++ b/src/underworld3/systems/solvers.py
@@ -1693,13 +1693,21 @@ def bodyforce(self, value):
         if isinstance(value, uw.function.expressions.UWexpression):
             self._bodyforce.sym = value.sym
         else:
-            # Convert UWQuantity objects to SymPy expressions before Matrix creation
+            # Convert UWQuantity components to their NON-DIMENSIONAL value before
+            # Matrix creation — the body force feeds the non-dimensional solver
+            # (see the ND<->units boundary contract / #282). The component may be
+            # symbolic (e.g. a buoyancy ρα g (T − T_ref) carries the T field), so
+            # non_dimensionalise yields a (1×1) array/Matrix whose element we
+            # extract. (The previous code called item._sympify_(), which does not
+            # exist on UWQuantity; _sympy_ raises on a dimensional quantity.)
             converted_value = []
             for item in value:
                 if isinstance(item, uw.function.quantities.UWQuantity):
-                    sympified = item._sympify_()
-                    # If UWQuantity contains a Matrix, extract the scalar element
-                    if hasattr(sympified, "shape") and sympified.shape == (1, 1):
+                    nd = uw.non_dimensionalise(item)
+                    sympified = nd.magnitude if hasattr(nd, "magnitude") else nd
+                    # If the non-dimensional value is a 1x1 array/Matrix, extract
+                    # the scalar element.
+                    if hasattr(sympified, "shape") and tuple(sympified.shape) == (1, 1):
                         converted_value.append(sympified[0, 0])
                     else:
                         converted_value.append(sympified)
diff --git a/tests/test_1011_bodyforce_units.py b/tests/test_1011_bodyforce_units.py
new file mode 100644
index 00000000..c512d78a
--- /dev/null
+++ b/tests/test_1011_bodyforce_units.py
@@ -0,0 +1,59 @@
+#!/usr/bin/env python3
+"""Regression: Stokes.bodyforce accepts a UWQuantity component (units-active).
+
+A units-active buoyancy force ``F = ρ₀ α g (T − T_ref) ẑ`` is a UWQuantity with a
+*symbolic* magnitude (it carries the T field). The bodyforce setter must
+non-dimensionalise it before building the (ND) body-force matrix. Previously it
+called ``item._sympify_()`` — which does not exist on UWQuantity (``_sympy_``
+also raises on a dimensional quantity) — so any units-active buoyancy assignment
+crashed (the second blocker, after #267/#282, for the thermal-convection units
+tutorial).
+"""
+import sympy
+import pytest
+
+import underworld3 as uw
+
+pytestmark = [pytest.mark.level_1, pytest.mark.tier_a]
+
+
+def test_bodyforce_accepts_uwquantity_buoyancy():
+    uw.reset_default_model()
+    model = uw.get_default_model()
+    model.set_reference_quantities(
+        domain_depth=uw.quantity(1000, "km"),
+        plate_velocity=uw.quantity(5, "cm/year"),
+        mantle_viscosity=uw.quantity(1e21, "Pa*s"),
+        temperature_difference=uw.quantity(1000, "K"),
+    )
+    mesh = uw.meshing.StructuredQuadBox(
+        elementRes=(4, 4), minCoords=(0.0, 0.0), maxCoords=(1000.0, 1000.0), units="km"
+    )
+    v = uw.discretisation.MeshVariable("v", mesh, mesh.dim, degree=2)
+    p = uw.discretisation.MeshVariable("p", mesh, 1, degree=1)
+    T = uw.discretisation.MeshVariable("T", mesh, 1, degree=2, units="K")
+
+    stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p)
+
+    rho = uw.quantity(3300, "kg/m**3")
+    alpha = uw.quantity(3e-5, "1/K")
+    g = uw.quantity(9.8, "m/s**2")
+    T_ref = uw.quantity(500, "K")
+    buoyancy = rho * alpha * g * (T - T_ref)  # symbolic UWQuantity
+
+    # Must not raise; the vertical component becomes a non-dimensional sympy expr.
+    stokes.bodyforce = [0, buoyancy]
+    fz = stokes.bodyforce.sym[1]
+    assert isinstance(fz, sympy.Expr)
+    assert fz.free_symbols, "body force should depend on the T field"
+
+
+def test_bodyforce_plain_vector_unchanged():
+    """A plain (non-quantity) body force is unaffected."""
+    uw.reset_default_model()
+    mesh = uw.meshing.StructuredQuadBox(elementRes=(4, 4), minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0))
+    v = uw.discretisation.MeshVariable("v2", mesh, mesh.dim, degree=2)
+    p = uw.discretisation.MeshVariable("p2", mesh, 1, degree=1)
+    stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p)
+    stokes.bodyforce = sympy.Matrix([0.0, -1.0])
+    assert stokes.bodyforce.sym[1] == -1.0