diff --git a/COMPAT.md b/COMPAT.md
index 11ee6262a..5075c1721 100644
--- a/COMPAT.md
+++ b/COMPAT.md
@@ -37,6 +37,7 @@ unchanged, but the new path is preferred going forward.
 |---|---|---|
 | `LoopStructural.modelling.features._feature_converters.add_fold_to_feature` (private module, imported directly by the QGIS plugin) | `GeologicalModel.add_fold_to_feature` | 2026-07-24 |
 | `LoopStructural.modelling.features._feature_converters.convert_feature_to_structural_frame` | `GeologicalModel.convert_feature_to_structural_frame` | 2026-07-24 |
+| `GeologicalModel.to_recipe_json()` / `to_recipe_dict()` / `save_recipe()` silently wrote a lossy recipe for any model: only `name`, `kind`, `faults` for each feature (build parameters dropped, kinds mislabelled), and a model with a fault, unconformity or custom `data=` rebuilt wrongly or not at all | Each foliation entry now stores the parameters it was requested with (`build`), `group_mapping` is stored, and `from_recipe_*` replays them, so a model built with non-default `nelements`, `interpolatortype` or `regularisation` rebuilds to the same field. **A model containing a fault, unconformity or any other non-foliation feature, a feature built with its own `data=`, or a parameter JSON cannot hold, now raises `LoopValueError` from `to_recipe_*` instead of writing a recipe that rebuilds wrongly.** Older recipes (no `build`) still load, with a warning, but an older entry whose kind is not `foliation` is refused. Faults and unconformities are not yet supported | 2026-10-29 |
 | `GeologicalModel.create_and_add_fault(..., faults=[])` | `GeologicalModel.create_and_add_fault(..., faults=None)` (list built internally, same effective default) | 2026-07-27 |
 | `GeologicalModel.create_and_add_intrusion(..., intrusion_frame_parameters={}, geometric_scaling_parameters={})` | `GeologicalModel.create_and_add_intrusion(..., intrusion_frame_parameters=None, geometric_scaling_parameters=None)` (dicts built internally, same effective default) | 2026-07-27 |
 | `GeologicalModel.get_fault_surfaces(faults=[])` | `GeologicalModel.get_fault_surfaces(faults=None)` (list built internally, same effective default) | 2026-07-27 |
