Thin viscous sheet (plane stress) solver — with a published analytic benchmark already available

[lmoresi]

What

A thin-viscous-sheet solver: the vertically averaged, plane-stress formulation used for continental deformation (England & McKenzie 1982, 1983; Houseman & England 1986; England & Houseman 1986 for India-Asia).

Why now

The analytic benchmark for it already exists, in the same appendix we have just implemented the plane-strain solution from — Barr & Houseman 1996 (GJI 125, 473-490), equations (A10)-(A13). So this is the rare case of a new solver that can arrive with an exact solution to check it against, including one containing a fault.

The plane-strain half is PR #550. The plane-stress half was deliberately left out there, because it is a different equation set and would have been misleading to present as a benchmark for incompressible Stokes.

The formulation

From the paper (their eqs 5, 6). Substituting p = sigma_zz - tau_zz into the force balance and averaging through the layer thickness:

d(tau_ij)/dx_j - d(tau_zz)/dx_i = -d(sigma_zz)/dx_i        (5)

with sigma_ij, tau_ij now the vertical averages, and the vertical deviatoric stress fixed by incompressibility plus the constitutive relation:

tau_zz = -(tau_xx + tau_yy)                                 (6)

Vertical strain is therefore permitted — the in-plane divergence is NOT zero, which is exactly what distinguishes this from the current Stokes solver. For instantaneous solutions sigma_zz is taken constant; for evolving ones it carries the crustal thickness, which is where the geodynamics lives (thickening, gravitational potential energy, the Argand number).

The benchmark that comes with it

Same structure as the plane-strain case in #550, with three changes:

• f_m = (3m + 8) / (3m - 2) (A10), replacing (m+2)/m
• pressure recovered from the velocity, p = 2 eta [ du_r/dr + (1/r)(du_theta/dtheta + u_r) ] (A11) — note this is 2 eta times the in-plane divergence, which is non-zero here
• fault constraint a_(m+2) = -3 m b_m / (3m - 2) (A12), replacing (A7c)

Test problem (A12a, A12b) and its exact solution (A13a, A13b, A13c), with the same normalisation: relative slip on the fault at r = R0, theta = 0 is 2 U0. Their Table A1 gives finite-element convergence for it — 1.35% / 0.68% / 0.34% at 662 / 2342 / 8774 nodes, going as N^(-1/2).

Transcription of (A13a) is confirmed by the slip anchor: the half-integer terms 7 cos(t/2) + 9 cos(3t/2) give +-16, times U0/16, so 2 U0. The same sign question flagged in #550 applies to (A12b) vs (A13b) — resolve it the same way, by requiring the field to satisfy the plane-stress balance rather than by reading the scan.

Suggested sequence

1. The solver — equations (5) and (6), instantaneous first with sigma_zz constant.
2. The analytic solution (A13) into uw.analytic, verified symbolically the same way Barr & Houseman faulted-medium analytic solution in uw.analytic #550 was: check it satisfies (5) and (6) and the fault conditions, rather than comparing to a stored answer.
3. Connect them, once faults can reach a boundary (2-D fault zones cannot reach a boundary — the outcrop machinery exists but is never wired into the ribbon path #549 — the analytic domain has the fault running from the tip at the centre out to the perimeter).

Steps 1 and 2 are independent of #549 and of each other.

Underworld development team with AI support from Claude Code