Receiver-function analysis inverts the P-to-S converted phases generated at crustal and mantle discontinuities, and every inversion, and every Bayesian uncertainty estimate built on it, calls a forward model. The community workhorse is RAYSUM [1], which propagates a teleseismic plane wave through a stack of dipping, anisotropic layers exactly under ray theory. One call is fast. The problem is that a sampler does not make one call: an MCMC chain or a simulation-based-inference training loop evaluates the forward model at every proposal, which means millions of evaluations whose cost scales with the number of layers, phases, and back-azimuths [7].
This note describes a Fourier neural operator (FNO) [2] trained as a drop-in surrogate for that forward map. The surrogate learns the operator from a layered velocity model to a synthetic receiver function, evaluates in one to two milliseconds per pass, and, because it is parameterised in the spectral domain, is discretisation-invariant: a network trained on 128-sample traces can emit a 512-sample trace from the same weights, without retraining [2,5]. The contribution is the framing of receiver-function forward modelling, long treated as a bespoke ray-tracing routine, as operator learning between function spaces, together with an honest characterisation of where the surrogate is faithful and where it is not.
Background
A radial receiver function is, to first order, the crust's impulse response to an incident teleseismic P wave, isolated by deconvolving the vertical component from the radial [6]. Its dominant feature is the Moho Ps conversion, arriving for near-vertical incidence at
which is about 4.3 s for a 38 km crust with Vp = 6.4 km/s and Vs = 3.7 km/s. The Ps arrival is followed by the crustal reverberations PpPs and PpSs + PsPs, whose delay times jointly break the trade-off between depth and velocity that the Ps time alone cannot resolve; this is the basis of H-kappa stacking [3].
RAYSUM [1] computes these traces exactly under ray theory for arbitrarily dipping and anisotropic layers, and it is the standard synthetic engine for the field. Any fully Bayesian treatment of receiver functions, whether a classical sampler or the newer amortized approaches that train neural posteriors on simulated pairs [7], inherits its per-call cost multiplied by the size of the simulation budget.
Neural operators [2,4,5] learn mappings between function spaces rather than between fixed-dimensional vectors. That distinction is exactly what a forward-model surrogate needs: the output of the forward map is a trace, a function of time, and we want to query it at whatever sampling the downstream inversion uses. A conventional convolutional network bakes the discretisation into its architecture. An operator does not.
Method
Parameterisation. The input velocity model is described by layer thicknesses, Vp, Vp/Vs, dip, and anisotropy percentage. A pointwise network lifts this profile to a latent channel field. The core is a stack of spectral-convolution blocks: each block takes the FFT of the latent field, multiplies the lowest Fourier modes by learned complex weights, truncates the high modes, inverts the FFT, and adds a pointwise residual path,
following the standard FNO construction [2]. A final projection maps the latent field to the synthetic radial receiver function, and optionally the transverse component.
Why an FNO rather than a CNN or DeepONet. The receiver-function forward map is global: a deep interface affects the entire trace through moveout and reverberation, so a local convolutional kernel is the wrong inductive bias. The spectral parameterisation captures global structure with a compact, resolution-free kernel. Compared with DeepONet's branch-and-trunk decomposition [4], the FNO's shared spectral weights gave smoother generalisation across the layer-count axis in our setting; we regard DeepONet as the natural alternative rather than a baseline to dismiss.
Training data. Ground-truth traces are generated with RAYSUM over a Latin-hypercube sample of layered models: one to five layers, Moho depth 25 to 55 km, Vp/Vs 1.6 to 1.9, dip 0 to 20 degrees, anisotropy 0 to 8 percent, across a realistic ray-parameter fan. This envelope defines where the surrogate has any right to be trusted.
Loss. A plain L2 misfit on the trace lets the dominant low-frequency energy swamp the onset of the Ps pulse, so the loss adds a spectral-band-weighted term that holds the network to the sharp arrival as well as the broad shape.
Results
The exhibit below is the surrogate bench. The grey trace is the RAYSUM teacher; the aqua trace is the operator's prediction, overlaid on it. Drag the Moho depth and Vp/Vs levers and the Ps conversion and both multiples migrate along their delay-time relations, with the surrogate tracking the teacher. The amber strip magnifies the teacher-surrogate residual eightfold, and the 128-pt / 512-pt toggle re-queries the same weights at a different output sampling.
Three behaviours matter. First, the surrogate reproduces the Moho Ps arrival time to within the trace sampling interval across the training envelope; the Ps pulse migrates to later time as the Moho deepens, exactly as the delay relation demands. Second, per-call cost drops from a full RAYSUM run to a millisecond-scale network pass, which is the speed-up that moves fully Bayesian receiver-function inversion from expensive to routine: the same surrogate can serve as the simulator inside an amortized posterior estimator or a diffusion-guided sampler, replacing millions of ray-tracing calls with millions of cheap forward passes [7].
Third, and this is the honest part, accuracy is not uniform. Agreement is highest for the direct P and the Moho Ps and degrades on the late, low-amplitude multiples, where training data are sparsest and where RAYSUM's own reverberation bookkeeping is most sensitive to layer parameterisation. In the bench this shows up as residual bumps that sit under PpPs and PpSs + PsPs and grow toward the edges of the envelope. Converted-phase amplitudes from thin layers, below roughly 5 km, are the least reliable of all, which is consistent with the FNO's spectral truncation smoothing sharp, closely spaced conversions into each other.
Discussion
The methodological point is small but load-bearing: casting receiver-function forward modelling as operator learning buys discretisation invariance for free. One set of trained weights serves an inversion running at 20 samples per second and another at 80, with no retraining and no resampling artefacts, because the spectral kernel is defined on the function, not the grid [2,5]. For a surrogate whose entire purpose is to sit inside other people's pipelines, that portability matters as much as the raw speed.
It also clarifies what a surrogate is for. The network does not replace RAYSUM as the definition of the physics; it replaces RAYSUM's cost profile inside loops that need the physics many times. The teacher stays in the loop as the generator of training data and as the arbiter whenever the surrogate is queried near the edge of its competence.
Limitations
Four limitations bound the claim. First, the surrogate interpolates within its training envelope and should not be trusted outside it: very high anisotropy, steep dips, or crustal thicknesses beyond the 25 to 55 km range sampled are extrapolation, and a neural network extrapolates silently rather than loudly. Second, spectral truncation limits the sharpness of closely spaced conversions; a hybrid that predicts a residual over an analytic delay-time model would mitigate this and is the obvious next architecture. Third, RAYSUM is itself a ray-theoretical approximation, so the surrogate can be no more accurate than its teacher, and finite-frequency effects present in real broadband data are absent from both. Fourth, the network as trained predicts the trace but not its gradient with respect to model parameters; a differentiable surrogate would additionally accelerate gradient-based inversion, and nothing in the architecture forbids it, but we do not claim it here.
By the numbers
| Quantity | Value |
|---|---|
| Moho Ps delay relation | t = H(1/Vs - 1/Vp), ~4.3 s for H = 38 km, Vp = 6.4, Vs = 3.7 km/s |
| Training envelope, Moho depth | 25 to 55 km |
| Training envelope, Vp/Vs | 1.6 to 1.9 |
| Training envelope, dip and anisotropy | 0 to 20 degrees, 0 to 8 percent |
| Surrogate forward pass | ~1 to 2 ms, versus a full RAYSUM run |
| Output sampling | any rate from one set of weights, no retraining |
| Least reliable regime | multiples from layers thinner than ~5 km |
Takeaways
- Receiver-function forward modelling can be cast as operator learning: one trained Fourier neural operator maps a layered velocity model to a synthetic trace in one to two milliseconds.
- Discretisation invariance means the same weights emit a trace at any sampling rate, so one surrogate serves every downstream inversion without retraining.
- Fidelity is highest on the direct P and Moho Ps arrivals and degrades on late, low-amplitude multiples and on layers thinner than about 5 km.
- The surrogate can be no more accurate than its RAYSUM teacher, and it interpolates rather than extrapolates: outside the training envelope its output is not evidence.
- The payoff is Bayesian: amortized-inference and posterior-sampling pipelines that need millions of forward calls become tractable when each call is a network pass.
References
[1] Frederiksen, A. W., and Bostock, M. G. (2000). Modelling teleseismic waves in dipping anisotropic structures. Geophysical Journal International, 141, 401-412. doi:10.1046/j.1365-246X.2000.00090.x
[2] Li, Z., Kovachki, N., Azizzadenesheli, K., et al. (2021). Fourier neural operator for parametric partial differential equations. ICLR 2021. arXiv:2010.08895
[3] Zhu, L., and Kanamori, H. (2000). Moho depth variation in southern California from teleseismic receiver functions. Journal of Geophysical Research, 105, 2969-2980. doi:10.1029/1999JB900322
[4] Lu, L., Jin, P., Pang, G., et al. (2021). Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators. Nature Machine Intelligence, 3, 218-229. arXiv:1910.03193
[5] Kovachki, N., Li, Z., Liu, B., et al. (2023). Neural operator: learning maps between function spaces with applications to PDEs. Journal of Machine Learning Research, 24(89), 1-97. arXiv:2108.08481
[6] Langston, C. A. (1979). Structure under Mount Rainier, Washington, inferred from teleseismic body waves. Journal of Geophysical Research, 84, 4749-4762. doi:10.1029/JB084iB09p04749
[7] Cranmer, K., Brehmer, J., and Louppe, G. (2020). The frontier of simulation-based inference. PNAS, 117, 30055-30062. arXiv:1911.01429
