diff --git a/docs/advanced/crossing-fault-zones.md b/docs/advanced/crossing-fault-zones.md
new file mode 100644
index 00000000..497cf0d5
--- /dev/null
+++ b/docs/advanced/crossing-fault-zones.md
@@ -0,0 +1,346 @@
+# Fault zones that cross
+
+Two fault zones that cross are the geometry the finite-width
+representation exists for. Thickening each trace and letting the CAD
+kernel fuse the results turns the junction into ordinary cells of a
+single region: nothing downstream is told which branch a cell came from,
+and no junction geometry is ever constructed. The zero-thickness
+representation cannot take this route at all — the split refuses
+overlapping cuts.
+
+That ignorance is deliberate, and it has a price. The rheology at the
+overlap comes entirely from the `Surface` lookup, which gives every point
+the **director of the nearest surface**. Where the two zones cross, the
+material is weak in two directions and a single director represents one
+of them. This page works out when that matters.
+
+Scripts: `~/+Simulations/ribbon_network_2d/`.
+
+## Building the crossing
+
+One call places the whole network. The two ribbons may overlap — that is
+the point — and `assembly="fuse"` returns their union as one region with
+no internal seam:
+
+```python
+from underworld3.utilities.place_surface import place_thin_volume
+
+dm, info = place_thin_volume(mesh0.dm, [trace_a, trace_b], width=0.02,
+                             label="Zone", label_value=31,
+                             size=0.005, assembly="fuse")
+```
+
+Two details in that call decide whether the experiment can measure
+anything.
+
+**`size` must be passed.** It defaults to `0.9 * width`, which puts about
+two elements across the band at every width. Two elements cannot resolve
+which half of an overlap a rule owns, and the ownership split is the
+whole question. `width/4` is used throughout here.
+
+**`fragment` is the wrong assembly for a crossing.** It keeps each
+overlap piece as its own region, so the mesh conforms to the boundaries
+of the overlap; where two zones converge at a shallow angle the overlap
+is a spike, and its fragmented tip meshes to arbitrarily bad angles. The
+fused union has no such tip.
+
+## What the overlap should be
+
+Before solving anything, the question has a closed-form answer in 2-D,
+and the closed form says where to look.
+
+In plane strain with an incompressible medium the deviatoric strain-rate
+space is **two-dimensional**, spanned by
+
+$$
+\mathbf{e}_1 = \tfrac{1}{\sqrt 2}\begin{pmatrix} 1 & 0\\ 0 & -1\end{pmatrix},
+\qquad
+\mathbf{e}_2 = \tfrac{1}{\sqrt 2}\begin{pmatrix} 0 & 1\\ 1 & 0\end{pmatrix},
+$$
+
+and rotating the material by $\theta$ rotates a state in that space by
+$2\theta$. A weak plane with normal at $\phi$ is compliant in exactly one
+mode — on-plane shear — so it occupies a single **line** in that space,
+at angle $2\phi + 90^\circ$.
+
+Everything follows from the doubling. Two planes separated by $\Delta$ in
+the rock are separated by $2\Delta$ there, so:
+
+- $\Delta = 0^\circ$ and $\Delta = 90^\circ$ put both compliances on the
+  **same line**. Shear on a plane with normal $\hat{\mathbf x}$ *is*
+  shear on a plane with normal $\hat{\mathbf y}$ — the same $\tau_{xy}$.
+  A 90-degree X is degenerate in 2-D.
+- $\Delta = 45^\circ$ makes the two compliances **orthogonal**. Together
+  they span the whole deviatoric space, and the overlap is isotropic.
+
+So the geometry the phrase "high angle" suggests — the square X — is the
+one geometry where the question cannot be asked. The tensors to compare
+differ most at 45 degrees.
+
+```{warning}
+**This degeneracy is two-dimensional and does not carry to 3-D.** It
+exists because the deviatoric space here has only two dimensions, so two
+lines in it must either coincide or be distinct in one way. In 3-D that
+space has five dimensions, two weak planes at any angle occupy different
+subspaces, and there is no null. Do not carry "a 90-degree crossing is
+safe" into a 3-D model — the 3-D question is open and nothing on this
+page answers it.
+```
+
+### The orthotropic tensor is a TI tensor
+
+Summing the two planes' compliances and inverting gives a tensor with two
+eigenvalues and one principal direction in the deviatoric plane, which is
+precisely the structure of
+{py:class}`~underworld3.constitutive_models.TransverseIsotropicFlowModel`.
+So the orthotropic treatment needs no new constitutive machinery in 2-D:
+read off an effective director and effective $\eta_0$, $\eta_1$ and use
+the shipped model.
+
+```python
+S = np.eye(2) / (2 * eta_0)
+s = 1 / (2 * eta_1) - 1 / (2 * eta_0)
+for phi, lam in zip(normal_angles, weights):
+    u = weak_mode(phi)                    # unit vector at 2*phi + 90
+    S = S + lam * s * np.outer(u, u)
+C = np.linalg.inv(S)
+eta_1_eff, eta_0_eff = sorted(np.linalg.eigvalsh(C) / 2)
+```
+
+### Two candidate rules, not one
+
+The weights are a modelling choice and the two ends of the family behave
+differently:
+
+| $\Delta$ | compliance **sum** ($\lambda = 1$ each) | compliance **average** ($\lambda = \tfrac12$ each) |
+|---|---|---|
+| $0^\circ$ | anisotropy 199 — twice as weak as one plane | exactly the single plane |
+| $45^\circ$ | isotropic, $\eta_1$ | isotropic, $2\eta_1$ |
+| $90^\circ$ | anisotropy 199 | exactly the single plane |
+
+The sum is what "add the compliances" means literally: two independent
+fabrics superposed. Read it at $\Delta = 0$ — two *coincident* weak
+planes come out twice as weak as one — and note that the excess does not
+switch off away from coincidence. The average treats the overlap as a
+50/50 mixture of the two fabrics and degenerates to the single plane
+exactly. Both are carried below; neither is obviously right.
+
+```{important}
+There is **no external reference solution** here, and there cannot be one
+without modelling the fabric itself. Both constructions are assumptions,
+so what is measured below is the distance from the shipped rule to a
+*choice*, not to a truth. What makes that worth reading is that the two
+choices differ by a factor of two in how weak they make the overlap and
+still bracket the answer tightly — 61.8% against 70.0% at the worst
+angle, 92.3% against 93.0% at $45^\circ$. The conclusion is insensitive
+to the part that had to be assumed.
+```
+
+## The experiment
+
+Two straight ribbons crossing at the centre of a unit box, fused into one
+zone, driven in horizontal extension with a free top. The crossing angle
+$\Delta$ is swept, and five rheologies are solved **on the same mesh**:
+
+| case | the zone | the overlap |
+|---|---|---|
+| `none` | uniform $\eta_0$ | — (the no-fault control) |
+| `isotropic` | weak in every mode | weak in every mode |
+| `nearest` | TI, nearest trace's director | TI, nearest trace's director |
+| `ortho_sum` | TI, nearest trace's director | the compliance-sum tensor |
+| `ortho_avg` | TI, nearest trace's director | the compliance-average tensor |
+
+The last three are identical everywhere except the overlap cells, so any
+difference between them is the overlap rheology and nothing else.
+`isotropic` weakens the arms too, so it is not a candidate rule — only a
+bracket on how compliant the network could be.
+
+The metric is the **dissipation weakening**
+
+$$
+W = 1 - \frac{D}{D_{\rm no\ fault}},
+$$
+
+referenced to the `none` control on the same mesh at the same angle. It
+is a single number with no probe geometry in it. (The control dissipates
+exactly 4.000000 in this box, which is the analytic value, so the
+integral is checked as well.)
+
+Note what $W$ is made of at a shallow crossing: at $\Delta = 10^\circ$
+the two dissipations are 3.9521 and 3.9225 against a control of 4.0000,
+so the whole signal is a difference of about 0.03 between numbers that
+agree to 1%. The exact control and the resolution study below are what
+make that difference readable; a single run at one resolution would not
+be.
+
+### What comes out
+
+| $\Delta$ | $W$ nearest | $W$ orthotropic | $W$ isotropic | nearest, as % of orthotropic |
+|---|---|---|---|---|
+| $10^\circ$ | 1.198% | 1.710% | 3.628% | **70.0%** |
+| $20^\circ$ | 3.687% | 4.520% | 6.256% | 81.6% |
+| $30^\circ$ | 6.618% | 7.679% | 9.240% | 86.2% |
+| $45^\circ$ | 11.566% | 12.435% | 13.646% | 93.0% |
+| $60^\circ$ | 15.789% | 16.253% | 17.110% | 97.1% |
+| $75^\circ$ | 18.322% | 18.476% | 19.121% | 99.2% |
+| $90^\circ$ | 18.918% | 18.918% | 19.465% | **100.0%** |
+
+```{figure} figures/crossing_weakening.png
+:width: 100%
+
+Left: what the two faults are worth, against the no-fault control on the
+same mesh. Right: the overlap rule on its own — the two runs share the
+mesh and the arms exactly, so the ratio is the overlap and nothing else.
+```
+
+Read the last row first. At $\Delta = 90^\circ$ the two dissipations
+agree to seven significant figures (3.243286 against 3.243286). That is
+not a tuned result: the algebra says the orthotropic overlap tensor there
+is *identically* the single plane's, and a 2-D director is only defined
+mod $90^\circ$, so the nearest-fault rule assigns that same tensor on
+both sides of the medial axis. The experiment reproduces it exactly, and
+that is what licenses the other six rows.
+
+Now read up the column. The nearest-fault rule loses a substantial part
+of the network's weakening at a shallow crossing, and the error closes
+monotonically as the crossing squares up. The percentages in that table
+are computed at a zone mesh of $w/4$ and are **bounds**, not estimates —
+refining the zone mesh moves the rule closer to the orthotropic one, and
+the converged values are about 76% at $10^\circ$ and 95% at $45^\circ$
+(see the resolution control below).
+
+```{figure} figures/crossing_map.png
+:width: 100%
+
+The same three answers in physical space at $\Delta = 45^\circ$, on one
+colour scale. They are indistinguishable by eye, and the region in
+dispute is 34 cells out of 1480 (right). A difference of this kind is
+only visible in an integral — which is why the metric above is a
+dissipation and not a picture.
+```
+
+### The error is largest where the issue expected none
+
+The per-cell tensors differ most at $\Delta = 45^\circ$, and the measured
+curve shows no feature there at all. Both facts are right, and the
+resolution is geometric: the overlap **area** goes as $w^2/\sin\Delta$.
+A shallow crossing has a long lens of shared material, a square one has
+almost none — 141 overlap cells at $10^\circ$ against 24 at $90^\circ$.
+The tensor error peaks at $45^\circ$, the area falls monotonically, and
+the product is dominated by the area.
+
+So the geometry to worry about is the **shallow** crossing, not the
+high-angle one — the opposite of what "high angle" suggests, and the
+opposite of the case a designer would think to check.
+
+## Measuring the slip is harder than it looks
+
+Slip is read as the tangential velocity jump between probes either side
+of the band. A **fixed** standoff — some fraction of the width — is wrong
+at a crossing: near a shallow one, a probe placed at $0.75\,w$ from one
+trace sits *inside the other ribbon*, in weak material, and reports
+almost no jump. At $\Delta = 10^\circ$ that read 0.006 against a peak of
+0.18 on the same fault.
+
+The probes have to be walked clear of **every** ribbon before they are
+read. The jump then means the same thing at every angle: the velocity
+difference across the whole weak region.
+
+Even walked out, the local probe **cannot** be trusted at the crossing
+itself: that is exactly where the standoff changes discontinuously, and
+the reading there dips sharply — at $\Delta = 90^\circ$ it goes negative
+once the no-fault baseline is subtracted, which is not a physical slip.
+This is the reason the metric above is a dissipation integral. Use the
+slip probe for the arms and for pictures, not for the number that decides
+the question.
+
+Two further cautions on reading absolute slip here:
+
+- The pair is symmetric about the loading axes, so changing $\Delta$
+  necessarily changes how well each plane is oriented for slip. At
+  $\Delta = 10^\circ$ both planes are nearly vertical — nearly principal
+  planes — and barely slip at all. Compare ratios at fixed angle, not
+  slip across angles.
+- The no-fault control read through the same probes gives the wall rock's
+  own contribution over the probe separation. Subtract it.
+
+## How much of this is the mesh
+
+The overlap is a few tens of cells, so the numbers above are worth
+nothing without a resolution control. At $\Delta = 45^\circ$:
+
+| host cell | zone $h$ | zone cells | overlap | $W$ nearest | $W$ orthotropic | ratio |
+|---|---|---|---|---|---|---|
+| 0.030 | 0.0100 | 380 | 12 | 10.9522% | 12.4148% | 88.2% |
+| 0.030 | 0.0050 | 1480 | 34 | 11.5662% | 12.4351% | 93.0% |
+| 0.030 | 0.0025 | 5342 | 130 | 11.7717% | 12.4565% | 94.5% |
+| 0.040 | 0.0050 | 1480 | 34 | 11.5593% | 12.4278% | 93.0% |
+| 0.020 | 0.0050 | 1480 | 34 | 11.5714% | 12.4404% | 93.0% |
+
+The **host** mesh does not enter: 0.04, 0.03 and 0.02 give 93.0% three
+times over, with $W$ agreeing to four significant figures. The embedded
+zone mesh is identical in all three, which is the embedding working as
+designed and worth having as a measurement rather than a claim.
+
+The **zone** mesh does. Halving it twice moves the ratio 88.2 → 93.0 →
+94.5%, and the increments fall by about a factor of three, so the
+converged value here is near 95%.
+
+The same check at the worst angle, $\Delta = 10^\circ$, where the overlap
+is largest and the resolution matters most:
+
+| host cell | zone $h$ | zone cells | overlap | $W$ nearest | $W$ orthotropic | ratio |
+|---|---|---|---|---|---|---|
+| 0.030 | 0.0050 | 1374 | 141 | 1.1980% | 1.7103% | 70.0% |
+| 0.030 | 0.0025 | 4910 | 522 | 1.2731% | 1.7252% | 73.8% |
+| 0.030 | 0.0013 | 19562 | 2051 | 1.3030% | 1.7361% | 75.1% |
+
+Increments of +3.8 then +1.3, the same threefold contraction, giving a
+converged value near 76%.
+
+Refinement always moves the nearest-fault rule *closer* to the
+orthotropic one, so the sweep table earlier on this page **overstates the
+error** everywhere. The corrected statement is that the rule loses about
+a quarter of the network's weakening at the worst crossing angle, not a
+third — which changes the size of the effect and none of its structure.
+
+## `clearance` is not monotone
+
+The carve deletes host vertices near the assembly's skin, and whether
+what is left is one simple hole is a topological accident of the host
+mesh against a plus-shaped cavity. It does not improve monotonically with
+`clearance`: a 30-degree crossing builds only at `clearance=1.3`, a
+45-degree crossing at 1.0 but *not* at 1.3. Walk a ladder and record
+which value was used.
+
+```python
+for clearance in (1.0, 0.9, 1.1, 0.8, 1.2, 0.7, 1.3):
+    try:
+        dm, info = place_thin_volume(..., clearance=clearance)
+    except RuntimeError:
+        continue
+    break
+```
+
+## Postscript: this test was too kind
+
+Everything above is measured in pure extension with the fault pair
+symmetric about the load — a drive with **no rotation**. Repeating the
+question in a shear box with an asymmetric branching geometry (spin in
+the kinematics, nothing mirror-symmetric) inverts the verdict: the
+difference between director rules stays immaterial, but the *entire
+TI-arms family* loses factors — not percent — against the isotropic
+reference wherever the network must hand slip between differently
+oriented members. The overlap tensor was never the dominant term; the
+arms' constraint is. See
+[Branching Faults and Junctions](fault-branching-junctions.md), which
+supersedes the recommendation implied here for any network that
+redistributes slip.
+
+## See also
+
+- [Fault networks](fault-networks.md) — the zero-thickness toolkit, and
+  the offset-junction preparation the ribbon does not need.
+- [Split-node faults](split-node-faults.md) — the contact representation.
+- [Transverse isotropy](vep-transverse-isotropy-faults.md) — the weak-plane
+  rheology and its parameters.
diff --git a/docs/advanced/fault-branching-junctions.md b/docs/advanced/fault-branching-junctions.md
new file mode 100644
index 00000000..6847acee
--- /dev/null
+++ b/docs/advanced/fault-branching-junctions.md
@@ -0,0 +1,187 @@
+# Branching faults, junctions, and what a crack model cannot do
+
+The [crossing-zones page](crossing-fault-zones.md) compared overlap
+rheologies for fused fault zones in pure extension and found the
+differences modest at ordinary angles. That test had a hidden weakness:
+the configuration is maximally symmetric about the load, and the drive
+carries **no rotation**. This page repeats the question with both
+removed — and the verdict inverts. It ends at the reason: two of our
+fault representations are crack models, one is not, and the difference
+is not a numerical detail but the mechanics of branching itself.
+
+Scripts and caches: `~/+Simulations/ribbon_network_2d/`
+(`branch_study.py`, `branch_width.py`, `branch_damage.py`).
+
+## The shear box
+
+Simple shear supplies what extension cannot: spin, $\omega =
+\dot\gamma/2$. The test frame is a driven box — top plate $v = (+0.5,
+0)$, bottom $v = (-0.5, 0)$, sides open (traction-free) — holding an
+asymmetric Y: a **sub-horizontal main fault** at $\theta = 20/30/40^\circ$
+and a **horizontal branch** running in from the side and fusing into it.
+Nothing is mirror-symmetric; the branch lies in the shear plane
+(optimally oriented) while the main degrades with $\theta$ (resolved
+drive $\propto \cos 2\theta$). Bands at the union-jack default
+resolution, structure verified per build; the isotropic weak band is the
+reference throughout, with transversely isotropic (TI) variants against
+it. Everything is measured from the P0 cell strain rate — local, no
+recovery step, no probes.
+
+```{figure} figures/branch_maps_z20.png
+:width: 100%
+
+Strain rate (P0 invariant, log scale, one range). Rows: $\theta$ =
+20/30/40$^\circ$; columns: the isotropic reference and three TI
+variants differing only at the junction. In the reference the junction
+is the hottest point of the network and intensifies with $\theta$; in
+every TI variant it cools, and the misoriented main fades until its
+lower arm is a ghost.
+```
+
+## The junction is a slip exchanger — and the crack constraint cannot run it
+
+Slip on each fault, at the junction, TI variants as a fraction of the
+reference:
+
+| $\theta$ | main fault | branch |
+|---|---|---|
+| $20^\circ$ | 77% | ~75% |
+| $30^\circ$ | **46%** | ~78% |
+| $40^\circ$ | **19%** | ~52% |
+
+Three facts organise this. First, the reference's main-fault slip at the
+junction **rises** as the fault becomes misoriented (0.30 → 0.36): the
+optimally-oriented branch delivers its offset *into* the main through
+the junction, the isotropic corner rotating the deformation between the
+two orientations. Second, the three TI overlap treatments (nearest
+director, harmonic average, isotropic overlap cells) are
+indistinguishable from each other while all collapsing against the
+reference — the overlap rule was never the dominant term. Third, the
+mechanism is general: **a gradient of slip along a fault is a
+volumetric exchange between band and surroundings**, and that is
+precisely the deformation mode the TI constraint (like the slip
+surface) removes. Crack kinematics carry slip *along* a fault
+faithfully and cannot hand it *between* faults.
+
+Imposed isotropic "damage halos" around the junction recover the
+exchange roughly linearly in their radius (42% at $3w$, 73% at $5w$)
+with no plateau: there is no compact junction patch, because the
+constraint taxes the mechanics wherever slip varies.
+
+```{figure} figures/halo_recovery.png
+:width: 100%
+
+Slip profiles and recovery against halo radius: freeing the constraint
+where the slip gradient lives, by hand.
+```
+
+## The exit throat, and a law with nothing fitted
+
+Where is the reference's peak? Not at the junction point. It sits on
+the main's up-dip arm at $r \approx 2.1\,w$ — a **self-similar**
+position (2.05–2.11 $w$ across a fourfold width sweep) — where the
+fused Y narrows back to a single band. The combined slip of both faults
+must exit through that neck, and flux conservation predicts its
+intensity outright:
+
+$$
+\gamma_{\rm peak} \;\approx\; \frac{s_{\rm main} + s_{\rm branch}}{w},
+$$
+
+measured at 0.96–1.09 of the prediction over the three resolved widths.
+
+```{figure} figures/junction_peak_zoom.png
+:width: 100%
+
+The peak (circled) at four band widths, one log scale: always at the
+exit throat, always at $r \approx 2.1\,w$, in healthy cells.
+```
+
+The limit is the point: as $w \to 0$ the throat slides onto the corner
+and its amplitude diverges. **A branch point is a corner singularity of
+the continuum problem, and the band width is its regularisation
+length.** Neither description converges there — the crack models
+*exclude* the corner physics; the continuum model *concentrates* it
+without bound. The junction's width is therefore not a numerical
+parameter to refine away but a physical one: the process-zone size. (It
+is the [gouge page's](gouge-zones.md) "width is physics", at the place
+where it matters most.)
+
+## What the country rock feels
+
+Deviatoric stress in the matrix (where $\tau = 2\dot\varepsilon$
+exactly), zones masked:
+
+```{figure} figures/branch_stress_matrix.png
+:width: 100%
+
+Matrix stress, one linear scale. Means agree to 2–4% across all six
+cases — the load is conserved; its **location** is not.
+```
+
+The locked TI junctions park a wall-rock stress concentration at the
+elbow (the tail rises 11%) that the reference does not possess; freeing
+the junction trades it back for **tip lobes**, which is where the
+reference — and a real fault — concentrates stress. Ranking these
+models by peak stress alone would order them backwards; the location is
+the diagnostic. A constrained-network model read at face value predicts
+wall-rock damage and seismicity at junctions where the physics puts
+none, and under-predicts tip process zones.
+
+## Offered plasticity, the crack model builds a bypass
+
+A deliberately primitive damage law — any matrix cell whose stress
+exceeds $\tau_c$ becomes isotropic $\eta_1$, instantly and irreversibly;
+re-solve; repeat — asks what a yielding background does with each
+representation. With a negative control:
+
+```{figure} figures/damage_control_compare.png
+:width: 100%
+
+Damaged cells coloured by the pass at which they yielded, same law and
+threshold. Left, the isotropic reference: growth at the free **tips**
+only — the physical pattern — with the junction silent throughout.
+Right, the fully-TI network: the same tip growth **plus** a
+junction-nucleated front marching along the shear plane — a
+through-going bypass that retires the misoriented fault.
+```
+
+Both models grow faults at tips; that part is physics. Only the
+constrained model manufactures growth at the junction — nucleated by
+the stress its own kinematic restriction parks there — and its response
+to that stress is not to heal the junction but to **re-route the
+network around it**. The law is a cartoon (no softening, no length
+scale, grid-width bands), so the nucleation sites and topology are the
+trustworthy content, not rates or widths. One further consequence of
+the corner singularity: at realistic (thinner) widths the *reference's*
+corner stress also exceeds any finite yield threshold — so in a
+yielding crust a branch point **must** develop a damage zone. The
+damage halo at a junction is not a modelling patch; it is the
+material's own regularisation of a corner that crack mechanics cannot
+represent and continuum mechanics cannot bound.
+
+## Choosing, finally
+
+One of these is a crack model twice over — the split surface imposes
+slip-only kinematics exactly; the TI band imposes the same as a
+$10^{3}$ constitutive penalty — and one is not: the isotropic band is a
+model of fault *rock*, not fault *kinematics*. The crack assumptions
+fight each other at branches. Hence:
+
+- **Large-scale fault networks**: build from surfaces (contacts) — the
+  crack idealisation at its best, where segments are distinct,
+  well-oriented, and slip-dominated.
+- **Zones with interior physics** (gouge, heating, damage evolution):
+  thin volumes/ribbons — with TI where the zone should stay
+  crack-like, isotropic where it must evolve internally.
+- **Crossings and joins**: make them kinematically compatible by
+  construction — bend faults in to meet tangentially rather than
+  abutting at high angle, and surround the joins with isotropic damage
+  zones sized as process zones. Do not expect any director rule, patch,
+  or refinement to buy back a high-angle locked junction: the failure
+  is the crack assumption itself.
+- **Evolution in time** — where a misoriented fault would dismember
+  into en-échelon segments and the network would rebuild itself through
+  its damage field — requires transporting faults as material chains
+  and re-deriving them with the evolving damage. That is a programme,
+  not a feature; this page records the mechanics that motivates it.
diff --git a/docs/advanced/figures/branch_maps_z20.png b/docs/advanced/figures/branch_maps_z20.png
new file mode 100644
index 00000000..455a574d
Binary files /dev/null and b/docs/advanced/figures/branch_maps_z20.png differ
diff --git a/docs/advanced/figures/branch_slip_z20.png b/docs/advanced/figures/branch_slip_z20.png
new file mode 100644
index 00000000..7775a5f2
Binary files /dev/null and b/docs/advanced/figures/branch_slip_z20.png differ
diff --git a/docs/advanced/figures/branch_stress_matrix.png b/docs/advanced/figures/branch_stress_matrix.png
new file mode 100644
index 00000000..204e77c0
Binary files /dev/null and b/docs/advanced/figures/branch_stress_matrix.png differ
diff --git a/docs/advanced/figures/crossing_map.png b/docs/advanced/figures/crossing_map.png
new file mode 100644
index 00000000..7e5f411b
Binary files /dev/null and b/docs/advanced/figures/crossing_map.png differ
diff --git a/docs/advanced/figures/crossing_weakening.png b/docs/advanced/figures/crossing_weakening.png
new file mode 100644
index 00000000..ee8dcd61
Binary files /dev/null and b/docs/advanced/figures/crossing_weakening.png differ
diff --git a/docs/advanced/figures/damage_control_compare.png b/docs/advanced/figures/damage_control_compare.png
new file mode 100644
index 00000000..10416087
Binary files /dev/null and b/docs/advanced/figures/damage_control_compare.png differ
diff --git a/docs/advanced/figures/gouge_insulation.png b/docs/advanced/figures/gouge_insulation.png
new file mode 100644
index 00000000..503683c9
Binary files /dev/null and b/docs/advanced/figures/gouge_insulation.png differ
diff --git a/docs/advanced/figures/gouge_mechanics.png b/docs/advanced/figures/gouge_mechanics.png
new file mode 100644
index 00000000..c931a042
Binary files /dev/null and b/docs/advanced/figures/gouge_mechanics.png differ
diff --git a/docs/advanced/figures/gouge_thermal.png b/docs/advanced/figures/gouge_thermal.png
new file mode 100644
index 00000000..8c59057b
Binary files /dev/null and b/docs/advanced/figures/gouge_thermal.png differ
diff --git a/docs/advanced/figures/halo_recovery.png b/docs/advanced/figures/halo_recovery.png
new file mode 100644
index 00000000..15f4e919
Binary files /dev/null and b/docs/advanced/figures/halo_recovery.png differ
diff --git a/docs/advanced/figures/junction_peak_zoom.png b/docs/advanced/figures/junction_peak_zoom.png
new file mode 100644
index 00000000..4934334a
Binary files /dev/null and b/docs/advanced/figures/junction_peak_zoom.png differ
diff --git a/docs/advanced/gouge-zones.md b/docs/advanced/gouge-zones.md
new file mode 100644
index 00000000..eea767c5
--- /dev/null
+++ b/docs/advanced/gouge-zones.md
@@ -0,0 +1,165 @@
+# Gouge zones
+
+A fault with a gouge zone can be built either way. As a **ribbon** it is a
+band of crushed rock of width $w$ with its own material viscosity
+$\eta_1$. As a **contact** it is a surface with an interface law
+$\eta_I$, and the two agree when
+
+$$
+\eta_I = \frac{\eta_1}{w},
+$$
+
+which is the matched collapse the campaign established and this page
+re-checks rather than assumes.
+
+The collapse is exact for the slip. What it is not exact for is anything
+that depends on the state *inside* the zone, and the reason is visible in
+the formula: it preserves the **ratio** $\eta_1/w$ and destroys the pair.
+Any quantity that needs $\eta_1$ and $w$ separately is gone.
+
+Scripts: `~/+Simulations/ribbon_network_2d/`.
+
+## The premise, checked
+
+One straight normal fault in an extending unit box, both ways, at four
+widths with the gouge viscosity held at $\eta_1 = 10^{-3}$ so the
+interface law changes with the width:
+
+| $w$ | $\eta_I = \eta_1/w$ | ribbon peak slip | contact peak slip | ratio |
+|---|---|---|---|---|
+| 0.020 | 0.05 | 0.50255 | 0.44947 | 111.8% |
+| 0.010 | 0.10 | 0.48200 | 0.45374 | 106.2% |
+| 0.005 | 0.20 | 0.46351 | 0.44800 | 103.5% |
+| 0.002 | 0.50 | 0.42800 | 0.42101 | **101.7%** |
+
+```{figure} figures/gouge_mechanics.png
+:width: 100%
+
+Left: slip along the fault, ribbon (solid) against the matched contact
+(dashed), at each width. Right: the two agree in the limit and not
+before it.
+```
+
+The collapse is a $w \to 0$ limit and behaves like one: it is worth 12% at
+$w = 0.02$ and 1.7% at $w = 0.002$. Quote it as an identity only when the
+zone is thin against the fault length.
+
+```{note}
+The contact runs through the rotated machinery, which builds its **own**
+prefixed KSP that `stokes.petsc_options` never reaches. Left on its
+iterative default the velocity sub-solve here hit its 200-iteration cap
+and warned, while the outer SNES still reported convergence. Set
+`stokes._rotated_use_lu = True` on problems this size.
+
+The verdict then has to be physical, because the outer SNES reason does
+not describe that solver either. Use the constraint: the no-opening
+condition is strong, so the **normal** part of the pair jump is machine
+zero when the increment was really solved. It reads $5.5\times10^{-17}$
+against a tangential jump of 0.47 here.
+```
+
+## What the collapse loses
+
+Dissipation per unit **length** of fault is $\eta_I V^2$ — preserved. The
+far-field temperature is therefore the same either way. Dissipation per
+unit **volume** in the gouge is
+
+$$
+\Phi = \eta_1 \left(\frac{V}{w}\right)^2 = \frac{\eta_I V^2}{w},
+$$
+
+which is not preserved: it rises as the gouge narrows at fixed interface
+law. Feed that into conduction across a thin band with its edge
+temperatures as boundary data and the excess at the centre is
+
+$$
+\Delta T = \frac{\Phi\, w^2}{8\kappa} = \frac{\eta_1 V^2}{8\kappa}
+         = \frac{\eta_I V^2}{8\kappa}\, w .
+$$
+
+At fixed $\eta_I$ that is **proportional to the width**. Two models with
+the same interface law and different widths have the same slip and the
+same far field, and gouge temperatures in the ratio of their widths. For
+a contact it is zero at every width, because a contact has no width.
+
+Peak gouge temperature is what controls thermal weakening and frictional
+melting, so this is not a diagnostic nicety.
+
+```{figure} figures/gouge_thermal.png
+:width: 100%
+
+Left: the gouge's own temperature excess, measured as the departure from
+the straight line through its two edges — the part a contact leaves out.
+Right: measured against predicted, with nothing fitted. $\eta_1$ and
+$\kappa$ are set, $V$ and $\Delta T$ are measured independently.
+```
+
+| $w$ | $V$ measured | $\eta_1V^2/8\kappa$ | $\Delta T$ measured | agreement |
+|---|---|---|---|---|
+| 0.020 | 0.50255 | 0.003157 | 0.004119 | 130.5% |
+| 0.010 | 0.48200 | 0.002904 | 0.002945 | **101.4%** |
+| 0.005 | 0.46351 | 0.002686 | 0.002715 | **101.1%** |
+| 0.002 | 0.42800 | 0.002290 | 0.002300 | **100.4%** |
+
+The $w = 0.02$ row is the thin-band reduction failing where $w/L = 0.033$
+is not small; the exact 1-D quadrature on the measured dissipation
+profile converges the same way (58.1%, 78.6%, 92.3%, 98.1%).
+
+```{warning}
+Do not take the dissipation profile through `uw.function.evaluate`. It
+L2-projects any derivative composite, and viscous dissipation is one, so
+the smoothing lands exactly on the across-band structure being measured.
+Project to a discontinuous degree-0 variable instead — on that space the
+L2 projection is the cell average and nothing crosses a cell boundary.
+
+The first attempt at an oracle here also failed for a physical reason
+worth recording: treating the band as a uniformly heated slab against a
+uniform *background*, and subtracting a far-field $\Phi$, misses by a
+factor of six and gets the sign wrong. There is no background to
+subtract — the wall rock enters through the measured edge temperatures,
+and the source inside the band is the gouge's own dissipation alone.
+```
+
+## A split mesh is not a home for a second equation
+
+The other failure is structural rather than quantitative, and it does not
+involve the mechanics at all. Solve plain conduction between a hot wall
+and a cold wall — uniform diffusivity, no fault properties anywhere — on
+a plain mesh and on the split-node mesh:
+
+```{figure} figures/gouge_insulation.png
+:width: 100%
+
+The same problem, twice. A split-node mesh is split for the **mechanics**;
+heat does not jump across a fault, but the doubled nodes let it.
+```
+
+| | flux through the cold wall | $T$ across the fault line |
+|---|---|---|
+| plain mesh | −0.010000 (analytic $-\kappa\Delta T/L$ = −0.01) | +0.008660, the smooth ramp over the probe offset |
+| split mesh | −0.008091, 19% of the heat stopped | **+0.471241**, a real discontinuity |
+
+The plain row reproduces the analytic flux exactly, so the split row is
+the mesh and not the solver. The fault behaves as a perfect insulator,
+and no choice of fault parameter changes that — there is no thermal
+interface condition on a split node to set.
+
+That is the sharpest form of the argument for the ribbon. Where a second
+equation has to see the fault, the finite-width zone gives it cells, and
+everything is continuous across them.
+
+## Choosing
+
+- **Ribbon** when the width is physics: a gouge with its own rheology, a
+  damage zone another equation has to see, anything where a quantity
+  *inside* the zone is wanted. It costs DOFs as $\mathrm{length}/w$.
+- **Contact** when the fault is genuinely a surface at the scale of the
+  problem and slip on it is the quantity of interest. It costs DOFs as
+  $\mathrm{length}$, with no $w$ in it at all.
+
+## See also
+
+- [Crossing fault zones](crossing-fault-zones.md) — networks, where the
+  ribbon's other advantage lives.
+- [Split-node faults](split-node-faults.md) — the contact representation.
+- [Fault networks](fault-networks.md) — the zero-thickness toolkit.
diff --git a/docs/advanced/index.md b/docs/advanced/index.md
index 6316d432..bc723aae 100644
--- a/docs/advanced/index.md
+++ b/docs/advanced/index.md
@@ -35,6 +35,24 @@ fault-plane yield for modelling fault zones.
 
 **[→ VEP + Transverse Isotropy](vep-transverse-isotropy-faults.md)**
 
+### Fault Zones That Cross
+Finite-width fault zones fused into one region, and what the rheology of
+the overlap is worth.
+
+**[→ Crossing Fault Zones](crossing-fault-zones.md)**
+
+### Gouge Zones
+What a fault zone's width is for, and what collapsing it to a surface
+throws away.
+
+**[→ Gouge Zones](gouge-zones.md)**
+
+### Branching Faults and Junctions
+What a crack model cannot do at a branch point, measured in a shear box
+— and how to choose representations because of it.
+
+**[→ Branching Faults and Junctions](fault-branching-junctions.md)**
+
 ### Custom Meshes
 Create complex geometries with gmsh for research problems.
 
@@ -105,6 +123,9 @@ complex-rheologies
 vep-transverse-isotropy-faults
 split-node-faults
 fault-networks
+crossing-fault-zones
+gouge-zones
+fault-branching-junctions
 fault-mechanics-examples
 custom-meshes
 curved-boundary-conditions