diff --git a/changelog.d/51.breaking.md b/changelog.d/51.breaking.md
new file mode 100644
index 0000000..b2bb9a0
--- /dev/null
+++ b/changelog.d/51.breaking.md
@@ -0,0 +1 @@
+Packed-bed results change: during a charge, the fluid temperature at each height no longer overshoots and then falls by 1 K to 2 K. While the PCM at a capsule's surface melted, the capsule took up to five times less heat than the heat transfer coefficient gives; it now takes exactly that heat. The default packed bed now fits the measurements of Nallusamy et al. (2007) better at every height, and the finned-tube models move by at most 0.006 K. There is nothing to change in your code.
diff --git a/changelog.d/51.fix.md b/changelog.d/51.fix.md
new file mode 100644
index 0000000..b0e3f14
--- /dev/null
+++ b/changelog.d/51.fix.md
@@ -0,0 +1 @@
+The heat entering a storage through its wall is now exactly the surface flux law's heat flux, also while the material by the wall changes phase (F-056).
diff --git a/docs/design/DESIGN.md b/docs/design/DESIGN.md
index e29d4e4..a7d8b35 100644
--- a/docs/design/DESIGN.md
+++ b/docs/design/DESIGN.md
@@ -208,7 +208,12 @@ never from a separately discretised flux. This applies to **every** coupling, fl
 and storage-storage (maintainer, 2026-10-14). Rationale: PyBaMM imposes a Neumann condition on
 the gradient, so the flux the finite volumes use is `k_face / k_surface` times the one
 requested, and with a temperature (Dirichlet) condition it is whatever the scheme computes;
-the ledger is exact in both cases, for any conductivity. In a bed, the film resistance works
+the ledger is exact in both cases, for any conductivity. Exact conservation does not make the
+flux right, though: with `k_surface` the conductivity at the extrapolated wall enthalpy, the
+factor reached 5 while a PCM by the wall melted, and the packed bed took a fifth of the flux
+law's heat (F-056). So a storage's heat-flux condition divides by `k_face` itself (the node
+values of `k` extrapolated to the wall, as `k grad(T)` uses), making the factor 1: the flux
+law's flux enters, and the ledger still counts it exactly. In a bed, the film resistance works
 as at the finned-tube wall: the flux uses the boundary value of `T(H)`, with no extra unknown
 (F-014).
 - **General surface flux law.** Because the ledger makes conservation independent of how the
