diff --git a/docs/advanced/figures/listric_compare.png b/docs/advanced/figures/listric_compare.png new file mode 100644 index 00000000..2a3e0dd2 Binary files /dev/null and b/docs/advanced/figures/listric_compare.png differ diff --git a/docs/advanced/figures/listric_construction.png b/docs/advanced/figures/listric_construction.png new file mode 100644 index 00000000..37f15751 Binary files /dev/null and b/docs/advanced/figures/listric_construction.png differ diff --git a/docs/advanced/figures/listric_slip.png b/docs/advanced/figures/listric_slip.png new file mode 100644 index 00000000..99007c99 Binary files /dev/null and b/docs/advanced/figures/listric_slip.png differ diff --git a/docs/advanced/index.md b/docs/advanced/index.md index 6316d432..3ab80737 100644 --- a/docs/advanced/index.md +++ b/docs/advanced/index.md @@ -102,6 +102,7 @@ performance multigrid-preconditioning solver-iteration-callbacks complex-rheologies +listric-fault-extension vep-transverse-isotropy-faults split-node-faults fault-networks diff --git a/docs/advanced/listric-fault-extension.md b/docs/advanced/listric-fault-extension.md new file mode 100644 index 00000000..79088739 --- /dev/null +++ b/docs/advanced/listric-fault-extension.md @@ -0,0 +1,246 @@ +# Listric faults over a common detachment + +The classic geometry of a continental rift: listric normal faults — steep +near the surface, flattening with depth — soling into a **common +detachment** in weak rock, the hanging-wall blocks stepping down between +them. + +Underworld3 can build this system two ways, and this page builds it both +ways from the *same traces, the same base mesh and the same drive* so the +two can be compared directly: + +- **finite-width ribbons** ([`place_thin_volume`](conforming-surfaces-and-fault-zones.md)) + — each fault is a thin band of weak material, of a width you choose; +- **zero-thickness faults** ([split-node](split-node-faults.md), assembled + by [`FaultNetwork`](fault-networks.md)) — each fault is a surface across + which the velocity jumps, with no width at all. + +Scripts: `~/+Simulations/listric_extension/`. + +## The geometry + +A listric trace flattening toward the detachment depth $y_d$, then a +straight tail continuing that dip *through* the detachment: + +```python +DETACH_Y, LAM, SOLE_DIP = 0.4, 0.18, 0.066 + +def listric_trace(x0, n=30): + y = np.linspace(1.0 - TOP_GAP, DETACH_Y + SOLE_DIP, n) + x = x0 + LAM * np.log((1.0 - DETACH_Y) / (y - DETACH_Y)) + ... # straight tail across the sole +``` + +Both representations consume this same polyline. Neither reaches the free +surface: 2-D faults cannot daylight in either formulation yet, so the tips +stop a little below the top. + +The drive is plain instantaneous extension — sides pulled at $\pm 0.5$, +free-slip base, free top — and the country rock has viscosity 1 +throughout. + +## Building it both ways + +**Ribbons.** One call places the whole network. The junctions are +resolved in the CAD stage and become ordinary cells of the union, so the +fault runs *into* the detachment as continuous weak material: + +```python +from underworld3.utilities.place_surface import place_thin_volume + +dm, info = place_thin_volume(mesh0.dm, [sole, f1, f2], width=0.02, + label="Zone", label_value=31) +``` + +The width is a **real parameter** — constitutive in this formulation +($V = 2\dot\varepsilon\,w$) — not something that tracks the mesh. The +zone's cells carry $\eta = 10^{-3}$, assigned per cell off the label (P0, +discontinuous: interpolating a nodal viscosity across a sharp interface +smears exactly what the conforming zone resolves). + +**Split nodes.** The split formulation refuses shared vertices at a +junction, so `FaultNetwork` converts each T abutment to the *offset* +form — the junior trace is pulled back and a small damage plug bridges +the gap: + +```python +net = uw.meshing.FaultNetwork(traces(), hierarchy=["Sole", "F1", "F2"]) +net.prepare(h=CELL, ligament=1.0) +net.mesh = base_mesh().add_fault(net.prepared) # same base as the ribbons +... +stokes.constitutive_model.Parameters.yield_stress = net.damage_yield(v, dial=0.05) +net.apply_contact(stokes, conds=ETA_F) +``` + +The interface law is **matched to the band**: a shear band of viscosity +$\eta_{\rm band}$ and width $w$ has zero-thickness limit +$\eta_f = \eta_{\rm band}/w = 10^{-3}/0.02 = 0.05$. + +```{figure} figures/listric_construction.png +:width: 100% + +The F1/detachment junction, at the same scale. **Left**: the ribbon union +is through-going — the fault arrives as weak cells and the detachment +continues through them. **Right**: the split mesh, with the prepared trace +(dark) stopping short of where the raw trace would have met the sole (thin +red), and the damage plug (peach) bridging the gap. +``` + +The pull-back is the thing to watch here. It is quoted as a *ligament* in +multiples of $h$, but that is a clearance measured **across** the senior +fault, so the retreat **along** the junior trace is $\ell h/\sin\theta$ at +an intersection angle $\theta$. A listric fault sole is a glancing +junction by definition — $\theta \approx 20°$ here — so one element of +clearance costs $0.10$ of trace, three times what the same setting costs +at an orthogonal crossing. The ribbon union pays nothing, but its +junction cells are correspondingly slivered (minimum angle $0.9°$; the +solve is untroubled, and steepening the terminal dip is the lever if a +stiffer rheology objects). + +## What the two give + +```{figure} figures/listric_compare.png +:width: 100% + +The same system, both representations. Horizontal velocity (top), vertical +velocity (middle), strain rate (bottom, log scale). +``` + +The **kinematics agree**. Both build the same two hanging-wall blocks with +the same sharp velocity contrast across each fault, and the vertical +velocity — the subsidence pattern, the geologically meaningful signal — +matches in shape and in magnitude ($-0.648$ against $-0.635$ peak). + +The **strain rate does not, and cannot**. This is the clearest statement +of what separates the two formulations: the ribbon spends the +deformation as strain *inside* a band, so the fault is a bright line in +$\dot\varepsilon_{II}$; the split mesh spends it as a *jump*, which no +strain measure can see. The only thing lit up on the split panel is the +junction glue. Neither picture is more correct — they are different +currencies, and a comparison has to be made in a quantity both can pay. + +```{figure} figures/listric_slip.png +:width: 100% + +Slip on each fault, measured in each representation's own way: the +velocity difference across the band edges (ribbon) against the contact +pair jump (split node). The split curves stop short at the detachment — +that is the ligament. +``` + +Peak slip: + +| fault | ribbon | split node | ratio | +|-------|--------|------------|-------| +| F1 | 0.498 | 0.401 | 0.80 | +| F2 | 0.417 | 0.390 | 0.93 | +| Sole | 0.317 | 0.154 | 0.49 | + +On the listric faults the two agree to within 7–20% and the profile +*shapes* are the same. On the detachment the split model slips half as +much. + +## Chasing the detachment deficit + +The obvious suspect is the junction, and it is a large part of the +answer — but not all of it, and the way that came out is worth showing, +because the tempting conclusion is the wrong one. Every knob was turned, +one at a time, reading peak slip on the sole: + +| control | sole slip | effect | +|---------|-----------|--------| +| ribbon reference | 0.317 | — | +| split baseline ($h$ = 0.035, ligament 1$h$, dial 0.05) | 0.154 | — | +| interface law removed ($\eta_f = 0$) | 0.161 | +5% | +| ligament $1h \to 2h$ | 0.130 | −16% | +| refine $h$: 0.035 → 0.0233 → 0.0175 | 0.154 → 0.181 → 0.187 | +21%, converging | +| glue dial 0.05 → 0.01 → 0.002 | 0.181 → 0.204 → 0.210 | +16%, saturating | +| glue plug made *larger* at fixed strength | 0.181 → 0.167 | −8% | +| best case (finest $h$, dial at the ceiling) | **0.216** | 68% of the ribbon | + +Reading down that table: + +- **It is not the law.** Removing the interface condition entirely moves + slip by 5%. Both representations are already weak, as intended. +- **The junction gap is the biggest single lever.** Doubling the ligament + costs 16%; refining the mesh — which shrinks the gap, since the + ligament is measured in $h$ — buys 21% back. This is the practical + rule: buy fidelity with elements, not with physical size. +- **But closing the gap does not close the difference.** The refinement + sequence is converging to ≈0.19, not to the ribbon's 0.317, and the + ribbon itself is essentially $h$-independent over the same range + (0.3173 → 0.3155) — as it should be, since its width is a physical + parameter, which makes it a clean control. +- **The glue helps, and then stops helping.** Weakening the plug raises + transfer to a ceiling near 0.21 — the inviscid-plug limit the network + documentation describes. Making the plug *bigger* at the same strength + makes things *worse*: a larger weak region lets the corner deform + instead of passing slip along the sole. A junction plug is a hinge, not + a bridge. + +Two explanations we tested and discarded, in case they occur to you too. +The ribbon's band, being weak and sub-horizontal in horizontal extension, +might simply be stretching along its own length, with our probe reading +that as slip — but shrinking $w$ by 4$\times$ at fixed +$\eta_f = \eta_{\rm band}/w$ moves the sole by only 4.5% (0.317 → 0.303), +so it is not that. (That the ribbon is nearly $w$-independent at fixed +$\eta_f$ is a good result on its own: it says $\eta_f$ really is the +similarity parameter linking the two formulations.) And the plug does +shrink away under refinement — but as the table shows, a bigger plug +transmits less, not more. + +So a residual of roughly 30% survives every control we have. The best the +offset form achieves here is 68% of the through-going detachment slip, +and no single parameter closes the rest. + +## What that means for choosing + +The two formulations agree on what a fault *is* — the block kinematics, +the subsidence, the slip on the individual listric faults — and disagree +on what a *junction* is, in a way that is structural rather than +mistuned. A through-going weak band and a network of zero-thickness +contacts joined at offset junctions are not the same mechanical object +where faults meet, and on this geometry the sole — the one fault that has +to *receive* slip through two junctions — is where the disagreement +shows. The listric faults themselves reach 86–97% of the ribbon's slip at +the glue ceiling; the sole reaches 68%. + +- **Ribbons** when the width is physics — a gouge zone with its own + rheology, a damage zone another equation has to see — and, on this + evidence, when the system is **junction-dominated**: a detachment fed + by splays is exactly that, and the union carries the connectivity + natively. +- **Split nodes** when slip on individual faults is the quantity of + interest — earthquake cycles, stress transfer, frictional and + rate-and-state laws — and when conditioning matters: there is no thin + feature and no viscosity contrast for the Schur complement to fight. + +```{warning} +The residual is not explained. If you are modelling a detachment system +in the offset form, calibrate the transfer rather than assuming it: the +numbers above are for one geometry at one intersection angle, and the +glancing junction is the hard case. +``` + +## Which to use + +- **Ribbons** when the width is physics — a gouge zone with its own + rheology, a damage zone that another equation (porous flow, thermal) + has to see, or a network whose junctions must be through-going without + further thought. +- **Split nodes** when slip is the quantity of interest — earthquake + cycles, stress transfer, frictional and rate-and-state laws — and when + conditioning matters: there is no thin feature and no viscosity + contrast for the Schur complement to fight. + +Junctions are the place to think hardest. The ribbon's are free but +slivered; the split's are clean but offset, and on a glancing +intersection that offset is expensive and only partly recoverable. + +```{note} +The split solves here report an inner velocity-block KSP reaching its +iteration cap. The answers were checked against the direct-LU route +(`solver._rotated_use_lu`) and agree to 0.2%, and are unchanged by a +1000x tighter tolerance — the outer Krylov is flexible, so the warning +reflects preconditioner quality on this mesh, not the solution. +```