diff --git a/INTRUSIONS.md b/INTRUSIONS.md
new file mode 100644
index 00000000..bbaf3ee9
--- /dev/null
+++ b/INTRUSIONS.md
@@ -0,0 +1,443 @@
+# Intrusions Module — Review, User Guide & Hardening Plan
+
+This document is the "dedicated discussion" input for **ROADMAP.md's Stage 6
+(outcome 8: intrusion workflow hardened, possibly rewritten)**. It has four
+parts:
+
+1. **How the module works today** — architecture summary.
+2. **User guide** — what data and parameters are actually required to build
+   an intrusion, worked from the one working example in the codebase.
+3. **Review findings** — concrete bugs, dead code, and design smells, each
+   with a file:line.
+4. **Simplification & hardening plan** — a phased proposal for Stage 6.
+
+Everything here was verified against the code as of commit `160c6e2a`
+(2026-08-05), not against docs or memory of older versions.
+
+**Status (2026-08-05): Phases A-C below are done.** The fixes, the data-
+contract validation, and new regression tests (including one that exercises
+`marginal_faults` end-to-end for the first time, and one that pins the newly-
+discovered `intrusion_steps` breakage — see finding 1) are all in
+`tests/unit/modelling/intrusions/test_intrusions.py`. Findings are left
+in place below as the historical record of what was found and fixed;
+each fixed item is marked **[FIXED]**. Phase D (the larger
+simplification/rewrite) is still the open item for the dedicated Stage 6
+discussion.
+
+---
+
+## 1. Architecture summary
+
+Files in `LoopStructural/modelling/intrusions/`:
+
+| File | Role |
+|---|---|
+| `intrusion_frame.py` | `IntrusionFrame(StructuralFrame)` — marker subclass, no added behaviour. |
+| `intrusion_frame_builder.py` | `IntrusionFrameBuilder(StructuralFrameBuilder)` — builds the curvilinear coordinate system (`c0`=growth away from the reference contact, `c1`=position along strike, `c2`=lateral/side distance), optionally conditioned by fault "steps" and "marginal faults". |
+| `intrusion_builder.py` | `IntrusionBuilder(BaseBuilder)` — given the frame, prepares per-side/per-contact data and fits RBF interpolators over `(c1, c2)` for lateral thresholds and `(c1, c2)` for vertical thresholds, constrained by a "conceptual model" function. |
+| `intrusion_feature.py` | `IntrusionFeature(BaseFeature)` — evaluates a signed-distance-like scalar field to the intrusion contact from the frame + the fitted thresholds. |
+| `geom_conceptual_models.py` | Three conceptual-geometry functions: `ellipse_function`, `constant_function`, `obliquecone_function`. |
+| `geometric_scaling_functions.py` | Empirical length→thickness scaling from the literature (pluton/laccolith/sill scaling laws) — see bugs, this is **not wired up**. |
+| `intrusion_support_functions.py` | Fast-marching "shortest path" helpers. **Dead code** — see §3. |
+
+Entry point: `GeologicalModel.create_and_add_intrusion(...)` →
+`_build_intrusion(...)` (`LoopStructural/modelling/core/geological_model.py:1373-1525`),
+registered in the `FeatureBuilderRegistry` under `"intrusion"`.
+
+Build sequence inside `_build_intrusion`:
+
+```
+IntrusionFrameBuilder.set_intrusion_frame_parameters(intrusion_data, params)
+IntrusionFrameBuilder.create_constraints_for_c0()      # synthesise c0=0 points/gradients
+IntrusionFrameBuilder.set_intrusion_frame_data(frame_data)
+IntrusionFrameBuilder.build(...)                        # -> IntrusionFrame
+IntrusionBuilder(frame, lateral_extent_model=..., vertical_extent_model=...)
+IntrusionBuilder.set_data_for_extent_calculation(intrusion_data)
+IntrusionBuilder.update_build_arguments({"geometric_scaling_parameters": ...})
+# lazy: IntrusionBuilder.build() runs on first evaluate_value() via up_to_date()
+```
+
+---
+
+## 2. User guide — what's required to build an intrusion
+
+As of 2026-08-05, `tests/unit/modelling/intrusions/test_intrusions.py`
+covers: a tabular intrusion built from a single roof-or-floor contact plus
+one stratigraphic conformable feature as the inflation-gradient proxy
+(`load_tabular_intrusion()`, the original path), and — new — one
+`marginal_faults` example built from scratch
+(`test_intrusion_marginal_faults`/`_build_marginal_fault_model`) showing a
+sill offset against a single bounding fault. `intrusion_steps` is *not*
+demonstrated as working, because it currently can't be — see §3, finding 1;
+`test_intrusion_steps_broken_with_current_stratigraphic_column` pins the
+failure instead. The shortest-path network method remains undemonstrated
+and is now dead code (§3, finding 5) rather than just untested.
+
+### 2.1 Two blocks of input data
+
+Both blocks live in `model.data` (one `pandas.DataFrame`, distinguished by
+`feature_name`), matching `intrusion_name` and `intrusion_frame_name` passed
+to `create_and_add_intrusion`.
+
+**A. Intrusion frame data** (`feature_name == intrusion_frame_name`) — this
+is an ordinary structural-frame dataset, the same schema used for faults and
+fold frames: rows tagged `coord` 0/1/2, each with either a `val` (point lies
+on that isovalue) or a gradient (`nx, ny, nz`, vector the coordinate's
+gradient should be parallel to). You do **not** need to supply `coord=0`
+data yourselves — `IntrusionFrameBuilder.set_intrusion_frame_data` (`intrusion_frame_builder.py:935-978`)
+synthesises coord-0 point and gradient constraints from the intrusion
+contact data below and appends them for you. You only need to supply
+`coord=1` and `coord=2` constraints (an origin point + two orthogonal
+direction vectors is enough — see the example).
+
+**B. Intrusion contact data** (`feature_name == intrusion_name`) — points
+on the intrusion margin:
+
+| Column | Required? | Meaning |
+|---|---|---|
+| `X, Y, Z` | yes | point location |
+| `intrusion_contact_type` | yes | `"roof"`/`"top"` or `"floor"`/`"base"` — which margin this point is on |
+| `intrusion_side` | yes for lateral extent | boolean; `True` marks points used to constrain the lateral (lengthwise) margins. Split into "min side"/"max side" by the sign of `c2` at that point. |
+| `intrusion_anisotropy` | only for steps/marginal-faults/shortest-path | name of the host-rock series or fault feature the point is associated with |
+
+None of these columns are validated: a missing `intrusion_contact_type`
+column raises a bare pandas `KeyError` deep inside `set_intrusion_frame_c0_data`
+with no indication of what's wrong (see §3, finding 3).
+
+### 2.2 `intrusion_frame_parameters` dict
+
+Passed to `create_and_add_intrusion(..., intrusion_frame_parameters={...})`.
+
+| Key | Default | Notes |
+|---|---|---|
+| `contact` | `"floor"` | which `intrusion_contact_type` value is the *reference* contact used to build `c0=0` |
+| `contact_anisotropies` | `None` | **de facto required** — a list of series-type features; `create_constraints_for_c0` unconditionally indexes `[0]` (`intrusion_frame_builder.py:890`) to get an inflation-gradient proxy. Omitting it, or passing `[]`, crashes with `TypeError`/`IndexError`, not a helpful error. |
+| `g_w` | `None` → 100 pts | weight/count controlling how many synthetic gradient constraints get added for `c0` |
+| `intrusion_steps` | `None` | **currently broken, don't use** — see §3, finding 1. Was meant to be a dict of step definitions, each needing `structure` (fault), `unit_from`/`unit_to` (stratigraphic unit names), `series_from`/`series_to` (series features) |
+| `marginal_faults` | `None` | dict of fault definitions bounding the intrusion; each needs `structure` (a `FaultSegment`/`FaultSegment`-like object, indexed as `structure[0]` for its coordinate-0 feature — same convention as `intrusion_steps`' `structure`), `block` (`"hanging wall"`/`"foot wall"`), `series` (a series feature). Working example: `test_intrusion_marginal_faults` in `test_intrusions.py`. |
+| `delta_c`, `delta_f` | `[1]*n` | multiplies the std-dev band used to detect points near a contact/fault when synthesising `c0` constraints |
+
+### 2.3 Conceptual model functions
+
+`intrusion_lateral_extent_model` and `intrusion_vertical_extent_model` are
+plain functions from `geom_conceptual_models.py` (or custom ones matching
+the same **dual-arity calling convention** — see §3, finding 8):
+`ellipse_function` (lateral) and `constant_function` or `obliquecone_function`
+(vertical) are the built-ins.
+
+### 2.4 Worked example (from the passing test)
+
+```python
+from LoopStructural import GeologicalModel
+from LoopStructural.datasets import load_tabular_intrusion
+from LoopStructural.modelling.intrusions import ellipse_function, constant_function
+
+data, bounding_box = load_tabular_intrusion()
+
+model = GeologicalModel(bounding_box[0, :], bounding_box[1, :])
+# NOT `model.data = data` -- this dataset has no nx/ny/nz/tx/ty/tz columns,
+# and `prepare_data` is what normalises those in; skipping it builds a model
+# that crashes with a bare KeyError the moment anything calls
+# `evaluate_value()` on it. See §3, finding 1b.
+model.data = model.prepare_data(data)
+
+stratigraphy = model.create_and_add_foliation("stratigraphy")
+
+intrusion = model.create_and_add_intrusion(
+    "tabular_intrusion",
+    "tabular_intrusion_frame",
+    intrusion_frame_parameters={
+        "contact": "roof",
+        "contact_anisotropies": [stratigraphy],
+    },
+    intrusion_lateral_extent_model=ellipse_function,
+    intrusion_vertical_extent_model=constant_function,
+)
+```
+
+`intrusion.evaluate_value(xyz)` then returns a signed pseudo-distance field
+to the intrusion contact (negated internally so it can double as an
+unconformity — `intrusion_feature.py:329-331`). Note `evaluate_gradient` is
+`NotImplementedError` (`intrusion_feature.py:233-235`) — anything in the
+model stack that needs gradients of a feature (fold axes, some fault
+throw calculations, vector-field visualisation) cannot be pointed at an
+`IntrusionFeature`.
+
+---
+
+## 3. Review findings
+
+Ordered roughly most → least impactful. **[FIXED]** markers were added
+2026-08-05 after implementing Phases A-C (see §4); the rest of the text is
+left as originally written for the record.
+
+1. **`intrusion_steps` is not just untested — it's provably broken against
+   the current `StratigraphicColumn`.** Discovered while writing the Phase C
+   regression test (below), this is worse than originally written here.
+   `IntrusionFrameBuilder.set_intrusion_steps_parameters` reads
+   `self.model.stratigraphic_column[series_from_name.name][unit_from_name].get("id"/"min"/"max")`
+   — i.e. it expects `model.stratigraphic_column` to still be the old
+   nested `{group_name: {unit_name: {...}}}` dict. It no longer is:
+   `GeologicalModel.stratigraphic_column` is now a `StratigraphicColumn`
+   object (`modelling/core/stratigraphic_column.py`) whose `__getitem__`
+   looks up a single element by `uuid`, not by group/series name, and
+   returns a `StratigraphicColumnElement` (no `.get()`). Worse,
+   `GeologicalModel.set_stratigraphic_column` — the old dict-based setter —
+   now unconditionally `raise DeprecationWarning(...)` after logging, so
+   there is no longer any way to put the model into the shape
+   `intrusion_steps` expects at all. This is a real regression from an
+   unrelated stratigraphic-column refactor that was never propagated to the
+   intrusions module. **Confirmed by** `test_intrusion_steps_broken_with_current_stratigraphic_column`
+   in `test_intrusions.py`, which pins the exact current failure
+   (`KeyError: 'No element found with uuid: ...'`). Fixing this for real
+   needs a design decision (what should `unit_from`/`series_from` resolve
+   against now?) — that's Stage 6 work, not something patched in Phase A-C.
+
+1b. **[FIXED, as a documentation/test fix]** The packaged
+   `datasets/data/tabular_intrusion.csv` — the *only* dataset this module
+   has — has no `nx`/`ny`/`nz`/`tx`/`ty`/`tz` columns, only `gx`/`gy`/`gz`.
+   `GeologicalModel.data = df` does not normalise columns (unlike
+   `GeologicalModel.prepare_data(df)`, which adds any of `all_heading()`
+   missing as `NaN`). `GeologicalFeatureBuilder.add_data_to_interpolator`
+   unconditionally indexes `normal_vec_names()`/`tangent_vec_names()` on a
+   builder's own data, so **any full end-to-end build+`evaluate_value()`
+   call through the real `create_and_add_intrusion` API** (as opposed to
+   the two pre-existing tests, which manually construct
+   `IntrusionFrameBuilder`/`IntrusionBuilder` and pass already-`prepare_data`-d
+   slices, or never call `evaluate_value` at all) crashes with a bare
+   pandas `KeyError` on the missing columns. This is a general
+   `GeologicalModel`/`GeologicalFeatureBuilder` sharp edge, not specific to
+   intrusions, but intrusions is where it surfaced, precisely because (per
+   finding 1) nothing had ever driven this dataset through the real,
+   public, end-to-end path before. The fix used throughout the new tests is
+   simply `model.data = model.prepare_data(df)` instead of `model.data = df`
+   — worth calling out prominently in any future user guide/docs for this
+   dataset, since the obvious `model.data = data` (as the original
+   `test_intrusion_frame_builder`/`test_intrusion_builder` imply, and as
+   this document's own §2.4 example originally showed) silently sets up a
+   model that will crash the moment anyone actually evaluates it.
+
+2. **[FIXED] `geometric_scaling_parameters` is fully non-functional.**
+   `IntrusionBuilder.create_geometry_using_geometric_scaling`
+   (`intrusion_builder.py:113-154`) always raised `NotImplementedError`,
+   even on the branch where `thickness` *is* provided — the
+   `raise NotImplementedError("Not implemented")` at line 146 was
+   unconditional, after the early-return-shaped `if estimated_thickness is
+   None` check. The only way to reach this function is if the *other*
+   contact has zero data points **and** `geometric_scaling_parameters` is
+   non-empty, so passing this parameter to `create_and_add_intrusion` for a
+   single-contact intrusion always crashed. `geometric_scaling_functions.py`'s
+   `contact_pts_using_geometric_scaling` (which this should call) is fully
+   implemented and unit-testable, but dead — its only caller was commented
+   out. **Fix applied:** the function now raises immediately, with a clear
+   message naming what's missing, instead of after several lines that look
+   like partial progress; the dead commented-out code and the now-pointless
+   wildcard import of `geometric_scaling_functions` were removed. Pinned by
+   `test_intrusion_geometric_scaling_not_implemented`.
+
+3. **[FIXED] A typo silently drops a user-supplied build weight.**
+   `GeologicalModel._build_intrusion` called
+   `intrusion_frame_builder.build(nelements=..., w2=weights[0],
+   w1=weights[1], gxygz=weights[2])`. `StructuralFrameBuilder.build`
+   only recognises the deprecated alias `gyxgz` for `w3` — `gxygz` (letters
+   transposed) fell into `**kwargs` and was silently ignored. Any caller
+   passing `gyxgz=...` to `create_and_add_intrusion` to weight the frame's
+   third-coordinate orthogonality constraint had no effect; `w3` was always
+   `1.0`. This was the only call site in the codebase using this kwarg name
+   for a frame build — faults/foliations don't have this bug, which is
+   consistent with the intrusion path being far less exercised. **Fix
+   applied:** `gxygz` → `gyxgz`. Pinned by
+   `test_intrusion_gyxgz_weight_reaches_frame_build`.
+
+4. **[FIXED] No data-contract validation.** Missing `intrusion_contact_type`,
+   `intrusion_side`, or `intrusion_anisotropy` columns failed with bare
+   `KeyError`/`AttributeError` deep in helper methods, not at the
+   `create_and_add_intrusion` boundary where the user could get a useful
+   message. Likewise `contact_anisotropies` being `None`/empty crashed with
+   no hint of what's missing. **Fix applied:**
+   `GeologicalModel._validate_intrusion_inputs` now runs at the top of
+   `_build_intrusion` and raises a clear `ValueError` naming the missing
+   feature/column/parameter. Pinned by
+   `test_intrusion_missing_data_raises_clear_error`,
+   `test_intrusion_missing_contact_type_column_raises_clear_error`,
+   `test_intrusion_missing_contact_anisotropies_raises_clear_error`.
+
+5. **Dead code — mostly [FIXED]:**
+   - **[FIXED, deleted]** `intrusion_support_functions.py` (shortest-path
+     grid/graph helpers, ~390 lines) — not imported anywhere outside itself.
+   - **[FIXED, deleted]** `IntrusionFeature.evaluate_value_test` — an entire
+     parallel implementation of `evaluate_value`, unused.
+   - **[FIXED, deleted]** `IntrusionFrameBuilder.update()` was an exact
+     duplicate of the parent `StructuralFrameBuilder.update()` — redundant
+     override.
+   - **[FIXED, deleted]** the two discarded
+     `self.marginal_faults[fault_i].get("emplacement_mechanism")` calls.
+   - **[FIXED]** `IntrusionFeature.evaluate_value`'s dead `if/else` with
+     identical bodies (the `marginal_faults is not None` branch) collapsed
+     to one unconditional assignment.
+   - **Left as-is (not removed):** `IntrusionFeature.add_assisting_faults`/
+     `self.assisting_faults` is still set but never populated by any caller
+     in the codebase, so the "asymmetry weight" branch in `evaluate_value`
+     stays unreachable in practice — removing a public method felt like a
+     bigger call than Phase A's "no behaviour change" scope; flagging again
+     for Phase D. `intrusion_builder.py`'s comment
+     `# intrusion_frame_builder.post_intrusion_faults = faults  # LG unused?`
+     (your own prior TODO) is also still there, unaddressed.
+
+6. **Reproducibility, corrected.** Originally written up here as "KMeans
+   hardcodes `random_state=0` instead of going through the shared,
+   seedable `rng`" — that framing turned out to be backwards.
+   `loop_common.utils.rng` (`from ...utils import rng`, used e.g. for
+   `rng.shuffle`) is `np.random.default_rng()` created **unseeded**, once,
+   at import time — it's not actually reproducible across process runs
+   either, and there's no `set_seed`-style entry point anywhere in the
+   codebase. Routing `KMeans` through it would have made cluster labels
+   (and therefore the synthesised `c0` constraints for steps/marginal
+   faults) vary run-to-run instead of being the one deterministic piece.
+   **Fix applied instead:** kept the fixed seed (it's the right call), but
+   replaced the two duplicated `random_state=0` literals and the stale
+   `# TODO create global loopstructural random state variable` comment
+   (which described a solution — the shared `rng` — that doesn't actually
+   solve this) with one named module constant, `_KMEANS_RANDOM_STATE`.
+
+7. **Scikit-learn dependency, corrected.** Originally flagged here as "a
+   hard dependency on scikit-learn just for the steps/marginal-faults code
+   path." On checking `pyproject.toml`, `scikit-learn` is already a
+   top-level `dependencies` entry for the whole package (also used in
+   `LoopStructural/utils/helper.py` and `_transformation.py`), so this
+   isn't an extra install cost specific to intrusions — retracted as a
+   finding. The module-scope `try/except ImportError: ... raise` in
+   `intrusion_frame_builder.py` is still slightly unusual style (failing at
+   import time rather than at first use), but not a real problem worth
+   spending Phase A/B effort on.
+
+8. **Conceptual-model functions overload call arity to mean two different
+   things.** `ellipse_function`/`constant_function`/`obliquecone_function`
+   are called once with **no arguments** (during
+   `set_conceptual_models_parameters`, to fetch bounds:
+   `intrusion_builder.py:277-278`) and later called **with data** to
+   evaluate the model at points. This works only because every branch
+   checks `if <data>.empty` and returns a different tuple shape. It's an
+   implicit, undocumented protocol — a custom conceptual model has to
+   reverse-engineer this from the three built-ins rather than from any
+   docstring or interface. A named two-method interface (e.g. `.bounds()`
+   / `.evaluate(data, ...)`, whether a `Protocol` or small ABC) would make
+   custom conceptual models actually writable by a user.
+
+9. **`IntrusionFrameBuilder` is a ~20-attribute god object** mixing frame
+   geometry, steps, marginal faults, deprecated shortest-path indicator
+   functions (`IFf`, `IFc`, single/double-letter names), and sill-splitting
+   in one class. `create_constraints_for_c0` alone is ~200 lines covering
+   four largely-independent concerns (steps / marginal faults / gradient
+   constraints / sill-splitting) in one method with deep nesting.
+
+10. **Copy-pasted min/max-side blocks.** `IntrusionBuilder.set_data_for_lateral_thresholds`
+    (`intrusion_builder.py:344-563`, ~220 lines) repeats the same
+    conceptual/residual computation twice — once for the `L<0` side, once
+    for `L>0` — with only sign/index flips between them. Same shape of
+    duplication in `interpolate_lateral_thresholds`
+    (`intrusion_feature.py:60-137`). Both are natural candidates for a
+    single side-parameterised helper.
+
+11. **Fragile cross-feature reach-through.** Sill-splitting
+    (`create_constraints_for_c0`, `intrusion_frame_builder.py:807-815`) does
+    `self.model.__getitem__(splits_from_sill_name).intrusion_frame.builder.intrusion_steps`
+    — reaching through a private-ish attribute chain on another feature
+    entirely, with no cycle detection if two sills reference each other.
+
+12. **[FIXED] Wildcard import.** `intrusion_builder.py:6` had
+    `from .geometric_scaling_functions import *`. Removed as part of the
+    finding-2 fix, once `create_geometry_using_geometric_scaling` no longer
+    called anything from that module.
+
+13. **Naming/typo debt** (low severity, but adds up for anyone trying to
+    read this code for the first time): `strigraphic` (`intrusion_frame_builder.py:306,474`),
+    `framce` (`:307,475`), `aftecting` (`:904`), `yo may increase` (`:822`),
+    `indentify` (multiple), single-letter/cryptic names (`If`, `Ic`, `IFf`,
+    `IFc`, `s`, `d`/`e`/`f` as boolean masks in `evaluate_value`).
+
+---
+
+## 4. Simplification & hardening plan
+
+Proposed as the concrete content for **ROADMAP.md Stage 6**. Phased so each
+step is independently shippable and testable, following the same
+small-stages philosophy as the rest of the v2 roadmap.
+
+**Phase A — stop the bleeding (no behaviour change to the working path) —
+DONE 2026-08-05:**
+- Fixed the `gxygz`/`gyxgz` typo (finding 3); regression test
+  `test_intrusion_gyxgz_weight_reaches_frame_build`.
+- Deleted confirmed-dead code: `intrusion_support_functions.py`,
+  `evaluate_value_test`, the unused `emplacement_mechanism` `.get()` calls,
+  the redundant `IntrusionFrameBuilder.update()` override, the now-pointless
+  wildcard import.
+- `geometric_scaling_parameters` now fails immediately with a clear message
+  instead of after several lines that looked like partial progress
+  (finding 2); regression test `test_intrusion_geometric_scaling_not_implemented`.
+  `add_assisting_faults` was deliberately **not** touched — removing a
+  public method felt out of scope for "no behaviour change"; still flagged
+  for Phase D.
+- Replaced the two hardcoded `KMeans(random_state=0)` with one named
+  constant and corrected the misleading comment about routing through the
+  shared `rng` (finding 6 — that would have made results *less*
+  reproducible, not more; see the corrected finding for why).
+
+**Phase B — make the failure modes legible — DONE 2026-08-05:**
+- `GeologicalModel._validate_intrusion_inputs` now runs at the top of
+  `_build_intrusion`: checks both `feature_name`s have data, required
+  columns (`intrusion_contact_type`, `intrusion_side`) are present, and
+  `contact_anisotropies` is non-empty — clear `ValueError`s instead of a
+  `KeyError` several calls deep. Regression tests
+  `test_intrusion_missing_data_raises_clear_error`,
+  `test_intrusion_missing_contact_type_column_raises_clear_error`,
+  `test_intrusion_missing_contact_anisotropies_raises_clear_error`.
+- The "make scikit-learn optional" item was **retracted**: scikit-learn is
+  already a top-level dependency of the whole package (finding 7,
+  corrected), so there was nothing to fix here.
+
+**Phase C — cover what's untested before touching it further — DONE
+2026-08-05, with a significant finding:**
+- Added `test_intrusion_marginal_faults` (`_build_marginal_fault_model`): a
+  from-scratch synthetic sill offset by one bounding fault, built and
+  evaluated end-to-end through the real `create_and_add_intrusion` public
+  API — the first time this code path has ever been exercised.
+- Attempting the equivalent for `intrusion_steps` surfaced finding 1:
+  it doesn't work at all against the current `StratigraphicColumn`, not
+  just "untested." `test_intrusion_steps_broken_with_current_stratigraphic_column`
+  pins the current failure instead of demonstrating success. **This means
+  Phase D's scope for `intrusion_steps` is bigger than originally framed
+  here**: it's not "simplify working code," it's "decide what this feature
+  should even look like against the new stratigraphic column API, or drop
+  it." That decision belongs in the dedicated Stage 6 discussion.
+- Also surfaced finding 1b (the packaged `tabular_intrusion.csv` dataset
+  needs `model.prepare_data()`, not a bare `model.data = df`, to survive an
+  end-to-end `evaluate_value()` call) — fixed in the tests and in this
+  document's own worked example (§2.4), which had the same bug.
+
+**Phase D — the actual simplification (larger, still open — needs the
+Stage 6 discussion):**
+- Split `IntrusionFrameBuilder` along its four concerns (frame geometry /
+  steps / marginal faults / deprecated shortest-path) into
+  composable pieces, so the common tabular case doesn't carry ~20 unused
+  attributes.
+- Replace the conceptual-model dual-arity calling convention (finding 8)
+  with an explicit two-method interface; migrate the three built-ins.
+- Collapse the min/max-side and lateral/vertical duplicated blocks
+  (findings 10) into single side-parameterised helpers.
+- Decide whether the shortest-path network method is kept (and given the
+  same test/doc treatment as the rest) or removed outright — it's now fully
+  dead code (finding 5/1), not half-present: `intrusion_network_type` is
+  initialised to `None` and never set to `"shortest path"` anywhere reachable.
+- Decide what `intrusion_steps` should even look like against the current
+  `StratigraphicColumn` object (finding 1) — the old nested-dict lookup it
+  was written against no longer exists, and there's no bridge back to it.
+  This is bigger than a refactor: it needs a real design decision before
+  any of the "simplify the god object" work below can safely touch it.
+
+Phases A-C landed 2026-08-05 (see the "DONE" notes above) — they didn't
+need to wait for the Stage-5 graph backend. Phase D is what still needs the
+"dedicated discussion" ROADMAP.md defers to post-Stage-5, since a graph
+backend may change how frame/feature dependencies (sill-splitting, step
+faults) are expressed anyway, and the `intrusion_steps`/`StratigraphicColumn`
+question needs a decision either way.
diff --git a/LoopStructural/modelling/core/_model_exporter.py b/LoopStructural/modelling/core/_model_exporter.py
new file mode 100644
index 00000000..7bdfc679
--- /dev/null
+++ b/LoopStructural/modelling/core/_model_exporter.py
@@ -0,0 +1,129 @@
+"""Surface/block-model export logic for GeologicalModel (see API.md).
+
+Extracted from GeologicalModel to separate export/visualization data prep
+from feature-container orchestration. GeologicalModel's @public_api methods
+(``get_fault_surfaces``, ``get_stratigraphic_surfaces``, ``get_block_model``,
+``save``) stay defined directly on the class -- their __qualname__ is part
+of the CI-checked stable API surface -- and delegate to the staticmethods
+here.
+"""
+
+import pathlib
+
+from ...geometry import StructuredGrid
+from ...utils import getLogger
+
+logger = getLogger(__name__)
+
+
+class ModelExporter:
+    @staticmethod
+    def get_fault_surfaces(model, faults=None):
+        if faults is None:
+            faults = []
+        surfaces = []
+        if len(faults) == 0:
+            faults = model.fault_names()
+
+        for f in faults:
+            surfaces.extend(model.get_feature_by_name(f).surfaces([0], model.bounding_box))
+        return surfaces
+
+    @staticmethod
+    def get_stratigraphic_surfaces(model, units=None, bottoms=True):
+        if units is None:
+            units = []
+        ## TODO change the stratigraphic column to its own class and have methods to get the relevant surfaces
+        surfaces = []
+        units = []
+        if model.stratigraphic_column is None:
+            return []
+        units = model.stratigraphic_column.get_isovalues()
+        units_for_group = {}
+        for name, u in units.items():
+            if u['group'] not in model:
+                logger.warning(f"Group {u['group']} not found in model")
+                continue
+            if u['group'] not in units_for_group:
+                units_for_group[u['group']] = []
+            u['name'] = name
+            units_for_group[u['group']].append(u)
+        for group, us in units_for_group.items():
+            feature = model.get_feature_by_name(group)
+            values = [u['value'] for u in us]
+            colours = [u['colour'] for u in us]
+            names = [u['name'] for u in us]
+            surfaces.extend(
+                feature.surfaces(values, model.bounding_box, name=names, colours=colours)
+            )
+
+        return surfaces
+
+    @staticmethod
+    def get_block_model(model, name='block model'):
+        # NOTE: bounding_box.structured_grid() returns loop_common's
+        # interpolation-support StructuredGrid (no properties dict); use
+        # LoopStructural's own geometry StructuredGrid for storing values.
+        grid = StructuredGrid(
+            origin=model.bounding_box.origin,
+            step_vector=model.bounding_box.step_vector,
+            nsteps=model.bounding_box.nsteps,
+            name=name,
+        )
+
+        grid.cell_properties['stratigraphy'] = model.evaluate_model(
+            model.rescale(model.bounding_box.cell_centres())
+        )
+        return grid, model.stratigraphic_ids()
+
+    @staticmethod
+    def save(
+        model,
+        filename: str,
+        block_model: bool = True,
+        stratigraphic_surfaces=True,
+        fault_surfaces=True,
+        stratigraphic_data=True,
+        fault_data=True,
+    ):
+        path = pathlib.Path(filename)
+        extension = path.suffix
+        parent = path.parent
+        name = path.stem
+        stratigraphic_surfaces = model.get_stratigraphic_surfaces()
+        if fault_surfaces:
+            for s in model.get_fault_surfaces():
+                ## geoh5 can save everything into the same file
+                if extension == ".geoh5" or extension == '.omf':
+                    s.save(filename)
+                else:
+                    s.save(f'{parent}/{name}_{s.name}{extension}')
+        if stratigraphic_surfaces:
+            for s in model.get_stratigraphic_surfaces():
+                if extension == ".geoh5" or extension == '.omf':
+                    s.save(filename)
+                else:
+                    s.save(f'{parent}/{name}_{s.name}{extension}')
+        if block_model:
+            grid, _ids = model.get_block_model()
+            if extension == ".geoh5" or extension == '.omf':
+                grid.save(filename)
+            else:
+                grid.save(f'{parent}/{name}_block_model{extension}')
+        if stratigraphic_data and model.stratigraphic_column is not None:
+            for group in model.stratigraphic_column:
+                if group == "faults":
+                    continue
+                for data in model.__getitem__(group).get_data():
+                    if extension == ".geoh5" or extension == '.omf':
+                        data.save(filename)
+                    else:
+                        data.save(f'{parent}/{name}_{group}_data{extension}')
+        if fault_data:
+            for f in model.fault_names():
+                for d in model.__getitem__(f).get_data():
+                    if extension == ".geoh5" or extension == '.omf':
+
+                        d.save(filename)
+                    else:
+                        d.save(f'{parent}/{name}_{group}{extension}')
diff --git a/LoopStructural/modelling/core/_model_feature_factory.py b/LoopStructural/modelling/core/_model_feature_factory.py
new file mode 100644
index 00000000..300562de
--- /dev/null
+++ b/LoopStructural/modelling/core/_model_feature_factory.py
@@ -0,0 +1,758 @@
+"""Feature-construction orchestration for GeologicalModel (see API.md).
+
+Extracted from GeologicalModel to separate the mechanics of building each
+feature type (foliation, fold frame, folded foliation, folded fold frame,
+intrusion, domain fault, fault) from feature-container state and the public
+``create_and_add_*`` API. GeologicalModel's @public_api ``create_and_add_*``
+methods stay defined directly on the class -- their __qualname__ is part of
+the CI-checked stable API surface -- and dispatch through
+``GeologicalModel.create_and_add_feature`` ->
+``FeatureBuilderRegistry.create`` -> the staticmethods here (wired up at the
+bottom of ``geological_model.py``).
+"""
+
+import numpy as np
+import pandas as pd
+
+from LoopStructural import LoopStructuralConfig
+
+from ...modelling.features import FeatureType, UnconformityFeature
+from ...modelling.features.builders import (
+    FaultBuilder,
+    FoldedFeatureBuilder,
+    GeologicalFeatureBuilder,
+    StructuralFrameBuilder,
+)
+from ...modelling.features.fold import FoldEvent, FoldFrame
+from ...modelling.intrusions import IntrusionBuilder, IntrusionFrameBuilder
+from ...utils import getLogger
+
+logger = getLogger(__name__)
+
+
+class ModelFeatureFactory:
+    @staticmethod
+    def build_foliation(
+        model,
+        series_surface_name: str,
+        *,
+        index: int | None = None,
+        data: pd.DataFrame | None = None,
+        interpolatortype: str = "FDI",
+        nelements: int = LoopStructuralConfig.nelements,
+        tol=None,
+        faults=None,
+        **kwargs,
+    ):
+        """
+        Parameters
+        ----------
+        series_surface_name : string
+            corresponding to the feature_name in the data
+        series_surface_data : pd.DataFrame, optional
+            data frame containing the surface data
+        interpolatortype : str
+            the type of interpolator to use, default is 'FDI'
+        nelements : int
+            the number of elements to use in the series surface
+        tol : float, optional
+            tolerance for the solver, if not specified uses the model default
+        faults : list, optional
+            list of faults to be used in the series surface, if not specified uses the model faults
+        kwargs
+
+        Returns
+        -------
+        feature : GeologicalFeature
+            the created geological feature
+
+        Notes
+        ------
+        This function creates an instance of a
+        :class:`LoopStructural.modelling.features.builders.GeologicalFeatureBuilder` and will return
+        a :class:`LoopStructural.modelling.features.builders.GeologicalFeature`
+        The feature is not interpolated until either
+        :meth:`LoopStructural.modelling.features.builders.GeologicalFeature.evaluate_value` is called or
+        :meth:`LoopStructural.modelling.core.GeologicalModel.update`
+
+        An interpolator will be chosen by calling :meth:`LoopStructural.GeologicalModel.get_interpolator`
+
+        """
+
+        # if tol is not specified use the model default
+        if tol is None:
+            tol = model.tol
+
+        series_builder = GeologicalFeatureBuilder(
+            bounding_box=model.bounding_box,
+            interpolatortype=interpolatortype,
+            nelements=nelements,
+            name=series_surface_name,
+            model=model,
+            **kwargs,
+        )
+        # add data
+        if data is None:
+            data = model.data.loc[model.data["feature_name"] == series_surface_name]
+
+        if data.shape[0] == 0:
+            logger.warning("No data for {series_surface_data}, skipping")
+            return
+        series_builder.add_data_from_data_frame(model.prepare_data(data, include_feature_name=False))
+        model._add_faults(series_builder, features=faults)
+
+        # build feature
+        series_feature = series_builder.feature
+        series_builder.update_build_arguments(kwargs | {"domain": True, 'tol': tol})
+        # this support is built for the entire model domain? Possibly would
+        # could just pass a regular grid of points - mask by any above unconformities??
+
+        series_feature.type = FeatureType.INTERPOLATED
+        model._add_feature(series_feature, index=index)
+        return series_feature
+
+    @staticmethod
+    def build_fold_frame(
+        model,
+        fold_frame_name: str,
+        *,
+        index: int | None = None,
+        data=None,
+        interpolatortype="FDI",
+        nelements=LoopStructuralConfig.nelements,
+        tol=None,
+        buffer=0.1,
+        **kwargs,
+    ):
+        """
+        Parameters
+        ----------
+        fold_frame_name : string
+            unique string in feature_name column
+        fold_frame_data : pandas data frame
+            if not specified uses the model data
+        interpolatortype : str
+            the type of interpolator to use, default is 'FDI'
+        nelements : int
+            the number of elements to use in the fold frame
+        tol : float, optional
+            tolerance for the solver
+        buffer : float
+            buffer to add to the bounding box of the fold frame
+        **kwargs : dict
+            additional parameters to be passed to the
+            :class:`LoopStructural.modelling.features.builders.StructuralFrameBuilder`
+            and :meth:`LoopStructural.modelling.features.builders.StructuralFrameBuilder.setup`
+            and the interpolator, such as `domain` or `tol`
+
+
+        Returns
+        -------
+        fold_frame : FoldFrame
+            the created fold frame
+        """
+
+        if tol is None:
+            tol = model.tol
+
+        # create fault frame
+        #
+        fold_frame_builder = StructuralFrameBuilder(
+            interpolatortype=interpolatortype,
+            bounding_box=model.bounding_box.with_buffer(buffer),
+            name=fold_frame_name,
+            frame=FoldFrame,
+            nelements=nelements,
+            model=model,
+            **kwargs,
+        )
+        # add data
+        if data is None:
+            data = model.data.loc[model.data["feature_name"] == fold_frame_name]
+        if data.shape[0] == 0:
+            logger.warning(f"No data for {fold_frame_name}, skipping")
+            return
+        fold_frame_builder.add_data_from_data_frame(
+            model.prepare_data(data, include_feature_name=False)
+        )
+        model._add_faults(fold_frame_builder[0])
+        model._add_faults(fold_frame_builder[1])
+        model._add_faults(fold_frame_builder[2])
+        kwargs["tol"] = tol
+        fold_frame_builder.build(**kwargs)
+        fold_frame = fold_frame_builder.frame
+
+        fold_frame.type = FeatureType.STRUCTURALFRAME
+        fold_frame.builder = fold_frame_builder
+        model._add_feature(fold_frame, index=index)
+
+        return fold_frame
+
+    @staticmethod
+    def build_folded_foliation(
+        model,
+        foliation_name,
+        *,
+        index: int | None = None,
+        data=None,
+        interpolatortype="DFI",
+        nelements=LoopStructuralConfig.nelements,
+        buffer=0.1,
+        fold_frame=None,
+        svario=True,
+        tol=None,
+        invert_fold_norm=False,
+        **kwargs,
+    ):
+        """
+        Create a folded foliation field from data and a fold frame
+
+        Parameters
+        ----------
+        foliation_data : str
+            unique string in type column of data frame
+        fold_frame :  FoldFrame
+        svario  : Boolean
+            whether to calculate svariograms, saves time if avoided
+        kwargs
+            additional kwargs to be passed through to other functions
+
+        Returns
+        -------
+        feature : GeologicalFeature
+            created geological feature
+
+        Notes
+        -----
+
+        - Building a folded foliation uses the fold interpolation code from Laurent et al., 2016
+        and fold profile fitting from Grose et al., 2017. For more information about the fold modelling
+        see :class:`LoopStructural.modelling.features.fold.FoldEvent`,
+        :class:`LoopStructural.modelling.features.builders.FoldedFeatureBuilder`
+
+        """
+
+        if tol is None:
+            tol = model.tol
+
+        if fold_frame is None:
+            logger.info("Using last feature as fold frame")
+            fold_frame = model.features[-1]
+        if not isinstance(fold_frame, FoldFrame):
+            raise TypeError("Please specify a FoldFrame")
+
+        fold = FoldEvent(fold_frame, name=f"Fold_{foliation_name}", invert_norm=invert_fold_norm)
+
+        if interpolatortype != "DFI":
+            logger.warning("Folded foliation only supports DFI interpolator, changing to DFI")
+            interpolatortype = "DFI"
+        series_builder = FoldedFeatureBuilder(
+            interpolatortype=interpolatortype,
+            bounding_box=model.bounding_box.with_buffer(buffer),
+            nelements=nelements,
+            fold=fold,
+            name=foliation_name,
+            svario=svario,
+            model=model,
+            **kwargs,
+        )
+        if data is None:
+            data = model.data.loc[model.data["feature_name"] == foliation_name]
+        if data.shape[0] == 0:
+            logger.warning(f"No data for {foliation_name}, skipping")
+            return
+        series_builder.add_data_from_data_frame(model.prepare_data(data, include_feature_name=False))
+
+        model._add_faults(series_builder)
+        # build feature
+
+        kwargs["tol"] = tol
+
+        series_feature = series_builder.feature
+        series_builder.update_build_arguments(kwargs)
+        series_feature.type = FeatureType.FOLDED
+        series_feature.fold = fold
+
+        model._add_feature(series_feature, index)
+        return series_feature
+
+    @staticmethod
+    def build_folded_fold_frame(
+        model,
+        fold_frame_name: str,
+        *,
+        index: int | None = None,
+        data: pd.DataFrame | None = None,
+        interpolatortype="FDI",
+        nelements=LoopStructuralConfig.nelements,
+        fold_frame=None,
+        tol=None,
+        **kwargs,
+    ):
+        """
+
+        Parameters
+        ----------
+        fold_frame_name : string
+            name of the feature to be added
+        fold_frame_data : pandas data frame, optional
+            data frame containing the fold frame data, if not specified uses the model data
+        interpolatortype : str
+            the type of interpolator to use, default is 'FDI' (unused) 5/6/2025
+        fold_frame : StructuralFrame, optional
+            the fold frame for the fold if not specified uses last feature added
+        nelements : int
+            the number of elements to use in the fold frame
+        tol : float, optional
+            tolerance for the solver, if not specified uses the model default
+        **kwargs : dict
+            additional parameters to be passed to the
+            :class:`LoopStructural.modelling.features.builders.StructuralFrameBuilder`
+            and :meth:`LoopStructural.modelling.features.builders.StructuralFrameBuilder.setup`
+
+        Returns
+        -------
+        fold_frame : FoldFrame
+            created fold frame
+
+        Notes
+        -----
+        This function build a structural frame where the first coordinate is constrained
+        with a fold interpolator.
+        Keyword arguments can be included to constrain
+
+        - :meth:`LoopStructural.GeologicalModel.get_interpolator`
+        - :class:`LoopStructural.StructuralFrameBuilder`
+        - :meth:`LoopStructural.StructuralFrameBuilder.setup`
+         - Building a folded foliation uses the fold interpolation code from Laurent et al., 2016
+        and fold profile fitting from Grose et al., 2017. For more information about the fold modelling
+        see :class:`LoopStructural.modelling.features.fold.FoldEvent`,
+        :class:`LoopStructural.modelling.features.builders.FoldedFeatureBuilder`
+        """
+
+        if tol is None:
+            tol = model.tol
+
+        if fold_frame is None:
+            logger.info("Using last feature as fold frame")
+            fold_frame = model.features[-1]
+        if not isinstance(fold_frame, FoldFrame):
+            raise TypeError("Please specify a FoldFrame")
+        fold = FoldEvent(fold_frame, name=f"Fold_{fold_frame_name}")
+
+        interpolatortypes = [
+            "DFI",
+            "FDI",
+            "FDI",
+        ]
+        fold_frame_builder = StructuralFrameBuilder(
+            interpolatortype=interpolatortypes,
+            bounding_box=model.bounding_box.with_buffer(kwargs.get("buffer", 0.1)),
+            nelements=[nelements, nelements, nelements],
+            name=fold_frame_name,
+            fold=fold,
+            frame=FoldFrame,
+            model=model,
+            **kwargs,
+        )
+        if data is None:
+            data = model.data[model.data["feature_name"] == fold_frame_name]
+        fold_frame_builder.add_data_from_data_frame(
+            model.prepare_data(data, include_feature_name=False)
+        )
+
+        for i in range(3):
+            model._add_faults(fold_frame_builder[i])
+        # build feature
+        kwargs["frame"] = FoldFrame
+        kwargs["tol"] = tol
+        fold_frame_builder.build(**kwargs)
+        folded_fold_frame = fold_frame_builder.frame
+        folded_fold_frame.builder = fold_frame_builder
+
+        folded_fold_frame.type = FeatureType.STRUCTURALFRAME
+
+        model._add_feature(folded_fold_frame, index=index)
+
+        return folded_fold_frame
+
+    @staticmethod
+    def validate_intrusion_inputs(
+        intrusion_name,
+        intrusion_frame_name,
+        intrusion_data,
+        intrusion_frame_data,
+        intrusion_frame_parameters,
+    ):
+        """Fail fast, at the `create_and_add_intrusion` boundary, with a clear
+        message naming the missing piece -- instead of a bare `KeyError`
+        several calls deep inside `IntrusionFrameBuilder`/`IntrusionBuilder`
+        once building has already started. See ``INTRUSIONS.md`` finding 4.
+        """
+        if intrusion_data.empty:
+            raise ValueError(
+                f"No data found for intrusion '{intrusion_name}': check that "
+                "model.data contains rows with feature_name == "
+                f"'{intrusion_name}'"
+            )
+        if intrusion_frame_data.empty:
+            raise ValueError(
+                f"No data found for intrusion frame '{intrusion_frame_name}': "
+                "check that model.data contains rows with feature_name == "
+                f"'{intrusion_frame_name}'"
+            )
+        required_columns = ["intrusion_contact_type", "intrusion_side"]
+        missing_columns = [c for c in required_columns if c not in intrusion_data.columns]
+        if missing_columns:
+            raise ValueError(
+                f"Intrusion data for '{intrusion_name}' is missing required "
+                f"column(s) {missing_columns}: 'intrusion_contact_type' marks "
+                "each point as 'roof'/'top' or 'floor'/'base', and "
+                "'intrusion_side' (boolean) marks points used to constrain "
+                "the lateral extent"
+            )
+        contact_anisotropies = intrusion_frame_parameters.get("contact_anisotropies")
+        if not contact_anisotropies:
+            raise ValueError(
+                "intrusion_frame_parameters['contact_anisotropies'] is "
+                "required: provide a non-empty list of series-type features "
+                "to use as the inflation-gradient proxy for the intrusion "
+                "frame's coordinate 0"
+            )
+
+    @staticmethod
+    def build_intrusion(
+        model,
+        intrusion_name,
+        intrusion_frame_name,
+        *,
+        intrusion_frame_parameters=None,
+        intrusion_lateral_extent_model=None,
+        intrusion_vertical_extent_model=None,
+        geometric_scaling_parameters=None,
+        **kwargs,
+    ):
+        """
+
+        Note
+        -----
+        An intrusion in built in two main steps:
+        (1) Intrusion builder: intrusion builder creates the intrusion structural frame.
+            This object is curvilinear coordinate system of the intrusion constrained with intrusion network points,
+            and flow and inflation measurements (provided by the user).
+            The intrusion network is a representation of the approximated location of roof or floor contact of the intrusion.
+            This object might be constrained using the anisotropies of the host rock if the roof (or floor) contact is not well constrained.
+
+        (2) Intrusion feature: simulation of lateral and vertical extent of intrusion within the model volume.
+            The simulations outcome consist in thresholds distances along the structural frame coordinates
+            that are used to constrained the extent of the intrusion.
+
+        Parameters
+        ----------
+        intrusion_name :  string,
+            name of intrusion feature in model data
+        intrusion_frame_name :  string,
+            name of intrusion frame in model data
+        intrusion_lateral_extent_model = function,
+            geometrical conceptual model for simulation of lateral extent
+        intrusion_vertical_extent_model = function,
+            geometrical conceptual model for simulation of vertical extent
+        intrusion_frame_parameters = dictionary
+
+        kwargs
+
+        Returns
+        -------
+        intrusion feature
+
+        """
+        if intrusion_frame_parameters is None:
+            intrusion_frame_parameters = {}
+        if geometric_scaling_parameters is None:
+            geometric_scaling_parameters = {}
+
+        intrusion_data = model.data[model.data["feature_name"] == intrusion_name].copy()
+        intrusion_frame_data = model.data[model.data["feature_name"] == intrusion_frame_name].copy()
+
+        ModelFeatureFactory.validate_intrusion_inputs(
+            intrusion_name,
+            intrusion_frame_name,
+            intrusion_data,
+            intrusion_frame_data,
+            intrusion_frame_parameters,
+        )
+
+        # -- get variables for intrusion frame interpolation
+        gxxgz = kwargs.get("gxxgz", 0)
+        gxxgy = kwargs.get("gxxgy", 0)
+        gyxgz = kwargs.get("gyxgz", 0)
+
+        interpolatortype = kwargs.get("interpolatortype", "PLI")
+        # buffer = kwargs.get("buffer", 0.1)
+        nelements = kwargs.get("nelements", LoopStructuralConfig.nelements)
+
+        weights = [gxxgz, gxxgy, gyxgz]
+
+        intrusion_frame_builder = IntrusionFrameBuilder(
+            interpolatortype=interpolatortype,
+            bounding_box=model.bounding_box.with_buffer(kwargs.get("buffer", 0.1)),
+            nelements=kwargs.get("nelements", LoopStructuralConfig.nelements),
+            name=intrusion_frame_name,
+            model=model,
+            **kwargs,
+        )
+
+        model._add_faults(intrusion_frame_builder)
+        # intrusion_frame_builder.post_intrusion_faults = faults  # LG unused?
+
+        # -- create intrusion frame using intrusion structures (steps and marginal faults) and flow/inflation measurements
+        if len(intrusion_frame_parameters) == 0:
+            logger.error("Please specify parameters to build intrusion frame")
+        intrusion_frame_builder.set_intrusion_frame_parameters(
+            intrusion_data, intrusion_frame_parameters
+        )
+        intrusion_frame_builder.create_constraints_for_c0()
+
+        intrusion_frame_builder.set_intrusion_frame_data(intrusion_frame_data)
+
+        ## -- create intrusion frame
+        intrusion_frame_builder.build(
+            nelements=nelements,
+            w2=weights[0],
+            w1=weights[1],
+            gyxgz=weights[2],
+        )
+
+        intrusion_frame = intrusion_frame_builder.frame
+
+        # -- create intrusion builder to compute distance thresholds along the frame coordinates
+        intrusion_builder = IntrusionBuilder(
+            intrusion_frame,
+            model=model,
+            # interpolator=interpolator,
+            name=f"{intrusion_name}_feature",
+            lateral_extent_model=intrusion_lateral_extent_model,
+            vertical_extent_model=intrusion_vertical_extent_model,
+            **kwargs,
+        )
+        intrusion_builder.set_data_for_extent_calculation(intrusion_data)
+
+        intrusion_builder.update_build_arguments(
+            {
+                "geometric_scaling_parameters": geometric_scaling_parameters,
+            }
+        )
+
+        intrusion_feature = intrusion_builder.feature
+        model._add_feature(intrusion_feature)
+
+        return intrusion_feature
+
+    @staticmethod
+    def build_domain_fault(
+        model,
+        fault_surface_data,
+        *,
+        nelements=LoopStructuralConfig.nelements,
+        interpolatortype="FDI",
+        index: int | None = None,
+        **kwargs,
+    ):
+        """
+        Parameters
+        ----------
+        fault_surface_data : string
+            name of the domain fault data in the data frame
+
+        Returns
+        -------
+        domain_Fault : GeologicalFeature
+            the created domain fault
+
+        Notes
+        -----
+        * :meth:`LoopStructural.GeologicalModel.get_interpolator`
+
+        """
+        domain_fault_feature_builder = GeologicalFeatureBuilder(
+            bounding_box=model.bounding_box,
+            interpolatortype=interpolatortype,
+            nelements=nelements,
+            name=fault_surface_data,
+            model=model,
+            **kwargs,
+        )
+
+        # add data
+        unconformity_data = model.data.loc[model.data["feature_name"] == fault_surface_data]
+
+        domain_fault_feature_builder.add_data_from_data_frame(unconformity_data)
+        # look through existing features if there is a fault before an
+        # unconformity
+        # then add to the feature, once we get to an unconformity stop
+        model._add_faults(domain_fault_feature_builder)
+
+        # build feature
+        domain_fault = domain_fault_feature_builder.feature
+        domain_fault_feature_builder.update_build_arguments(kwargs)
+        domain_fault.type = FeatureType.DOMAINFAULT
+        model._add_feature(domain_fault, index=index)
+        model._add_domain_fault_below(domain_fault)
+
+        domain_fault_uc = UnconformityFeature(domain_fault, 0)
+        # iterate over existing features and add the unconformity as a region
+        # so the feature is only evaluated where the unconformity is positive
+        return domain_fault_uc
+
+    @staticmethod
+    def build_fault(
+        model,
+        fault_name: str,
+        displacement: float,
+        *,
+        index: int | None = None,
+        data: pd.DataFrame | None = None,
+        interpolatortype="FDI",
+        tol=None,
+        fault_slip_vector=None,
+        fault_normal_vector=None,
+        fault_center=None,
+        major_axis=None,
+        minor_axis=None,
+        intermediate_axis=None,
+        faultfunction="BaseFault",
+        faults=None,
+        force_mesh_geometry: bool = False,
+        points: bool = False,
+        fault_buffer=0.2,
+        fault_trace_anisotropy=0.0,
+        fault_dip=90,
+        fault_dip_anisotropy=0.0,
+        fault_pitch=None,
+        **kwargs,
+    ):
+        """
+        Parameters
+        ----------
+        fault_name : string
+            name of the fault surface data in the dataframe
+        displacement : displacement magnitude
+            displacement magnitude of the fault, in model units
+        fault_data : pd.DataFrame, optional
+            data frame containing the fault data, if not specified uses the model data
+        major_axis : [type], optional
+            [description], by default None
+        minor_axis : [type], optional
+            [description], by default None
+        intermediate_axis : [type], optional
+            [description], by default None
+        kwargs : additional kwargs for Fault and interpolators
+
+        Returns
+        -------
+        fault : FaultSegment
+            created fault
+
+        Notes
+        -----
+        * :meth:`LoopStructural.GeologicalModel.get_interpolator`
+        * :class:`LoopStructural.modelling.features.builders.FaultBuilder`
+        * :meth:`LoopStructural.modelling.features.builders.FaultBuilder.setup`
+        """
+        if faults is None:
+            faults = []
+        if "fault_extent" in kwargs and major_axis is None:
+            major_axis = kwargs["fault_extent"]
+        if "fault_influence" in kwargs and minor_axis is None:
+            minor_axis = kwargs["fault_influence"]
+        if "fault_vectical_radius" in kwargs and intermediate_axis is None:
+            intermediate_axis = kwargs["fault_vectical_radius"]
+
+        logger.info(f'Creating fault "{fault_name}"')
+        logger.info(f"Displacement: {displacement}")
+        logger.info(f"Tolerance: {tol}")
+        logger.info(f"Fault function: {faultfunction}")
+        logger.info(f"Fault slip vector: {fault_slip_vector}")
+        logger.info(f"Fault center: {fault_center}")
+        logger.info(f"Major axis: {major_axis}")
+        logger.info(f"Minor axis: {minor_axis}")
+        logger.info(f"Intermediate axis: {intermediate_axis}")
+        if fault_slip_vector is not None:
+            fault_slip_vector = np.array(fault_slip_vector, dtype="float")
+        if fault_center is not None:
+            fault_center = np.array(fault_center, dtype="float")
+
+        for k, v in kwargs.items():
+            logger.info(f"{k}: {v}")
+
+        if tol is None:
+            tol = model.tol
+            # divide the tolerance by half of the minor axis, as this is the equivalent of the distance
+            # of the unit vector
+            # if minor_axis:
+            # tol *= 0.1*minor_axis
+
+        if displacement == 0:
+            logger.warning(f"{fault_name} displacement is 0")
+
+        if "data_region" in kwargs:
+            kwargs.pop("data_region")
+            logger.error("kwarg data_region currently not supported, disabling")
+        displacement_scaled = displacement
+        fault_frame_builder = FaultBuilder(
+            interpolatortype,
+            bounding_box=model.bounding_box,
+            nelements=kwargs.pop("nelements", LoopStructuralConfig.nelements),
+            name=fault_name,
+            model=model,
+            **kwargs,
+        )
+        if data is None:
+            data = model.data.loc[model.data["feature_name"] == fault_name]
+        if data.shape[0] == 0:
+            logger.warning(f"No data for {fault_name}, skipping")
+            return
+
+        model._add_faults(fault_frame_builder, features=faults)
+        # add data
+
+        if fault_center is not None and ~np.isnan(fault_center).any():
+            fault_center = model.scale(fault_center, inplace=False)
+        # Keep the supplied fault-axis values unchanged; the previous self-assignment
+        # was only present to satisfy a linter and did not affect behavior.
+        fault_frame_builder.create_data_from_geometry(
+            fault_frame_data=model.prepare_data(data, include_feature_name=False),
+            fault_center=fault_center,
+            fault_normal_vector=fault_normal_vector,
+            fault_slip_vector=fault_slip_vector,
+            minor_axis=minor_axis,
+            major_axis=major_axis,
+            intermediate_axis=intermediate_axis,
+            points=points,
+            force_mesh_geometry=force_mesh_geometry,
+            fault_buffer=fault_buffer,
+            fault_trace_anisotropy=fault_trace_anisotropy,
+            fault_dip=fault_dip,
+            fault_dip_anisotropy=fault_dip_anisotropy,
+            fault_pitch=fault_pitch,
+        )
+        if "force_mesh_geometry" not in kwargs:
+            fault_frame_builder.set_mesh_geometry(kwargs.get("fault_buffer", 0.2), 0)
+        if "splay" in kwargs and "splayregion" in kwargs:
+            fault_frame_builder.add_splay(kwargs["splay"], kwargs["splayregion"])
+
+        kwargs["tol"] = tol
+        fault_frame_builder.build(**kwargs)
+        fault = fault_frame_builder.frame
+        fault.displacement = displacement_scaled
+        fault.faultfunction = faultfunction
+
+        for f in reversed(model.features):
+            if f.type == FeatureType.UNCONFORMITY:
+                fault.add_region(f)
+                break
+        if displacement == 0:
+            fault.type = FeatureType.INACTIVEFAULT
+        model._add_feature(fault, index=index)
+
+        return fault
diff --git a/LoopStructural/modelling/core/_model_relationships.py b/LoopStructural/modelling/core/_model_relationships.py
new file mode 100644
index 00000000..dcb59c53
--- /dev/null
+++ b/LoopStructural/modelling/core/_model_relationships.py
@@ -0,0 +1,177 @@
+"""Fault/unconformity relationship bookkeeping for GeologicalModel (see API.md).
+
+Extracted from GeologicalModel to separate feature-stack relationship logic
+(domain faults, unconformities) from feature-container state and
+construction. GeologicalModel's @public_api methods (``add_unconformity``,
+``add_onlap_unconformity``) stay defined directly on the class -- their
+__qualname__ is part of the CI-checked stable API surface -- and delegate to
+the staticmethods here. The private ``_add_faults``/``_add_domain_fault_*``
+helpers also delegate, since they're called throughout GeologicalModel's
+feature-construction methods.
+"""
+
+from ...modelling.features import FeatureType, UnconformityFeature
+from ...utils import getLogger
+
+logger = getLogger(__name__)
+
+
+class FeatureRelationshipManager:
+    @staticmethod
+    def add_faults(model, feature_builder, features=None):
+        """Adds all existing faults to a geological feature builder
+
+        Parameters
+        ----------
+        model : GeologicalModel
+        feature_builder : GeologicalFeatureBuilder/StructuralFrameBuilder
+            The feature buider to add the faults to
+        features : list, optional
+            A specific list of features rather than all features in the model
+        """
+        if features is None:
+            features = model.features
+        for f in reversed(features):
+            if isinstance(f, str):
+                f = model.__getitem__(f)
+            if f.type == FeatureType.FAULT:
+                feature_builder.add_fault(f)
+
+    @staticmethod
+    def add_domain_fault_above(model, feature):
+        """
+        Looks through the feature list and adds any domain faults to the feature. The domain fault masks everything
+        where the fault scalar field is < 0 as being active when added to feature.
+
+        Parameters
+        ----------
+        model : GeologicalModel
+        feature : GeologicalFeatureBuilder
+            the feature being added to the model where domain faults should be added
+        """
+        for f in reversed(model.features):
+            if f.name == feature.name:
+                continue
+            if f.type == "domain_fault":
+                feature.add_region(lambda pos, fault=f: fault.evaluate_value(pos) < 0)
+                break
+
+    @staticmethod
+    def add_domain_fault_below(model, domain_fault):
+        """
+        Looks through the feature list and adds any the domain_fault to the features
+        that already exist in the stack until an unconformity is reached. domain faults
+        to the feature. The domain fault masks everything where the fault scalar field
+         is < 0 as being active when added to feature.
+
+        Parameters
+        ----------
+        model : GeologicalModel
+        domain_fault : GeologicalFeatureBuilder
+            the feature being added to the model where domain faults should be added
+        """
+        for f in reversed(model.features):
+            if f.name == domain_fault.name:
+                continue
+            f.add_region(lambda pos: domain_fault.evaluate_value(pos) > 0)
+            if f.type == FeatureType.UNCONFORMITY:
+                break
+
+    @staticmethod
+    def add_unconformity_above(model, feature):
+        """
+        Adds a region to the feature to prevent the value from being
+        interpolated where the unconformities exists above e.g.
+        if there is another feature above and the unconformity is at 0
+        then the features added below (after) will only be visible where the
+        uncomformity is <0
+
+        Parameters
+        ----------
+        model : GeologicalModel
+        feature - GeologicalFeature
+        """
+
+        if feature.type == FeatureType.FAULT:
+            return
+        for f in reversed(model.features):
+            if f.type == FeatureType.UNCONFORMITY and f.name != feature.name:
+                logger.info(f"Adding {f.name} as unconformity to {feature.name}")
+                feature.add_region(f)
+            if f.type == FeatureType.ONLAPUNCONFORMITY and f.name != feature.name:
+                feature.add_region(f)
+                break
+
+    @staticmethod
+    def add_unconformity(model, feature, value, index=None):
+        """
+        Use an existing feature to add an unconformity to the model.
+
+        Parameters
+        ----------
+        model : GeologicalModel
+        feature : GeologicalFeature
+            existing geological feature
+        value : float
+            scalar value of isosurface that represents
+
+        Returns
+        -------
+        unconformity : GeologicalFeature
+            unconformity feature
+        """
+        logger.debug(f"Adding {feature.name} as unconformity at {value}")
+        if feature is None:
+            logger.warning("Cannot add unconformtiy, base feature is None")
+            return
+        # look backwards through features and add the unconformity as a region until
+        # we get to an unconformity
+        uc_feature = UnconformityFeature(feature, value)
+        feature.add_region(uc_feature.inverse())
+        for f in reversed(model.features):
+            if f.type == FeatureType.UNCONFORMITY:
+                logger.debug(f"Reached unconformity {f.name}")
+                break
+            logger.debug(f"Adding {uc_feature.name} as unconformity to {f.name}")
+            if f.type == FeatureType.FAULT or f.type == FeatureType.INACTIVEFAULT:
+                continue
+            if f == feature:
+                continue
+            else:
+                f.add_region(uc_feature)
+        # now add the unconformity to the feature list
+        model._add_feature(uc_feature, index=index)
+        return uc_feature
+
+    @staticmethod
+    def add_onlap_unconformity(model, feature, value, index=None):
+        """
+        Use an existing feature to add an unconformity to the model.
+
+        Parameters
+        ----------
+        model : GeologicalModel
+        feature : GeologicalFeature
+            existing geological feature
+        value : float
+            scalar value of isosurface that represents
+
+        Returns
+        -------
+        unconformity_feature : GeologicalFeature
+            the created unconformity
+        """
+        feature.regions = []
+        uc_feature = UnconformityFeature(feature, value, False, onlap=True)
+        feature.add_region(uc_feature.inverse())
+        for f in reversed(model.features):
+            if f.type in (FeatureType.UNCONFORMITY, FeatureType.ONLAPUNCONFORMITY):
+                logger.debug(f"Reached unconformity {f.name}")
+                break
+            if f.type == FeatureType.FAULT or f.type == FeatureType.INACTIVEFAULT:
+                continue
+            if f != feature:
+                f.add_region(uc_feature)
+        model._add_feature(uc_feature.inverse(), index=index)
+
+        return uc_feature
diff --git a/LoopStructural/modelling/core/_model_serializer.py b/LoopStructural/modelling/core/_model_serializer.py
new file mode 100644
index 00000000..cbf4d112
--- /dev/null
+++ b/LoopStructural/modelling/core/_model_serializer.py
@@ -0,0 +1,240 @@
+"""Recipe (JSON) and pickle serialization logic for GeologicalModel.
+
+Extracted from GeologicalModel to separate serialization concerns from
+feature-container orchestration (see API.md). GeologicalModel's
+``@public_api``-decorated methods (``to_dict``, ``to_recipe_dict``,
+``from_recipe_dict``, ``to_recipe_json``, ``from_recipe_json``,
+``save_recipe``, ``load_recipe``, ``to_file``, ``from_file``) stay defined
+directly on the class -- their ``__qualname__`` is part of the CI-checked
+stable API surface (``tests/unit/test_public_api_contract.py``) -- and just
+delegate to the staticmethods here.
+"""
+
+import json
+import pathlib
+
+import pandas as pd
+
+from ...geometry import BoundingBox
+from ...utils import LoopValueError, getLogger
+from ..features import GeologicalFeature, StructuralFrame, UnconformityFeature
+from ..features.fault import FaultSegment
+from .stratigraphic_column import StratigraphicColumn
+
+logger = getLogger(__name__)
+
+
+class ModelSerializer:
+    @staticmethod
+    def feature_recipe_kind(feature):
+        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 to_dict(model):
+        result = {}
+        result["model"] = {}
+        result["model"]["features"] = [f.name for f in model.features]
+        result['model']['bounding_box'] = model.bounding_box.to_dict()
+        result["model"]["stratigraphic_column"] = model.stratigraphic_column
+        return result
+
+    @staticmethod
+    def to_recipe_dict(model, data_reference=None):
+        recipe = {
+            "schema": "LoopStructural.GeologicalModelRecipe",
+            "version": 1,
+            "model": {
+                "bounding_box": model.bounding_box.to_dict(),
+                "stratigraphic_column": model.stratigraphic_column.to_dict(),
+                "data_source": None,
+                "features": [],
+            },
+        }
+        for feature in model.features:
+            feature_entry = {
+                "name": feature.name,
+                "kind": ModelSerializer.feature_recipe_kind(feature),
+                "faults": [
+                    fault.name
+                    for fault in getattr(feature, "faults", [])
+                    if getattr(fault, "name", None)
+                ],
+                "regions": [],
+            }
+            recipe["model"]["features"].append(feature_entry)
+        if data_reference is not None:
+            recipe["model"]["data_source"] = {
+                "kind": "reference",
+                "path": str(pathlib.Path(data_reference)),
+            }
+        elif not model.data.empty:
+            recipe["model"]["data_source"] = {
+                "kind": "inline",
+                "dataframe": model.data.to_dict(orient="split"),
+            }
+        return recipe
+
+    @staticmethod
+    def from_recipe_dict(cls, recipe):
+        if not isinstance(recipe, dict):
+            raise TypeError("recipe must be a dictionary")
+
+        model_data = recipe.get("model", recipe)
+        bounding_box = model_data.get("bounding_box")
+        if isinstance(bounding_box, dict):
+            bounding_box = BoundingBox.from_dict(bounding_box)
+        if not isinstance(bounding_box, BoundingBox):
+            raise TypeError("recipe must include a bounding_box dictionary")
+
+        model = cls(bounding_box)
+
+        data_source = model_data.get("data_source")
+        if isinstance(data_source, dict):
+            kind = data_source.get("kind")
+            if kind == "reference":
+                model.data = pd.read_csv(pathlib.Path(data_source["path"]))
+            elif kind == "inline":
+                dataframe = data_source.get("dataframe")
+                if not isinstance(dataframe, dict):
+                    raise TypeError("inline data_source must include a dataframe dictionary")
+                model.data = pd.DataFrame(**dataframe)
+            elif kind is not None:
+                raise ValueError(f"Unsupported data_source kind: {kind}")
+        elif isinstance(data_source, str):
+            model.data = pd.read_csv(pathlib.Path(data_source))
+        elif data_source is not None:
+            raise TypeError("data_source must be a dictionary, string path, or None")
+
+        stratigraphic_column = model_data.get("stratigraphic_column")
+        if isinstance(stratigraphic_column, dict):
+            model.stratigraphic_column = StratigraphicColumn.from_dict(stratigraphic_column)
+        elif stratigraphic_column is not None:
+            raise TypeError("stratigraphic_column must be a dictionary or None")
+
+        features = model_data.get("features", [])
+        if features is None:
+            features = []
+        if not isinstance(features, list):
+            raise TypeError("features must be a list")
+
+        feature_map = {}
+        for feature_entry in features:
+            if not isinstance(feature_entry, dict):
+                raise TypeError("each feature entry must be a dictionary")
+            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)
+            if feature is None:
+                raise ValueError(f"Could not recreate feature '{feature_name}' from recipe")
+            feature_map[feature_name] = feature
+
+        for feature_entry in features:
+            feature_name = feature_entry.get("name")
+            fault_names = feature_entry.get("faults", [])
+            if not isinstance(fault_names, list):
+                raise TypeError("faults must be a list")
+            if fault_names:
+                feature = feature_map[feature_name]
+                feature.faults = [feature_map[name] for name in fault_names if name in feature_map]
+
+        return model
+
+    @staticmethod
+    def to_recipe_json(model, data_reference=None, indent=2):
+        recipe = ModelSerializer.to_recipe_dict(model, data_reference=data_reference)
+        return json.dumps(recipe, indent=indent)
+
+    @staticmethod
+    def from_recipe_json(cls, json_str):
+        if not isinstance(json_str, str):
+            raise TypeError("json_str must be a string")
+        try:
+            recipe = json.loads(json_str)
+        except json.JSONDecodeError as e:
+            raise TypeError(f"json_str is not valid JSON: {e}")
+        return ModelSerializer.from_recipe_dict(cls, recipe)
+
+    @staticmethod
+    def save_recipe(model, filename, data_reference=None):
+        filename = pathlib.Path(filename)
+        recipe = ModelSerializer.to_recipe_dict(model, data_reference=data_reference)
+        with open(filename, "w") as f:
+            json.dump(recipe, f, indent=2)
+        logger.info(f"Recipe saved to {filename}")
+
+    @staticmethod
+    def load_recipe(cls, filename):
+        filename = pathlib.Path(filename)
+        if not filename.exists():
+            raise FileNotFoundError(f"Recipe file not found: {filename}")
+        with open(filename, "r") as f:
+            recipe = json.load(f)
+        logger.info(f"Recipe loaded from {filename}")
+        return ModelSerializer.from_recipe_dict(cls, recipe)
+
+    @staticmethod
+    def to_file(model, file):
+        try:
+            import dill as pickle
+        except ImportError:
+            logger.error("Cannot write to file, dill not installed \n" "pip install dill")
+            return
+        try:
+            logger.info(f"Writing GeologicalModel to: {file}")
+            with open(file, "wb") as handle:
+                pickle.dump(model, handle)
+        except pickle.PicklingError:
+            logger.error("Error saving file")
+
+    @staticmethod
+    def from_file(cls, file, allow_pickle: bool = True):
+        if not allow_pickle:
+            raise LoopValueError(
+                "Pickle-based loading is disabled (allow_pickle=False). "
+                f"Refusing to unpickle '{file}' because deserialising untrusted "
+                "pickle/dill data can execute arbitrary code. If you generated "
+                "this file yourself and trust its contents, call "
+                "GeologicalModel.from_file(file, allow_pickle=True). Otherwise, "
+                "use the JSON-based GeologicalModel.from_recipe_dict "
+                "(paired with GeologicalModel.to_recipe_dict) as a safe "
+                "alternative serialisation format."
+            )
+        logger.warning(
+            f"Loading GeologicalModel from '{file}' using dill/pickle. "
+            "Only load model files from trusted sources: deserialising a "
+            "pickle file can execute arbitrary code. Pass allow_pickle=False "
+            "to refuse pickle-based loading, or use "
+            "GeologicalModel.from_recipe_dict for untrusted/JSON-based input."
+        )
+        try:
+            import dill as pickle
+        except ImportError:
+            logger.error("Cannot import from file, dill not installed")
+            return None
+        path = pathlib.Path(file)
+        if not path.is_file():
+            raise LoopValueError(f"Cannot load model, file does not exist: {file}")
+        try:
+            with open(path, "rb") as f:
+                model = pickle.load(f)
+        except Exception as e:
+            logger.error(f"Failed to load model from {file}: {e}")
+            raise LoopValueError(f"Failed to load model from {file}: {e}") from e
+        if isinstance(model, cls):
+            logger.info("GeologicalModel initialised from file")
+            return model
+        else:
+            logger.error(f"{file} does not contain a geological model")
+            return None
diff --git a/LoopStructural/modelling/core/geological_model.py b/LoopStructural/modelling/core/geological_model.py
index 72e66120..36c6fb00 100644
--- a/LoopStructural/modelling/core/geological_model.py
+++ b/LoopStructural/modelling/core/geological_model.py
@@ -3,15 +3,14 @@
 """
 from __future__ import annotations
 
-import json
-import pathlib
+import warnings
 
 import numpy as np
 import pandas as pd
 
 from LoopStructural import LoopStructuralConfig
 
-from ...geometry import BoundingBox, StructuredGrid
+from ...geometry import BoundingBox
 from ...modelling.features import (
     BaseFeature,
     FeatureType,
@@ -19,18 +18,8 @@
     StructuralFrame,
     UnconformityFeature,
 )
-from ...modelling.features.builders import (
-    FaultBuilder,
-    FoldedFeatureBuilder,
-    GeologicalFeatureBuilder,
-    StructuralFrameBuilder,
-)
 from ...modelling.features.fault import FaultSegment
-from ...modelling.features.fold import (
-    FoldEvent,
-    FoldFrame,
-)
-from ...modelling.intrusions import IntrusionBuilder, IntrusionFrameBuilder
+from ...modelling.features.fold import FoldFrame
 from ...utils import LoopValueError, getLogger, public_api, strikedip2vector, timed_stage
 from ...utils.helper import (
     all_heading,
@@ -43,6 +32,10 @@
     convert_feature_to_structural_frame as _convert_feature_to_structural_frame,
 )
 from ._feature_registry import FeatureBuilderRegistry
+from ._model_exporter import ModelExporter
+from ._model_feature_factory import ModelFeatureFactory
+from ._model_relationships import FeatureRelationshipManager
+from ._model_serializer import ModelSerializer
 from .stratigraphic_column import StratigraphicColumn
 
 logger = getLogger(__name__)
@@ -153,25 +146,7 @@ def to_dict(self):
         json : str
             json string of the geological model
         """
-        json = {}
-        json["model"] = {}
-        json["model"]["features"] = [f.name for f in self.features]
-        json['model']['bounding_box'] = self.bounding_box.to_dict()
-        json["model"]["stratigraphic_column"] = self.stratigraphic_column
-        # json["features"] = [f.to_json() for f in self.features]
-        return json
-
-    @staticmethod
-    def _feature_recipe_kind(feature):
-        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()
+        return ModelSerializer.to_dict(self)
 
     @public_api(tier="provisional")
     def to_recipe_dict(self, data_reference=None):
@@ -181,110 +156,13 @@ def to_recipe_dict(self, data_reference=None):
         box, stratigraphic column, and either inline model data or a file
         reference to it.
         """
-        recipe = {
-            "schema": "LoopStructural.GeologicalModelRecipe",
-            "version": 1,
-            "model": {
-                "bounding_box": self.bounding_box.to_dict(),
-                "stratigraphic_column": self.stratigraphic_column.to_dict(),
-                "data_source": None,
-                "features": [],
-            },
-        }
-        for feature in self.features:
-            feature_entry = {
-                "name": feature.name,
-                "kind": self._feature_recipe_kind(feature),
-                "faults": [
-                    fault.name
-                    for fault in getattr(feature, "faults", [])
-                    if getattr(fault, "name", None)
-                ],
-                "regions": [],
-            }
-            recipe["model"]["features"].append(feature_entry)
-        if data_reference is not None:
-            recipe["model"]["data_source"] = {
-                "kind": "reference",
-                "path": str(pathlib.Path(data_reference)),
-            }
-        elif not self.data.empty:
-            recipe["model"]["data_source"] = {
-                "kind": "inline",
-                "dataframe": self.data.to_dict(orient="split"),
-            }
-        return recipe
+        return ModelSerializer.to_recipe_dict(self, data_reference=data_reference)
 
     @classmethod
     @public_api(tier="provisional")
     def from_recipe_dict(cls, recipe):
         """Rebuild a geological model from a recipe dictionary."""
-        if not isinstance(recipe, dict):
-            raise TypeError("recipe must be a dictionary")
-
-        model_data = recipe.get("model", recipe)
-        bounding_box = model_data.get("bounding_box")
-        if isinstance(bounding_box, dict):
-            bounding_box = BoundingBox.from_dict(bounding_box)
-        if not isinstance(bounding_box, BoundingBox):
-            raise TypeError("recipe must include a bounding_box dictionary")
-
-        model = cls(bounding_box)
-
-        data_source = model_data.get("data_source")
-        if isinstance(data_source, dict):
-            kind = data_source.get("kind")
-            if kind == "reference":
-                model.data = pd.read_csv(pathlib.Path(data_source["path"]))
-            elif kind == "inline":
-                dataframe = data_source.get("dataframe")
-                if not isinstance(dataframe, dict):
-                    raise TypeError("inline data_source must include a dataframe dictionary")
-                model.data = pd.DataFrame(**dataframe)
-            elif kind is not None:
-                raise ValueError(f"Unsupported data_source kind: {kind}")
-        elif isinstance(data_source, str):
-            model.data = pd.read_csv(pathlib.Path(data_source))
-        elif data_source is not None:
-            raise TypeError("data_source must be a dictionary, string path, or None")
-
-        stratigraphic_column = model_data.get("stratigraphic_column")
-        if isinstance(stratigraphic_column, dict):
-            model.stratigraphic_column = StratigraphicColumn.from_dict(stratigraphic_column)
-        elif stratigraphic_column is not None:
-            raise TypeError("stratigraphic_column must be a dictionary or None")
-
-        features = model_data.get("features", [])
-        if features is None:
-            features = []
-        if not isinstance(features, list):
-            raise TypeError("features must be a list")
-
-        feature_map = {}
-        for feature_entry in features:
-            if not isinstance(feature_entry, dict):
-                raise TypeError("each feature entry must be a dictionary")
-            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)
-            if feature is None:
-                raise ValueError(f"Could not recreate feature '{feature_name}' from recipe")
-            feature_map[feature_name] = feature
-
-        for feature_entry in features:
-            feature_name = feature_entry.get("name")
-            fault_names = feature_entry.get("faults", [])
-            if not isinstance(fault_names, list):
-                raise TypeError("faults must be a list")
-            if fault_names:
-                feature = feature_map[feature_name]
-                feature.faults = [feature_map[name] for name in fault_names if name in feature_map]
-
-        return model
+        return ModelSerializer.from_recipe_dict(cls, recipe)
 
     @public_api(tier="provisional")
     def to_recipe_json(self, data_reference=None, indent=2):
@@ -303,8 +181,7 @@ def to_recipe_json(self, data_reference=None, indent=2):
         str
             JSON-formatted recipe string.
         """
-        recipe = self.to_recipe_dict(data_reference=data_reference)
-        return json.dumps(recipe, indent=indent)
+        return ModelSerializer.to_recipe_json(self, data_reference=data_reference, indent=indent)
 
     @classmethod
     @public_api(tier="provisional")
@@ -326,13 +203,7 @@ def from_recipe_json(cls, json_str):
         TypeError
             If json_str is not a string or does not parse as valid JSON.
         """
-        if not isinstance(json_str, str):
-            raise TypeError("json_str must be a string")
-        try:
-            recipe = json.loads(json_str)
-        except json.JSONDecodeError as e:
-            raise TypeError(f"json_str is not valid JSON: {e}")
-        return cls.from_recipe_dict(recipe)
+        return ModelSerializer.from_recipe_json(cls, json_str)
 
     @public_api(tier="provisional")
     def save_recipe(self, filename, data_reference=None):
@@ -346,11 +217,7 @@ def save_recipe(self, filename, data_reference=None):
             Path to an external CSV file to reference instead of embedding
             data inline in the JSON. If None, data is embedded.
         """
-        filename = pathlib.Path(filename)
-        recipe = self.to_recipe_dict(data_reference=data_reference)
-        with open(filename, "w") as f:
-            json.dump(recipe, f, indent=2)
-        logger.info(f"Recipe saved to {filename}")
+        ModelSerializer.save_recipe(self, filename, data_reference=data_reference)
 
     @classmethod
     @public_api(tier="provisional")
@@ -367,13 +234,7 @@ def load_recipe(cls, filename):
         GeologicalModel
             The reconstructed geological model.
         """
-        filename = pathlib.Path(filename)
-        if not filename.exists():
-            raise FileNotFoundError(f"Recipe file not found: {filename}")
-        with open(filename, "r") as f:
-            recipe = json.load(f)
-        logger.info(f"Recipe loaded from {filename}")
-        return cls.from_recipe_dict(recipe)
+        return ModelSerializer.load_recipe(cls, filename)
 
     def __str__(self):
         return f"GeologicalModel with {len(self.features)} features"
@@ -529,44 +390,7 @@ def from_file(cls, file, allow_pickle: bool = True):
         GeologicalModel
             the geological model object
         """
-        if not allow_pickle:
-            raise LoopValueError(
-                "Pickle-based loading is disabled (allow_pickle=False). "
-                f"Refusing to unpickle '{file}' because deserialising untrusted "
-                "pickle/dill data can execute arbitrary code. If you generated "
-                "this file yourself and trust its contents, call "
-                "GeologicalModel.from_file(file, allow_pickle=True). Otherwise, "
-                "use the JSON-based GeologicalModel.from_recipe_dict "
-                "(paired with GeologicalModel.to_recipe_dict) as a safe "
-                "alternative serialisation format."
-            )
-        logger.warning(
-            f"Loading GeologicalModel from '{file}' using dill/pickle. "
-            "Only load model files from trusted sources: deserialising a "
-            "pickle file can execute arbitrary code. Pass allow_pickle=False "
-            "to refuse pickle-based loading, or use "
-            "GeologicalModel.from_recipe_dict for untrusted/JSON-based input."
-        )
-        try:
-            import dill as pickle
-        except ImportError:
-            logger.error("Cannot import from file, dill not installed")
-            return None
-        path = pathlib.Path(file)
-        if not path.is_file():
-            raise LoopValueError(f"Cannot load model, file does not exist: {file}")
-        try:
-            with open(path, "rb") as f:
-                model = pickle.load(f)
-        except Exception as e:
-            logger.error(f"Failed to load model from {file}: {e}")
-            raise LoopValueError(f"Failed to load model from {file}: {e}") from e
-        if isinstance(model, GeologicalModel):
-            logger.info("GeologicalModel initialised from file")
-            return model
-        else:
-            logger.error(f"{file} does not contain a geological model")
-            return None
+        return ModelSerializer.from_file(cls, file, allow_pickle=allow_pickle)
 
     def __getitem__(self, feature_name):
         """Accessor for feature in features using feature_name_index
@@ -680,17 +504,7 @@ def to_file(self, file):
         file : string
             path to file location
         """
-        try:
-            import dill as pickle
-        except ImportError:
-            logger.error("Cannot write to file, dill not installed \n" "pip install dill")
-            return
-        try:
-            logger.info(f"Writing GeologicalModel to: {file}")
-            with open(file, "wb") as handle:
-                pickle.dump(self, handle)
-        except pickle.PicklingError:
-            logger.error("Error saving file")
+        ModelSerializer.to_file(self, file)
 
     def _add_feature(self, feature, index: int | None = None):
         """
@@ -721,9 +535,9 @@ def _add_feature(self, feature, index: int | None = None):
                 self.features.append(feature)
                 self.feature_name_index[feature.name] = len(self.features) - 1
                 logger.info(f"Adding {feature.name} to model at location {len(self.features)}")
-        self._add_domain_fault_above(feature)
+        FeatureRelationshipManager.add_domain_fault_above(self, feature)
         if feature.type == FeatureType.INTERPOLATED:
-            self._add_unconformity_above(feature)
+            FeatureRelationshipManager.add_unconformity_above(self, feature)
         feature.model = self
 
     def data_for_feature(self, feature_name: str) -> pd.DataFrame:
@@ -830,13 +644,15 @@ def set_stratigraphic_column(self, stratigraphic_column, cmap="tab20"):
         }
 
         """
+        warnings.warn(
+            "set_stratigraphic_column is deprecated, use model.stratigraphic_column.add_units instead",
+            DeprecationWarning,
+            stacklevel=2,
+        )
         self.stratigraphic_column.clear(basement=False)
         # if the colour for a unit hasn't been specified we can just sample from
         # a colour map e.g. tab20
         logger.info("Adding stratigraphic column to model")
-        raise DeprecationWarning(
-            "set_stratigraphic_column is deprecated, use model.stratigraphic_column.add_units instead"
-        )
         for i, g in enumerate(stratigraphic_column.keys()):
             if g == 'faults':
                 logger.info('Not adding faults to stratigraphic column')
@@ -909,7 +725,8 @@ def create_and_add_foliation(
     ):
         """Create a foliation feature and add it to the model.
 
-        See :meth:`_build_foliation` for parameter documentation. Thin
+        See :meth:`~._model_feature_factory.ModelFeatureFactory.build_foliation`
+        for parameter documentation. Thin
         wrapper around :meth:`create_and_add_feature` (see ``API.md``);
         kept as a stable, unchanged entry point.
         """
@@ -925,86 +742,6 @@ def create_and_add_foliation(
             **kwargs,
         )
 
-    def _build_foliation(
-        self,
-        series_surface_name: str,
-        *,
-        index: int | None = None,
-        data: pd.DataFrame | None = None,
-        interpolatortype: str = "FDI",
-        nelements: int = LoopStructuralConfig.nelements,
-        tol=None,
-        faults=None,
-        **kwargs,
-    ):
-        """
-        Parameters
-        ----------
-        series_surface_name : string
-            corresponding to the feature_name in the data
-        series_surface_data : pd.DataFrame, optional
-            data frame containing the surface data
-        interpolatortype : str
-            the type of interpolator to use, default is 'FDI'
-        nelements : int
-            the number of elements to use in the series surface
-        tol : float, optional
-            tolerance for the solver, if not specified uses the model default
-        faults : list, optional
-            list of faults to be used in the series surface, if not specified uses the model faults
-        kwargs
-
-        Returns
-        -------
-        feature : GeologicalFeature
-            the created geological feature
-
-        Notes
-        ------
-        This function creates an instance of a
-        :class:`LoopStructural.modelling.features.builders.GeologicalFeatureBuilder` and will return
-        a :class:`LoopStructural.modelling.features.builders.GeologicalFeature`
-        The feature is not interpolated until either
-        :meth:`LoopStructural.modelling.features.builders.GeologicalFeature.evaluate_value` is called or
-        :meth:`LoopStructural.modelling.core.GeologicalModel.update`
-
-        An interpolator will be chosen by calling :meth:`LoopStructural.GeologicalModel.get_interpolator`
-
-        """
-
-        # if tol is not specified use the model default
-        if tol is None:
-            tol = self.tol
-
-        series_builder = GeologicalFeatureBuilder(
-            bounding_box=self.bounding_box,
-            interpolatortype=interpolatortype,
-            nelements=nelements,
-            name=series_surface_name,
-            model=self,
-            **kwargs,
-        )
-        # add data
-        if data is None:
-            data = self.data.loc[self.data["feature_name"] == series_surface_name]
-
-        if data.shape[0] == 0:
-            logger.warning("No data for {series_surface_data}, skipping")
-            return
-        series_builder.add_data_from_data_frame(self.prepare_data(data, include_feature_name=False))
-        self._add_faults(series_builder, features=faults)
-
-        # build feature
-        # series_feature = series_builder.build(**kwargs)
-        series_feature = series_builder.feature
-        series_builder.update_build_arguments(kwargs | {"domain": True, 'tol': tol})
-        # this support is built for the entire model domain? Possibly would
-        # could just pass a regular grid of points - mask by any above unconformities??
-
-        series_feature.type = FeatureType.INTERPOLATED
-        self._add_feature(series_feature, index=index)
-        return series_feature
-
     @public_api(tier="stable")
     def create_and_add_fold_frame(
         self,
@@ -1020,7 +757,8 @@ def create_and_add_fold_frame(
     ):
         """Create a fold frame and add it to the model.
 
-        See :meth:`_build_fold_frame` for parameter documentation. Thin
+        See :meth:`~._model_feature_factory.ModelFeatureFactory.build_fold_frame`
+        for parameter documentation. Thin
         wrapper around :meth:`create_and_add_feature` (see ``API.md``);
         kept as a stable, unchanged entry point.
         """
@@ -1036,82 +774,6 @@ def create_and_add_fold_frame(
             **kwargs,
         )
 
-    def _build_fold_frame(
-        self,
-        fold_frame_name: str,
-        *,
-        index: int | None = None,
-        data=None,
-        interpolatortype="FDI",
-        nelements=LoopStructuralConfig.nelements,
-        tol=None,
-        buffer=0.1,
-        **kwargs,
-    ):
-        """
-        Parameters
-        ----------
-        fold_frame_name : string
-            unique string in feature_name column
-        fold_frame_data : pandas data frame
-            if not specified uses the model data
-        interpolatortype : str
-            the type of interpolator to use, default is 'FDI'
-        nelements : int
-            the number of elements to use in the fold frame
-        tol : float, optional
-            tolerance for the solver
-        buffer : float
-            buffer to add to the bounding box of the fold frame
-        **kwargs : dict
-            additional parameters to be passed to the
-            :class:`LoopStructural.modelling.features.builders.StructuralFrameBuilder`
-            and :meth:`LoopStructural.modelling.features.builders.StructuralFrameBuilder.setup`
-            and the interpolator, such as `domain` or `tol`
-
-
-        Returns
-        -------
-        fold_frame : FoldFrame
-            the created fold frame
-        """
-
-        if tol is None:
-            tol = self.tol
-
-        # create fault frame
-        #
-        fold_frame_builder = StructuralFrameBuilder(
-            interpolatortype=interpolatortype,
-            bounding_box=self.bounding_box.with_buffer(buffer),
-            name=fold_frame_name,
-            frame=FoldFrame,
-            nelements=nelements,
-            model=self,
-            **kwargs,
-        )
-        # add data
-        if data is None:
-            data = self.data.loc[self.data["feature_name"] == fold_frame_name]
-        if data.shape[0] == 0:
-            logger.warning(f"No data for {fold_frame_name}, skipping")
-            return
-        fold_frame_builder.add_data_from_data_frame(
-            self.prepare_data(data, include_feature_name=False)
-        )
-        self._add_faults(fold_frame_builder[0])
-        self._add_faults(fold_frame_builder[1])
-        self._add_faults(fold_frame_builder[2])
-        kwargs["tol"] = tol
-        fold_frame_builder.build(**kwargs)
-        fold_frame = fold_frame_builder.frame
-
-        fold_frame.type = FeatureType.STRUCTURALFRAME
-        fold_frame.builder = fold_frame_builder
-        self._add_feature(fold_frame, index=index)
-
-        return fold_frame
-
     @public_api(tier="stable")
     def create_and_add_folded_foliation(
         self,
@@ -1130,7 +792,8 @@ def create_and_add_folded_foliation(
     ):
         """Create a folded foliation and add it to the model.
 
-        See :meth:`_build_folded_foliation` for parameter documentation.
+        See :meth:`~._model_feature_factory.ModelFeatureFactory.build_folded_foliation`
+        for parameter documentation.
         Thin wrapper around :meth:`create_and_add_feature` (see
         ``API.md``); kept as a stable, unchanged entry point.
         """
@@ -1149,95 +812,6 @@ def create_and_add_folded_foliation(
             **kwargs,
         )
 
-    def _build_folded_foliation(
-        self,
-        foliation_name,
-        *,
-        index: int | None = None,
-        data=None,
-        interpolatortype="DFI",
-        nelements=LoopStructuralConfig.nelements,
-        buffer=0.1,
-        fold_frame=None,
-        svario=True,
-        tol=None,
-        invert_fold_norm=False,
-        **kwargs,
-    ):
-        """
-        Create a folded foliation field from data and a fold frame
-
-        Parameters
-        ----------
-        foliation_data : str
-            unique string in type column of data frame
-        fold_frame :  FoldFrame
-        svario  : Boolean
-            whether to calculate svariograms, saves time if avoided
-        kwargs
-            additional kwargs to be passed through to other functions
-
-        Returns
-        -------
-        feature : GeologicalFeature
-            created geological feature
-
-        Notes
-        -----
-
-        - Building a folded foliation uses the fold interpolation code from Laurent et al., 2016
-        and fold profile fitting from Grose et al., 2017. For more information about the fold modelling
-        see :class:`LoopStructural.modelling.features.fold.FoldEvent`,
-        :class:`LoopStructural.modelling.features.builders.FoldedFeatureBuilder`
-
-        """
-
-        if tol is None:
-            tol = self.tol
-
-        if fold_frame is None:
-            logger.info("Using last feature as fold frame")
-            fold_frame = self.features[-1]
-        if not isinstance(fold_frame, FoldFrame):
-            raise TypeError("Please specify a FoldFrame")
-
-        fold = FoldEvent(fold_frame, name=f"Fold_{foliation_name}", invert_norm=invert_fold_norm)
-
-        if interpolatortype != "DFI":
-            logger.warning("Folded foliation only supports DFI interpolator, changing to DFI")
-            interpolatortype = "DFI"
-        series_builder = FoldedFeatureBuilder(
-            interpolatortype=interpolatortype,
-            bounding_box=self.bounding_box.with_buffer(buffer),
-            nelements=nelements,
-            fold=fold,
-            name=foliation_name,
-            svario=svario,
-            model=self,
-            **kwargs,
-        )
-        if data is None:
-            data = self.data.loc[self.data["feature_name"] == foliation_name]
-        if data.shape[0] == 0:
-            logger.warning(f"No data for {foliation_name}, skipping")
-            return
-        series_builder.add_data_from_data_frame(self.prepare_data(data, include_feature_name=False))
-
-        self._add_faults(series_builder)
-        # series_builder.add_data_to_interpolator(True)
-        # build feature
-
-        kwargs["tol"] = tol
-
-        # series_feature = series_builder.build(**kwargs)
-        series_feature = series_builder.feature
-        series_builder.update_build_arguments(kwargs)
-        series_feature.type = FeatureType.FOLDED
-        series_feature.fold = fold
-
-        self._add_feature(series_feature, index)
-        return series_feature
-
     @public_api(tier="stable")
     def create_and_add_folded_fold_frame(
         self,
@@ -1253,7 +827,8 @@ def create_and_add_folded_fold_frame(
     ):
         """Create a folded fold frame and add it to the model.
 
-        See :meth:`_build_folded_fold_frame` for parameter documentation.
+        See :meth:`~._model_feature_factory.ModelFeatureFactory.build_folded_fold_frame`
+        for parameter documentation.
         Thin wrapper around :meth:`create_and_add_feature` (see
         ``API.md``); kept as a stable, unchanged entry point.
         """
@@ -1269,106 +844,6 @@ def create_and_add_folded_fold_frame(
             **kwargs,
         )
 
-    def _build_folded_fold_frame(
-        self,
-        fold_frame_name: str,
-        *,
-        index: int | None = None,
-        data: pd.DataFrame | None = None,
-        interpolatortype="FDI",
-        nelements=LoopStructuralConfig.nelements,
-        fold_frame=None,
-        tol=None,
-        **kwargs,
-    ):
-        """
-
-        Parameters
-        ----------
-        fold_frame_name : string
-            name of the feature to be added
-        fold_frame_data : pandas data frame, optional
-            data frame containing the fold frame data, if not specified uses the model data
-        interpolatortype : str
-            the type of interpolator to use, default is 'FDI' (unused) 5/6/2025
-        fold_frame : StructuralFrame, optional
-            the fold frame for the fold if not specified uses last feature added
-        nelements : int
-            the number of elements to use in the fold frame
-        tol : float, optional
-            tolerance for the solver, if not specified uses the model default
-        **kwargs : dict
-            additional parameters to be passed to the
-            :class:`LoopStructural.modelling.features.builders.StructuralFrameBuilder`
-            and :meth:`LoopStructural.modelling.features.builders.StructuralFrameBuilder.setup`
-
-        Returns
-        -------
-        fold_frame : FoldFrame
-            created fold frame
-
-        Notes
-        -----
-        This function build a structural frame where the first coordinate is constrained
-        with a fold interpolator.
-        Keyword arguments can be included to constrain
-
-        - :meth:`LoopStructural.GeologicalModel.get_interpolator`
-        - :class:`LoopStructural.StructuralFrameBuilder`
-        - :meth:`LoopStructural.StructuralFrameBuilder.setup`
-         - Building a folded foliation uses the fold interpolation code from Laurent et al., 2016
-        and fold profile fitting from Grose et al., 2017. For more information about the fold modelling
-        see :class:`LoopStructural.modelling.features.fold.FoldEvent`,
-        :class:`LoopStructural.modelling.features.builders.FoldedFeatureBuilder`
-        """
-
-        if tol is None:
-            tol = self.tol
-
-        if fold_frame is None:
-            logger.info("Using last feature as fold frame")
-            fold_frame = self.features[-1]
-        if not isinstance(fold_frame, FoldFrame):
-            raise TypeError("Please specify a FoldFrame")
-        fold = FoldEvent(fold_frame, name=f"Fold_{fold_frame_name}")
-
-        interpolatortypes = [
-            "DFI",
-            "FDI",
-            "FDI",
-        ]
-        fold_frame_builder = StructuralFrameBuilder(
-            interpolatortype=interpolatortypes,
-            bounding_box=self.bounding_box.with_buffer(kwargs.get("buffer", 0.1)),
-            nelements=[nelements, nelements, nelements],
-            name=fold_frame_name,
-            fold=fold,
-            frame=FoldFrame,
-            model=self,
-            **kwargs,
-        )
-        if data is None:
-            data = self.data[self.data["feature_name"] == fold_frame_name]
-        fold_frame_builder.add_data_from_data_frame(
-            self.prepare_data(data, include_feature_name=False)
-        )
-
-        for i in range(3):
-            self._add_faults(fold_frame_builder[i])
-        # build feature
-        kwargs["frame"] = FoldFrame
-        kwargs["tol"] = tol
-        fold_frame_builder.build(**kwargs)
-        # fold_frame_builder.build_arguments = kwargs
-        folded_fold_frame = fold_frame_builder.frame
-        folded_fold_frame.builder = fold_frame_builder
-
-        folded_fold_frame.type = FeatureType.STRUCTURALFRAME
-
-        self._add_feature(folded_fold_frame, index=index)
-
-        return folded_fold_frame
-
     @public_api(tier="stable")
     def create_and_add_intrusion(
         self,
@@ -1383,7 +858,8 @@ def create_and_add_intrusion(
     ):
         """Create an intrusion and add it to the model.
 
-        See :meth:`_build_intrusion` for parameter documentation. Thin
+        See :meth:`~._model_feature_factory.ModelFeatureFactory.build_intrusion`
+        for parameter documentation. Thin
         wrapper around :meth:`create_and_add_feature` (see ``API.md``);
         kept as a stable, unchanged entry point.
         """
@@ -1402,128 +878,6 @@ def create_and_add_intrusion(
             **kwargs,
         )
 
-    def _build_intrusion(
-        self,
-        intrusion_name,
-        intrusion_frame_name,
-        *,
-        intrusion_frame_parameters=None,
-        intrusion_lateral_extent_model=None,
-        intrusion_vertical_extent_model=None,
-        geometric_scaling_parameters=None,
-        **kwargs,
-    ):
-        """
-
-        Note
-        -----
-        An intrusion in built in two main steps:
-        (1) Intrusion builder: intrusion builder creates the intrusion structural frame.
-            This object is curvilinear coordinate system of the intrusion constrained with intrusion network points,
-            and flow and inflation measurements (provided by the user).
-            The intrusion network is a representation of the approximated location of roof or floor contact of the intrusion.
-            This object might be constrained using the anisotropies of the host rock if the roof (or floor) contact is not well constrained.
-
-        (2) Intrusion feature: simulation of lateral and vertical extent of intrusion within the model volume.
-            The simulations outcome consist in thresholds distances along the structural frame coordinates
-            that are used to constrained the extent of the intrusion.
-
-        Parameters
-        ----------
-        intrusion_name :  string,
-            name of intrusion feature in model data
-        intrusion_frame_name :  string,
-            name of intrusion frame in model data
-        intrusion_lateral_extent_model = function,
-            geometrical conceptual model for simulation of lateral extent
-        intrusion_vertical_extent_model = function,
-            geometrical conceptual model for simulation of vertical extent
-        intrusion_frame_parameters = dictionary
-
-        kwargs
-
-        Returns
-        -------
-        intrusion feature
-
-        """
-        if intrusion_frame_parameters is None:
-            intrusion_frame_parameters = {}
-        if geometric_scaling_parameters is None:
-            geometric_scaling_parameters = {}
-        # if intrusions is False:
-        #     logger.error("Libraries not installed")
-        #     raise Exception("Libraries not installed")
-
-        intrusion_data = self.data[self.data["feature_name"] == intrusion_name].copy()
-        intrusion_frame_data = self.data[self.data["feature_name"] == intrusion_frame_name].copy()
-
-        # -- get variables for intrusion frame interpolation
-        gxxgz = kwargs.get("gxxgz", 0)
-        gxxgy = kwargs.get("gxxgy", 0)
-        gyxgz = kwargs.get("gyxgz", 0)
-
-        interpolatortype = kwargs.get("interpolatortype", "PLI")
-        # buffer = kwargs.get("buffer", 0.1)
-        nelements = kwargs.get("nelements", LoopStructuralConfig.nelements)
-
-        weights = [gxxgz, gxxgy, gyxgz]
-
-        intrusion_frame_builder = IntrusionFrameBuilder(
-            interpolatortype=interpolatortype,
-            bounding_box=self.bounding_box.with_buffer(kwargs.get("buffer", 0.1)),
-            nelements=kwargs.get("nelements", LoopStructuralConfig.nelements),
-            name=intrusion_frame_name,
-            model=self,
-            **kwargs,
-        )
-
-        self._add_faults(intrusion_frame_builder)
-        # intrusion_frame_builder.post_intrusion_faults = faults  # LG unused?
-
-        # -- create intrusion frame using intrusion structures (steps and marginal faults) and flow/inflation measurements
-        if len(intrusion_frame_parameters) == 0:
-            logger.error("Please specify parameters to build intrusion frame")
-        intrusion_frame_builder.set_intrusion_frame_parameters(
-            intrusion_data, intrusion_frame_parameters
-        )
-        intrusion_frame_builder.create_constraints_for_c0()
-
-        intrusion_frame_builder.set_intrusion_frame_data(intrusion_frame_data)
-
-        ## -- create intrusion frame
-        intrusion_frame_builder.build(
-            nelements=nelements,
-            w2=weights[0],
-            w1=weights[1],
-            gxygz=weights[2],
-        )
-
-        intrusion_frame = intrusion_frame_builder.frame
-
-        # -- create intrusion builder to compute distance thresholds along the frame coordinates
-        intrusion_builder = IntrusionBuilder(
-            intrusion_frame,
-            model=self,
-            # interpolator=interpolator,
-            name=f"{intrusion_name}_feature",
-            lateral_extent_model=intrusion_lateral_extent_model,
-            vertical_extent_model=intrusion_vertical_extent_model,
-            **kwargs,
-        )
-        intrusion_builder.set_data_for_extent_calculation(intrusion_data)
-
-        intrusion_builder.update_build_arguments(
-            {
-                "geometric_scaling_parameters": geometric_scaling_parameters,
-            }
-        )
-
-        intrusion_feature = intrusion_builder.feature
-        self._add_feature(intrusion_feature)
-
-        return intrusion_feature
-
     def _add_faults(self, feature_builder, features=None):
         """Adds all existing faults to a geological feature builder
 
@@ -1537,34 +891,7 @@ def _add_faults(self, feature_builder, features=None):
         -------
 
         """
-        if features is None:
-            features = self.features
-        for f in reversed(features):
-            if isinstance(f, str):
-                f = self.__getitem__(f)
-            if f.type == FeatureType.FAULT:
-                feature_builder.add_fault(f)
-
-    def _add_domain_fault_above(self, feature):
-        """
-        Looks through the feature list and adds any domain faults to the feature. The domain fault masks everything
-        where the fault scalar field is < 0 as being active when added to feature.
-
-        Parameters
-        ----------
-        feature : GeologicalFeatureBuilder
-            the feature being added to the model where domain faults should be added
-
-        Returns
-        -------
-
-        """
-        for f in reversed(self.features):
-            if f.name == feature.name:
-                continue
-            if f.type == "domain_fault":
-                feature.add_region(lambda pos, fault=f: fault.evaluate_value(pos) < 0)
-                break
+        FeatureRelationshipManager.add_faults(self, feature_builder, features=features)
 
     def _add_domain_fault_below(self, domain_fault):
         """
@@ -1582,40 +909,7 @@ def _add_domain_fault_below(self, domain_fault):
         -------
 
         """
-        for f in reversed(self.features):
-            if f.name == domain_fault.name:
-                continue
-            f.add_region(lambda pos: domain_fault.evaluate_value(pos) > 0)
-            if f.type == FeatureType.UNCONFORMITY:
-                break
-
-    def _add_unconformity_above(self, feature):
-        """
-
-        Adds a region to the feature to prevent the value from being
-        interpolated where the unconformities exists above e.g.
-        if there is another feature above and the unconformity is at 0
-        then the features added below (after) will only be visible where the
-        uncomformity is <0
-
-        Parameters
-        ----------
-        feature - GeologicalFeature
-
-        Returns
-        -------
-
-        """
-
-        if feature.type == FeatureType.FAULT:
-            return
-        for f in reversed(self.features):
-            if f.type == FeatureType.UNCONFORMITY and f.name != feature.name:
-                logger.info(f"Adding {f.name} as unconformity to {feature.name}")
-                feature.add_region(f)
-            if f.type == FeatureType.ONLAPUNCONFORMITY and f.name != feature.name:
-                feature.add_region(f)
-                break
+        FeatureRelationshipManager.add_domain_fault_below(self, domain_fault)
 
     @public_api(tier="stable")
     def add_unconformity(
@@ -1637,28 +931,7 @@ def add_unconformity(
             unconformity feature
 
         """
-        logger.debug(f"Adding {feature.name} as unconformity at {value}")
-        if feature is None:
-            logger.warning("Cannot add unconformtiy, base feature is None")
-            return
-        # look backwards through features and add the unconformity as a region until
-        # we get to an unconformity
-        uc_feature = UnconformityFeature(feature, value)
-        feature.add_region(uc_feature.inverse())
-        for f in reversed(self.features):
-            if f.type == FeatureType.UNCONFORMITY:
-                logger.debug(f"Reached unconformity {f.name}")
-                break
-            logger.debug(f"Adding {uc_feature.name} as unconformity to {f.name}")
-            if f.type == FeatureType.FAULT or f.type == FeatureType.INACTIVEFAULT:
-                continue
-            if f == feature:
-                continue
-            else:
-                f.add_region(uc_feature)
-        # now add the unconformity to the feature list
-        self._add_feature(uc_feature, index=index)
-        return uc_feature
+        return FeatureRelationshipManager.add_unconformity(self, feature, value, index=index)
 
     @public_api(tier="stable")
     def add_onlap_unconformity(
@@ -1680,20 +953,7 @@ def add_onlap_unconformity(
             the created unconformity
 
         """
-        feature.regions = []
-        uc_feature = UnconformityFeature(feature, value, False, onlap=True)
-        feature.add_region(uc_feature.inverse())
-        for f in reversed(self.features):
-            if f.type in (FeatureType.UNCONFORMITY, FeatureType.ONLAPUNCONFORMITY):
-                logger.debug(f"Reached unconformity {f.name}")
-                break
-            if f.type == FeatureType.FAULT or f.type == FeatureType.INACTIVEFAULT:
-                continue
-            if f != feature:
-                f.add_region(uc_feature)
-        self._add_feature(uc_feature.inverse(), index=index)
-
-        return uc_feature
+        return FeatureRelationshipManager.add_onlap_unconformity(self, feature, value, index=index)
 
     @public_api(tier="provisional")
     def add_fold_to_feature(
@@ -1762,7 +1022,8 @@ def create_and_add_domain_fault(
     ):
         """Create a domain fault and add it to the model.
 
-        See :meth:`_build_domain_fault` for parameter documentation. Thin
+        See :meth:`~._model_feature_factory.ModelFeatureFactory.build_domain_fault`
+        for parameter documentation. Thin
         wrapper around :meth:`create_and_add_feature` (see ``API.md``);
         kept as a stable, unchanged entry point.
         """
@@ -1775,62 +1036,6 @@ def create_and_add_domain_fault(
             **kwargs,
         )
 
-    def _build_domain_fault(
-        self,
-        fault_surface_data,
-        *,
-        nelements=LoopStructuralConfig.nelements,
-        interpolatortype="FDI",
-        index: int | None = None,
-        **kwargs,
-    ):
-        """
-        Parameters
-        ----------
-        fault_surface_data : string
-            name of the domain fault data in the data frame
-
-        Returns
-        -------
-        domain_Fault : GeologicalFeature
-            the created domain fault
-
-        Notes
-        -----
-        * :meth:`LoopStructural.GeologicalModel.get_interpolator`
-
-        """
-        domain_fault_feature_builder = GeologicalFeatureBuilder(
-            bounding_box=self.bounding_box,
-            interpolatortype=interpolatortype,
-            nelements=nelements,
-            name=fault_surface_data,
-            model=self,
-            **kwargs,
-        )
-
-        # add data
-        unconformity_data = self.data.loc[self.data["feature_name"] == fault_surface_data]
-
-        domain_fault_feature_builder.add_data_from_data_frame(unconformity_data)
-        # look through existing features if there is a fault before an
-        # unconformity
-        # then add to the feature, once we get to an unconformity stop
-        self._add_faults(domain_fault_feature_builder)
-
-        # build feature
-        # domain_fault = domain_fault_feature_builder.build(**kwargs)
-        domain_fault = domain_fault_feature_builder.feature
-        domain_fault_feature_builder.update_build_arguments(kwargs)
-        domain_fault.type = FeatureType.DOMAINFAULT
-        self._add_feature(domain_fault, index=index)
-        self._add_domain_fault_below(domain_fault)
-
-        domain_fault_uc = UnconformityFeature(domain_fault, 0)
-        # iterate over existing features and add the unconformity as a region
-        # so the feature is only evaluated where the unconformity is positive
-        return domain_fault_uc
-
     @public_api(tier="stable")
     def create_and_add_fault(
         self,
@@ -1860,7 +1065,8 @@ def create_and_add_fault(
     ):
         """Create a fault and add it to the model.
 
-        See :meth:`_build_fault` for parameter documentation. Thin
+        See :meth:`~._model_feature_factory.ModelFeatureFactory.build_fault`
+        for parameter documentation. Thin
         wrapper around :meth:`create_and_add_feature` (see ``API.md``);
         kept as a stable, unchanged entry point.
         """
@@ -1892,158 +1098,6 @@ def create_and_add_fault(
             **kwargs,
         )
 
-    def _build_fault(
-        self,
-        fault_name: str,
-        displacement: float,
-        *,
-        index: int | None = None,
-        data: pd.DataFrame | None = None,
-        interpolatortype="FDI",
-        tol=None,
-        fault_slip_vector=None,
-        fault_normal_vector=None,
-        fault_center=None,
-        major_axis=None,
-        minor_axis=None,
-        intermediate_axis=None,
-        faultfunction="BaseFault",
-        faults=None,
-        force_mesh_geometry: bool = False,
-        points: bool = False,
-        fault_buffer=0.2,
-        fault_trace_anisotropy=0.0,
-        fault_dip=90,
-        fault_dip_anisotropy=0.0,
-        fault_pitch=None,
-        **kwargs,
-    ):
-        """
-        Parameters
-        ----------
-        fault_name : string
-            name of the fault surface data in the dataframe
-        displacement : displacement magnitude
-            displacement magnitude of the fault, in model units
-        fault_data : pd.DataFrame, optional
-            data frame containing the fault data, if not specified uses the model data
-        major_axis : [type], optional
-            [description], by default None
-        minor_axis : [type], optional
-            [description], by default None
-        intermediate_axis : [type], optional
-            [description], by default None
-        kwargs : additional kwargs for Fault and interpolators
-
-        Returns
-        -------
-        fault : FaultSegment
-            created fault
-
-        Notes
-        -----
-        * :meth:`LoopStructural.GeologicalModel.get_interpolator`
-        * :class:`LoopStructural.modelling.features.builders.FaultBuilder`
-        * :meth:`LoopStructural.modelling.features.builders.FaultBuilder.setup`
-        """
-        if faults is None:
-            faults = []
-        if "fault_extent" in kwargs and major_axis is None:
-            major_axis = kwargs["fault_extent"]
-        if "fault_influence" in kwargs and minor_axis is None:
-            minor_axis = kwargs["fault_influence"]
-        if "fault_vectical_radius" in kwargs and intermediate_axis is None:
-            intermediate_axis = kwargs["fault_vectical_radius"]
-
-        logger.info(f'Creating fault "{fault_name}"')
-        logger.info(f"Displacement: {displacement}")
-        logger.info(f"Tolerance: {tol}")
-        logger.info(f"Fault function: {faultfunction}")
-        logger.info(f"Fault slip vector: {fault_slip_vector}")
-        logger.info(f"Fault center: {fault_center}")
-        logger.info(f"Major axis: {major_axis}")
-        logger.info(f"Minor axis: {minor_axis}")
-        logger.info(f"Intermediate axis: {intermediate_axis}")
-        if fault_slip_vector is not None:
-            fault_slip_vector = np.array(fault_slip_vector, dtype="float")
-        if fault_center is not None:
-            fault_center = np.array(fault_center, dtype="float")
-
-        for k, v in kwargs.items():
-            logger.info(f"{k}: {v}")
-
-        if tol is None:
-            tol = self.tol
-            # divide the tolerance by half of the minor axis, as this is the equivalent of the distance
-            # of the unit vector
-            # if minor_axis:
-            # tol *= 0.1*minor_axis
-
-        if displacement == 0:
-            logger.warning(f"{fault_name} displacement is 0")
-
-        if "data_region" in kwargs:
-            kwargs.pop("data_region")
-            logger.error("kwarg data_region currently not supported, disabling")
-        displacement_scaled = displacement
-        fault_frame_builder = FaultBuilder(
-            interpolatortype,
-            bounding_box=self.bounding_box,
-            nelements=kwargs.pop("nelements", LoopStructuralConfig.nelements),
-            name=fault_name,
-            model=self,
-            **kwargs,
-        )
-        if data is None:
-            data = self.data.loc[self.data["feature_name"] == fault_name]
-        if data.shape[0] == 0:
-            logger.warning(f"No data for {fault_name}, skipping")
-            return
-
-        self._add_faults(fault_frame_builder, features=faults)
-        # add data
-
-        if fault_center is not None and ~np.isnan(fault_center).any():
-            fault_center = self.scale(fault_center, inplace=False)
-        # Keep the supplied fault-axis values unchanged; the previous self-assignment
-        # was only present to satisfy a linter and did not affect behavior.
-        fault_frame_builder.create_data_from_geometry(
-            fault_frame_data=self.prepare_data(data, include_feature_name=False),
-            fault_center=fault_center,
-            fault_normal_vector=fault_normal_vector,
-            fault_slip_vector=fault_slip_vector,
-            minor_axis=minor_axis,
-            major_axis=major_axis,
-            intermediate_axis=intermediate_axis,
-            points=points,
-            force_mesh_geometry=force_mesh_geometry,
-            fault_buffer=fault_buffer,
-            fault_trace_anisotropy=fault_trace_anisotropy,
-            fault_dip=fault_dip,
-            fault_dip_anisotropy=fault_dip_anisotropy,
-            fault_pitch=fault_pitch,
-        )
-        if "force_mesh_geometry" not in kwargs:
-            fault_frame_builder.set_mesh_geometry(kwargs.get("fault_buffer", 0.2), 0)
-        if "splay" in kwargs and "splayregion" in kwargs:
-            fault_frame_builder.add_splay(kwargs["splay"], kwargs["splayregion"])
-
-        kwargs["tol"] = tol
-        fault_frame_builder.build(**kwargs)
-        fault = fault_frame_builder.frame
-        fault.displacement = displacement_scaled
-        fault.faultfunction = faultfunction
-
-        for f in reversed(self.features):
-            if f.type == FeatureType.UNCONFORMITY:
-                fault.add_region(f)
-                break
-        if displacement == 0:
-            fault.type = FeatureType.INACTIVEFAULT
-        self._add_feature(fault, index=index)
-
-        return fault
-
     # TODO move rescale to bounding box/transformer
     @public_api(tier="stable")
     def rescale(self, points: np.ndarray, *, inplace: bool = False) -> np.ndarray:
@@ -2383,62 +1437,15 @@ def stratigraphic_ids(self):
 
     @public_api(tier="stable")
     def get_fault_surfaces(self, faults: list[str] | None = None):
-        if faults is None:
-            faults = []
-        surfaces = []
-        if len(faults) == 0:
-            faults = self.fault_names()
-
-        for f in faults:
-            surfaces.extend(self.get_feature_by_name(f).surfaces([0], self.bounding_box))
-        return surfaces
+        return ModelExporter.get_fault_surfaces(self, faults=faults)
 
     @public_api(tier="stable")
     def get_stratigraphic_surfaces(self, units: list[str] | None = None, bottoms: bool = True):
-        if units is None:
-            units = []
-        ## TODO change the stratigraphic column to its own class and have methods to get the relevant surfaces
-        surfaces = []
-        units = []
-        if self.stratigraphic_column is None:
-            return []
-        units = self.stratigraphic_column.get_isovalues()
-        units_for_group = {}
-        for name, u in units.items():
-            if u['group'] not in self:
-                logger.warning(f"Group {u['group']} not found in model")
-                continue
-            if u['group'] not in units_for_group:
-                units_for_group[u['group']] = []
-            u['name'] = name
-            units_for_group[u['group']].append(u)
-        for group, us in units_for_group.items():
-            feature = self.get_feature_by_name(group)
-            values = [u['value'] for u in us]
-            colours = [u['colour'] for u in us]
-            names = [u['name'] for u in us]
-            surfaces.extend(
-                feature.surfaces(values, self.bounding_box, name=names, colours=colours)
-            )
-
-        return surfaces
+        return ModelExporter.get_stratigraphic_surfaces(self, units=units, bottoms=bottoms)
 
     @public_api(tier="stable")
     def get_block_model(self, name='block model'):
-        # NOTE: bounding_box.structured_grid() returns loop_common's
-        # interpolation-support StructuredGrid (no properties dict); use
-        # LoopStructural's own geometry StructuredGrid for storing values.
-        grid = StructuredGrid(
-            origin=self.bounding_box.origin,
-            step_vector=self.bounding_box.step_vector,
-            nsteps=self.bounding_box.nsteps,
-            name=name,
-        )
-
-        grid.cell_properties['stratigraphy'] = self.evaluate_model(
-            self.rescale(self.bounding_box.cell_centres())
-        )
-        return grid, self.stratigraphic_ids()
+        return ModelExporter.get_block_model(self, name=name)
 
     @public_api(tier="stable")
     def save(
@@ -2450,78 +1457,49 @@ def save(
         stratigraphic_data=True,
         fault_data=True,
     ):
-        path = pathlib.Path(filename)
-        extension = path.suffix
-        parent = path.parent
-        name = path.stem
-        stratigraphic_surfaces = self.get_stratigraphic_surfaces()
-        if fault_surfaces:
-            for s in self.get_fault_surfaces():
-                ## geoh5 can save everything into the same file
-                if extension == ".geoh5" or extension == '.omf':
-                    s.save(filename)
-                else:
-                    s.save(f'{parent}/{name}_{s.name}{extension}')
-        if stratigraphic_surfaces:
-            for s in self.get_stratigraphic_surfaces():
-                if extension == ".geoh5" or extension == '.omf':
-                    s.save(filename)
-                else:
-                    s.save(f'{parent}/{name}_{s.name}{extension}')
-        if block_model:
-            grid, _ids = self.get_block_model()
-            if extension == ".geoh5" or extension == '.omf':
-                grid.save(filename)
-            else:
-                grid.save(f'{parent}/{name}_block_model{extension}')
-        if stratigraphic_data and self.stratigraphic_column is not None:
-            for group in self.stratigraphic_column:
-                if group == "faults":
-                    continue
-                for data in self.__getitem__(group).get_data():
-                    if extension == ".geoh5" or extension == '.omf':
-                        data.save(filename)
-                    else:
-                        data.save(f'{parent}/{name}_{group}_data{extension}')
-        if fault_data:
-            for f in self.fault_names():
-                for d in self.__getitem__(f).get_data():
-                    if extension == ".geoh5" or extension == '.omf':
-
-                        d.save(filename)
-                    else:
-                        d.save(f'{parent}/{name}_{group}{extension}')
+        ModelExporter.save(
+            self,
+            filename,
+            block_model=block_model,
+            stratigraphic_surfaces=stratigraphic_surfaces,
+            fault_surfaces=fault_surfaces,
+            stratigraphic_data=stratigraphic_data,
+            fault_data=fault_data,
+        )
 
 
 # Wire the built-in feature types up to GeologicalModel.create_and_add_feature
-# (see FeatureBuilderRegistry / API.md). Each factory reuses the existing
-# _build_* method unchanged; new feature types register here without
-# modifying GeologicalModel's source.
+# (see FeatureBuilderRegistry / API.md). Each factory calls the corresponding
+# ModelFeatureFactory.build_* staticmethod; new feature types register here
+# without modifying GeologicalModel's source.
 FeatureBuilderRegistry.register(
-    "foliation", lambda model, name, **params: model._build_foliation(name, **params)
+    "foliation", lambda model, name, **params: ModelFeatureFactory.build_foliation(model, name, **params)
 )
 FeatureBuilderRegistry.register(
-    "fold_frame", lambda model, name, **params: model._build_fold_frame(name, **params)
+    "fold_frame",
+    lambda model, name, **params: ModelFeatureFactory.build_fold_frame(model, name, **params),
 )
 FeatureBuilderRegistry.register(
     "folded_foliation",
-    lambda model, name, **params: model._build_folded_foliation(name, **params),
+    lambda model, name, **params: ModelFeatureFactory.build_folded_foliation(model, name, **params),
 )
 FeatureBuilderRegistry.register(
     "folded_fold_frame",
-    lambda model, name, **params: model._build_folded_fold_frame(name, **params),
+    lambda model, name, **params: ModelFeatureFactory.build_folded_fold_frame(model, name, **params),
 )
 FeatureBuilderRegistry.register(
     "intrusion",
-    lambda model, name, **params: model._build_intrusion(
-        name, params.pop("intrusion_frame_name"), **params
+    lambda model, name, **params: ModelFeatureFactory.build_intrusion(
+        model, name, params.pop("intrusion_frame_name"), **params
     ),
 )
 FeatureBuilderRegistry.register(
     "domain_fault",
-    lambda model, name, **params: model._build_domain_fault(name, **params),
+    lambda model, name, **params: ModelFeatureFactory.build_domain_fault(model, name, **params),
 )
 FeatureBuilderRegistry.register(
     "fault",
-    lambda model, name, **params: model._build_fault(name, params.pop("displacement"), **params),
+    lambda model, name, **params: ModelFeatureFactory.build_fault(
+        model, name, params.pop("displacement"), **params
+    ),
 )
diff --git a/LoopStructural/modelling/features/_geological_feature.py b/LoopStructural/modelling/features/_geological_feature.py
index 7c47ac0c..4660591b 100644
--- a/LoopStructural/modelling/features/_geological_feature.py
+++ b/LoopStructural/modelling/features/_geological_feature.py
@@ -202,6 +202,7 @@ def evaluate_gradient(
             tetrahedron = regular_tetraherdron_for_points(pos, element_scale_parameter)
 
             while not resolved:
+                resolved = True
                 for f in self.faults:
                     v = (
                         f[0]
@@ -215,8 +216,7 @@ def evaluate_gradient(
                         )
                         element_scale_parameter *= 0.5
                         tetrahedron = regular_tetraherdron_for_points(pos, element_scale_parameter)
-
-                resolved = True
+                        resolved = False
 
             tetrahedron_faulted = self._apply_faults(np.array(tetrahedron.reshape(-1, 3))).reshape(
                 tetrahedron.shape
diff --git a/LoopStructural/modelling/intrusions/intrusion_builder.py b/LoopStructural/modelling/intrusions/intrusion_builder.py
index 703cb929..c4d80609 100644
--- a/LoopStructural/modelling/intrusions/intrusion_builder.py
+++ b/LoopStructural/modelling/intrusions/intrusion_builder.py
@@ -3,7 +3,6 @@
 
 from ...utils import getLogger, rng
 from ..features.builders import BaseBuilder
-from .geometric_scaling_functions import *
 from .intrusion_feature import IntrusionFeature
 
 logger = getLogger(__name__)
@@ -113,45 +112,25 @@ def set_data_for_extent_calculation(self, intrusion_data: pd.DataFrame):
     def create_geometry_using_geometric_scaling(
         self, geometric_scaling_parameters, reference_contact_data
     ):
-
-        geometric_scaling_parameters.get("intrusion_type", None)
-        intrusion_length = geometric_scaling_parameters.get("intrusion_length", None)
-        geometric_scaling_parameters.get("inflation_vector", np.array([[0, 0, 1]]))
-        thickness = geometric_scaling_parameters.get("thickness", None)
-
-        if (
-            self.intrusion_frame.builder.intrusion_network_contact == "floor"
-            or self.intrusion_frame.builder.intrusion_network_contact == "base"
-        ):
-            geometric_scaling_parameters.get("inflation_vector", np.array([[0, 0, 1]]))
-        else:
-            geometric_scaling_parameters.get("inflation_vector", np.array([[0, 0, -1]]))
-
-        if intrusion_length is None and thickness is None:
-            raise ValueError(
-                f"No {self.intrusion_frame.builder.intrusion_other_contact} data. Add intrusion_type and intrusion_length (or thickness) to geometric_scaling_parameters dictionary"
-            )
-
-        else:  # -- create synthetic data to constrain interpolation using geometric scaling
-            estimated_thickness = thickness
-            if estimated_thickness is None:
-                raise NotImplementedError("Not implemented")
-                # estimated_thickness = thickness_from_geometric_scaling(
-                #     intrusion_length, intrusion_type
-                # )
-
-            logger.info(
-                f"Building tabular intrusion using geometric scaling parameters: estimated thicknes = {round(estimated_thickness)} meters"
-            )
-            raise NotImplementedError("Not implemented")
-            # (
-            #     other_contact_data_temp,
-            #     other_contact_data_xyz_temp,
-            # ) = contact_pts_using_geometric_scaling(
-            #     estimated_thickness, reference_contact_data, inflation_vector
-            # )
-
-            # return other_contact_data_temp
+        """Not currently implemented.
+
+        This is meant to synthesise the missing contact (roof or floor) from
+        an estimated thickness (either given directly or derived from
+        empirical length/thickness scaling laws, see
+        ``geometric_scaling_functions.thickness_from_geometric_scaling``) and
+        an inflation vector, via
+        ``geometric_scaling_functions.contact_pts_using_geometric_scaling``.
+        That wiring was never completed, so every call path here always
+        raised ``NotImplementedError`` regardless of what was passed in
+        (see ``INTRUSIONS.md`` finding 2). Raising immediately, rather than
+        after partially validating parameters, makes that unambiguous.
+        """
+        raise NotImplementedError(
+            "geometric_scaling_parameters is not currently supported: "
+            f"'{self.intrusion_frame.builder.intrusion_other_contact}' contact "
+            "has no data, and synthesising it from geometric scaling is not "
+            "implemented. Provide explicit data for both contacts instead."
+        )
 
     def prepare_data(self, geometric_scaling_parameters):
         """Prepare the data to compute distance thresholds along the frame coordinates.
diff --git a/LoopStructural/modelling/intrusions/intrusion_feature.py b/LoopStructural/modelling/intrusions/intrusion_feature.py
index bc76c74c..d4f7ab0c 100644
--- a/LoopStructural/modelling/intrusions/intrusion_feature.py
+++ b/LoopStructural/modelling/intrusions/intrusion_feature.py
@@ -265,13 +265,8 @@ def evaluate_value(self, pos):
             intrusion_coord1_pts
         )
 
-        if self.intrusion_frame.builder.marginal_faults is not None:
-            c2_minside_threshold = thresholds[0]  # np.zeros_like(intrusion_coord2_pts)
-            c2_maxside_threshold = thresholds[1]
-
-        else:
-            c2_minside_threshold = thresholds[0]
-            c2_maxside_threshold = thresholds[1]
+        c2_minside_threshold = thresholds[0]
+        c2_maxside_threshold = thresholds[1]
 
         thresholds, _residuals, _conceptual = self.interpolate_vertical_thresholds(
             intrusion_coord1_pts, intrusion_coord2_pts
@@ -332,81 +327,6 @@ def evaluate_value(self, pos):
 
         return intrusion_sf
 
-    def evaluate_value_test(self, points):
-        """
-        Computes a distance scalar field to the intrusion contact (isovalue = 0).
-
-        Parameters
-        ------------
-        points : numpy array (x,y,z),  points where the IntrusionFeature is evaluated.
-
-        Returns
-        ------------
-        intrusion_sf : numpy array, contains distance to intrusion contact
-
-        """
-        self.builder.up_to_date()
-
-        # compute coordinates values for each evaluated point
-        intrusion_coord0_pts = self.intrusion_frame[0].evaluate_value(points)
-        intrusion_coord1_pts = self.intrusion_frame[1].evaluate_value(points)
-        intrusion_coord2_pts = self.intrusion_frame[2].evaluate_value(points)
-
-        self.evaluated_points = [
-            points,
-            intrusion_coord0_pts,
-            intrusion_coord1_pts,
-            intrusion_coord2_pts,
-        ]
-
-        thresholds, _residuals, _conceptual = self.interpolate_lateral_thresholds(
-            intrusion_coord1_pts
-        )
-
-        if self.intrusion_frame.builder.marginal_faults is not None:
-            c2_minside_threshold = np.zeros_like(intrusion_coord2_pts)
-            c2_maxside_threshold = thresholds[1]
-
-        else:
-            c2_minside_threshold = thresholds[0]
-            c2_maxside_threshold = thresholds[1]
-
-        thresholds, _residuals, _conceptual = self.interpolate_vertical_thresholds(
-            intrusion_coord1_pts, intrusion_coord2_pts
-        )
-        c0_minside_threshold = thresholds[1]
-        c0_maxside_threshold = thresholds[0]
-
-        mid_point = c0_minside_threshold + ((c0_maxside_threshold - c0_minside_threshold) / 2)
-
-        mod_intrusion_coord0_pts = intrusion_coord0_pts - mid_point
-        mod_c0_minside_threshold = c0_minside_threshold - mid_point
-        mod_c0_maxside_threshold = c0_maxside_threshold + mid_point
-
-        a = (
-            (mod_intrusion_coord0_pts >= mid_point)
-            * (c2_minside_threshold < intrusion_coord2_pts)
-            * (intrusion_coord2_pts < c2_maxside_threshold)
-        )
-        b = (
-            (mod_intrusion_coord0_pts <= mid_point)
-            * (c2_minside_threshold < intrusion_coord2_pts)
-            * (intrusion_coord2_pts < c2_maxside_threshold)
-        )
-        c = (
-            (mod_intrusion_coord0_pts <= mid_point)
-            * (mod_intrusion_coord0_pts >= mod_c0_minside_threshold)
-            * (c2_minside_threshold < intrusion_coord2_pts)
-            * (intrusion_coord2_pts < c2_maxside_threshold)
-        )
-
-        intrusion_sf = mod_intrusion_coord0_pts
-        intrusion_sf[a] = mod_intrusion_coord0_pts[a] - mod_c0_maxside_threshold[a]
-        intrusion_sf[b] = abs(mod_c0_minside_threshold[b] + mod_intrusion_coord0_pts[b])
-        intrusion_sf[c] = mod_intrusion_coord0_pts[c] - mod_c0_minside_threshold[c]
-
-        return intrusion_sf
-
     def get_data(self, value_map: dict | None = None):
         pass
 
diff --git a/LoopStructural/modelling/intrusions/intrusion_frame_builder.py b/LoopStructural/modelling/intrusions/intrusion_frame_builder.py
index 1cd7f7e4..9e100d6a 100644
--- a/LoopStructural/modelling/intrusions/intrusion_frame_builder.py
+++ b/LoopStructural/modelling/intrusions/intrusion_frame_builder.py
@@ -17,6 +17,12 @@
     logger.error('Scikitlearn cannot be imported')
     raise
 
+# Fixed (not derived from the shared `rng`) so that repeated builds of the
+# same intrusion produce the same contact/fault clustering: `loop_common`'s
+# shared `rng` is a fresh, unseeded `np.random.default_rng()` per process, so
+# routing this through it would make cluster labels vary run-to-run instead.
+_KMEANS_RANDOM_STATE = 0
+
 
 class IntrusionFrameBuilder(StructuralFrameBuilder):
     def __init__(
@@ -228,10 +234,9 @@ def add_contact_anisotropies(self, series_list: list | None = None, **kwargs):
 
             # -- use scalar field values to find different contacts
             series_i_vals_mod = series_i_vals.reshape(len(series_i_vals), 1)
-            # TODO create global loopstructural random state variable
-            contact_clustering = KMeans(n_clusters=n_contacts, random_state=0).fit(
-                series_i_vals_mod
-            )
+            contact_clustering = KMeans(
+                n_clusters=n_contacts, random_state=_KMEANS_RANDOM_STATE
+            ).fit(series_i_vals_mod)
 
             for j in range(n_contacts):
                 z = np.ma.masked_not_equal(contact_clustering.labels_, j)
@@ -358,7 +363,9 @@ def set_intrusion_steps_parameters(self):
                 )
                 series_values = series_from_name.evaluate_value(data_points_xyz)
                 series_values_mod = series_values.reshape(len(series_values), 1)
-                contact_clustering = KMeans(n_clusters=2, random_state=0).fit(series_values_mod)
+                contact_clustering = KMeans(
+                    n_clusters=2, random_state=_KMEANS_RANDOM_STATE
+                ).fit(series_values_mod)
 
                 # contact 0
                 z = np.ma.masked_not_equal(contact_clustering.labels_, 0)
@@ -491,7 +498,6 @@ def set_marginal_faults_parameters(self):
         for fault_i in self.marginal_faults:
             marginal_fault = self.marginal_faults[fault_i].get("structure")
             block = self.marginal_faults[fault_i].get("block")  # hanging wall or foot wall
-            self.marginal_faults[fault_i].get("emplacement_mechanism")
             series_name = self.marginal_faults[fault_i].get("series")
 
             series_values_temp = series_name.evaluate_value(intrusion_frame_c0_data_xyz)
@@ -828,7 +834,6 @@ def create_constraints_for_c0(self, **kwargs):
                 delta_contact = self.marginal_faults[fault_i].get("delta_c", 1)
                 marginal_fault = self.marginal_faults[fault_i].get("structure")
                 block = self.marginal_faults[fault_i].get("block")  # hanging wall or foot wall
-                self.marginal_faults[fault_i].get("emplacement_mechanism")
                 series_name = self.marginal_faults[fault_i].get("series")
 
                 fault_gridpoints_vals = marginal_fault[0].evaluate_value(grid_points)
@@ -976,7 +981,3 @@ def set_intrusion_frame_data(self, intrusion_frame_data):  # , intrusion_network
 
         self.add_data_from_data_frame(intrusion_frame_data_complete)
         self.update_geometry(intrusion_frame_data_complete[["X", "Y", "Z"]].to_numpy())
-
-    def update(self):
-        for i in range(3):
-            self.builders[i].update()
diff --git a/LoopStructural/modelling/intrusions/intrusion_support_functions.py b/LoopStructural/modelling/intrusions/intrusion_support_functions.py
deleted file mode 100644
index 086ef247..00000000
--- a/LoopStructural/modelling/intrusions/intrusion_support_functions.py
+++ /dev/null
@@ -1,389 +0,0 @@
-## Support Functions for intrusion network simulated as the shortest path, and for simulations in general
-import numpy as np
-
-from ...utils import getLogger
-
-logger = getLogger(__name__)
-
-
-def sort_2_arrays(main_array, array):
-    """
-    Sort two arrays, considering values of only the main array
-
-    Parameters
-    ----------
-    main array: numpy array, considered to sort secondary array
-    array: numpy aray, array to be sorted
-
-    Returns
-    -------
-    sorted arrays
-    """
-    # function to sort 2 arrays, considering values of only the main array
-
-    for i in range(len(main_array)):
-        swap = i + np.argmin(main_array[i:])
-        (main_array[i], main_array[swap]) = (main_array[swap], main_array[i])
-        (array[i], array[swap]) = (array[swap], array[i])
-
-    return main_array, array
-
-
-def findMinDiff(arr, n):
-    """
-    Find the min diff by comparing difference of all possible pairs in given array
-
-    Parameters
-    ----------
-    arr: numpy array with values to compare and fin the minimum difference
-
-    Returns
-    -------
-    minimum difference between values in arr
-
-    """
-    # Initialize difference as infinite
-    diff = 10**20
-
-    if n < 2:
-        return diff
-
-    values = np.asarray(arr[:n], dtype=float)
-    pairwise_diff = np.abs(values[:, None] - values[None, :])
-    np.fill_diagonal(pairwise_diff, np.inf)
-    min_diff = pairwise_diff.min()
-    diff = min(diff, min_diff)
-
-    return diff
-
-
-def array_from_coords(df, section_axis, df_axis):
-    """
-    Create numpy array representing a section of the model
-    from a dataframe containing coordinates and values
-
-    Parameters
-    ----------
-    df: pandas dataframe, should have at least ['X', 'Y', 'Z', 'val_1'] columns
-    section_axis: string 'X' or 'Y', the cross section represented by the arrays is along the section_axis
-    df_axis: number of the column where the value on interest is
-
-    Returns
-    -------
-    array: numpy array, contains values (e.g. velocities at each point, scalar field value at each point, etc)
-
-    """
-
-    if section_axis == "X":
-        other_axis = "Y"
-    elif section_axis == "Y":
-        other_axis = "X"
-
-    col = len(df.columns)
-    if col < df_axis:
-        logger.error("Finding shortest path, dataframe axis out of range")
-
-    else:
-        df.sort_values([other_axis, "Z"], ascending=[True, False], inplace=True)
-        xys = df[other_axis].unique()
-        zs = df["Z"].unique()
-        rows = len(zs)
-        columns = len(xys)
-        # values are laid out column-major (column j occupies rows
-        # n:n+rows of the sorted dataframe, n increasing by rows each column)
-        values = df.iloc[:, df_axis].to_numpy()
-        array = values.reshape(columns, rows).T
-    return array
-
-
-def find_inout_points(velocity_field_array, velocity_parameters):
-    """
-    Looks for the indexes of the inlet and outle in an array.
-    Velocity parameters of anisotropies are used to find the indexes of inlet and outlet.
-    It is assumed that velocity_parameter[0] correspond to the inlet anisotropy and
-    velocity_parameter[len(velocity_parameter)-1] corresponds to the outlet anisotropy
-
-    Parameters
-    ----------
-    velocity_field_array: numpy array, containing values of the velocity field used to find the shortest path
-    velocity_parameters: list of numbers, each value correspond to a velocity assign to an anisotropy involved in intrusion emplacement
-
-    Returns
-    -------
-    inlet: list of indexes, [row index in array, column index in array]
-    outlet: list of indexes, [row index in array, column index in array]
-
-    """
-    inlet_point = [0, 0]
-    outlet_point = [0, 0]
-
-    inlet_velocity = velocity_parameters[0] + 0.1
-    outlet_velocity = velocity_parameters[len(velocity_parameters) - 1] + 0.1
-
-    # inlet: leftmost column containing inlet_velocity, take its last (deepest) row match
-    inlet_mask = velocity_field_array == inlet_velocity
-    col_has_inlet = inlet_mask.any(axis=0)
-    if col_has_inlet.any():
-        col = int(np.argmax(col_has_inlet))
-        rows_matching = np.nonzero(inlet_mask[:, col])[0]
-        inlet_point[0] = rows_matching[-1]
-        inlet_point[1] = col
-
-    # outlet: rightmost column containing outlet_velocity, take its first row match
-    outlet_mask = velocity_field_array == outlet_velocity
-    col_has_outlet = outlet_mask.any(axis=0)
-    if col_has_outlet.any():
-        col = len(col_has_outlet) - 1 - int(np.argmax(col_has_outlet[::-1]))
-        rows_matching = np.nonzero(outlet_mask[:, col])[0]
-        outlet_point[0] = rows_matching[0]
-        outlet_point[1] = col
-
-    return inlet_point, outlet_point
-
-
-def shortest_path(inlet, outlet, time_map):
-    """
-    Look for the shortest path between inlet and outlet, using a time map.
-    In practice, for a given point looks for the neighbour elements which is closer in time.
-    The search starts in the inlet, until it reaches the outlet.
-
-    Parameters
-    ----------
-    inlet: list of indexes (row and column) of inlet point.
-    outlet: list of indexes (row and column) of outlet point.
-    time_map: numpy array, contains values of time, computed using the fast-marching method.
-
-    Returns
-    -------
-    inet: (intrusion network) numpy array of same shape of time_map array.
-        0 where the intrusion network is found, 1 above it, and -1 below it
-
-    """
-    inet = np.ones_like(time_map)  # array to save shortest path with zeros
-    temp_inlet = inlet  # temporary inlet
-    inet[temp_inlet[0], temp_inlet[1]] = 0
-    i = 0
-
-    while True:
-        i = i + 1
-
-        neighbors = element_neighbour(
-            temp_inlet, time_map, inet
-        )  # identify neighbours elements of temporary outlet
-        direction = index_min(neighbors)  # obtain the location (index min) of minimun difference
-
-        if direction == 10:
-            break
-
-        temp_inlet = new_inlet(temp_inlet, direction)
-        row = temp_inlet[0]
-        col = temp_inlet[1]
-
-        # if row >= n_rows or col >= n_cols:
-        #     break
-        # else:
-        inet[row, col] = 0
-
-        if temp_inlet[0] == outlet[0] and temp_inlet[1] == outlet[1]:
-            break
-        else:
-            continue
-
-    # Assign -1 to points below intrusion network.
-    # For each column, find the first row where inet == 0 and set everything
-    # below it to -1. Columns with no zero are left untouched (matches the
-    # original loop, where h would reach the last row without breaking and
-    # inet[(h + 1):, j] = -1 is then a no-op empty slice).
-    mask_zero = inet == 0
-    has_zero = mask_zero.any(axis=0)
-    first_zero_row = np.argmax(mask_zero, axis=0)
-    row_idx = np.arange(inet.shape[0])[:, None]
-    below_mask = (row_idx > first_zero_row[None, :]) & has_zero[None, :]
-    inet[below_mask] = -1
-
-    return inet
-
-
-def element_neighbour(index, array, inet):
-    """
-    Identify the value of neighbours elements for a given element.
-
-    Parameters
-    ----------
-    index: list of indexes (row and column) of elements.
-    array: numpy array, containing values
-    inet: numpy array, temporal intrusion network. If one of the elements is already inet=0, the assign -1.
-
-    Returns
-    -------
-    values: numpy array, 1x5 with the values of the neighbours of a particular element
-
-    """
-
-    rows = len(array) - 1  # max index of rows of time_map array
-    cols = len(array[0]) - 1  # max index of columns of time_map arrays
-
-    # fixed offsets of the 8 neighbours, in the same order as the original
-    # k/h indices (0: above-left, 1: above, 2: above-right, 3: left, 4: right,
-    # 5: below-left, 6: below, 7: below-right)
-    offsets = np.array(
-        [
-            [-1, -1],
-            [-1, 0],
-            [-1, 1],
-            [0, -1],
-            [0, 1],
-            [1, -1],
-            [1, 0],
-            [1, 1],
-        ]
-    )
-    neighbour_idx = np.asarray(index) + offsets
-    valid = (
-        (neighbour_idx[:, 0] >= 0)
-        & (neighbour_idx[:, 0] <= rows)
-        & (neighbour_idx[:, 1] >= 0)
-        & (neighbour_idx[:, 1] <= cols)
-    )
-
-    values = np.full(8, -1.0)
-    if valid.any():
-        valid_rows = neighbour_idx[valid, 0]
-        valid_cols = neighbour_idx[valid, 1]
-        values[valid] = array[valid_rows, valid_cols]
-
-        # check if some of the neighbours is already part of the intrusion network
-        already_in_network = inet[valid_rows, valid_cols] == 0
-        valid_positions = np.nonzero(valid)[0]
-        values[valid_positions[already_in_network]] = -2
-
-    return values
-
-
-def index_min(array):
-    """
-    Given an array of 1x8, dentify the index of the minimum value within the array.
-
-    Parameters
-    ----------
-    array: numpy array, 1x8
-
-    Returns
-    -------
-    index_min: integer, index of minimum value in array
-
-    """
-
-    # return the index value of the minimum value in an array of 1x8
-    array = np.asarray(array)
-    mask = array >= 0
-
-    if mask.any():
-        masked = np.where(mask, array, np.inf)
-        minimum_val = masked.min()
-        # original loop keeps overwriting index_min for every matching key
-        # in increasing order, so ties resolve to the LAST (highest) index
-        matches = np.nonzero(masked == minimum_val)[0]
-        index_min = int(matches[-1])
-    else:
-        index_min = 10
-
-    return index_min
-
-
-def new_inlet(inlet, direction):
-    """
-    Determine new inlet indexes, given current inlet and direction of minimum difference in time map.
-
-    Parameters
-    ----------
-    inlet: list of indexes [row, column]
-    direction: integers e[0,7] (0: above-left, 1: above, 2: above right, 3: left, 4: right, 5: below left, 6: below, 7:below right)
-
-    Returns
-    -------
-    new_inlet: list of indexes [row, column]
-
-    """
-    pot_new_inlets = {}
-
-    pot_new_inlets.update({"0": np.array([inlet[0] - 1, inlet[1] - 1])})
-    pot_new_inlets.update({"1": np.array([inlet[0] - 1, inlet[1]])})
-    pot_new_inlets.update({"2": np.array([inlet[0] - 1, inlet[1] + 1])})
-    pot_new_inlets.update({"3": np.array([inlet[0], inlet[1] - 1])})
-    pot_new_inlets.update({"4": np.array([inlet[0], inlet[1] + 1])})
-    pot_new_inlets.update({"5": np.array([inlet[0] + 1, inlet[1] - 1])})
-    pot_new_inlets.update({"6": np.array([inlet[0] + 1, inlet[1]])})
-    pot_new_inlets.update({"7": np.array([inlet[0] + 1, inlet[1] + 1])})
-    new_inlet = np.zeros(2)
-
-    for key, value in pot_new_inlets.items():
-        if key == str(direction):
-            new_inlet = value
-    return new_inlet
-
-
-def grid_from_array(array, fixed_coord, lower_extent, upper_extent):
-    """
-    Create an numpy matrix of [i,j,x,y,z,values in array], given an array of 2 dimensions (any combination between x, y an z)
-
-    Parameters
-    ----------
-    array: numpy array, two dimension. Represents a cross section of the model, and its values could be any property
-    fixed_coord: list, containing coordinate and value,
-            ie, [0,2] means section is in x=2, or [1, .45] means sections is in y= 0.45
-            the cross section is along this coordinate
-    lower_extent: numpy array 1x3, lower extent of the model
-    upper_extent: numpy array 1x3, upper extent of the model
-
-    Returns
-    -------
-    values: numpy matrix of [i,j,x,y,z,values in array]
-            (i,j) indexed in array
-            (x,y,z) coordinates considering lower and upper extent of model
-            values, from array
-
-    """
-
-    array = np.asarray(array)
-    spacing_i = len(array)  # number of rows
-    spacing_j = len(array[0])  # number of columns
-    values = np.zeros([spacing_i * spacing_j, 6])
-
-    # original loops iterate outer j, inner i, with l incrementing each
-    # inner step, so i is the fast-varying axis and j the slow-varying axis
-    i_flat = np.tile(np.arange(spacing_i), spacing_j)
-    j_flat = np.repeat(np.arange(spacing_j), spacing_i)
-    array_vals = array[spacing_i - 1 - i_flat, j_flat]
-
-    if fixed_coord[0] == "X":
-        y = np.linspace(lower_extent[1], upper_extent[1], spacing_j)
-        z = np.linspace(lower_extent[2], upper_extent[2], spacing_i)
-        values[:, 0] = i_flat
-        values[:, 1] = j_flat
-        values[:, 2] = fixed_coord[1]
-        values[:, 3] = y[j_flat]
-        values[:, 4] = z[i_flat]
-        values[:, 5] = array_vals
-
-    if fixed_coord[0] == "Y":
-        x = np.linspace(lower_extent[0], upper_extent[0], spacing_j)
-        z = np.linspace(lower_extent[2], upper_extent[2], spacing_i)
-        values[:, 0] = i_flat
-        values[:, 1] = j_flat
-        values[:, 2] = x[j_flat]
-        values[:, 3] = fixed_coord[1]
-        values[:, 4] = z[i_flat]
-        values[:, 5] = array_vals
-
-    if fixed_coord[0] == "Z":
-        x = np.linspace(lower_extent[0], upper_extent[0], spacing_j)
-        y = np.linspace(lower_extent[1], upper_extent[1], spacing_i)
-        values[:, 0] = spacing_i - 1 - i_flat
-        values[:, 1] = spacing_j - 1 - j_flat
-        values[:, 2] = x[j_flat]
-        values[:, 3] = y[i_flat]
-        values[:, 4] = fixed_coord[1]
-        values[:, 5] = array_vals
-
-    return values
diff --git a/ROADMAP.md b/ROADMAP.md
index 3ee00c1d..ab97f260 100644
--- a/ROADMAP.md
+++ b/ROADMAP.md
@@ -457,7 +457,28 @@ just at release time.
   using the Stage 3 YAML schema as the serialization contract and the
   `GeologicalModel` API as a compat facade.
 - [ ] **Stage 6 — Intrusion workflow (outcome 8).** Dedicated design
-  discussion once the graph backend lands.
+  discussion once the graph backend lands. `INTRUSIONS.md` (added
+  2026-08-05) is the input for that discussion: a full review of the
+  current module (bugs, dead code, design smells), a user guide for the
+  data/parameters it actually requires, and a phased hardening plan
+  (A: fix/delete dead code, B: validate the data contract, C: cover
+  `intrusion_steps`/`marginal_faults` with tests before touching them,
+  D: the larger simplification — split the god-object builder, fix the
+  conceptual-model calling convention, decide the fate of the shortest-path
+  method). **Phases A-C done 2026-08-05** (didn't need to wait for Stage 5):
+  fixed a silently-dropped build weight (`gxygz`/`gyxgz` typo), deleted
+  confirmed dead code, made `geometric_scaling_parameters` fail fast with a
+  clear message, added data-contract validation at the
+  `create_and_add_intrusion` boundary, and added the first-ever test/example
+  for `marginal_faults`. Attempting the same for `intrusion_steps` surfaced
+  a real regression: it no longer works at all against the current
+  `StratigraphicColumn` object (an unrelated earlier refactor moved that API
+  from a nested dict to an object with different lookup semantics, and
+  nothing updated the intrusions module to match) — pinned by a test rather
+  than fixed, since fixing it for real is a Phase D design decision, not a
+  mechanical patch. Full detail in `INTRUSIONS.md`'s finding 1 and its
+  Phase A-C status notes. Phase D remains the open item for the dedicated
+  discussion.
 
 ## Status log
 
diff --git a/pyproject.toml b/pyproject.toml
index 06eff0af..e2bd24cb 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -282,7 +282,6 @@ allow-dict-calls-with-keyword-arguments = true
 "LoopStructural/modelling/intrusions/intrusion_feature.py" = ["D", "ANN"]
 "LoopStructural/modelling/intrusions/intrusion_frame.py" = ["D", "ANN"]
 "LoopStructural/modelling/intrusions/intrusion_frame_builder.py" = ["D", "ANN"]
-"LoopStructural/modelling/intrusions/intrusion_support_functions.py" = ["D", "ANN"]
 "LoopStructural/utils/__init__.py" = ["D", "ANN"]
 "LoopStructural/utils/_api_registry.py" = ["D", "ANN"]
 "LoopStructural/utils/_log_sinks.py" = ["D", "ANN"]
diff --git a/tests/unit/modelling/intrusions/test_intrusions.py b/tests/unit/modelling/intrusions/test_intrusions.py
index e700114f..7d3b1868 100644
--- a/tests/unit/modelling/intrusions/test_intrusions.py
+++ b/tests/unit/modelling/intrusions/test_intrusions.py
@@ -1,3 +1,7 @@
+import numpy as np
+import pandas as pd
+import pytest
+
 # Loop library
 from LoopStructural import GeologicalModel
 from LoopStructural.datasets import load_tabular_intrusion
@@ -150,6 +154,301 @@ def counting_prepare_data(*args, **kwargs):
     assert intrusion_builder._up_to_date is True
 
 
+def test_intrusion_gyxgz_weight_reaches_frame_build():
+    """Regression test for a `gxygz`/`gyxgz` typo in
+    `GeologicalModel._build_intrusion` that silently dropped the caller's
+    coordinate-2 orthogonality weight: it was passed to
+    `IntrusionFrameBuilder.build` under the wrong keyword (`gxygz`), which
+    `StructuralFrameBuilder.build` doesn't recognise, so it fell into
+    `**kwargs` and `w3` silently stayed at its default. See INTRUSIONS.md
+    finding 3.
+    """
+    model = GeologicalModel(boundary_points[0, :], boundary_points[1, :])
+    model.data = data
+    model.nsteps = [10, 10, 10]
+
+    conformable_feature = model.create_and_add_foliation("stratigraphy")
+
+    captured_kwargs = {}
+    original_build = IntrusionFrameBuilder.build
+
+    def capturing_build(self, *args, **kwargs):
+        captured_kwargs.update(kwargs)
+        return original_build(self, *args, **kwargs)
+
+    IntrusionFrameBuilder.build = capturing_build
+    try:
+        model.create_and_add_intrusion(
+            "tabular_intrusion",
+            "tabular_intrusion_frame",
+            intrusion_frame_parameters={
+                "contact": "roof",
+                "contact_anisotropies": [conformable_feature],
+            },
+            intrusion_lateral_extent_model=ellipse_function,
+            intrusion_vertical_extent_model=constant_function,
+            gyxgz=5,
+        )
+    finally:
+        IntrusionFrameBuilder.build = original_build
+
+    assert captured_kwargs.get("gyxgz") == 5
+
+
+def test_intrusion_geometric_scaling_not_implemented():
+    """`geometric_scaling_parameters` is only reached when one of the two
+    contacts (roof/floor) has no data. Before this test, that path always
+    raised `NotImplementedError` regardless of what was passed in, several
+    calls deep inside a helper that looked like it partially worked (see
+    INTRUSIONS.md finding 2). This pins the simplified, immediate failure
+    so the behaviour stays an explicit "not supported" rather than silently
+    regressing to something that looks like it works.
+    """
+    model = GeologicalModel(boundary_points[0, :], boundary_points[1, :])
+    single_contact_data = data[
+        ~(
+            (data["feature_name"] == "tabular_intrusion")
+            & (data["intrusion_contact_type"] == "floor")
+        )
+    ].copy()
+    # `tabular_intrusion.csv` doesn't have nx/ny/nz/tx/ty/tz columns, so it
+    # must be routed through `prepare_data` before evaluating any feature
+    # end-to-end (`model.data = ...` alone does not normalise columns) -- see
+    # INTRUSIONS.md finding 1b.
+    model.data = model.prepare_data(single_contact_data)
+    model.nsteps = [10, 10, 10]
+
+    conformable_feature = model.create_and_add_foliation("stratigraphy")
+
+    intrusion_feature = model.create_and_add_intrusion(
+        "tabular_intrusion",
+        "tabular_intrusion_frame",
+        intrusion_frame_parameters={
+            "contact": "roof",
+            "contact_anisotropies": [conformable_feature],
+        },
+        intrusion_lateral_extent_model=ellipse_function,
+        intrusion_vertical_extent_model=constant_function,
+        geometric_scaling_parameters={"thickness": 5.0},
+    )
+
+    with pytest.raises(
+        NotImplementedError, match="geometric_scaling_parameters is not currently supported"
+    ):
+        intrusion_feature.evaluate_value(np.array([[2.5, 2.5, 1.5]]))
+
+
+def test_intrusion_missing_data_raises_clear_error():
+    """Before this test, an intrusion_name/intrusion_frame_name typo (or a
+    feature_name with no matching rows) would fail deep inside
+    `IntrusionFrameBuilder`/`IntrusionBuilder` with a bare `KeyError`, not at
+    the `create_and_add_intrusion` boundary. See INTRUSIONS.md finding 4.
+    """
+    model = GeologicalModel(boundary_points[0, :], boundary_points[1, :])
+    model.data = model.prepare_data(data)
+    model.nsteps = [10, 10, 10]
+
+    conformable_feature = model.create_and_add_foliation("stratigraphy")
+
+    with pytest.raises(ValueError, match="No data found for intrusion 'not_a_real_feature'"):
+        model.create_and_add_intrusion(
+            "not_a_real_feature",
+            "tabular_intrusion_frame",
+            intrusion_frame_parameters={
+                "contact": "roof",
+                "contact_anisotropies": [conformable_feature],
+            },
+            intrusion_lateral_extent_model=ellipse_function,
+            intrusion_vertical_extent_model=constant_function,
+        )
+
+    with pytest.raises(ValueError, match="No data found for intrusion frame 'not_a_real_frame'"):
+        model.create_and_add_intrusion(
+            "tabular_intrusion",
+            "not_a_real_frame",
+            intrusion_frame_parameters={
+                "contact": "roof",
+                "contact_anisotropies": [conformable_feature],
+            },
+            intrusion_lateral_extent_model=ellipse_function,
+            intrusion_vertical_extent_model=constant_function,
+        )
+
+
+def test_intrusion_missing_contact_type_column_raises_clear_error():
+    model = GeologicalModel(boundary_points[0, :], boundary_points[1, :])
+    incomplete_data = model.prepare_data(data).drop(columns=["intrusion_contact_type"])
+    model.data = incomplete_data
+    model.nsteps = [10, 10, 10]
+
+    conformable_feature = model.create_and_add_foliation("stratigraphy")
+
+    with pytest.raises(ValueError, match="missing required column"):
+        model.create_and_add_intrusion(
+            "tabular_intrusion",
+            "tabular_intrusion_frame",
+            intrusion_frame_parameters={
+                "contact": "roof",
+                "contact_anisotropies": [conformable_feature],
+            },
+            intrusion_lateral_extent_model=ellipse_function,
+            intrusion_vertical_extent_model=constant_function,
+        )
+
+
+def test_intrusion_missing_contact_anisotropies_raises_clear_error():
+    model = GeologicalModel(boundary_points[0, :], boundary_points[1, :])
+    model.data = model.prepare_data(data)
+    model.nsteps = [10, 10, 10]
+
+    with pytest.raises(ValueError, match="contact_anisotropies"):
+        model.create_and_add_intrusion(
+            "tabular_intrusion",
+            "tabular_intrusion_frame",
+            intrusion_frame_parameters={"contact": "roof"},
+            intrusion_lateral_extent_model=ellipse_function,
+            intrusion_vertical_extent_model=constant_function,
+        )
+
+
+def _build_marginal_fault_model():
+    """A minimal sill offset by a single bounding fault, to exercise the
+    `marginal_faults` code path in `IntrusionFrameBuilder`
+    (`set_marginal_faults_parameters`/`create_constraints_for_c0`), which
+    -- unlike the plain tabular-intrusion path -- has no test or example
+    anywhere in the codebase (see INTRUSIONS.md finding 1). Reuses the
+    `stratigraphy` feature and the `tabular_intrusion`/`tabular_intrusion_frame`
+    geometry from `load_tabular_intrusion()` (proven to build end-to-end),
+    just renamed, plus one fault.
+    """
+    stratigraphy_rows = data[data["feature_name"] == "stratigraphy"].copy()
+
+    fault_rows = pd.DataFrame(
+        [
+            [2.5, 2.5, 2.5, 0, 1, 0, 0, "marginal_fault", 0],
+            [2.5, 2.5, 2.5, 1, 0, 0, 1, "marginal_fault", 0],
+            [2.5, 2.5, 2.5, 0, 0, 1, 2, "marginal_fault", 0],
+        ],
+        columns=["X", "Y", "Z", "nx", "ny", "nz", "coord", "feature_name", "val"],
+    )
+
+    sill_frame_rows = pd.DataFrame(
+        [
+            [2.00, 2.00, 2.00, 0, np.nan, 0, 0, -1, "sill_frame"],
+            [3.00, 1.00, 2.00, 0, np.nan, 0, 0, -1, "sill_frame"],
+            [1.00, 3.00, 2.00, 0, np.nan, 0, 0, -1, "sill_frame"],
+            [3.00, 2.00, 1.00, 1, 0, np.nan, np.nan, np.nan, "sill_frame"],
+            [3.00, 2.00, 1.00, 1, np.nan, 0, 1, 0, "sill_frame"],
+            [2.50, 1.00, 1.00, 2, 0, np.nan, np.nan, np.nan, "sill_frame"],
+            [2.50, 2.00, 1.00, 2, 0, np.nan, np.nan, np.nan, "sill_frame"],
+            [2.50, 2.00, 1.00, 2, np.nan, 1, 0, 0, "sill_frame"],
+        ],
+        columns=["X", "Y", "Z", "coord", "val", "gx", "gy", "gz", "feature_name"],
+    )
+
+    # Roof/floor contact points spanning both sides of the Y=2.5 fault, and
+    # a handful flagged for the lateral (side) extent -- same shape as
+    # `tabular_intrusion`'s own contact data, just relabelled.
+    sill_contact_rows = pd.DataFrame(
+        [
+            [3.04, 2.12, 2.18, "roof", False],
+            [4.02, 3.85, 2.14, "roof", True],
+            [2.02, 3.16, 2.02, "roof", False],
+            [2.12, 2.16, 2.02, "roof", False],
+            [1.14, 3.50, 2.04, "roof", True],
+            [4.08, 3.02, 1.18, "floor", True],
+            [1.14, 1.06, 1.16, "floor", True],
+            [4.20, 2.18, 1.12, "floor", True],
+            [1.02, 3.02, 1.12, "floor", True],
+        ],
+        columns=["X", "Y", "Z", "intrusion_contact_type", "intrusion_side"],
+    )
+    sill_contact_rows["feature_name"] = "sill"
+
+    model_data = pd.concat(
+        [stratigraphy_rows, fault_rows, sill_frame_rows, sill_contact_rows],
+        ignore_index=True,
+    )
+
+    model = GeologicalModel(boundary_points[0, :], boundary_points[1, :])
+    model.data = model.prepare_data(model_data)
+    model.nsteps = [10, 10, 10]
+
+    stratigraphy = model.create_and_add_foliation("stratigraphy")
+    fault = model.create_and_add_fault("marginal_fault", 0.5, nelements=1e4)
+
+    return model, stratigraphy, fault
+
+
+def test_intrusion_marginal_faults():
+    model, stratigraphy, fault = _build_marginal_fault_model()
+
+    intrusion_feature = model.create_and_add_intrusion(
+        "sill",
+        "sill_frame",
+        intrusion_frame_parameters={
+            "contact": "roof",
+            "contact_anisotropies": [stratigraphy],
+            "marginal_faults": {
+                0: {
+                    "structure": fault,
+                    "block": "hanging wall",
+                    "series": stratigraphy,
+                }
+            },
+        },
+        intrusion_lateral_extent_model=ellipse_function,
+        intrusion_vertical_extent_model=constant_function,
+    )
+
+    value = intrusion_feature.evaluate_value(np.array([[2.5, 2.5, 1.5]]))
+    assert np.isfinite(value).all()
+
+
+def test_intrusion_steps_broken_with_current_stratigraphic_column():
+    """`intrusion_steps` is currently unusable, not just untested.
+
+    `IntrusionFrameBuilder.set_intrusion_steps_parameters`
+    (`intrusion_frame_builder.py`) reads
+    `self.model.stratigraphic_column[series_from_name.name][unit_from_name]`,
+    i.e. it expects `model.stratigraphic_column` to be the old nested
+    `{group_name: {unit_name: {...}}}` dict. `GeologicalModel.stratigraphic_column`
+    is now a `StratigraphicColumn` object whose `__getitem__` looks up a
+    single element by `uuid` (not by group/series name), and
+    `GeologicalModel.set_stratigraphic_column` (the old dict-based setter)
+    unconditionally raises `DeprecationWarning` -- so there is no longer any
+    way to put the model into the shape `intrusion_steps` expects. This is a
+    real regression from an unrelated stratigraphic-column refactor that was
+    never propagated to the intrusions module (see INTRUSIONS.md finding 1,
+    revised). This test pins the current failure so a future fix (Stage 6)
+    has a clear target and so this doesn't regress further/silently.
+    """
+    model, stratigraphy, fault = _build_marginal_fault_model()
+    model.stratigraphic_column.add_unit("unitA", thickness=1)
+    model.stratigraphic_column.add_unit("unitB", thickness=1)
+
+    with pytest.raises(KeyError, match="No element found with uuid"):
+        model.create_and_add_intrusion(
+            "sill",
+            "sill_frame",
+            intrusion_frame_parameters={
+                "contact": "roof",
+                "contact_anisotropies": [stratigraphy],
+                "intrusion_steps": {
+                    0: {
+                        "structure": fault,
+                        "unit_from": "unitA",
+                        "series_from": stratigraphy,
+                        "unit_to": "unitB",
+                        "series_to": stratigraphy,
+                    }
+                },
+            },
+            intrusion_lateral_extent_model=ellipse_function,
+            intrusion_vertical_extent_model=constant_function,
+        )
+
+
 # if __name__ == "__main__":
 #     test_intrusion_freame_builder()
 #     test_intrusion_builder()