diff --git a/LoopStructural/modelling/core/_model_serializer.py b/LoopStructural/modelling/core/_model_serializer.py
index cbf4d1120..5696c5abc 100644
--- a/LoopStructural/modelling/core/_model_serializer.py
+++ b/LoopStructural/modelling/core/_model_serializer.py
@@ -10,9 +10,11 @@
 delegate to the staticmethods here.
 """
 
+import copy
 import json
 import pathlib
 
+import numpy as np
 import pandas as pd
 
 from ...geometry import BoundingBox
@@ -24,19 +26,99 @@
 logger = getLogger(__name__)
 
 
+# Parameters of create_and_add_feature that are not stored with a feature's creation record:
+# ``index`` is the position (the recipe keeps the order), ``faults`` are other features (the recipe
+# keeps their names), and ``data`` is a data frame (see ``custom_data``).
+_NOT_RECORDED = ("index", "faults", "data")
+
+# Feature types a recipe can rebuild from the parameters stored with it.
+_SUPPORTED_FEATURE_TYPES = ("foliation",)
+
+
+def _json_safe(value: object, where: str) -> object:
+    """Return ``value`` as something ``json`` can write, or raise ``LoopValueError``."""
+    if value is None or isinstance(value, (bool, int, float, str)):
+        return value
+    if isinstance(value, np.generic):
+        return value.item()
+    if isinstance(value, np.ndarray):
+        return value.tolist()
+    if isinstance(value, (list, tuple)):
+        return [_json_safe(v, where) for v in value]
+    if isinstance(value, dict):
+        return {str(k): _json_safe(v, f"{where}[{k!r}]") for k, v in value.items()}
+    raise LoopValueError(
+        f"{where} is a {type(value).__name__}, which cannot be written to a recipe. "
+        "Only numbers, strings, booleans, None, lists and dictionaries of these can be stored."
+    )
+
+
 class ModelSerializer:
     @staticmethod
     def feature_recipe_kind(feature):
+        # subclasses first: FaultSegment is a StructuralFrame and UnconformityFeature a
+        # GeologicalFeature
+        if isinstance(feature, FaultSegment):
+            return "fault"
+        if isinstance(feature, UnconformityFeature):
+            return "unconformity"
         if isinstance(feature, GeologicalFeature):
             return "foliation"
         if isinstance(feature, StructuralFrame):
             return "structural_frame"
-        if isinstance(feature, UnconformityFeature):
-            return "unconformity"
-        if isinstance(feature, FaultSegment):
-            return "fault"
         return feature.__class__.__name__.lower()
 
+    @staticmethod
+    def record_creation(feature: object, feature_type: str, params: dict) -> None:
+        """Remember how ``feature`` was created, for the recipe.
+
+        Called by ``GeologicalModel.create_and_add_feature`` with the parameters that were
+        *requested* (for example ``nelements``, not the element count that was realised).
+        """
+        feature._recipe_creation = {
+            "feature_type": feature_type,
+            # deep copy: a list/dict/array the caller changes later must not change the recipe
+            "parameters": {
+                k: copy.deepcopy(v) for k, v in params.items() if k not in _NOT_RECORDED
+            },
+            "custom_data": params.get("data") is not None,
+        }
+
+    @staticmethod
+    def feature_build_entry(feature: object) -> dict:
+        """Return the ``build`` part of a feature's recipe entry.
+
+        Raises
+        ------
+        LoopValueError
+            if the feature cannot be rebuilt from a recipe yet
+        """
+        kind = ModelSerializer.feature_recipe_kind(feature)
+        creation = getattr(feature, "_recipe_creation", None)
+        if creation is None:
+            raise LoopValueError(
+                f"Feature '{feature.name}' ({kind}) cannot be stored in a recipe: there is no "
+                "record of how it was created. Recipes store features created with "
+                "GeologicalModel.create_and_add_feature (or a create_and_add_* method)."
+            )
+        feature_type = creation["feature_type"]
+        if kind != "foliation" or feature_type not in _SUPPORTED_FEATURE_TYPES:
+            raise LoopValueError(
+                f"Feature '{feature.name}' ({kind}, created as '{feature_type}') cannot be "
+                "stored in a recipe yet: recipes currently support foliations only."
+            )
+        if creation["custom_data"]:
+            raise LoopValueError(
+                f"Feature '{feature.name}' was created with its own data=, which a recipe cannot "
+                "store: a recipe rebuilds a feature from the rows of the model data whose "
+                "feature_name matches. Put the data in model.data instead."
+            )
+        parameters = {
+            key: _json_safe(value, f"parameter '{key}' of feature '{feature.name}'")
+            for key, value in creation["parameters"].items()
+        }
+        return {"feature_type": feature_type, "parameters": parameters}
+
     @staticmethod
     def to_dict(model):
         result = {}
@@ -68,6 +150,7 @@ def to_recipe_dict(model, data_reference=None):
                     if getattr(fault, "name", None)
                 ],
                 "regions": [],
+                "build": ModelSerializer.feature_build_entry(feature),
             }
             recipe["model"]["features"].append(feature_entry)
         if data_reference is not None:
@@ -132,10 +215,27 @@ def from_recipe_dict(cls, recipe):
             feature_name = feature_entry.get("name")
             if not isinstance(feature_name, str):
                 raise TypeError("each feature entry must include a string name")
-            feature_data = model.data.loc[model.data["feature_name"] == feature_name].copy()
-            if feature_data.empty:
-                feature_data = None
-            feature = model.create_and_add_foliation(feature_name, data=feature_data)
+            build = feature_entry.get("build")
+            if build is None:
+                # a recipe written before build parameters were stored
+                kind = feature_entry.get("kind", "foliation")
+                if kind != "foliation":
+                    raise LoopValueError(
+                        f"Cannot rebuild feature '{feature_name}' of kind '{kind}' from this "
+                        "recipe: recipes currently support foliations only."
+                    )
+                logger.warning(
+                    f"The recipe has no build parameters for feature '{feature_name}' (it was "
+                    "written by an older version); building it with the default parameters"
+                )
+                build = {"feature_type": "foliation", "parameters": {}}
+            if not isinstance(build, dict) or not isinstance(build.get("parameters", {}), dict):
+                raise TypeError("a feature's build entry must be a dictionary with parameters")
+            feature = model.create_and_add_feature(
+                build.get("feature_type", "foliation"),
+                feature_name,
+                **build.get("parameters", {}),
+            )
             if feature is None:
                 raise ValueError(f"Could not recreate feature '{feature_name}' from recipe")
             feature_map[feature_name] = feature
diff --git a/LoopStructural/modelling/core/geological_model.py b/LoopStructural/modelling/core/geological_model.py
index cabaff0cf..bf844e121 100644
--- a/LoopStructural/modelling/core/geological_model.py
+++ b/LoopStructural/modelling/core/geological_model.py
@@ -708,7 +708,10 @@ def create_and_add_feature(self, feature_type: str, name: str, **params):
         feature : BaseFeature
             the created geological feature, or None if it could not be built
         """
