An H-κ stack hands you two numbers and a picture of a peak. The numbers are usually fine. The picture is where the trouble starts, because that peak is the maximum of a stacking amplitude, not the maximum of a likelihood, and the spread around it carries no probabilistic meaning whatsoever. Reading an error bar off that spread is a habit, not an inference.
The alternative is to train a density estimator on simulations and let it return the entire posterior over crustal thickness H and bulk Vp/Vs ratio κ. What comes back is not a tighter number. It is a differently shaped object, and the shape is the part that changes decisions.
The peak is a score, not a probability
H-κ stacking [1] sweeps a grid of candidate values, sums receiver-function amplitude along the predicted arrival times of the direct Ps conversion and its crustal multiples, and reports where the sum is largest. Fast, standard, and perfectly defensible as an estimator.
What it is not is a probability statement. The quantity being maximised is a stacking amplitude. Nothing in the construction ties that amplitude to the probability of the observed data under a crustal model, so the width of the peak is not an uncertainty, and two stations with equally broad peaks are not equally uncertain in any sense you can carry into a siting decision. Noise gets handled by the same route, which is to say by hand.
The answer has a shape, and the shape is a ridge
Simulation-based inference [2] is built for the case where simulating data is easy and writing down a likelihood is not. Receiver functions qualify. RAYSUM [3] will produce a synthetic for any layered, dipping, anisotropic model you hand it, and nobody can give you the matching likelihood in closed form. Neural posterior estimation [4] trains a conditional density estimator, here a normalizing flow [5], on pairs of subsurface model and simulated receiver function. At inference it takes an observed receiver function and returns q(H, κ | d) directly.
The first thing that falls out is geometry. H and κ are not separately pinned down by a Ps delay. The delay fixes a combination of the two [1], so the models compatible with an observed Ps time lie along a curve in the (H, κ) plane, and the posterior is a curved ridge rather than a compact cloud.
That has a consequence worth stating precisely. A point estimate plus two marginal error bars, one on H and one on κ, describes a rectangle. The ridge runs diagonally across that rectangle, so most of the box holds models the data already rejected, and its corners are exactly where a screening call goes wrong without anyone noticing. The joint region is the object. The marginals are shadows of it.
What breaks the trade-off is a different arrival
Once the credible region is drawn as a ridge, the obvious question is what shortens it, and the classical method already answers. The Ps delay constrains one combination of H and κ. The crustal multiples arrive on a different combination, so bringing them in cuts across the ridge instead of sliding along it [1].
More events of the same kind move you up and down the ridge and leave its length alone. A second arrival is what collapses it. That is a survey-design conclusion rather than a statistical one, and a point estimate cannot express it, because the information sits in the orientation of a region that a point estimate never draws.
Calibration is what turns an interval into a claim
A flow will always return a posterior. It will return one just as confidently for data unlike anything it was trained on, which is the failure mode that makes practitioners rightly suspicious of learned uncertainty.
Simulation-based calibration [6] is the check. Over many simulated observations, an interval quoted at credibility α should contain the true value α of the time. Run the test and the intervals either mean what they say or they do not.
This is a different discipline from reading a spread off a stack. An interval becomes a testable claim, and the test costs little once a simulator exists. An estimator that fails coverage is not a slightly worse estimator. It is one whose numbers should not be quoted at all.
One network, every station on the line
The training cost is paid once. After that, a station's posterior is a single forward pass, which is what makes the profile-scale version practical. The source figure behind this note is a two-dimensional P-receiver-function profile through a craton, from Dr Tannistha Maiti's 2018 PhD thesis on receiver-function imaging of the Moho and the lithosphere-asthenosphere boundary. Turning a transect like that into per-station posteriors becomes a batch of forward passes rather than a queue of independent samplers.
On clean data the maximum of the learned posterior lands where the H-κ peak lands. The method does not overturn the classical result; it reproduces it and attaches a calibrated distribution to it. That is the only honest argument for putting it into a workflow that already functions.
Limitations
- The posterior is only as good as the simulator and the prior behind it. Structure the simulator cannot represent, such as dipping interfaces, anisotropy or sediment layers, biases the answer, and a coverage test run against the same simulator will not catch it.
- Normalizing flows can be confidently wrong out of distribution. Calibration has to be measured, never assumed.
- The parameterisation here is small: H, κ, and a few nuisance layer parameters. Richer models raise the simulation budget rather than coming for free.
- Summary statistics matter. Compressing a receiver function before the estimator sees it discards information, and the loss reappears as an inflated posterior.
- The instrument is a teaching schematic. The ridge follows the standard Ps moveout relation between H and κ, but the region widths and their growth with picking noise illustrate the mechanism and are not the output of a trained estimator.
By the numbers
Parameters the classical stack searches
Inference cost per station, once trained
Shape of the joint credible region
Schematic anchor in the source figure
Key takeaways
- An H-κ stack maximises a stacking amplitude, not a likelihood, so the spread around its peak is not an uncertainty and should not be quoted as one.
- Simulation-based inference trains a conditional density estimator on simulated receiver functions and returns a full posterior over crustal thickness and Vp/Vs from one observation.
- H and κ are jointly constrained by the Ps delay, so the credible region is a curved ridge; a point plus two marginal error bars describes a box that mostly holds rejected models.
- The ridge shortens when a different arrival is added, not when more of the same events are stacked, which makes the region orientation a survey-design signal.
- Coverage testing is what turns a quoted interval into a testable claim, and an estimator that fails it should not have its numbers used.
References
[1] L. Zhu, H. Kanamori. Moho depth variation in southern California from teleseismic receiver functions. J. Geophys. Res., 2000. doi:10.1029/1999JB900322
[2] K. Cranmer, J. Brehmer, G. Louppe. The frontier of simulation-based inference. PNAS, 2020. arXiv:1911.01429
[3] A. W. Frederiksen, M. G. Bostock. Modelling teleseismic waves in dipping anisotropic structures. Geophys. J. Int., 2000. doi:10.1046/j.1365-246X.2000.00090.x
[4] G. Papamakarios, I. Murray. Fast ε-free Inference of Simulation Models with Bayesian Conditional Density Estimation. NeurIPS, 2016. arXiv:1605.06376
[5] G. Papamakarios, E. Nalisnick, D. J. Rezende, et al. Normalizing Flows for Probabilistic Modeling and Inference. JMLR, 2021. arXiv:1912.02762
[6] S. Talts, M. Betancourt, D. Simpson, et al. Validating Bayesian Inference Algorithms with Simulation-Based Calibration. 2018. arXiv:1804.06788
[7] Source: 2018 PhD thesis by Dr Tannistha Maiti, receiver-function imaging of the Moho and LAB, Fig 5.7 (two-dimensional P-receiver-function profile through a craton). The posterior geometry, credible-interval widths and H and κ readouts in the instrument are schematic illustrations of amortized inference, not outputs of a trained estimator.
