diff --git a/packages/loop_interpolation/src/loop_interpolation/_discrete_interpolator.py b/packages/loop_interpolation/src/loop_interpolation/_discrete_interpolator.py
index a51c1a3e3..699fcd17f 100644
--- a/packages/loop_interpolation/src/loop_interpolation/_discrete_interpolator.py
+++ b/packages/loop_interpolation/src/loop_interpolation/_discrete_interpolator.py
@@ -1008,6 +1008,9 @@ def _run_constant_norm_polish(
             polish_solver_kwargs = dict(stable_solver_kwargs)
             polish_solver_kwargs["x0"] = lambda _support, x0=x0_seed: np.array(x0, copy=True)
             polish_solver_kwargs["constant_norm_iterations"] = 0
+            # warn about non-convergence once per solve: only the final polish iteration reports
+            # (the earlier ones are intermediate, warm-started steps)
+            polish_solver_kwargs["warn_on_nonconvergence"] = i == int(iterations) - 1
 
             solved = self.solve_system("admm", solver_kwargs=polish_solver_kwargs)
             if not solved:
@@ -1232,6 +1235,11 @@ def _solve_with_admm(
             logger.warning("No inequality constraints, using lsmr")
             return self.solve_system("lsmr", solver_kwargs=solver_kwargs)
 
+        if constant_norm_iterations > 0 and constant_norm_weight > 0.0:
+            # the polish re-solves follow and the last one reports non-convergence for the
+            # whole solve, so this first solve stays quiet
+            solver_kwargs = {**solver_kwargs, "warn_on_nonconvergence": False}
+
         try:
             c, history, ok = _solver_strategy.solve_with_admm(
                 A,
@@ -1396,6 +1404,11 @@ def solve_system(
             if x0 is not None and scaling_matrix is not None:
                 # the solver works on y where c = S @ y
                 x0 = x0 / scaling_matrix.diagonal()
+            if scaling_matrix is not None and Q.shape[0] > 0:
+                # the solver works on y where c = S @ y, so the inequality rows must act on y
+                # too: Q @ c = (Q @ S) @ y. Without this the bounds are applied to y, not to the
+                # field.
+                Q = (Q @ scaling_matrix).tocsr()
 
             if callable(solver_choice):
                 self.up_to_date = self._solve_with_callable(solver_choice, A, b, timing)
diff --git a/packages/loop_interpolation/src/loop_interpolation/_solver_strategy.py b/packages/loop_interpolation/src/loop_interpolation/_solver_strategy.py
index a9f447f17..e34dfa786 100644
--- a/packages/loop_interpolation/src/loop_interpolation/_solver_strategy.py
+++ b/packages/loop_interpolation/src/loop_interpolation/_solver_strategy.py
@@ -183,6 +183,7 @@ def extract_admm_kwargs(solver_kwargs: dict, admm_solve) -> tuple[str, dict]:
         "inner_maxiter_schedule": solver_kwargs.pop("inner_maxiter_schedule", "linear"),
         "model_update_tol": solver_kwargs.pop("model_update_tol", 0.0),
         "return_history": solver_kwargs.pop("return_history", False),
+        "warn_on_nonconvergence": solver_kwargs.pop("warn_on_nonconvergence", True),
     }
     supported_optional = set(inspect.signature(admm_solve).parameters.keys())
     admm_kwargs = {k: v for k, v in admm_kwargs.items() if k in supported_optional}
diff --git a/packages/loop_interpolation/src/loop_interpolation/loopsolver/admm_solver.py b/packages/loop_interpolation/src/loop_interpolation/loopsolver/admm_solver.py
index e033ea852..0ab484597 100644
--- a/packages/loop_interpolation/src/loop_interpolation/loopsolver/admm_solver.py
+++ b/packages/loop_interpolation/src/loop_interpolation/loopsolver/admm_solver.py
@@ -5,11 +5,14 @@
 from dataclasses import dataclass
 
 import numpy as np
+from loop_common.logging import get_logger
 from scipy.sparse import csr_matrix, diags, vstack
 from scipy.sparse.linalg import LinearOperator, cg, lsmr, lsqr
 
 from .admm_method import ADMM
 
+logger = get_logger(__name__)
+
 
 @dataclass
 class Config:
@@ -95,6 +98,7 @@ def admm_solve(
     model_update_tol: float = 0.0,
     min_iterations: int = 5,
     return_history: bool = False,
+    warn_on_nonconvergence: bool = True,
 ):
     if A.shape[1] != x0.shape[0]:
         raise ValueError("Number of columns in interpolation matrix does not match x0")
@@ -320,6 +324,11 @@ def _solve_linsys(system_matrix, rhs, kwargs):
         )
 
     history = []
+    # how the major loop ended: it converged, it stopped because the model stopped changing (the
+    # caller asked for that with model_update_tol), or it used all of nmajor iterations
+    converged = False
+    stopped_on_model_update = False
+    last_residuals = None
     inner_x0 = None
     active_rows_prev = None
     use_cached_normal_eq = solver_name == "cg" and not active_set
@@ -415,6 +424,8 @@ def _solve_linsys(system_matrix, rhs, kwargs):
             eps_pri = float(n_ie_sqrt * admm_abs_tol + admm_rel_tol * max(qx_norm, z_norm))
             eps_dual = float(n_ie_sqrt * admm_abs_tol + admm_rel_tol * u_norm)
 
+            last_residuals = (primal_norm, eps_pri, dual_norm, eps_dual)
+
             if return_history:
                 history.append(
                     {
@@ -431,6 +442,7 @@ def _solve_linsys(system_matrix, rhs, kwargs):
                 )
 
             if (_i + 1) >= min_iterations and primal_norm <= eps_pri and dual_norm <= eps_dual:
+                converged = True
                 break
 
             if model_update_tol > 0.0 and (_i + 1) >= min_iterations:
@@ -440,6 +452,7 @@ def _solve_linsys(system_matrix, rhs, kwargs):
                         history[-1]["stopped_on_model_update"] = True
                         history[-1]["dx_norm"] = dx_norm
                         history[-1]["model_norm"] = model_norm_before
+                    stopped_on_model_update = True
                     break
 
             if adaptive_rho_enabled and (_i + 1) >= min_iterations:
@@ -458,6 +471,32 @@ def _solve_linsys(system_matrix, rhs, kwargs):
             if dx_norm <= model_update_tol * model_scale:
                 break
 
+    if (
+        warn_on_nonconvergence
+        and Q.shape[0] > 0
+        and not converged
+        and not stopped_on_model_update
+        and last_residuals
+    ):
+        primal_norm, eps_pri, dual_norm, eps_dual = last_residuals
+        if primal_norm > eps_pri:
+            consequence = "the inequality constraints may not be fully satisfied"
+            advice = (
+                f"Try a larger admm_weight (the final value used was {rho:g}) or a larger nmajor."
+            )
+        else:
+            # the bounds hold to tolerance; only the dual residual (optimality) missed
+            consequence = (
+                "the inequality constraints are satisfied to tolerance but the solution may "
+                "not be fully converged"
+            )
+            advice = f"Try a larger nmajor (the final admm_weight was {rho:g})."
+        logger.warning(
+            f"ADMM used all {nmajor} iterations without converging: primal residual "
+            f"{primal_norm:.3g} (tolerance {eps_pri:.3g}), dual residual {dual_norm:.3g} "
+            f"(tolerance {eps_dual:.3g}); {consequence}. {advice}"
+        )
+
     if return_history:
         return model, history
     return model
diff --git a/packages/loop_interpolation/tests/test_inequality_enforcement.py b/packages/loop_interpolation/tests/test_inequality_enforcement.py
new file mode 100644
index 000000000..1b6c1f80a
--- /dev/null
+++ b/packages/loop_interpolation/tests/test_inequality_enforcement.py
@@ -0,0 +1,256 @@
+"""Inequality constraints must be enforced, or at least not silently ignored.
+
+``InterpolatorBuilder.add_inequality_constraints`` takes rows ``x, y, z, lower, upper`` and the
+documented meaning is ``lower <= f(x, y, z) <= upper``. Inequality constraints are kept apart
+from the equality system and are only used by the ADMM solver.
+
+The scenario is a sparse model on a gaussian dome whose true field is known: one borehole (three
+contact points on the 0.3, 0.5 and 0.7 horizons) and two orientations. Thirty random points whose
+true value is clearly above or below the 0.5 horizon become one-sided bounds on it:
+``0.5 <= f`` above, ``f <= 0.5`` below. Without the bounds this data leaves 4 of the 24 bounded
+points on the wrong side of 0.5.
+
+Diagnosis:
+
+* ``solve_system`` column-scales the equality system (``apply_scaling_matrix``, on by default) and
+  solves for ``y`` where ``c = S @ y``. The inequality matrix ``Q`` was passed to ADMM unscaled,
+  so the bounds were applied to ``y`` instead of to the field ``Q @ c``. With scaling on, even a
+  large penalty left bounds violated (``admm_weight=10``: 10 of 24). Fixed by passing ``Q @ S``.
+* With that fixed ADMM behaves like ordinary ADMM. Its default penalty ``admm_weight=0.01`` is too
+  weak for this problem: after the default 200 outer iterations the primal residual is 0.056
+  against a tolerance of 0.0025 (still falling after 5000 iterations: 0.0145), so 4 of 24 bounds
+  are still violated. ``admm_weight=0.1`` stops at iteration 150, ``admm_weight=1`` at iteration
+  90 with every bound satisfied. The solver reaches ``nmajor`` without meeting its tolerance and
+  returns without a warning.
+* ``nmajor=10, maxiter=50`` (the arguments of the existing builder test) enforces the bounds
+  only by accident: it is an unconverged, truncated solve, and ``maxiter=100`` with the same
+  ``nmajor`` violates them again.
+* With the ``cg`` solver the inequality data have no effect, and nothing is logged.
+
+The default-penalty and ``cg`` cases are strict xfails. When they start to pass they XPASS, which
+fails the run: remove the markers then.
+"""
+import logging
+
+import numpy as np
+import pytest
+from loop_common.geometry import BoundingBox
+from loop_interpolation import InterpolatorBuilder
+from scipy.optimize import brentq
+
+TOLERANCE = 1e-3
+LEVEL = 0.5
+
+
+def truth(xyz):
+    r2 = (xyz[:, 0] - 0.5) ** 2 + (xyz[:, 1] - 0.5) ** 2
+    return xyz[:, 2] - 0.2 * np.exp(-r2 / (2 * 0.2**2))
+
+
+def truth_gradient(xyz):
+    r2 = (xyz[:, 0] - 0.5) ** 2 + (xyz[:, 1] - 0.5) ** 2
+    e = 0.2 * np.exp(-r2 / (2 * 0.2**2)) / 0.2**2
+    return np.c_[e * (xyz[:, 0] - 0.5), e * (xyz[:, 1] - 0.5), np.ones(len(xyz))]
+
+
+def sparse_data():
+    contacts = []
+    for level in (0.3, 0.5, 0.7):
+        z = brentq(lambda z, level=level: truth(np.array([[0.5, 0.5, z]]))[0] - level, 0.0, 1.0)
+        contacts.append([0.5, 0.5, z, level])
+    orientation_xyz = np.array([[0.4, 0.4, 0.5], [0.6, 0.6, 0.5]])
+    normals = np.c_[orientation_xyz, truth_gradient(orientation_xyz)]
+
+    rng = np.random.default_rng(9)
+    candidates = rng.uniform(0.1, 0.9, (30, 3))
+    rows = []
+    for point, value in zip(candidates, truth(candidates)):
+        if value > LEVEL + 0.05:
+            rows.append([*point, LEVEL, 10.0])  # f >= 0.5
+        elif value < LEVEL - 0.05:
+            rows.append([*point, -10.0, LEVEL])  # f <= 0.5
+    return np.array(contacts), normals, np.array(rows)
+
+
+def builder(with_inequality=True):
+    contacts, normals, inequality = sparse_data()
+    bounding_box = BoundingBox(np.array([0.0, 0.0, 0.0]), np.array([1.0, 1.0, 1.0]))
+    built = InterpolatorBuilder(
+        "FDI", bounding_box=bounding_box, nelements=4000, buffer=0.1
+    )
+    built.add_value_constraints(contacts).add_normal_constraints(normals)
+    if with_inequality:
+        built.add_inequality_constraints(inequality)
+    return built
+
+
+def violations(interpolator):
+    """How far each bounded point is outside its bounds (0 where satisfied)."""
+    _, _, inequality = sparse_data()
+    value = interpolator.evaluate_value(inequality[:, :3])
+    below = np.maximum(inequality[:, 3] - value, 0.0)
+    above = np.maximum(value - inequality[:, 4], 0.0)
+    return below + above
+
+
+def test_the_unconstrained_solution_violates_some_bounds():
+    """The scenario is only meaningful if the data alone get some bounded points wrong."""
+    interpolator = builder(with_inequality=False).setup_interpolator().solve().build()
+    assert np.count_nonzero(violations(interpolator) > TOLERANCE) >= 3
+
+
+def test_admm_with_the_arguments_of_the_existing_builder_test_enforces_the_bounds():
+    interpolator = (
+        builder()
+        .setup_interpolator()
+        .use_solver("admm", nmajor=10, admm_weight=0.01, maxiter=50)
+        .solve()
+        .build()
+    )
+    assert violations(interpolator).max() <= TOLERANCE
+
+
+@pytest.mark.xfail(
+    strict=True,
+    reason="the default penalty admm_weight=0.01 is too weak here: after 200 iterations the primal "
+    "residual (0.056) is far above the tolerance (0.0025) and 4 of 24 bounds are violated (the "
+    "solver now warns about this, see the non-convergence tests below)",
+)
+def test_admm_with_the_default_penalty_enforces_the_bounds():
+    interpolator = builder().setup_interpolator().use_solver("admm").solve().build()
+    assert violations(interpolator).max() <= TOLERANCE
+
+
+@pytest.mark.parametrize("apply_scaling_matrix", [True, False])
+@pytest.mark.parametrize("admm_weight", [1.0, 10.0])
+def test_admm_enforces_the_bounds_with_and_without_column_scaling(apply_scaling_matrix, admm_weight):
+    """The inequality rows must be scaled with the equality system (they were not)."""
+    built = builder()
+    built.interpolator.apply_scaling_matrix = apply_scaling_matrix
+    interpolator = built.setup_interpolator().use_solver("admm", admm_weight=admm_weight).solve().build()
+    assert violations(interpolator).max() <= TOLERANCE
+
+
+@pytest.mark.xfail(
+    strict=True,
+    reason="expectation, not a documented contract: with solver 'cg' the inequality data are "
+    "dropped without a warning",
+)
+def test_cg_with_inequality_constraints_warns_that_they_are_not_used(caplog):
+    with caplog.at_level(logging.WARNING):
+        builder().setup_interpolator().use_solver("cg").solve().build()
+    assert any("inequality" in record.getMessage().lower() for record in caplog.records)
+
+
+# ------------------------------------------------------------------ non-convergence is reported
+def not_converged_warnings(caplog):
+    return [r for r in caplog.records if "without converging" in r.getMessage()]
+
+
+def test_the_default_penalty_that_does_not_converge_says_so(caplog):
+    with caplog.at_level(logging.WARNING):
+        builder().setup_interpolator().use_solver("admm").solve().build()
+    (warning,) = not_converged_warnings(caplog)
+    message = warning.getMessage()
+    assert "ADMM used all 200 iterations" in message
+    assert "primal residual" in message and "tolerance" in message
+    assert "admm_weight" in message
+    # the primal residual missed its tolerance here, so the bounds really may be violated
+    assert "may not be fully satisfied" in message
+
+
+def test_the_constant_norm_polish_warns_once_not_once_per_iteration(caplog):
+    """The polish re-solves with ADMM each iteration; the warning is for the whole solve."""
+    with caplog.at_level(logging.WARNING):
+        (
+            builder()
+            .setup_interpolator()
+            .use_solver("admm", constant_norm_iterations=3, constant_norm_weight=0.1)
+            .solve()
+            .build()
+        )
+    assert len(not_converged_warnings(caplog)) == 1
+
+
+def test_a_penalty_that_converges_does_not_warn(caplog):
+    with caplog.at_level(logging.WARNING):
+        interpolator = (
+            builder().setup_interpolator().use_solver("admm", admm_weight=1.0).solve().build()
+        )
+    assert not not_converged_warnings(caplog)
+    assert violations(interpolator).max() <= TOLERANCE
+
+
+def small_problem():
+    """Data say the unknowns sum to 3; the first two are bounded to [0, 0.2]."""
+    from scipy import sparse
+
+    n = 6
+    a = sparse.vstack(
+        [sparse.identity(n, format="csr") * 0.2, sparse.csr_matrix(np.ones((1, n)))]
+    ).tocsr()
+    b = np.r_[np.zeros(n), 3.0]
+    q = sparse.csr_matrix(np.eye(n)[:2])
+    bounds = np.array([[0.0, 0.2, 1.0], [0.0, 0.2, 1.0]])
+    return a, b, q, bounds, np.zeros(n)
+
+
+def solve_small(**kwargs):
+    from loop_interpolation.loopsolver import admm_solve
+
+    a, b, q, bounds, x0 = small_problem()
+    return admm_solve(a, b, q, bounds, x0, linsys_solver="lsmr", **kwargs)
+
+
+def test_solver_converged_in_time_does_not_warn(caplog):
+    with caplog.at_level(logging.WARNING):
+        x = solve_small(admm_weight=1.0, nmajor=200)
+    assert not not_converged_warnings(caplog)
+    np.testing.assert_allclose(x[:2], 0.2, atol=1e-3)
+
+
+def test_solver_out_of_iterations_warns_and_reports_the_residuals(caplog):
+    with caplog.at_level(logging.WARNING):
+        x = solve_small(admm_weight=1e-4, nmajor=20)
+    (warning,) = not_converged_warnings(caplog)
+    assert "all 20 iterations" in warning.getMessage()
+    assert "the final value used was 0.0001" in warning.getMessage()
+    # and the bounds really are violated, so the warning is true
+    assert x[:2].max() > 0.2 + 0.1
+
+
+def test_dual_only_miss_does_not_claim_the_bounds_are_violated(caplog):
+    """admm_weight=1 with 20 iterations meets the primal tolerance but not the dual one."""
+    with caplog.at_level(logging.WARNING):
+        interpolator = (
+            builder()
+            .setup_interpolator()
+            .use_solver("admm", admm_weight=1.0, nmajor=20)
+            .solve()
+            .build()
+        )
+    (warning,) = not_converged_warnings(caplog)
+    message = warning.getMessage()
+    assert "satisfied to tolerance" in message
+    assert "may not be fully satisfied" not in message
+    assert violations(interpolator).max() <= TOLERANCE  # and the message is true
+
+
+def test_warning_can_be_silenced_per_solve(caplog):
+    with caplog.at_level(logging.WARNING):
+        solve_small(admm_weight=1e-4, nmajor=20, warn_on_nonconvergence=False)
+    assert not not_converged_warnings(caplog)
+
+
+def test_solver_stopped_by_model_update_tolerance_does_not_warn(caplog):
+    """Stopping because the model stopped changing was asked for; it is not a failure."""
+    with caplog.at_level(logging.WARNING):
+        solve_small(admm_weight=1e-4, nmajor=20, model_update_tol=10.0)
+    assert not not_converged_warnings(caplog)
+
+
+def test_a_solver_with_fewer_iterations_than_the_minimum_warns(caplog):
+    """It cannot be declared converged before min_iterations, so it used all of its iterations."""
+    with caplog.at_level(logging.WARNING):
+        solve_small(admm_weight=1.0, nmajor=3)
+    assert len(not_converged_warnings(caplog)) == 1
diff --git a/tests/unit/modelling/test_recipe_reproduces_field.py b/tests/unit/modelling/test_recipe_reproduces_field.py
new file mode 100644
index 000000000..27ecac309
--- /dev/null
+++ b/tests/unit/modelling/test_recipe_reproduces_field.py
@@ -0,0 +1,64 @@
+"""A model recipe should rebuild a model that gives the same scalar field.
+
+``GeologicalModel.to_recipe_json`` / ``from_recipe_json`` are documented as a way to describe and
+rebuild a model. The existing recipe tests check the names of features, the data and the
+stratigraphic column, but not the field, and the field depends on how each feature was built:
+the number of elements, the interpolator type and the regularisation weight.
+
+Those build parameters are not stored: a feature entry in the recipe has only ``name``, ``kind``,
+``faults`` and ``regions``. A model built with the defaults rebuilds exactly, because the rebuild
+uses the defaults. A model built with anything else rebuilds with the defaults and gives a
+different field (on the Claudius data, differences of 12 to 62 on a field of range about 560).
+
+The non-default cases are strict xfails. When the recipe stores the build parameters they will
+XPASS, which fails the run: remove the markers then.
+"""
+import numpy as np
+import pytest
+
+from LoopStructural import GeologicalModel
+from LoopStructural.datasets import load_claudius
+
+RECIPE_DROPS_BUILD_PARAMETERS = pytest.mark.xfail(
+    strict=True,
+    reason="to_recipe_json stores only name, kind, faults and regions for a feature, not the "
+    "build parameters (nelements, interpolatortype, regularisation, ...), so the rebuilt "
+    "model uses the defaults",
+)
+
+CASES = [
+    pytest.param({}, id="defaults"),
+    pytest.param(
+        {"nelements": 2000}, id="nelements", marks=RECIPE_DROPS_BUILD_PARAMETERS
+    ),
+    pytest.param(
+        {"interpolatortype": "PLI"}, id="interpolatortype", marks=RECIPE_DROPS_BUILD_PARAMETERS
+    ),
+    pytest.param(
+        {"regularisation": 10.0}, id="regularisation", marks=RECIPE_DROPS_BUILD_PARAMETERS
+    ),
+]
+
+
+def build_model(**build_parameters):
+    data, bb = load_claudius()
+    model = GeologicalModel(bb[0, :], bb[1, :])
+    model.data = data
+    model.stratigraphic_column.clear(basement=False)
+    model.create_and_add_foliation("strati", **build_parameters)
+    model.update()
+    return model, bb
+
+
+@pytest.mark.parametrize("build_parameters", CASES)
+def test_rebuilt_model_gives_the_same_field(build_parameters):
+    model, bb = build_model(**build_parameters)
+
+    rebuilt = GeologicalModel.from_recipe_json(model.to_recipe_json())
+    rebuilt.update()
+
+    points = np.random.default_rng(0).uniform(bb[0], bb[1], (300, 3))
+    original = model.evaluate_feature_value("strati", points)
+    again = rebuilt.evaluate_feature_value("strati", points)
+    np.testing.assert_allclose(again, original, atol=1e-6)
+