Improving moment tensor solutions under Earth structure uncertainty with simulation-based inference
Authors: A. A. Saoulis, T. -S. Pham, A. M. G. Ferreira
Organizations: Department of Physics & Astronomy, University College London, Gower Street, London, WC1E 6BT, United Kingdom · Department of Earth Sciences, University College London, 5 Gower Place, London, WC1E 6BS, United Kingdom · Research School of Earth Sciences, The Australian National University, Canberra, ACT, Australia
Bayesian inference represents a principled way to incorporate Earth structure uncertainty in full-waveform moment tensor inversions, but traditional approaches generally require significant approximations that risk biasing the resulting solutions. We introduce a robust method for handling theory errors using simulation-based inference (SBI), a machine learning approach that empirically models their impact on the observations. This framework retains the rigour of Bayesian inference while avoiding restrictive assumptions about the functional form of the uncertainties. We begin by demonstrating that the common Gaussian parametrisation of theory errors breaks down under minor (1−3%) 1-D Earth model uncertainty. To address this issue, we develop two formalisms for utilising SBI to improve the quality of the moment tensor solutions: one using physics-based insights into the theory errors, and another utilising an end-to-end deep learning algorithm. We then compare the results of moment tensor inversion with the standard Gaussian approach and SBI, and demonstrate that Gaussian assumptions induce bias and significantly under-report moment tensor uncertainties. We also show that these effects are particularly problematic when inverting short period data and for shallow, isotropic events. On the other hand, SBI produces more reliable, better calibrated posteriors of the earthquake source mechanism. Finally, we successfully apply our methodology to two well studied moderate magnitude earthquakes: one from the 1997 Long Valley Caldera volcanic earthquake sequence, and the 2020 Zagreb earthquake.
Figures & tables
Figure 1: Panel (a) shows the effect of κ=3% fractional perturbations to the layer depths, Vp and Vs of the Southern California 1-D velocity model ( Dreger & Helmberger, 1990 ) . Panel (b) shows the resulting uncertainty in the observations for a single station component, as well as the variance in the seismograms (i.e. the diagonal of Ct in (c)). Panel (c) shows the three independent contributions to the Gaussian covariance C . Ct(m) is computed for an event and receiver configuration similar to the LV2 earthquake studied in this work (see Sections 6 & 7.1 ). Note the differences in the colourbar scales, where the theory errors dominate for this moderate magnitude event.
Figure 2: Illustration of the two SBI frameworks introduced in this work. Panel (a) shows optimal score compression, which projects the residuals D−μ∗ for stations ST01, ST02, ST03, etc. onto error weighted sensitivity kernels G to produce compressed observations t . Panel (b) instead uses a deep learning algorithm to learn a flexible, non-linear compression operation. A shared CNN processes full-waveforms into dense features with T time-steps and C feature channels. These station-level features are passed to an N-block axial transformer to aggregate information across the station array before producing a final compressed representation. Panel (c) visualises a normalizing flow architecture, used by both compression frameworks to perform empirical density modelling. The normalizing flow is trained to model the posterior over source parameters qϕ(m∣t) by transforming samples from a simple latent distribution p0(z) using learnable transformations fϕi .
Figure 3: Earthquake location (gold star), focal mechanism (red beachball), and station configuration (brown triangles) for the two events studied in this manuscript. Panel a) shows the LV2 volcanic event in California originally studied in Dreger et al. (2000) , with focal mechanism solution from Phạm & Tkalčić (2021) . We also use this source-receiver configuration for the synthetic experiments in the main text. Panel b) shows the 2020 tectonic event near Zagreb, Croatia with focal mechanism solution from Hu et al. (2025) .
Figure 4: A comparison between the analytically expected χred2 statistic distribution against their empirical distributions under varying levels of Earth structure uncertainty ( κ∈[0.1,1,3,5] ). The degree of agreement between the analytic and empirical distributions probes how well the Gaussian likelihood models the observed variability in the data observations D . We quantify this with the KS statistic in Equation 26 . We find that even under very minor Earth structure perturbations ( ∼1% ) a Gaussian likelihood becomes a very poor approximation.
Figure 5: Examples of the per-trace χred2 statistic under Earth structure perturbations for κ=3% . Higher values of χred2 are denoted by a brighter-copper colour, with significant departure from those expected under a Gaussian likelihood. Even under modest perturbations in the Earth structure, there are significant non-Gaussian effects in the observations.
Figure 6: Posterior inference on an artificial example comparing the results using a Gaussian likelihood (red) against the two SBI frameworks: score compression with shallow density estimation (blue) and deep learning based compression and density estimation (purple). Panel (a) shows the observation and posterior predictive checks for the two SBI approaches. The focal mechanism parameters γ,δ,MW and best fitting double-couple plane orientation posterior are shown in (b), as well as the location on the lune plot in (c) (the top left pair plot in (b) shows the same lune coordinates). The ±[1,2]σ contours are shown for each posterior, and the true artificial solution is marked by the black dashed line and black beachball in (c).
Figure 7: Posterior inference on a more problematic artificial example, where the Gaussian likelihood fails to handle the randomly perturbed 1-D velocity structure. Panel (a) shows the observation and posterior predictive checks for the two SBI approaches. The focal mechanism posteriors are shown in (b), as well as the location on the lune plot in (c) (the top left pair plot in (b) shows the same lune coordinates), with the artificial solution marked in black. Plot details same as in Figure 6.
Figure 8: Calibration statistics for each of the three inference methods over 400 artificial inversions with randomly sampled input source models and perturbed 1-D Earth models. The ideal calibration is denoted by the black dashed line on the diagonal. The Gaussian likelihood approach yields overconfident ( ∼50% ) and slightly biased posterior contours, resulting from inaccuracies in the approximate likelihood. The score compression SBI approach yields much more faithful posterior contours, with a minor degree of bias likely caused by the prior truncation procedure (discussed in the main text). The ML-based compression SBI approach is slightly conservative ( ∼10% ) but nonetheless trustworthy and unbiased. The ±[1,2]σ uncertainty regions are estimated from bootstrapping the coverage results 50 times.
Method
Compression
Metric
γ
δ
Mw(10−2)
strike
dip
rake
SBI
Optimal score
σ
1.51
2.10
0.96
2.56
1.91
2.94
Bias
-0.17
0.49
-0.23
0.05
-0.24
0.29
Deep learning
σ
0.71
0.84
0.37
1.25
0.95
1.52
Bias
-0.06
0.06
-0.01
-0.23
0.06
-0.26
Gaussian
-
σ
0.66
0.87
0.45
0.79
0.71
1.04
Bias
-0.15
0.15
-0.32
0.07
-0.18
0.53
Table 1: Average posterior standard deviations ( σ ) and biases for each inference method over the 400 artificial inversions with Earth structure uncertainty of κ=5% (see also Figure 6 ). Bias is computed per inversion by computing the difference between the mean of the modelled posterior and the true source parameters. Mw values are reported in units of 10−2 , all other values in units of degrees.
Figure 9: The calibration performance of the Gaussian likelihood inversions under varied experimental configurations. Panel (a) shows the results under varying levels of Earth structure uncertainty, parametrised by fractional perturbations κ . Panel (b) explores two different experimental configurations: one with a balanced network with full azimuthal coverage, and another that includes shorter periods in the inversion.
Figure 10: Inference results with each approach on 10 random artificial shallow isotropic sources. We show the isotropic half of the lune plot for each inversion, with the true solution marked by an orange diamond. The focal mechanism decomposition of the true source into its fractional isotropic, CLVD and double-couple components is shown to the left of each lune. The posterior solutions from the Gaussian likelihood approach (red), the score-compression SBI approach (blue) and the ML-based compression SBI approach (purple) are shown as contours denoting the ±[1,2] sigma regions.
Figure 11: Posterior inference on an the LV2 volcanic event in Southern California. Panels (a), (b) and (c) are as before in Figure 6 , with the gold coloured beachball in (c) showing the LV2 focal mechanism solution from Phạm & Tkalčić (2021) .
Figure 12: Posterior inference on the Zagreb 2020 event. Panel (a) shows the posterior predictive checks for the two SBI approaches. For visualisation purposes, for the ML-based compression approach we plot a random subset of the 10% lowest χ2 synthetics (to avoid showing the broad prior over station time shifts p(t) ). Panels (b) and (c) are as before, with the gold coloured beachball in (c) showing the Zagreb 2020 mainshock (deviatoric) solution from Hu et al. (2025) .
Appendix figures & tables12 assets
Supplementary material from the paper’s appendix.
Appendix
Figure S1: Full posterior inference on an artificial short period inversion comparing the results using a Gaussian likelihood (red) against the two SBI frameworks: score compression with shallow density estimation (blue) and deep learning based compression and density estimation (purple). Panel (a) shows a zoom in on the observation and posterior predictive checks for the two SBI approaches. The focal mechanism parameters γ,δ,MW and best fitting double-couple plane orientation posterior are shown in (b), as well as the location on the lune plot in (c) (the top left pair plot in (b) shows the same lune coordinates). The ±[1,2]σ contours are shown for each posterior, and the true artificial solution is marked by the black dashed line and black beachball in (c).
Figure S2: The (artificial) balanced network configuration for one of the synthetic tests in the main text. Stations US.MNV and dummy receivers with network code DM were added to this example, and were not used for the rest of the main text.
Figure S3: The calibration performance of the SBI with optimal score compression under varied experimental configurations. Panel (a) shows the results under varying levels of Earth structure uncertainty, parametrised by κ . Panel (b) explores two different experimental configurations: one with a balanced network with full azimuthal coverage, and another that includes shorter periods in the inversion.
Figure S4: The calibration performance of the ML-based compression SBI approach under varied experimental configurations. Panel (a) shows the results under varying levels of Earth structure uncertainty, parametrised by κ and differing filtering bands. Panel (b) explores model performance when exposed to greater or less Earth structure uncertainty than during initial training.
Figure S5: A detailed look at the posterior solutions for the 300 artificial inversions used to perform coverage testing, for the κ=3% example. Panel (a) shows the per-parameter distributions of the bias between the posterior mean and the true solution. Panel (b) shows the same for the z -score distribution. Panel (c) shows the ordered per-parameter z -scores for all inversions against the expected distribution for well-calibrated (Gaussian) posteriors.
Figure S6: A detailed look at the posterior solutions for the 300 artificial inversions used to perform coverage testing, for the κ=5% example. Panel (a) shows the per-parameter distributions of the bias between the posterior mean and the true solution. Panel (b) shows the same for the z -score distribution. Panel (c) shows the ordered per-parameter z -scores for all inversions against the expected distribution for well-calibrated (Gaussian) posteriors.
Figure S7: A detailed look at the posterior solutions for the 300 artificial inversions used to perform coverage testing, for the short period example. Panel (a) shows the per-parameter distributions of the bias between the posterior mean and the true solution. Panel (b) shows the same for the z -score distribution. Panel (c) shows the ordered per-parameter z -scores for all inversions against the expected distribution for well-calibrated (Gaussian) posteriors.
Figure S8: The 2-D composite Earth model adapted from Hu et al. (2025) to include 1-D velocity model perturbations. Panel (a) shows the network configuration and which model was used in each region. Panel (b) shows the model velocities as a function of depth and the corresponding κ=3% perturbations used in the main text. We keep model numbering consistent with Hu et al. (2025) .
Method
Compression
Metric
γ
δ
Mw(10−2)
strike
dip
rake
SBI
Optimal score
σ
0.42
0.77
0.28
0.64
0.49
0.78
Bias
-0.06
-0.07
-0.02
-0.09
-0.03
-0.02
Deep learning
σ
0.53
0.65
0.26
0.69
0.56
0.87
Bias
-0.01
-0.04
-0.01
-0.02
0.05
-0.01
Gaussian
-
σ
0.26
0.42
0.17
0.38
0.30
0.47
Bias
-0.05
0.01
-0.08
-0.07
-0.11
-0.09
Appendix
Table S1: Average posterior standard deviations ( σ ) and biases for each inference method over the 300 artificial inversions with Earth structure uncertainty of κ=3% . Bias is computed per inversion by computing the difference between the mean of the modelled posterior and the true source parameters. Mw values are reported in units of 10−2 , all other values in units of degrees.
Method
Compression
Metric
γ
δ
Mw(10−2)
strike
dip
rake
SBI
Optimal score
σ
1.49
1.77
1.12
3.02
1.95
3.04
Bias
-0.01
-0.57
0.34
-0.43
-0.14
0.16
Deep learning
σ
0.75
1.01
0.40
1.12
0.99
1.50
Bias
-0.21
-0.02
-0.05
-0.09
0.10
-0.07
Gaussian
-
σ
0.52
0.68
0.37
0.92
0.58
0.90
Bias
0.03
-0.29
0.50
-0.33
0.11
-0.26
Appendix
Table S2: Average posterior standard deviations ( σ ) and biases for each inference method over the 300 artificial inversions with Earth structure uncertainty of κ=5% and shorter period ( 6−50 s) data. Units as before.
Method
χred2
Correlation misfit
Envelope misfit
Power misfit
MSE
Score compression SBI
0.563
0.052
0.421
0.001
3.72e-13
ML compression SBI
0.590
0.095
0.455
0.001
5.04e-13
Gaussian likelihood
0.606
0.110
0.419
0.001
7.09e-13
Appendix
Table S3: LV2: Best-subset posterior predictive performance. Each column reports the mean value of a given metric computed on the top 50 posterior predictive samples selected according to that metric. Lower values indicate closer fits for all metrics.
Method
χred2
Correlation misfit
Envelope misfit
Power misfit
MSE
Score compression SBI
2.648
0.055
0.482
0.008
1.53e-11
ML compression SBI
2.808
0.083
0.458
0.005
1.90e-11
Gaussian likelihood
3.575
0.182
0.438
0.002
3.97e-11
Appendix
Table S4: Croatia: Best-subset posterior predictive performance. Each column reports the mean value of a given metric computed on the top 50 posterior predictive samples selected according to that metric. Lower values indicate closer fits for all metrics.
Simulation-based inference (SBI) with machine learning is an increasingly important tool for solving inverse problems in science and engineering, including parameter inference and the inversion of detector effects. We provide an overview of the Bayesian and frequentist statistical frameworks, describe how machine-learning-based SBI methods, such as neural posterior estimation and neural likelihood estimation, can be used for parameter estimation within these frameworks, and show that the same methods can also be applied to Empirical Bayes or unfolding tasks. We also discuss how to validate inference results and the limitations of SBI with machine learning.
Maximilian Dax, Theo Heimel, Gilles Louppe
ELLIS Institute Tübingen, Max Planck Institute for Intelligent Systems, Tübingen AI Center, Tübingen, Germany · CP3, Université catholique de Louvain, Louvain-la-Neuve, Belgium · Montefiore Institute, University of Liège, Belgium
Accurate calibration of particle-based simulators is crucial for robotic earthwork simulation, but analytical calibration is challenging due to this task's highly nonlinear particle dynamics and the black-box nature of conventional simulators. Although simulation-based inference (SBI) can estimate posterior distributions over simulation parameters solely from forward simulations, applying SBI directly to high-fidelity (HF) particle simulators is often computationally prohibitive. Low-fidelity (LF) simulators with coarser particles can reduce this cost, but changes in particle size and particle count shift the parameter values needed to reproduce the same observation, producing biased LF posteriors. We propose Bridged SBI, which leverages a biased but informative LF posterior to guide HF inference. This method first uses inexpensive LF simulations to identify a coarse high-density parameter region, and then it learns a local residual bridge to transport LF posterior samples toward HF-consistent regions by correcting the LF--HF discrepancy. We analyze how sequential multi-fidelity SBI (Naive-MF) can suffer from LF-induced posterior miscoverage when it directly relies on the LF posterior without discrepancy correction. We then show that Bridged SBI is designed to alleviate this issue by explicitly modeling the LF--HF discrepancy through residual correction. Experiments on both sim-to-sim particle-parameter calibration and real-to-sim calibration with real soil observation show that Bridged SBI produces more accurate and reliable HF posteriors than HF-only SBI or the Naive-MF baseline, especially under limited HF simulation costs.
Gahee Kim, Yuki Kadokawa, Sandro M. Alcantara Tacora +5
NARA INSTITUTE OF SCIENCE AND TECHNOLOGY, Nara, Japan · PUBLIC WORKS RESEARCH INSTITUTE, Ibaraki, Japan
Accurate 3D seismic velocity modeling through seismic travel-time tomography using both active- and passive-source data provides critical underpinning models for seismicity monitoring and hazard assessment. Because travel-time tomography is an inherently ill-posed inverse problem, UQ of the estimated models using Bayesian methods is also important for reliable downstream interpretations and analyses. However, Bayesian inference for 3D tomography based on conventional grid-based representations faces the ``curse of dimensionality'' and severe computational bottlenecks. Consequently, rigorous Bayesian UQ for margin-wide 3D travel-time tomography has remained largely unexplored. In this study, we propose a meshless 3D Bayesian travel-time tomography method that combines PINNs with a neural representation of the velocity structure, enabling tractable and data-efficient Bayesian inference through function-space particle-based variational inference. To efficiently integrate passive-source data into the Bayesian estimation of the velocity structure, we conduct analytical marginalization treating uncertain source parameters as nuisance parameters, with passive-source relocation carried out in post-processing. We validated the capability of our approach for 3D problems through synthetic experiments. Furthermore, we applied the method to a real-world dataset from marine active-source surveys and natural earthquakes off the Kii Peninsula, Nankai Trough. Our probabilistic 3D ensemble successfully resolves key geological features and provides data-consistent uncertainty maps. The posterior mean hypocenters shifted mainly in the vertical direction by 10-15 km, consistent with a previous relocation result. Finally, the neural representation drastically reduces storage requirements for the entire ensemble velocity model, highlighting the scalability and data efficiency of the proposed framework.
Ryoichiro Agata, Kazuya Shiraishi, Gou Fujie +1
Japan Agency for Marine-Earth Science and Technology · Earth Science New Zealand