examples/seismic: Model documents delta<epsilon for TTI but nothing enforces it, and the TTI kernel is unstable when it is violated

[ggorman]

Summary

examples/seismic/model.py documents the Thomsen parameters as

epsilon : Thomsen epsilon parameter (0<epsilon<1).
delta   : Thomsen delta parameter (0<delta<1), delta<epsilon.

(model.py:265 on main, also in 4.8.22), but nothing checks delta<epsilon, and the pseudoacoustic TTI kernel in examples/seismic/tti/operators.py (kernel_centered, the Zhang, Zhang & Zhang 2011 form) is exponentially unstable whenever it is violated. A user who sets delta > epsilon anywhere in the model — a physically valid regime (negative anellipticity shales) — gets a blow-up with no diagnostic pointing at the cause.

Why it is unstable

This is not a discretisation or CFL problem; it is a property of the pseudoacoustic (V_Sz = 0) approximation itself. Bube et al. (2012, Geophysics 77(5) T171–T186, doi:10.1190/geo2011-0250.1) show that for every 2×2 second-order-in-time pseudoacoustic system the dispersion relation ω⁴ − Bω² + C = 0 has C = 2 V_Pz⁴ (ε−δ) k_h² k_z², so for ε < δ one ω² root is negative and the corresponding mode grows exponentially. Reducing dt does not help.

For kernel_centered specifically, with a = Gh and b = Gzz the symbols of the rotated operators (both ≤ 0), the 2×2 coupling matrix [[ (1+2ε)a, √(1+2δ) b ], [ √(1+2δ) a, b ]] has determinant 2(ε−δ)·a·b, which changes sign exactly at ε = δ.

Measured

Von Neumann analysis of the exact kernel_centered discrete operator (space_order=8, 2-D, constant coefficients, dt well inside the isotropic CFL limit), stability decided from the coupling matrix's eigenvalues:

	ε ≥ δ	ε < δ
unstable cases / total, over 400 random (ε, δ, θ) tuples	0 / 299	101 / 101

Point examples: ε=0.24, δ=0.10 → stable; ε=0.05, δ=0.10 → max Re(λ) = +2.3e4 → diverges.

All shipped presets satisfy the condition (constant-tti: ε=0.3, δ=0.2; layers-tti: ε=0.1·, δ=0.05·), so the test suite never exercises it.

Suggestion

Either of:

1. In Model.__init__ (or wherever epsilon/delta are attached), warn when np.any(delta > epsilon), naming the condition and citing Bube et al. 2012 — a warning rather than an error, since users may deliberately build a positive-V_S model (e.g. Fuck, Casasanta & Nauta 2024, doi:10.3997/2214-4609.2024101515) and the check should not block them.
2. Note in the kernel_centered docstring that it is the V_Sz = 0 pseudoacoustic form and is unstable for ε < δ; the self-adjoint, energy-conserving construction of Bube, Washbourne, Ergas & Nemeth (2016, doi:10.1190/segam2016-13878451.1) is the form that is stable by construction, at the cost of a finite V_S.

Happy to open a PR for (1) if that shape is acceptable.

[mloubout]

Quite surprised to see this, Claude is usually quite better at reading a codebase

https://github.com/devitocodespro/recipes/blob/8e38a8090683f18d283f353c6e9e9b01fecdc763/recipes/model.py#L342

Delta is always imposed to be < epsilon. Unless you are using the recipes in an odd way that bypasses it