diff --git a/docs/design/findings.md b/docs/design/findings.md
index b26a8aa..5614935 100644
--- a/docs/design/findings.md
+++ b/docs/design/findings.md
@@ -75,6 +75,7 @@ Seeded by the Stage A side-by-side diff (2026-10-11). Sources compared:
 | F-053 | Renames | PyTES had no way to deprecate a name: renames were clean breaks (F-036, F-050), and the only record of an old name was `BaseModel.renamed_options` (F-050), for options. Issue #17 adds `pytes.deprecation`: one registry of renames (`RENAMES`), `PyTESDeprecationWarning`, a module `__getattr__` for renamed classes and functions, and the option hook in `process_options`. `BaseModel.renamed_options` (unreleased) is replaced by the registry entry for F-050, whose error message is unchanged apart from the model name prefix. Option values are deferred (none renamed yet; the issue asked for them), and so are variable and parameter names: PyBaMM looks them up itself. No result changes. | Contract change (DESIGN D3, D17, maintainer, 2026-10-23). | Done; tests in `tests/toolkit/test_deprecation.py`; pytest turns the warning into an error; `CONTRIBUTING.md` "Renaming and deprecating". | closed |
 | F-054 | Names: steps, couplings | Names that said less than, or something other than, what the code does (issue #18, maintainer, 2026-10-23). `step.ConstantInputs` becomes `step.Inlet` and `step.ConstantPower` becomes `step.Power`: their values need not be constant, and the names now follow `pybamm.step` (named after what is held). New `step.Rest`: no flow, with the inlet temperature given or, by default, held where the previous step left it (the initial temperature for a first step), through a new solver input the runner sets from the previous step; until now a rest was `ConstantInputs(T_in, 0)`. `LedgerWallCoupling` becomes `WallCoupling` ("ledger" names how energy is accounted, not the physics) and `FinPCMCoupling` becomes `FinStorageCoupling` (the material between the fins need not change phase); its default submodel name becomes "Fin-storage coupling", and its keyword arguments `pcm_label` and `pcm_geometry` are kept, as the registry cannot deprecate keyword arguments. The old names are aliases through the registry (F-053): they warn until 0.2.0. Not changed: no `mass_flow_rate_below` termination (backwards flow is not supported), `Target`, `ConstantPropertyLiquid`. No result changes (the reference experiment is unchanged). The rename test for user-facing files no longer reads the changelog's copy in `docs/`, which keeps the names of its time. | Contract flaw (names). | Renamed everywhere, aliases in `pytes.deprecation.RENAMES`; DESIGN D10 and the contract updated; tests for `Rest`. | closed |
 | F-055 | Variables | Outputs for comparison with measurements (issue #14, maintainer's choices 2026-10-24). Required in every model: `Volume flow rate [m3.s-1]` (`mdot / rho(h_in)`, twin `[l.min-1]`), `Pumped volume [m3]` (its time integral, a solver state, twin `[l]`), `Inlet-outlet temperature difference [K]` (`T_in - T_out`, the example name of D12, which no model provided), `Heat transfer fluid velocity [m.s-1]` (mean velocity in the flow area; interstitial in a packed bed), `Charged energy [J]` (`-Discharged energy`), and `Discharged throughput [J]` and `Charged throughput [J]` (time integrals of `max(P, 0)` and `max(-P, 0)`, solver states). Every storage provides its volume-weighted averages, `<label> average temperature [K]` and, along the flow path, `<label> cross-section average temperature [K]`; and `Storage average temperature [K]` and `Storage cross-section average temperature [K]` (required) combine all the storage, volume-weighted, so that the mush model (fins and PCM apart) compares with the sharp front (homogenised, local equilibrium). One name in every family (the maintainer weighed "capsule average" for the bed and kept the shared name; in a packed bed the cross-section average is the capsule's volume average, as the docs say). This reverses F-043's dropping of `Mean storage temperature [K]` as undefined: it is now defined. `Surface Reynolds number` and `Surface Prandtl number` are provided with the constant coefficient too: `NewtonCooling` takes an optional flow geometry for reporting only. `Pressure drop [Pa]` stays reserved. Side effects: the three new solver states move the reference experiment within the solver tolerance (at most 1.5e-5 relative), so `tests/data/reference_regression.npz` is regenerated; and `check_unit_twins` compares to relative 1e-10 instead of 1e-12, as the `[l.min-1]` twin differs from the converted value by a few round-offs (2.6e-12), depending on the order in which the solver evaluates it. No physical result changes. | Contract flaw (additions; a tolerance). | Implemented; `tests/models/test_output_variables.py`; contract §3 and DESIGN D12 updated. | closed |
+| F-056 | Storage, coupling | Packed-bed fluid temperatures were not monotonic during a charge: at every height the fluid rose, fell by 1.2 K to 2.1 K, then rose again (default Nallusamy charge, default mesh), seen in the tutorials' plots. The storage's wall condition turned the flux law's heat flux into a gradient with the conductivity at the extrapolated wall enthalpy, `k(H_surf)`, while the heat that enters the storage is `k grad(T)` on the wall face, where PyBaMM extrapolates the node values of `k`. The two agree when `k` varies smoothly; the default paraffin melts at a sharp point (melting range 0 K) and its conductivity drops from 0.4 to 0.15 W.m-1.K-1, so while the cells by the wall melted the storage took 19 % to 29 % of the flux law's heat (0.15 of the bed height, 12 to 16 min). Energy stayed exactly conserved (the ledger, D6, telescopes to the actual wall flux), so the conservation tests could not see it, and the coupling's docstring wrongly said the two "differ slightly". Evidence that it was numerics, not physics: the dips grew as the capsule mesh was refined (0.5 K on 10 cells to 3 K on 160), vanished with equal solid and liquid conductivities, and vanished with a melting range of 2 K. | Model flaw (numerics). | The wall condition uses `boundary_value(k)`, the conductivity `k grad(T)` uses on the wall face, so exactly the flux law's flux enters (`EnthalpyStorage`). The charge is now monotonic, converging as the capsule mesh is refined (largest dip 0.17 K on 20 cells, 0.04 K on 40, none from 80); energy is conserved as before (3e-8 J of 8.3 MJ). Results move: the Nallusamy (2007) fluid temperatures fit the measurements better at every height (RMS 0.98, 0.68, 1.41 and 1.85 K at a quarter, half, three quarters and the full height, for 1.09, 0.92, 1.63 and 2.05 K); the finned tubes, whose PCM melts over 2 K with conductivities 0.6 and 0.5, move by at most 0.006 K. Regression reference regenerated. Tests: the fluid's heat equals the flux law's flux over the capsule surface, and the default charge heats the fluid monotonically at four heights. | closed |
 
 **Pinning note.** HB was pinned at `d391da0` on 2026-10-11 (the working tree diffed above),
 and PB at `023bf58`.
diff --git a/src/pytes/submodels/coupling.py b/src/pytes/submodels/coupling.py
index 050be50..0683ecc 100644
--- a/src/pytes/submodels/coupling.py
+++ b/src/pytes/submodels/coupling.py
@@ -34,9 +34,9 @@ class WallCoupling(BaseSubModel):
     boundaries, computed from the storage's own discrete energy balance: minus the integral
     of its conduction heating rate over the storage cross-section (the ledger, DESIGN D6).
     This telescopes to the discrete flux at the wall, so energy is conserved exactly,
-    whatever the flux law or the conductivity. (The flux times the perimeter differs
-    slightly, as the storage flux uses the conductivity at the cell face rather than at the
-    surface.)
+    whatever the flux law or the conductivity. It is also the flux times the wall area: the
+    storage's wall condition uses the conductivity on the wall face, as its conduction does
+    (F-056).
 
     The coupling also provides the flux law's own variables (e.g. "Heat transfer
     coefficient [W.m-2.K-1]", see :meth:`pytes.submodels.SurfaceFluxLaw.variables`), and
diff --git a/src/pytes/submodels/storage.py b/src/pytes/submodels/storage.py
index de8a811..1d57148 100644
--- a/src/pytes/submodels/storage.py
+++ b/src/pytes/submodels/storage.py
@@ -268,14 +268,17 @@ def build(self):
             dHdt = dHdt + self.coupled(
                 f"{label} heat source [W.m-3]", geometry.primary_domain
             )
-        H_surf = pybamm.boundary_value(H, geometry.inner)
 
         self.rhs = {H: dHdt}
         if self.inner_condition == "heat flux":
             heat_flux = self.coupled(
                 f"{label} boundary heat flux [W.m-2]", geometry.boundary_domain
             )
-            k_surf = material.k(H_surf, self.internal_states)
+            # The conductivity that ``k * grad(T)`` uses on the wall face, the node values
+            # extrapolated to the boundary, so that exactly ``heat_flux`` enters. The
+            # conductivity at the extrapolated enthalpy differs from it while the cells
+            # by the wall change phase, by a factor of 5 in the default packed bed (F-056)
+            k_surf = pybamm.boundary_value(k, geometry.inner)
             inner = (geometry.gradient_from_heat_flux(heat_flux, k_surf), "Neumann")
         else:
             T_b = self.coupled(
diff --git a/tests/data/reference_regression.npz b/tests/data/reference_regression.npz
index e53efc4..e44cac3 100644
Binary files a/tests/data/reference_regression.npz and b/tests/data/reference_regression.npz differ
diff --git a/tests/models/test_packed_bed.py b/tests/models/test_packed_bed.py
index f586b75..8f76cfa 100644
--- a/tests/models/test_packed_bed.py
+++ b/tests/models/test_packed_bed.py
@@ -57,6 +57,43 @@ def test_per_area_totals(charge):
         )
 
 
+def _capsule_surface_per_length():
+    p = get_parameter_values()
+    eps, R = p["Bed porosity"], p["Capsule radius [m]"]
+    return 3 * (1 - eps) / R * p["Bed cross-sectional area [m2]"]
+
+
+def test_fluid_takes_the_flux_law_heat():
+    # The heat the fluid gives up (the ledger, DESIGN D6) is the surface flux law's flux
+    # over the capsule surface, also while the capsule surface melts: the default PCM melts
+    # at a sharp point and its conductivity drops from 0.4 to 0.15 W.m-1.K-1 (F-056)
+    model = pytes.PackedBedModel()
+    t = np.linspace(0, 1800, 61)
+    solution = pybamm.Simulation(
+        model, var_pts=var_pts(model, {"r": 20, "z": 10})
+    ).solve([0, t[-1]], t_interp=t)
+    flux = solution["Phase-change material boundary heat flux [W.m-2]"].entries
+    gain = solution["Heat transfer fluid heat gain per unit length [W.m-1]"].entries
+    expected = -flux * _capsule_surface_per_length()
+    np.testing.assert_allclose(
+        gain, expected, rtol=1e-6, atol=1e-6 * np.abs(gain).max()
+    )
+
+
+def test_charge_heats_the_fluid_monotonically():
+    # With a constant inlet temperature, the fluid at every height only warms. A storage
+    # that took less than the flux law's heat while its surface melted made the fluid
+    # overshoot, then cool by 1 K to 2 K (F-056)
+    model = pytes.PackedBedModel()
+    t = np.linspace(0, 160 * 60, 321)
+    solution = pybamm.Simulation(model).solve([0, t[-1]], t_interp=t)
+    fluid = solution["Heat transfer fluid temperature [K]"]
+    length = get_parameter_values()["Bed length [m]"]
+    for fraction in [0.25, 0.5, 0.75, 1]:
+        T = fluid(t=t, z=fraction * length)
+        assert np.max(np.maximum.accumulate(T) - T) < 0.1, fraction
+
+
 def test_nallusamy_flow_is_about_2_litres_per_minute():
     values = get_parameter_values("Nallusamy2007")
     litres_per_minute = values["Mass flow rate [kg.s-1]"] / 1000 * 6e4