diff --git a/cpp/core/coupled_rates.cpp b/cpp/core/coupled_rates.cpp
index a70efd1..60abf1f 100644
--- a/cpp/core/coupled_rates.cpp
+++ b/cpp/core/coupled_rates.cpp
@@ -96,6 +96,7 @@ void CoupledRatePlan::validate() const {
       case RateOp::growth_rate:
       case RateOp::cell_type:
       case RateOp::cell_volume:
+      case RateOp::cell_volume_change_rate:
       case RateOp::cell_surface_area:
         break;
       case RateOp::negate:
diff --git a/cpp/core/mechanics_integration.cpp b/cpp/core/mechanics_integration.cpp
index 25e4a6e..bd9e310 100644
--- a/cpp/core/mechanics_integration.cpp
+++ b/cpp/core/mechanics_integration.cpp
@@ -72,12 +72,14 @@ void integrate_mechanics_result(WorldState& state, const MechanicsSolveResult& r
       throw std::invalid_argument("desired length increments must be finite and non-negative");
     }
     const auto applied_length_increment =
-        cell.fixed ? desired_increment : std::max(0.0F, desired_increment + correction.length);
+        desired_length_increments.empty()
+            ? 0.0F
+            : (cell.fixed ? desired_increment
+                          : std::max(0.0F, desired_increment + correction.length));
     const auto new_position = cell.fixed ? cell.position : cell.position + correction.translation;
-    const auto new_direction =
-        cell.fixed ? cell.direction
-                   : rotate_axis_angle(cell.direction, correction.rotation,
-                                       parameters.max_rotation_radians);
+    const auto new_direction = cell.fixed ? cell.direction
+                                          : rotate_axis_angle(cell.direction, correction.rotation,
+                                                              parameters.max_rotation_radians);
     const auto new_length = cell.length + applied_length_increment;
     if (!finite(new_position) || !finite(new_direction) || !std::isfinite(new_length)) {
       throw std::overflow_error("mechanics integration produced non-finite geometry");
diff --git a/cpp/core/species.cpp b/cpp/core/species.cpp
index 3753489..3be757c 100644
--- a/cpp/core/species.cpp
+++ b/cpp/core/species.cpp
@@ -85,6 +85,7 @@ void SpeciesRatePlan::validate() const {
       case RateOp::growth_rate:
       case RateOp::cell_type:
       case RateOp::cell_volume:
+      case RateOp::cell_volume_change_rate:
       case RateOp::cell_surface_area:
         break;
       case RateOp::negate:
@@ -125,6 +126,8 @@ void SpeciesRatePlan::validate() const {
 }
 
 float effective_cell_volume(float length, float radius) noexcept {
+  // Conserved biomass volume for the endpoint-preserving division rule.
+  // This is a biochemical measure, not the geometric capsule volume.
   constexpr float pi = 3.14159265358979323846F;
   return pi * radius * radius * (length + 2.0F * radius);
 }
diff --git a/cpp/cpu/cpu_coupled.cpp b/cpp/cpu/cpu_coupled.cpp
index 0b73526..a8efeb2 100644
--- a/cpp/cpu/cpu_coupled.cpp
+++ b/cpp/cpu/cpu_coupled.cpp
@@ -16,7 +16,7 @@ namespace {
 float evaluate_instruction(const RateInstruction& instruction, std::span<const float> workspace,
                            std::span<const float> species, std::span<const float> signals,
                            const CellGeometryView& geometry, const CellAttributeView& attributes,
-                           std::size_t cell) {
+                           std::size_t cell, float volume_change_rate) {
   switch (instruction.operation) {
     case RateOp::constant:
       return instruction.value;
@@ -38,6 +38,8 @@ float evaluate_instruction(const RateInstruction& instruction, std::span<const f
       return attributes.growth_rates[cell];
     case RateOp::cell_type:
       return static_cast<float>(attributes.cell_types[cell]);
+    case RateOp::cell_volume_change_rate:
+      return volume_change_rate;
     case RateOp::cell_volume:
       return effective_cell_volume(geometry.lengths[cell], geometry.radii[cell]);
     case RateOp::cell_surface_area:
@@ -161,8 +163,13 @@ SignalSolveReport advance_coupled_cpu(WorldState& state, SignalGrid& grid,
     const auto cell_signals =
         std::span<const float>(sampled).subspan(cell * plan.signal_count(), plan.signal_count());
     for (std::size_t index = 0; index < plan.instructions().size(); ++index) {
-      workspace[index] = evaluate_instruction(plan.instructions()[index], workspace, cell_species,
-                                              cell_signals, geometry, attributes, cell);
+      workspace[index] = evaluate_instruction(
+          plan.instructions()[index], workspace, cell_species, cell_signals, geometry, attributes,
+          cell,
+          dt == 0.0F ? 0.0F
+                     : (effective_cell_volume(geometry.lengths[cell], geometry.radii[cell]) -
+                        effective_cell_volume(previous_lengths[cell], geometry.radii[cell])) /
+                           dt);
       if (!std::isfinite(workspace[index])) {
         throw std::domain_error("coupled rate instruction " + std::to_string(index) +
                                 " produced a non-finite value");
diff --git a/cpp/cpu/cpu_species.cpp b/cpp/cpu/cpu_species.cpp
index f6ad98b..89d7490 100644
--- a/cpp/cpu/cpu_species.cpp
+++ b/cpp/cpu/cpu_species.cpp
@@ -13,7 +13,8 @@ namespace {
 
 float evaluate_instruction(const RateInstruction& instruction, std::span<const float> workspace,
                            std::span<const float> species, const CellGeometryView& geometry,
-                           const CellAttributeView& attributes, std::size_t cell) {
+                           const CellAttributeView& attributes, std::size_t cell,
+                           float volume_change_rate) {
   switch (instruction.operation) {
     case RateOp::constant:
       return instruction.value;
@@ -35,6 +36,8 @@ float evaluate_instruction(const RateInstruction& instruction, std::span<const f
       return attributes.growth_rates[cell];
     case RateOp::cell_type:
       return static_cast<float>(attributes.cell_types[cell]);
+    case RateOp::cell_volume_change_rate:
+      return volume_change_rate;
     case RateOp::cell_volume:
       return effective_cell_volume(geometry.lengths[cell], geometry.radii[cell]);
     case RateOp::cell_surface_area:
@@ -119,8 +122,12 @@ void advance_species_cpu(WorldState& state, const SpeciesRatePlan& plan,
     const auto cell_species =
         std::span<const float>(next_levels).subspan(offset, state.species_count());
     for (std::size_t index = 0; index < plan.instructions().size(); ++index) {
-      workspace[index] = evaluate_instruction(plan.instructions()[index], workspace, cell_species,
-                                              geometry, attributes, cell);
+      workspace[index] = evaluate_instruction(
+          plan.instructions()[index], workspace, cell_species, geometry, attributes, cell,
+          dt == 0.0F ? 0.0F
+                     : (effective_cell_volume(geometry.lengths[cell], geometry.radii[cell]) -
+                        effective_cell_volume(previous_lengths[cell], geometry.radii[cell])) /
+                           dt);
       if (!std::isfinite(workspace[index])) {
         throw std::domain_error("species rate instruction " + std::to_string(index) +
                                 " produced a non-finite value");
diff --git a/cpp/cuda/kernels/coupled_rates.cu b/cpp/cuda/kernels/coupled_rates.cu
index 5f8b67b..55f98ce 100644
--- a/cpp/cuda/kernels/coupled_rates.cu
+++ b/cpp/cuda/kernels/coupled_rates.cu
@@ -92,7 +92,8 @@ __device__ float sample_signal(const float* levels, SignalGridShapeGpu shape, fl
 
 __device__ float evaluate_instruction(const RateInstructionGpu& instruction, const float* workspace,
                                       const float* species, const float* signals, float4 center,
-                                      float4 geometry, float growth_rate, std::int32_t cell_type) {
+                                      float4 geometry, float growth_rate, std::int32_t cell_type,
+                                      float volume_change_rate) {
   switch (instruction.operation) {
     case 0:
       return instruction.value;
@@ -112,6 +113,8 @@ __device__ float evaluate_instruction(const RateInstructionGpu& instruction, con
       return growth_rate;
     case 8:
       return static_cast<float>(cell_type);
+    case 28:
+      return volume_change_rate;
     case 9:
       return effective_volume(geometry.x, geometry.y);
     case 10:
@@ -190,7 +193,11 @@ __global__ void advance_coupled_cells(
   for (std::uint32_t index = 0; index < instruction_count; ++index) {
     const auto value =
         evaluate_instruction(instructions[index], cell_workspace, cell_species, cell_signals,
-                             centers[cell], geometry[cell], growth_rates[cell], cell_types[cell]);
+                             centers[cell], geometry[cell], growth_rates[cell], cell_types[cell],
+                             dt == 0.0f ? 0.0f
+                                        : (effective_volume(geometry[cell].x, radius) -
+                                           effective_volume(previous_lengths[cell], radius)) /
+                                              dt);
     cell_workspace[index] = value;
     if (!isfinite(value)) {
       atomicOr(error, 1U);
diff --git a/cpp/cuda/kernels/species.cu b/cpp/cuda/kernels/species.cu
index ae27b1b..afbce91 100644
--- a/cpp/cuda/kernels/species.cu
+++ b/cpp/cuda/kernels/species.cu
@@ -16,7 +16,8 @@ __device__ float effective_surface_area(float length, float radius) {
 
 __device__ float evaluate_instruction(const RateInstructionGpu& instruction, const float* workspace,
                                       const float* species, float4 center, float4 geometry,
-                                      float growth_rate, std::int32_t cell_type) {
+                                      float growth_rate, std::int32_t cell_type,
+                                      float volume_change_rate) {
   switch (instruction.operation) {
     case 0:
       return instruction.value;
@@ -36,6 +37,8 @@ __device__ float evaluate_instruction(const RateInstructionGpu& instruction, con
       return growth_rate;
     case 8:
       return static_cast<float>(cell_type);
+    case 28:
+      return volume_change_rate;
     case 9:
       return effective_volume(geometry.x, geometry.y);
     case 10:
@@ -104,7 +107,11 @@ __global__ void advance_species(float* levels, const float* previous_lengths, co
   for (std::uint32_t index = 0; index < instruction_count; ++index) {
     const auto value =
         evaluate_instruction(instructions[index], cell_workspace, cell_species, centers[cell],
-                             geometry[cell], growth_rates[cell], cell_types[cell]);
+                             geometry[cell], growth_rates[cell], cell_types[cell],
+                             dt == 0.0f ? 0.0f
+                                        : (effective_volume(geometry[cell].x, radius) -
+                                           effective_volume(previous_lengths[cell], radius)) /
+                                              dt);
     cell_workspace[index] = value;
     if (!isfinite(value)) {
       atomicOr(error, 1U);
diff --git a/cpp/include/cm/species.hpp b/cpp/include/cm/species.hpp
index ee11509..1dc4a86 100644
--- a/cpp/include/cm/species.hpp
+++ b/cpp/include/cm/species.hpp
@@ -38,6 +38,7 @@ enum class RateOp : std::uint8_t {
   equal = 25,
   select = 26,
   signal = 27,
+  cell_volume_change_rate = 28,
 };
 
 struct RateInstruction {
diff --git a/cpp/metal/kernels/coupled_rates.metal b/cpp/metal/kernels/coupled_rates.metal
index 72093db..04b6244 100644
--- a/cpp/metal/kernels/coupled_rates.metal
+++ b/cpp/metal/kernels/coupled_rates.metal
@@ -99,7 +99,7 @@ float sample_signal(device const float* levels, GridShape shape, float4 origin,
 
 float evaluate_instruction(const RateInstruction instruction, device const float* workspace,
                            device const float* species, device const float* signals, float4 center,
-                           float4 geometry, float growth_rate, int cell_type) {
+                           float4 geometry, float growth_rate, int cell_type, float volume_change_rate) {
   switch (instruction.operation) {
     case 0:
       return instruction.value;
@@ -119,6 +119,8 @@ float evaluate_instruction(const RateInstruction instruction, device const float
       return growth_rate;
     case 8:
       return float(cell_type);
+    case 28:
+      return volume_change_rate;
     case 9:
       return effective_volume(geometry.x, geometry.y);
     case 10:
@@ -201,7 +203,9 @@ kernel void advance_coupled_cells(
   for (uint index = 0; index < instruction_count; ++index) {
     float value =
         evaluate_instruction(instructions[index], cell_workspace, cell_species, cell_signals,
-                             centers[cell], geometry[cell], growth_rates[cell], cell_types[cell]);
+                             centers[cell], geometry[cell], growth_rates[cell], cell_types[cell],
+            dt == 0.0f ? 0.0f : (effective_volume(geometry[cell].x, radius) -
+                                 effective_volume(previous_lengths[cell], radius)) / dt);
     cell_workspace[index] = value;
     if (!isfinite(value)) {
       atomic_fetch_or_explicit(error, 1u, memory_order_relaxed);
diff --git a/cpp/metal/kernels/species.metal b/cpp/metal/kernels/species.metal
index 15b84ba..8e07760 100644
--- a/cpp/metal/kernels/species.metal
+++ b/cpp/metal/kernels/species.metal
@@ -26,7 +26,7 @@ float evaluate_instruction(const RateInstruction instruction,
                            float4 center,
                            float4 geometry,
                            float growth_rate,
-                           int cell_type) {
+                           int cell_type, float volume_change_rate) {
   switch (instruction.operation) {
     case 0:
       return instruction.value;
@@ -46,6 +46,8 @@ float evaluate_instruction(const RateInstruction instruction,
       return growth_rate;
     case 8:
       return float(cell_type);
+    case 28:
+      return volume_change_rate;
     case 9:
       return effective_volume(geometry.x, geometry.y);
     case 10:
@@ -122,7 +124,9 @@ kernel void advance_species(
   for (uint index = 0; index < instruction_count; ++index) {
     float value = evaluate_instruction(instructions[index], cell_workspace, cell_species,
                                        centers[cell], geometry[cell], growth_rates[cell],
-                                       cell_types[cell]);
+                                       cell_types[cell],
+            dt == 0.0f ? 0.0f : (effective_volume(geometry[cell].x, radius) -
+                                 effective_volume(previous_lengths[cell], radius)) / dt);
     cell_workspace[index] = value;
     if (!isfinite(value)) {
       atomic_fetch_or_explicit(error, 1u, memory_order_relaxed);
diff --git a/cpp/python/bindings.cpp b/cpp/python/bindings.cpp
index ac17539..e6ef22e 100644
--- a/cpp/python/bindings.cpp
+++ b/cpp/python/bindings.cpp
@@ -49,6 +49,7 @@ NB_MODULE(_core, module) {
       .value("GROWTH_RATE", cm::RateOp::growth_rate)
       .value("CELL_TYPE", cm::RateOp::cell_type)
       .value("CELL_VOLUME", cm::RateOp::cell_volume)
+      .value("CELL_VOLUME_CHANGE_RATE", cm::RateOp::cell_volume_change_rate)
       .value("CELL_SURFACE_AREA", cm::RateOp::cell_surface_area)
       .value("ADD", cm::RateOp::add)
       .value("SUBTRACT", cm::RateOp::subtract)
@@ -215,7 +216,7 @@ NB_MODULE(_core, module) {
       .def_prop_ro("instructions",
                    [](const cm::SpeciesRatePlan& plan) {
                      return std::vector<cm::RateInstruction>(plan.instructions().begin(),
-                                                              plan.instructions().end());
+                                                             plan.instructions().end());
                    })
       .def_prop_ro("outputs",
                    [](const cm::SpeciesRatePlan& plan) {
@@ -234,7 +235,7 @@ NB_MODULE(_core, module) {
       .def_prop_ro("instructions",
                    [](const cm::CoupledRatePlan& plan) {
                      return std::vector<cm::RateInstruction>(plan.instructions().begin(),
-                                                              plan.instructions().end());
+                                                             plan.instructions().end());
                    })
       .def_prop_ro("species_outputs",
                    [](const cm::CoupledRatePlan& plan) {
@@ -271,7 +272,7 @@ NB_MODULE(_core, module) {
       .def_prop_ro("contacts",
                    [](const cm::ContactGraph& graph) {
                      return std::vector<cm::CellContact>(graph.contacts().begin(),
-                                                          graph.contacts().end());
+                                                         graph.contacts().end());
                    })
       .def("__len__", &cm::ContactGraph::size)
       .def(
@@ -357,7 +358,7 @@ NB_MODULE(_core, module) {
       .def_prop_ro("contacts",
                    [](const cm::ExternalContactGraph& graph) {
                      return std::vector<cm::ExternalContact>(graph.contacts().begin(),
-                                                              graph.contacts().end());
+                                                             graph.contacts().end());
                    })
       .def("__len__", &cm::ExternalContactGraph::size)
       .def(
@@ -400,8 +401,8 @@ NB_MODULE(_core, module) {
       .def(nb::init<cm::BackendKind, std::size_t, std::size_t, std::uint32_t>(),
            "backend"_a = cm::BackendKind::cpu, "reserved_capacity"_a = 0, "species_count"_a = 0,
            "device_index"_a = 0)
-      .def(nb::init<cm::BackendKind, const cm::SimulationCheckpoint&, std::uint32_t>(),
-           "backend"_a, "checkpoint"_a, "device_index"_a = 0)
+      .def(nb::init<cm::BackendKind, const cm::SimulationCheckpoint&, std::uint32_t>(), "backend"_a,
+           "checkpoint"_a, "device_index"_a = 0)
       .def_prop_ro("backend_info", &cm::Simulation::backend_info)
       .def("supports", &cm::Simulation::supports, "feature"_a)
       .def_prop_ro("time", &cm::Simulation::time)
diff --git a/docs/architecture/0024-biomass-accounting.md b/docs/architecture/0024-biomass-accounting.md
new file mode 100644
index 0000000..e0adddf
--- /dev/null
+++ b/docs/architecture/0024-biomass-accounting.md
@@ -0,0 +1,11 @@
+# Biomass accounting
+
+The conserved biochemical volume is B = pi r²(l + 2r), with rod centerline length l and radius r. At constant biomass density rho, biomass mass is rho B and an intracellular concentration c represents amount c B. This effective volume is distinct from the geometric capsule volume pi r²l + 4pi r³/3. The latter changes when a parent capsule is replaced by two rounded daughters, so it must not be used interchangeably with B in biomass-dependent feedback.
+
+Native division preserves the parent's outer endpoints and satisfies l1 + l2 = l - 2r. It therefore conserves B exactly in real arithmetic, including unequal division. The division fraction partitions the available cylindrical length, not total biomass. Both daughters inherit concentrations, so intracellular amount is conserved too. Geometry remains a coarse representation of division rather than a volume-preserving model of septum formation.
+
+Growth retains the elongation law dl/dt = g l. Consequently dB/dt = pi r² g l, not g B. `RatePlanBuilder.cell_volume_change_rate()` returns the realized change (B_after - B_before)/dt during the current growth step, with zero at dt=0. A nutrient sink `-rates.cell_volume_change_rate()/yield` therefore consumes exactly the biomass increment divided by yield, up to floating-point error. `rates.cell_volume()` and `microsimulator.biomass.biomass_volume` expose B. A negative concentration is a rejected step, never an invitation to clip away a mass-balance error.
+
+Contact relaxation without explicitly prescribed length increments preserves each rod's length. Numerical overlap corrections therefore cannot manufacture biomass after the biological update. Direct geometry edits and explicitly prescribed mechanical length increments are externally imposed changes; callers must account for their biomass and intracellular amounts.
+
+Colony resistance uses a calibrated biomass-density closure derived from B. Its density must be deposited conservatively over a physical averaging scale independent of voxel size. It is not an exact solid-occupancy reconstruction of individual capsules.
diff --git a/docs/architecture/README.md b/docs/architecture/README.md
index efdf323..219d0c4 100644
--- a/docs/architecture/README.md
+++ b/docs/architecture/README.md
@@ -29,3 +29,5 @@ The [numerical contract](numerical-contract.md) is the best starting point for w
 - [Versioned columnar analysis datasets](0013-analysis-datasets.md)
 
 Current backend status and the tests required to support it are documented in [testing and validation](../development/validation.md).
+
+- [Biomass accounting](0024-biomass-accounting.md)
diff --git a/python/src/microsimulator/_core.pyi b/python/src/microsimulator/_core.pyi
index f097eff..927af77 100644
--- a/python/src/microsimulator/_core.pyi
+++ b/python/src/microsimulator/_core.pyi
@@ -35,6 +35,7 @@ class RateOp(Enum):
     GROWTH_RATE: RateOp
     CELL_TYPE: RateOp
     CELL_VOLUME: RateOp
+    CELL_VOLUME_CHANGE_RATE: RateOp
     CELL_SURFACE_AREA: RateOp
     ADD: RateOp
     SUBTRACT: RateOp
diff --git a/python/src/microsimulator/biomass.py b/python/src/microsimulator/biomass.py
new file mode 100644
index 0000000..523d055
--- /dev/null
+++ b/python/src/microsimulator/biomass.py
@@ -0,0 +1,19 @@
+"""The conserved biochemical volume and its distinct geometric counterpart."""
+
+import math
+
+
+def biomass_volume(length: float, radius: float) -> float:
+    """Return pi*r**2*(length + 2*r), conserved by native cell division.
+
+    At constant biomass density, multiplying this measure by density gives
+    biomass mass. It is the same volume used by native concentration dilution.
+    """
+    if not math.isfinite(length) or length < 0 or not math.isfinite(radius) or radius <= 0:
+        raise ValueError("length must be finite and nonnegative; radius finite and positive")
+    return math.pi * radius * radius * (length + 2 * radius)
+
+
+def capsule_volume(length: float, radius: float) -> float:
+    """Return the geometric capsule volume, which is not the biomass measure."""
+    return biomass_volume(length, radius) - (2 / 3) * math.pi * radius**3
diff --git a/python/src/microsimulator/checkpoint.py b/python/src/microsimulator/checkpoint.py
index d2c2c4a..d0e219c 100644
--- a/python/src/microsimulator/checkpoint.py
+++ b/python/src/microsimulator/checkpoint.py
@@ -93,6 +93,7 @@ class CheckpointBundle:
     RateOp.GROWTH_RATE: "growth_rate",
     RateOp.CELL_TYPE: "cell_type",
     RateOp.CELL_VOLUME: "cell_volume",
+    RateOp.CELL_VOLUME_CHANGE_RATE: "cell_volume_change_rate",
     RateOp.CELL_SURFACE_AREA: "cell_surface_area",
     RateOp.ADD: "add",
     RateOp.SUBTRACT: "subtract",
diff --git a/python/src/microsimulator/rates.py b/python/src/microsimulator/rates.py
index cca57ae..90eed4a 100644
--- a/python/src/microsimulator/rates.py
+++ b/python/src/microsimulator/rates.py
@@ -165,8 +165,18 @@ def cell_type(self) -> RateExpression:
         return self._source(RateOp.CELL_TYPE)
 
     def cell_volume(self) -> RateExpression:
+        """Effective biomass volume pi*r**2*(length + 2*r)."""
         return self._source(RateOp.CELL_VOLUME)
 
+    def cell_volume_change_rate(self) -> RateExpression:
+        """Realized biomass-volume increase / dt in this step (zero at dt=0).
+
+        Use ``-cell_volume_change_rate() / yield`` for growth-linked uptake.
+        This includes the discrete growth increment instead of approximating
+        it with the elongation rate times the post-growth cell volume.
+        """
+        return self._source(RateOp.CELL_VOLUME_CHANGE_RATE)
+
     def cell_surface_area(self) -> RateExpression:
         return self._source(RateOp.CELL_SURFACE_AREA)
 
diff --git a/python/tests/test_biomass.py b/python/tests/test_biomass.py
new file mode 100644
index 0000000..256d2d7
--- /dev/null
+++ b/python/tests/test_biomass.py
@@ -0,0 +1,76 @@
+from __future__ import annotations
+
+import math
+from pathlib import Path
+
+import pytest
+from microsimulator import (
+    BackendKind,
+    CellInit,
+    GridShape,
+    RatePlanBuilder,
+    SignalGridSpec,
+    Simulation,
+    Vec3,
+    backend_available,
+    load_checkpoint,
+    save_checkpoint,
+)
+from microsimulator.biomass import biomass_volume, capsule_volume
+
+
+@pytest.mark.parametrize("fraction", [0.2, 0.5, 0.8])
+def test_division_preserves_biomass_and_intracellular_amount(fraction: float) -> None:
+    simulation = Simulation(species_count=1)
+    cell = CellInit()
+    cell.length, cell.radius, cell.species = 6, 0.5, [3]
+    parent = simulation.add_cell(cell)
+    initial = biomass_volume(cell.length, cell.radius)
+    daughters = [simulation.cell(i) for i in simulation.divide(parent, fraction)]
+    assert math.isclose(
+        sum(biomass_volume(c.length, c.radius) for c in daughters), initial, rel_tol=1e-6
+    )
+    assert math.isclose(
+        sum(c.species[0] * biomass_volume(c.length, c.radius) for c in daughters),
+        3 * initial,
+        rel_tol=1e-6,
+    )
+    assert sum(capsule_volume(c.length, c.radius) for c in daughters) < capsule_volume(
+        cell.length, cell.radius
+    )
+
+
+@pytest.mark.parametrize("backend", list(BackendKind))
+def test_realized_growth_consumes_exactly_its_yield_on_every_backend(
+    backend: BackendKind, tmp_path: Path
+) -> None:
+    if not backend_available(backend):
+        pytest.skip("backend unavailable")
+    simulation = Simulation(backend, species_count=1)
+    spec = SignalGridSpec()
+    shape = GridShape()
+    shape.x = shape.y = shape.z = 1
+    spec.shape = shape
+    spec.spacing = Vec3(4, 4, 4)
+    spec.signal_count, spec.diffusion, spec.advection = 1, [0], [Vec3()]
+    simulation.configure_signal_grid(spec, [10])
+    rates = RatePlanBuilder()
+    simulation.set_coupled_rate_plan(
+        rates.coupled_plan(1, 1, (rates.constant(0),), (-rates.cell_volume_change_rate() / 0.4,))
+    )
+    cell = CellInit()
+    cell.length, cell.radius, cell.growth_rate, cell.species = 2, 0.5, 0.7, [3]
+    cid = simulation.add_cell(cell)
+    initial = biomass_volume(cell.length, cell.radius)
+    for dt in [0, 0.03, 0.1, 0.2]:
+        simulation.step(dt)
+        current = simulation.cell(cid)
+        volume = biomass_volume(current.length, current.radius)
+        consumed = (10 - simulation.signal_levels[0]) * spec.voxel_volume
+        assert math.isclose(0.4 * consumed, volume - initial, abs_tol=3e-5)
+        assert math.isclose(current.species[0] * volume, 3 * initial, rel_tol=2e-6)
+    save_checkpoint(simulation, tmp_path / "biomass.json")
+    restored = load_checkpoint(tmp_path / "biomass.json", backend=backend)
+    simulation.step(0.02)
+    restored.step(0.02)
+    assert restored.signal_levels == simulation.signal_levels
diff --git a/tests/cpp/mechanics_integration_test.cpp b/tests/cpp/mechanics_integration_test.cpp
index 0ce2741..1d693f6 100644
--- a/tests/cpp/mechanics_integration_test.cpp
+++ b/tests/cpp/mechanics_integration_test.cpp
@@ -75,6 +75,16 @@ void test_validation_is_atomic_and_requires_convergence() {
   assert(close(state.cell(id).position.x, 0.0F));
 }
 
+void test_relaxation_cannot_create_biomass() {
+  cm::WorldState state;
+  const auto id = state.add_cell(cm::CellInit{});
+  const auto before = state.cell(id).length;
+  cm::MechanicsSolveResult result;
+  result.corrections = {cm::CellCorrection{{}, {}, 0.5F}};
+  cm::integrate_mechanics_result(state, result);
+  assert(state.cell(id).length == before);
+}
+
 void test_fixed_cell_integration_only_applies_declared_growth() {
   cm::WorldState state;
   cm::CellInit cell;
@@ -129,6 +139,7 @@ void test_simulation_relaxation_reduces_penetration() {
 }  // namespace
 
 int main() {
+  test_relaxation_cannot_create_biomass();
   test_integration_applies_declared_geometry_semantics();
   test_validation_is_atomic_and_requires_convergence();
   test_fixed_cell_integration_only_applies_declared_growth();