-        return FeatureBuilderRegistry.create(feature_type, self, name, **params)
+        feature = FeatureBuilderRegistry.create(feature_type, self, name, **params)
+        if feature is not None:
+            ModelSerializer.record_creation(feature, feature_type, params)
+        return feature
 
     @public_api(tier="stable")
     def create_and_add_foliation(
diff --git a/LoopStructural/modelling/core/stratigraphic_column.py b/LoopStructural/modelling/core/stratigraphic_column.py
index a2d1971f1..b1f83c049 100644
--- a/LoopStructural/modelling/core/stratigraphic_column.py
+++ b/LoopStructural/modelling/core/stratigraphic_column.py
@@ -641,6 +641,7 @@ def to_dict(self):
         """
         return {
             "elements": [element.to_dict() for element in self.order],
+            "group_mapping": dict(self.group_mapping),
         }
 
     def update_from_dict(self, data):
@@ -658,6 +659,7 @@ def update_from_dict(self, data):
                 else:
                     element = StratigraphicUnit.from_dict(element_data)
                 self.add_element(element)
+            self.group_mapping = dict(data.get("group_mapping", {}))
 
     @classmethod
     def from_dict(cls, data):
@@ -675,6 +677,7 @@ def from_dict(cls, data):
             else:
                 element = StratigraphicUnit.from_dict(element_data)
             column.add_element(element)
+        column.group_mapping = dict(data.get("group_mapping", {}))
         return column
 
     def get_isovalues(self, where: str = 'bottom') -> dict[str, float]:
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..b0c4a018b
--- /dev/null
+++ b/packages/loop_interpolation/tests/test_inequality_enforcement.py
@@ -0,0 +1,122 @@
+"""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.
+
+What was observed (not diagnosed):
+
+* ADMM with ``nmajor=10, admm_weight=0.01, maxiter=50`` (the arguments the existing builder test
+  uses) enforces every bound, and the error against the truth falls.
+* ADMM with its default arguments, or with only some of those arguments, returns essentially the
+  unconstrained solution: 4 of 24 bounds still violated, worst by 0.056.
+* With the ``cg`` solver the inequality data have no effect, and nothing is logged.
+
+The last two 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="ADMM with its default arguments returns essentially the unconstrained solution "
+    "here (4 of 24 bounds violated, worst 0.056); cause not diagnosed",
+)
+def test_admm_with_default_arguments_enforces_the_bounds():
+    interpolator = builder().setup_interpolator().use_solver("admm").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)
diff --git a/tests/unit/modelling/test_recipe_limits.py b/tests/unit/modelling/test_recipe_limits.py
new file mode 100644
index 000000000..fe9c673b6
--- /dev/null
+++ b/tests/unit/modelling/test_recipe_limits.py
@@ -0,0 +1,175 @@
+"""What a recipe can and cannot store, and that it says so instead of rebuilding something wrong."""
+import json
+
+import numpy as np
+import pandas as pd
+import pytest
+
+from LoopStructural import GeologicalModel
+from LoopStructural.datasets import load_claudius
+from LoopStructural.modelling.core._model_serializer import ModelSerializer
+from LoopStructural.utils import LoopValueError
+
+
+def claudius_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", **{"nelements": 1000, **build_parameters})
+    return model, data, bb
+
+
+def fault_model():
+    data = pd.DataFrame(
+        [
+            [5, 5, 5, 0, 0.70710678, 0.0, 0.70710678, 0, "fault"],
+            [5, 5, 5, 0, -0.70710678, 0.0, 0.70710678, 1, "fault"],
+            [8, 5, 5, 0, 0, 0, 1, np.nan, "strati"],
+        ],
+        columns=["X", "Y", "Z", "val", "nx", "ny", "nz", "coord", "feature_name"],
+    )
+    model = GeologicalModel([0, 0, 0], [10, 10, 10])
+    model.data = data
+    model.create_and_add_fault(
+        "fault", 1, nelements=1000, buffer=0.5, major_axis=10, minor_axis=3, intermediate_axis=10
+    )
+    model.create_and_add_foliation("strati", nelements=1000, faults=[model["fault"]])
+    return model
+
+
+# ------------------------------------------------------------------ what is stored
+def test_the_requested_parameters_are_stored_not_the_realised_ones():
+    model, _, _ = claudius_model(nelements=2000, regularisation=10.0)
+    entry = json.loads(model.to_recipe_json())["model"]["features"][0]
+    parameters = entry["build"]["parameters"]
+    assert entry["build"]["feature_type"] == "foliation"
+    assert parameters["nelements"] == 2000
+    assert parameters["regularisation"] == 10.0
+    assert parameters["interpolatortype"] == "FDI"
+    # the realised element count is different, and storing it would not reproduce the grid
+    model.update()
+    assert model["strati"].builder.build_arguments["nelements"] != 2000
+
+
+def test_data_and_other_features_are_not_stored_as_parameters():
+    model, _, _ = claudius_model()
+    parameters = json.loads(model.to_recipe_json())["model"]["features"][0]["build"]["parameters"]
+    assert not {"data", "faults", "index"} & set(parameters)
+
+
+def test_group_mapping_survives_a_round_trip():
+    model, _, _ = claudius_model()
+    model.stratigraphic_column.add_unit("unit_0", thickness=60, id=0)
+    model.stratigraphic_column.group_mapping["Group_0"] = "strati"
+    rebuilt = GeologicalModel.from_recipe_json(model.to_recipe_json())
+    assert rebuilt.stratigraphic_column.group_mapping == {"Group_0": "strati"}
+
+
+def test_a_column_without_a_group_mapping_still_loads():
+    """Recipes written before group_mapping was stored have no such key."""
+    model, _, _ = claudius_model()
+    model.stratigraphic_column.add_unit("unit_0", thickness=60, id=0)
+    recipe = model.to_recipe_dict()
+    del recipe["model"]["stratigraphic_column"]["group_mapping"]
+    rebuilt = GeologicalModel.from_recipe_dict(recipe)
+    assert rebuilt.stratigraphic_column.group_mapping == {}
+
+
+# ------------------------------------------------------------------ kind detection
+def test_a_fault_is_recognised_as_a_fault():
+    model = fault_model()
+    assert ModelSerializer.feature_recipe_kind(model["fault"]) == "fault"
+    assert ModelSerializer.feature_recipe_kind(model["strati"]) == "foliation"
+
+
+def test_an_unconformity_is_recognised_as_an_unconformity():
+    model, _, _ = claudius_model()
+    model.add_unconformity(model["strati"], 0)
+    kinds = {f.name: ModelSerializer.feature_recipe_kind(f) for f in model.features}
+    assert kinds["__strati_unconformity"] == "unconformity"
+
+
+# ------------------------------------------------------------------ what is refused
+def test_a_model_with_a_fault_is_refused_not_rebuilt_without_the_fault():
+    with pytest.raises(LoopValueError, match=r"fault.*foliations only"):
+        fault_model().to_recipe_json()
+
+
+def test_a_model_with_an_unconformity_is_refused():
+    model, _, _ = claudius_model()
+    model.add_unconformity(model["strati"], 0)
+    with pytest.raises(LoopValueError, match=r"__strati_unconformity.*unconformity"):
+        model.to_recipe_json()
+
+
+def test_a_feature_built_with_its_own_data_is_refused():
+    model, data, _ = claudius_model()
+    model.create_and_add_foliation(
+        "strati_2", data=data[data["feature_name"] == "strati"].iloc[::2], nelements=1000
+    )
+    with pytest.raises(LoopValueError, match=r"strati_2.*its own data"):
+        model.to_recipe_json()
+
+
+def test_a_parameter_json_cannot_hold_is_refused():
+    model, _, _ = claudius_model(solver=lambda a, b: None)
+    with pytest.raises(LoopValueError, match=r"parameter 'solver' of feature 'strati'"):
+        model.to_recipe_json()
+
+
+def test_a_feature_with_no_creation_record_is_refused():
+    model, _, _ = claudius_model()
+    del model["strati"]._recipe_creation
+    with pytest.raises(LoopValueError, match="no record of how it was created"):
+        model.to_recipe_json()
+
+
+# ------------------------------------------------------------------ older recipes
+def test_an_older_recipe_without_build_parameters_loads_with_a_warning(caplog):
+    import logging
+
+    model, _, _ = claudius_model()
+    recipe = model.to_recipe_dict()
+    del recipe["model"]["features"][0]["build"]
+    with caplog.at_level(logging.WARNING):
+        rebuilt = GeologicalModel.from_recipe_dict(recipe)
+    assert [f.name for f in rebuilt.features] == ["strati"]
+    assert any("no build parameters" in record.getMessage() for record in caplog.records)
+
+
+def test_an_older_recipe_entry_that_is_not_a_foliation_is_refused():
+    model, _, _ = claudius_model()
+    recipe = model.to_recipe_dict()
+    entry = recipe["model"]["features"][0]
+    del entry["build"]
+    entry["kind"] = "structural_frame"
+    with pytest.raises(LoopValueError, match="foliations only"):
+        GeologicalModel.from_recipe_dict(recipe)
+
+
+# ------------------------------------------------------- parameters are copied, not aliased
+def test_changing_a_parameter_after_creation_does_not_change_the_recipe():
+    model, _, _ = claudius_model()
+    values = [1, 2, 3]
+    weights = {"a": [1.0, 2.0]}
+    ModelSerializer.record_creation(
+        model["strati"], "foliation", {"nelements": 1000, "v": values, "w": weights}
+    )
+    before = model.to_recipe_json()
+    values.append(4)
+    weights["a"].append(3.0)
+    weights["b"] = 1
+    assert model.to_recipe_json() == before
+    parameters = json.loads(before)["model"]["features"][0]["build"]["parameters"]
+    assert parameters["v"] == [1, 2, 3]
+    assert parameters["w"] == {"a": [1.0, 2.0]}
+
+
+def test_a_numpy_array_parameter_is_stored_as_a_list_and_copied():
+    model, _, _ = claudius_model()
+    vector = np.array([0.0, 1.0, 0.5])
+    ModelSerializer.record_creation(model["strati"], "foliation", {"nelements": 1000, "vec": vector})
+    vector[0] = 99.0
+    parameters = json.loads(model.to_recipe_json())["model"]["features"][0]["build"]["parameters"]
+    assert parameters["vec"] == [0.0, 1.0, 0.5]
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..a6ac2bb7e
--- /dev/null
+++ b/tests/unit/modelling/test_recipe_reproduces_field.py
@@ -0,0 +1,46 @@
+"""A model recipe should rebuild a model that gives the same scalar field.
+
+The field depends on how each feature was built: the number of elements, the interpolator type
+and the regularisation weight. The recipe stores the parameters each feature was *requested*
+with (not the element count that was realised, which would not reproduce the grid), and the
+rebuilt model repeats the request.
+"""
+import numpy as np
+import pytest
+
+from LoopStructural import GeologicalModel
+from LoopStructural.datasets import load_claudius
+
+CASES = [
+    pytest.param({}, id="defaults"),
+    pytest.param({"nelements": 2000}, id="nelements"),
+    pytest.param({"interpolatortype": "PLI"}, id="interpolatortype"),
+    pytest.param({"regularisation": 10.0}, id="regularisation"),
+    pytest.param(
+        {"nelements": 2000, "regularisation": 10.0, "damp": False}, id="several-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)
+