What should the rheology be where two fault zones overlap at a high angle?

[lmoresi]

Rescoped 2026-08-16 after review. This was originally filed as three defects; two of them do not survive scrutiny and are withdrawn below. What remains is one open modelling question.

The question

Where two fault zones overlap, the fused zone carries a single label by design — that ignorance is what makes fuse (#541) remove junction detection, and it is deliberate. So the rheology at an overlap comes entirely from the Surface lookup, which today gives each node the director of the nearest surface.

At a low-angle crossing or a tangential merge the two directors nearly agree, so the choice does not matter. At a high-angle crossing the overlap is weak in two directions, and a single director represents only one of them.

The open question is whether that matters in practice. The ownership boundary is a measure-zero surface with cells either side taking opposite directors, so the net response of the overlap region may already approximate the orthotropic answer without anyone constructing it. That is plausible and unmeasured.

What would settle it

A high-angle crossing — the tangential geometries used so far cannot discriminate, since both directors agree there by construction. Compare, on one mesh:

• nearest-fault director (current behaviour)
• an explicit orthotropic treatment, e.g. summing the two weak planes' compliance tensors and inverting — which represents both weaknesses, is continuous, and degenerates correctly where only one surface is near

against the gmsh-union control. If they agree, the current rule is fine and this issue closes. If not, the gap is the error being carried today.

Note the connection to #539: an additive-compliance treatment needs the TI compliance tensor and the eta -> 1/(4 eta) structure that issue records. The two are related if this route is taken.

Withdrawn from the original report

• "Nearest vertex, not nearest surface." Real as stated — a densely-sampled trace annexes territory from a sparsely-sampled one — but this is a modelling-discipline requirement, not a defect: a fault set that builds one model should be consistently sampled. Generalising the KDTree to exact distance-to-segment would make a coarse triangulation work well enough to keep, which is the wrong service. Refine or retriangulate the surface instead, which fixes the distance query and everything else that reads it.
• "The director flips discontinuously across the partition boundary." Poor wording on my part: "partition" here meant the geometric medial axis between the two traces, not the MPI decomposition. Nothing is decomposition-dependent — the KDTree is built from replicated Surface vertices and queried with rank-local coordinates, so a given point gets the same answer at any rank count. The parallel contract is intact.

Underworld development team with AI support from Claude Code

[lmoresi]

Expanding on why this sits on the critical path for ribbon fault NETWORKS, and a possible route that connects it to #539.

fuse moved this problem rather than removing it

The point of fusing overlapping zones (#541) is that the zone mesh is deliberately ignorant of which branch a cell came from — that ignorance is what makes the connection region self-selecting and removes the geometric special-casing at junctions.

The consequence is that all branch information now has to come from the Surface lookup. Before fuse that lookup was incidental; after it, it carries the whole representation at a junction. So its defects matter more than they did.

Three problems, increasing in seriousness

1. Nearest VERTEX, not nearest SURFACE (the original report). Reproducibility bug with a clear fix: query segment distance per surface and take the minimum, rather than pooling all surfaces' vertices into one tree. Surface already has the polyline distance routine for the single-surface case.

2. The partition boundary carries no physics. Even with (1) fixed, ownership flips across the medial axis between traces and the director flips discontinuously with it. Benign for a near-tangential merge, where both directors nearly agree — which is why the listric campaign never tripped on it. Not benign at a high-angle junction.

3. At a genuine crossing no single director is correct. Two faults at a high angle leave the overlap weak in TWO directions. That is orthotropic, and picking a nearest fault discards one weakness entirely rather than partially.

A route for (3): compliances add

Rather than choosing an owner, let overlapping weak planes contribute their compliance tensors additively and take the effective stiffness as the inverse of the sum. That:

• represents both weaknesses instead of discarding one;
• is continuous across the overlap, so there is no medial-axis flip and no discontinuous director;
• degenerates correctly to a single TI fault wherever only one surface is near;
• needs no notion of ownership at all, which suits a fused zone that was designed not to know.

This is what #539 is asking for from the other end. That issue wants the TI compliance tensor and records that its inverse has the same structure with eta -> 1/(4 eta). That map is exactly what an additive-compliance combination needs. #539 and #544 look like one problem seen from two ends, and neither issue currently says so.

What would settle whether it matters

A high-angle crossing — not a tangential merge — comparing nearest-fault-director against additive-compliance, both against the gmsh-union control on the same mesh. The tangential geometries used so far cannot discriminate, because both directors agree there by construction.

Underworld development team with AI support from Claude Code

[lmoresi]

Rescoped in the description above (2026-08-16). Recording what changed so the earlier comments are not left standing unqualified:

The original report and my expansion comment both over-claimed. Two of the three points are withdrawn — the sampling-density complaint is a modelling requirement rather than a defect, and the "partition boundary" wording wrongly suggested MPI decomposition dependence when the medial axis is purely geometric and the lookup is decomposition-independent.

What remains is one unmeasured question: whether the nearest-fault director already approximates the orthotropic response of a high-angle overlap, given that cells either side of the ownership boundary take opposite directors. The tangential geometries the campaign has run cannot discriminate, so this needs a high-angle crossing specifically.

[lmoresi]

Draft comment for #544 — fill the table from sweep.log before posting

Measured. The rule is not a good approximation, and the geometry to worry about is the opposite of the one this issue names.

Two straight ribbons crossing at the centre of a unit box, fused into one
zone, driven in horizontal extension. Four rheologies on the SAME mesh at
each crossing angle, differing only in the overlap cells, against a
no-fault control on that same mesh. Scripts:
~/+Simulations/ribbon_network_2d/.

The algebra says the acceptance test as written would have measured nothing

In plane strain the deviatoric strain-rate space is two-dimensional,
and rotating the material by θ rotates a state in it by 2θ. A weak plane
is compliant in exactly one mode, so it occupies a single line in
that space, at 2 phi + 90 for a normal at phi. Two planes separated
by Delta in the rock are separated by 2 Delta there.

So Delta = 0 and Delta = 90 are the same case: both compliances land
on the same line. Shear on a plane with normal x̂ is shear on a plane
with normal ŷ — the same tau_xy. Under the compliance-average rule the
overlap at 90 degrees is identically the single plane, so the
nearest-fault rule and the orthotropic one coincide there by
construction. The 90-degree X this issue proposes as the test is the
one geometry where the question cannot be asked. The tensors differ
most at Delta = 45, where the two compliances are orthogonal and the
overlap is isotropic.

Checked numerically against the shipped Mühlhaus–Moresi tensor
(compliance_algebra.py). Two by-products:

• TI: expose the compliance (inverse) tensor — same structure, reciprocal scale #539 confirmed. The eta -> 1/(4 eta) substitution reproduces the
  compliance to round-off, and the director round-trips through the
  tensor exactly (mod 90 degrees, which is the 2-D degeneracy above).
• The orthotropic treatment is free in 2-D. The additive-compliance
  tensor has two eigenvalues and one principal direction, which is
  exactly the structure of TransverseIsotropicFlowModel. No new
  constitutive machinery is needed — read off an effective director and
  effective eta_0, eta_1.

"Sum the compliances" needs pinning down

Taken literally, summing makes two coincident weak planes twice as
weak as one, and the excess does not switch off away from coincidence: at
Delta = 15 the sum has already made the off-plane mode 14x more
compliant than the wall rock. That is a modelling choice (two independent
fabrics superposed), not a consequence. The compliance average (Reuss;
the overlap as a 50/50 mixture) degenerates to the single plane exactly.
Both were run.

The measurement

Metric: dissipation weakening W = 1 - D/D_none, referenced to the
no-fault control at the same angle, with no probe geometry in it. (The
control dissipates exactly 4.000000, the analytic value for this box, so
the integral is checked too.)

Delta	W nearest	W orthotropic	W isotropic	nearest, as % of orthotropic
10	1.198%	1.710%	3.628%	70.0%
20	3.687%	4.520%	6.256%	81.6%
30	6.618%	7.679%	9.240%	86.2%
45	11.566%	12.435%	13.646%	93.0%
60	15.789%	16.253%	17.110%	97.1%
75	18.322%	18.476%	19.121%	99.2%
90	18.918%	18.918%	19.465%	100.0%

(orthotropic = compliance average; the sum column runs 61.8 / 78.5 / 84.7
/ 92.3 / 96.7 / 98.9 / 99.7, the small residual at 90 being the sum rule
double-counting one mode in 24 cells.)

The 90-degree row is the negative control and it passes exactly: the
two dissipations agree to seven significant figures, 3.243286 against
3.243286. Nothing was tuned to make that happen — the algebra above
requires it — and it is what licenses the other six rows.

The error is largest where this issue expected none

The per-cell tensors differ most at 45 degrees, and the measured curve
shows no feature there at all. Both are right, and the resolution is
geometric: the overlap AREA goes as w^2 / sin(Delta), so there are 141
overlap cells at 10 degrees against 24 at 90. The tensor error peaks at
45, the area falls monotonically, and the product is dominated by the
area.

What this means

The nearest-fault rule is not a good approximation to an orthotropic
overlap, and the case to worry about is the shallow crossing. At a
10-degree crossing it loses about a quarter of the network's weakening;
the error closes monotonically as the crossing squares up and is exactly
zero at 90 degrees.

The table above is computed at a zone mesh of w/4 and its percentages
are BOUNDS: refining the zone mesh always moves the rule closer to the
orthotropic one. At Delta = 10 the ratio runs 70.0 -> 73.8 -> 75.1% as
the zone mesh halves twice (overlap 141 -> 522 -> 2051 cells), with the
increments contracting threefold, so the converged value is near 76%. At
Delta = 45 it runs 88.2 -> 93.0 -> 94.5%, converging near 95%. The HOST
mesh does not enter at all — three host sizes give the same ratio to four
significant figures, which is the embedding working as designed.

That is the opposite of the geometry a designer would think to check, and
it is the opposite of what this issue proposed as the acceptance test.

There are two defensible ways to close it, and this issue should pick one
rather than stay open:

1. Document the regime and keep the rule. The error is bounded by
   about 30% and only reaches that at crossings so shallow that the two
   traces are nearly the same fault. A model with a genuine X at 60-90
   degrees is inside a few percent today.
2. Implement the orthotropic overlap. In 2-D this is nearly free —
   the tensor is a TI tensor, so the shipped constitutive model already
   expresses it. What is missing is in the Surface lookup, which would
   have to return the two nearest surfaces and their distances rather
   than only the nearest. The 3-D case is not free and is not answered
   here.

My recommendation is (1) plus this page in the docs, because the regime
where the rule fails is also the regime where "two fault zones" is a
questionable description of the geometry — but the choice is a modelling
call, not a technical one.

One thing this issue asks for that cannot be delivered as stated

"Summing the two weak planes' compliance tensors" makes two coincident
weak planes twice as weak as one, and the excess does not switch off away
from coincidence: at Delta = 15 the sum has already made the off-plane
mode 14x more compliant than the wall rock. If the additive route is
taken, the weighting has to be decided explicitly — the sum and the
average are the two ends of a one-parameter family and they differ by a
factor of two in the overlap's weakness.

Caveats, stated

• 2-D only. The 90-degree degeneracy is a property of the two-dimensional
  deviatoric space and does not carry to 3-D, where two weak planes
  at any angle occupy different subspaces of a five-dimensional deviatoric
  space. The 3-D question is open and is not answered here.
• The overlap is a few tens of cells; crossing_resolution.py varies the
  zone mesh size and the host mesh size to show the answer has stopped
  moving.
• The angle sweep necessarily changes the loading as well as the crossing
  (the pair is symmetric about the loading axes), so absolute slip is not
  comparable across angles. Every ratio quoted is between runs at the
  same angle on the same mesh.

Underworld development team with AI support from Claude Code

[lmoresi]

Final recommendation, after the shear-box test

The extension-frame measurement above stands, but a further test (@lmoresi's design: a shear box with spin in the drive, holding an asymmetric Y — sub-horizontal main fault plus a horizontal branch fusing into it) shows it was answering a smaller question than it appeared to. In that frame the three overlap treatments — nearest director, harmonic average, isotropic overlap cells — are indistinguishable from each other while all losing factors (77% / 46% / 19% of the isotropic reference at junction angles 20/30/40°) wherever the network hands slip between differently oriented members. The overlap tensor was never the dominant term; the arms' TI constraint is. Mechanism: a slip gradient along a fault is a volumetric exchange between band and surroundings, and that is the mode the weak-plane constraint removes — the same reason the split-node representation needs junction glue.

So the resolution of this issue as filed:

1. The nearest-director rule needs no replacement. Its error relative to any orthotropic construction is real but subordinate (≤ 25% of the network weakening at the worst 2-D crossing angle, resolution-bounded), symmetric alternatives are one small change if ever wanted (the harmonic average is exactly a TI tensor in 2-D), and no overlap rule addresses the term that matters.
2. What matters is documented instead: docs/advanced/crossing-fault-zones.md (the extension measurement, with its scope stated) and docs/advanced/fault-branching-junctions.md (the shear-box study: the slip-exchanger junction, the exit-throat law γ ≈ (Σ slip)/w at r ≈ 2.1 w, the corner-singularity argument that junction width is a process-zone size, the matrix-stress pattern flip, and the damage experiment with its negative control). The choosing rules there — surfaces for distinct well-oriented segments, isotropic zones wherever slip is exchanged, joins made kinematically compatible with damage zones around them — are the practical answer to the question this issue asked.

Recommend closing on that basis; the 3-D overlap question noted earlier remains open and is a different (and smaller) matter than the branching mechanics documented here.

Underworld development team with AI support from Claude Code