Paper deep dive
Improving moment tensor solutions under Earth structure uncertainty with simulation-based inference
A. A. Saoulis, T. -S. Pham, A. M. G. Ferreira
Intelligence
Status: succeeded | Model: google/gemini-3.1-flash-lite-preview | Prompt: intel-v1 | Confidence: 94%
Last extracted: 3/22/2026, 6:11:08 AM
Summary
The paper introduces a robust method for moment tensor inversion using simulation-based inference (SBI) to account for Earth structure uncertainty. By replacing traditional Gaussian likelihood approximations with machine learning-based density estimation, the authors demonstrate that SBI provides more reliable, better-calibrated posteriors for earthquake source mechanisms, particularly for short-period data and shallow, isotropic events.
Entities (5)
Relation Signals (3)
Simulation-based inference → appliedto → Zagreb earthquake
confidence 95% · 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.
Simulation-based inference → improves → Moment tensor inversion
confidence 95% · we develop two formalisms for utilising SBI to improve the quality of the moment tensor solutions
Gaussian likelihood → inducedbiasin → Moment tensor inversion
confidence 90% · demonstrate that Gaussian assumptions induce bias and significantly under-report moment tensor uncertainties.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract: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.
Tags
Links
- Source: https://arxiv.org/abs/2603.18925v1
- Canonical: https://arxiv.org/abs/2603.18925v1
Trouble viewing inline? Open PDF directly →
Full Text
140,974 characters extracted from source content.
Expand or collapse full text
(2026) 000, 1–23 Improving moment tensor solutions under Earth structure uncertainty with simulation-based inference A. A. Saoulis 1,2⋆ , T.-S. Pha . m 3 , A. M. G. Ferreira 2 1 Department of Physics & Astronomy, University College London, Gower Street, London, WC1E 6BT, United Kingdom 2 Department of Earth Sciences, University College London, 5 Gower Place, London, WC1E 6BS, United Kingdom 3 Research School of Earth Sciences, The Australian National University, Canberra, ACT, Australia 20 March 2026 SUMMARY Bayesian inference represents a principled way to incorporate Earth structure uncertainty in full-waveform moment tensor inversions, but traditional approaches generally require signifi- cant 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 the- ory 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 ten- sor 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 in- duce 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 earth- quake 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. Key words: Earthquake source observations – Waveform inversion – Bayesian inference – Inverse theory – Machine learning 1 INTRODUCTION Focal mechanism solutions of full-waveform moment tensor inversion have been shown to rely on high-quality forward mod- elling, which primarily depends on accurate Earth structure mod- els (Hardebeck & Shearer, 2002; Hj ̈ orleifsd ́ ottir & Ekstr ̈ om, 2010; Takemura et al., 2020; R ̈ osler et al., 2021; Song et al., 2022; Simut ̇ e et al., 2023, see Pha . m et al., 2025 for a recent review). These effects are particularly pronounced in regions with strongly heterogeneous or poorly constrained crustal structure (Abercrombie & Ekstr ̈ om, 2001; K ˇ r ́ ı ˇ zov ́ a et al., 2013; Petersen et al., 2021). Reliable moment tensor solutions are also critical for downstream applications such as seismic tomography (Blom et al., 2023). Beyond characterising individual events, heterogeneities can lead to systematic biases in 1-D Earth model derived source cata- logs (Hingee et al., 2011; Hejrani et al., 2017; R ̈ osler et al., 2024a). These inaccuracies can induce biases in, and specious trade-offs between, moment tensor components (Dziewonski & Woodhouse, ⋆ a.saoulis@ucl.ac.uk 1983; Ferreira & Woodhouse, 2006; Weston et al., 2011; Ferreira et al., 2011). Such effects have been shown to cause, for instance, systematic spurious non-double-couple components in retrieved seismic source solutions for tectonic events (Sawade et al., 2022; R ̈ osler et al., 2024b). Moreover, unresolved Earth structure uncer- tainty continues to pose significant challenges for interpreting seis- mic source solutions (e.g., Akuhara et al., 2025). An appealing so- lution is Bayesian inference, which provides a framework for treat- ing theory or mismodelling errors directly in the uncertainty esti- mation procedure (Tarantola, 2005). We assume the uncertain 1-D Earth model to be the primary source of theory errors in the context of this work. Including the effect of theory errors in source inversions can significantly improve inversion quality (e.g., Yagi & Fukahata, 2011; Minson et al., 2013, 2014; Duputel et al., 2015; Pha . m & Tkal ˇ ci ́ c, 2021). This is often achieved by treating theory errors as an additional term in a Gaussian likelihood function, following Tarantola (2005). Several studies approximate the theory errors by sampling over the modelled Earth structure uncertainty empirically and estimating the resulting variability in the observations (Du- arXiv:2603.18925v1 [physics.geo-ph] 19 Mar 2026 2A. A. Saoulis et al. putel et al., 2012; Vasyura-Bathke et al., 2021; Pha . m & Tkal ˇ ci ́ c, 2021; Poppeliers & Preston, 2021; Pha . m et al., 2024), while oth- ers rely on theoretical approximations for the effect of theory er- rors on the observation covariance structure (Duputel et al., 2014; Hallo & Gallovi ˇ c, 2016). All these prior works show that incorpo- rating theory error contributions yields improved solutions over the na ̈ ıve approach. However, they perform significant approximations to make the inverse problem tractable by traditional means. The question of whether the resulting solutions are genuinely accurate and trustworthy thus still remains. There are many other hybrid or empirically motivated ap- proaches. Grid-search strategies have proven effective for trac- ing likelihood contours for uncertainty quantification (Sokos & Zahradn ́ ık, 2013; Vack ́ a ˇ r et al., 2017; Zahradn ́ ık & Sokos, 2025; Thurin et al., 2025), though likelihood contours often do not have a clear probabilistic interpretation. Additionally, making this ap- proach tractable for holistically treating theory errors may prove challenging. Another family of approaches relies on building solu- tion ensembles through hybrid techniques, including bootstrapping (W ́ eber, 2006; Valentine & Trampert, 2012) and by injecting em- pirically motivated theory errors (Poppeliers & Preston, 2022). A relevant example to this study is St ̈ ahler & Sigloch (2014, 2016), which used earthquake source catalogues to build empirical likeli- hood functions that better account for mismodelling errors. While these approaches remain valuable, particularly as a means for deal- ing with unknown (or challenging to model) theory uncertainties, it is desirable to continue to develop a complementary approach that directly integrates our ability to model theory errors. Simulation-based inference (SBI) is a promising alternative to fully treat the effects of theory errors without compromising (or abandoning) a rigorous Bayesian approach to moment tensor in- version. SBI uses machine learning (ML) models to learn the ob- servational uncertainties empirically through direct simulation, in principle allowing arbitrarily complex sources of uncertainty to be modelled without making parametric assumptions about their form (Alsing et al., 2019; Cranmer et al., 2020; Zammit-Mangion et al., 2025). Saoulis et al. (2025) demonstrated that standard Gaussian likelihood assumptions about full-waveform data errors yielded in- accurate moment tensor solutions. SBI, on the other hand, produced accurate results at a fraction of the computational cost. This work applies SBI to full-waveform moment tensor inversion while incor- porating theory errors. SBI has received significant attention across the physical sci- ences for enabling fast and accurate Bayesian inference (Gonc ̧alves et al., 2020; Stockman et al., 2024; von Wietersheim-Kramsta et al., 2025; Atlas Collaboration, 2025; Dax et al., 2025). A key practi- cal challenge in adapting SBI for a given domain is preparing the observational data to make it amenable to ML-based probabilistic modelling (see e.g. Deistler et al., 2025, for a practical perspec- tive). This is typically straightforward in low dimensions, but high- dimensional observations, such as full-waveform multi-station seis- mic data, can pose problems. The choice of data representation therefore plays a central role in determining both the accuracy and efficiency of SBI-based inversions (Gerardi et al., 2024; Park et al., 2025). Within this context, two broad strategies have emerged for handling high-dimensional observations in SBI. One approach uses carefully selected lower dimensional statistics or generic linear compression algorithms to reduce the dimensionality (e.g., Als- ing et al., 2019). An alternative approach relies on data-driven ML methods to learn compressed representations directly from the ob- servations. Increasingly, SBI has turned to deep learning to extract complex, non-linear information from high-dimensional data, im- proving the constraining power, accuracy, and efficiency of the in- verse problem (e.g., Gloeckler et al., 2024; Lanzieri et al., 2025; Jeffrey et al., 2025; Dax et al., 2025). Our work fuses this strand of SBI with a broad trend in seismology utilising deep learning for processing full-waveform data (Zhu & Beroza 2019; Zhu et al. 2019; Mousavi et al. 2020; M ̈ unchmeyer et al. 2021; Sun et al. 2023; McBrearty & Beroza 2023; Si et al. 2024b; see e.g. Mousavi & Beroza 2022; Kubo et al. 2024 for reviews). We propose two SBI frameworks for full-waveform moment tensor inversion that differ in how the observations are compressed. We assess a physics-motivated linear compression method, adapted from Saoulis et al. (2025) to address theory errors, and a data-driven deep learning approach. We explore the advantages of each ap- proach, and compare both SBI frameworks with the most recent attempts using explicit-likelihood Bayesian inversion technique for treating theory errors (Pha . m & Tkal ˇ ci ́ c, 2021; Pha . m et al., 2024). The manuscript is laid out as follows. Section 2 sets out the theoretical background, elucidating the approximations often made in parametric likelihood-based treatments of the theory errors, and explaining how SBI can bypass these limitations. Section 3 presents the Gaussian likelihood approach for full-waveform moment tensor inversions, which serves as a benchmark for SBI. Section 4 presents the technical and architectural details for our two SBI frameworks. Section 5 introduces methods to evaluate the accuracy of Gaussian likelihood assumptions for a given problem, and how we evaluate the moment tensor solution quality of each approach. We then explore a range of applications. Section 6 presents a systematic evaluation of each of the approaches across three case- studies: how each approach fares under moderate Earth structure uncertainty; how Gaussian likelihood assumptions produce inaccu- rate results across a wide range of experimental configurations; and how SBI more accurately recovers moment tensors, including with short-period data and shallow, isotropic sources. We present the re- sults of two real data inversions in Section 7, before discussing the broader limitations and implications of this work in Section 8 and Section 9. 2 THEORETICAL BACKGROUND 2.1 Posterior inference under Earth structure uncertainty For a given observation D, the joint posterior distribution over the source parameters m and the uncertain Earth model parameters Ω is given by p(m, Ω| D) = p(D| m, Ω)p(m)p(Ω) p(D) ,(1) where p(D | m, Ω) is the likelihood (e.g. obtained via the deter- ministic forward model g(Ω, m)), p(m) and p(Ω) are the prior distributions over the source and Earth model parameters respec- tively, and p(D) is the marginal likelihood (evidence). Our ultimate goal, however, is to infer the posterior over the source parameters m while properly accounting for unresolved un- certainty in the Earth model. Analytically, this requires marginalis- ing over Ω: p(m| D) = Z p(m, Ω| D)dΩ,(2) = p(m) R p(D| m, Ω)p(Ω)dΩ p(D) .(3) Theory errors in moment tensor inversion with SBI3 A conceptually straightforward route would be to sample from the joint posterior in Equation (1) directly using Markov Chain Monte Carlo (MCMC) or other joint-sampling algorithms, and then marginalise the empirical sample set over Ω as in Equation (2). This approach is exact in the limit of many posterior samples, and it ex- tends to hierarchical formulations (e.g. by introducing hyperparam- eters in the prior for Ω and sampling them jointly). The principal drawback for seismological problems is cost: each proposed Ω (i) typically requires recomputing Green’s functions, and the posterior dimensionality dimm, Ω may become very large, so joint sam- pling can be extremely expensive. The more popular approach is to attempt to approximate the marginalisation of the likelihood over all prior Earth structures in Equation (3) (Tarantola, 2005): p(D| m) = Z p(D| m, Ω)p(Ω)dΩ,(4) which is typically analytically intractable. Most prior work approx- imates Equation (4) by assuming that the variability introduced by Ω can be represented as an additional stochastic term in the like- lihood (Tarantola et al., 1982). In particular, assuming a Gaussian likelihood with data covariance matrix C d , we may write: p(D| m) = Z N D g(m, Ω), C d p(Ω)dΩ,(5) whereafter g(m, Ω) can be expanded perturbatively around some fiducial model Ω ∗ to linear order. This yields g(m, Ω)≈ g(m, Ω ∗ ) + K Ω | Ω ∗ δΩ + ...(6) where K Ω = ∇ Ω g(m, Ω) are the gradients (i.e., sensitivity ker- nels) with respect to the Earth model parameters. This can be sub- stituted into Equation (5): p(D| m) = Z N D g(m, Ω ∗ ) + K Ω δΩ, C d p(Ω)dΩ. (7) To make progress, one can then assume that the Earth model pertur- bations are also Gaussian distributed following δΩ ∼ N (0, Σ Ω ). Applying simple Gaussian convolution rules produces: p(D| m)≈N D g(m, Ω ∗ ), C d + C t (m) ,(8) where C t (m) = K Ω (m)Σ Ω K T Ω (m) is the theory covariance ma- trix. Note that kernels K Ω may in general be expensive to com- pute. Equation (8) can also be reached without assuming Gaussian- ity of δΩ (e.g. as in Duputel et al., 2014). In that case, however, the induced distribution p(D | m) is not strictly Gaussian, and a Gaussian approximation is being introduced at the level of the like- lihood. In practice, C t (m) can be estimated by Monte Carlo sam- pling. One can draw N samples from the prior p(Ω), evaluate the forward model, and compute the empirical covariance of the mod- elled waveforms: Ω (n) ∼ p(Ω), g (n) (m)≡ g(Ω (n) , m), n = 1,...,N, (9) b C t (m)≡ 1 N − 1 N X n=1 g (n) (m)−g(m) g (n) (m)−g(m) T , (10) where g(m) is a mean estimate of the synthetic samples g (n) (m). The Gaussian likelihood in Equation (8) can then be used for max- imum likelihood estimation and posterior inference. At this stage, it is worth taking stock of the significant simpli- fications required to make previous approaches tractable. We made the following approximations: (i) Perturbations to structure lead to small and additive pertur- bations on the forward model operator. (i) Both the prior p(Ω) and the data noise model p(D g(m, Ω), C d ) must be Gaussian. (i) The theory error contribution C t (m) is only valid locally for a given set of source parameters m, and is estimated empirically to avoid direct computation of Earth model kernels K Ω . (iv) The data and theory error contributions are independent and additive. We note that these are not novel observations; (i) and (i) are dis- cussed in Tarantola (2005, Sec. 5.8.6), for instance. Of these approximations, (i)-(i) all pose significant con- straints and may introduce inaccuracies during posterior inference. In particular, (i) may be violated for realistic Earth structure per- turbations that induce non-linear and multiplicative effects on the full waveforms, such as phase offsets and amplitude anomalies. (i) is highly restrictive since in general p(Ω) may not be Gaussian (for instance, if p(Ω) is taken from a tomographic ensemble), and the data noise errors may not be Gaussian (see e.g. Saoulis et al., 2025). For (i), the fact that the covariance varies over parameter space m poses several practical problems that increase the com- putational costs of inversions. In addition, estimating and using the dense matrix b C is often computationally challenging or intractable, so typically practitioners make simplifications (e.g. block sparse, per-station covariance matrices). 2.2 Simulation-based inference for theory errors Simulation-based inference (SBI) offers an alternative to such approximate analytic or sampling-based marginalisation. SBI draws simulations from the joint prior predictive model: m∼ p(m),Ω∼ p(Ω), D∼ p(D| m, Ω).(11) This builds a realistic data set of observationsD: D =m, D(12) which implicitly encodes the key quantities required for posterior inference, like the marginalised posterior p(m | D) from Equa- tion (2) or marginalised likelihood p(D| m) from Equation (4). Within this framework, a variety of probabilistic quantities may be modelled using machine learning, including the joint dis- tribution, the likelihood (or ratios between likelihoods) or the pos- terior itself (Lueckmann et al., 2019, 2021; Durkan et al., 2020b; Hermans et al., 2020; Cranmer et al., 2020). For example, an em- pirical likelihood can be combined with classical sampling-based approaches such as MCMC to perform posterior inference (Papa- makarios et al., 2019). In this work, however, we focus on neu- ral posterior estimation (NPE) in which the posterior density is modelled directly, enabling rapid and amortised posterior sampling once the model has been trained (Papamakarios & Murray, 2016; Greenberg et al., 2019; Deistler et al., 2022). SBI proceeds by usingD to train a machine learning model q φ with learnable parameters φ to approximate one of the key proba- bilistic quantities. This is achieved by training a neural network to minimise some loss function quantifying the divergence between the true and modelled probability distributions (p ∗ (x) and q φ (x), respectively). A typical choice of loss function is the Kullback- Liebler (KL) divergence (Papamakarios et al., 2021): L = KL p ∗ (x)∥q φ (x) =E p ∗ (x) [logp ∗ (x)− logq φ (x)] = const.−E p ∗ (x) [logq φ (x)], (13) 4A. A. Saoulis et al. Figure 1. Panel (a) shows the effect of κ = 3% fractional perturbations to the layer depths, V p and V s 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 C t in (c)). Panel (c) shows the three independent contributions to the Gaussian covariance C. C t (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. where the constant term can be ignored in an optimisation routine. The resulting objective function to be minimised for, say, the pos- terior p(m| D) is then: L =−E (m i ,D i )∼D [logq φ (m i | D i )].(14) This has a simple interpretation: given observation D i , the model q φ is trained to maximise the posterior density assigned to the each “true” source parameter m i . Since Ω is sampled from its prior dur- ing dataset generation, the learned conditional density q φ (m | D) implicitly marginalises over Earth model uncertainty, avoiding the need for costly joint MCMC sampling or simplifying assumptions about the likelihood. With this framework at hand, applying SBI hinges on the prac- tical question of how to perform robust empirical density estima- tion. Much of the success of SBI over the last decade has been driven by the development of specialised ML models, referred to as neural density estimators (NDEs), which have proven capable of reliable and flexible density modelling (see e.g., Alsing et al., 2019; Cranmer et al., 2020; Deistler et al., 2025). This work makes use of several popular classes of NDEs to empirically model the required density functions, with details presented in Section 4. 3 THE GAUSSIAN LIKELIHOOD 3.1 Approximating the likelihood We intend to make comparisons between moment tensor in- versions performed with a Gaussian likelihood, parametrised by some covariance b C, and inversions using SBI. We follow prior work by splitting the covariance into a theory error C t (m) and data error C d contributions (Tarantola, 2005): b C = C d + C t (m).(15) We estimate C t (m) using Monte Carlo sampling outlined in Sec- tion 2, Equation (10). Pha . m et al. (2024) demonstrated that es- timating g(m) as a mean of the Monte Carlo samples in Equa- tion (10) leads to bias. Instead, one should use the fiducial observa- tiong(m) = g(m, Ω ∗ ) for a better estimate. There is a wide range of approaches for specifying C d . Prior work has demonstrated that several of these, such as constant diag- onal or exponentially tapered covariances, are poor representations of the measurement noise that can bias the inversions (Vasyura- Bathke et al., 2021; Saoulis et al., 2025). Here, we use the empir- ically motivated exponential-tapered cosine covariance from Kolb & Leki ́ c (2014): C ij d = σ 2 e −λτ ij cos λω 0 τ ij ,(16) where τ ij = |t j − t i | is the time lag between samples and σ 2 is the noise variance. We use the empirical choice from Kolb & Theory errors in moment tensor inversion with SBI5 Leki ́ c (2014) for ω 0 = 4.4 and use λ = 0.05 to match maxi- mum frequency of the filtering band for our data. The noise level σ 2 is estimated from the filtered noise prior to an event for each re- ceiver–component. For the high signal-to-noise events studied here, theory errors dominate and the data covariance parametrisation C d does not impact the quality of the inversions. Note that, following prior work, we only model the per station- component covariances in the Gaussian likelihood covariance b C, leading to a block sparse form that assumes all cross-component and cross-station covariances are zero. This hugely reduces the computational complexity of estimating (and inverting) b C, mak- ing the problem tractable. However, velocity model perturbations will produce non-zero cross-component, cross-station covariances (as will realistic noise, to a lesser extent). All of these effects are therefore neglected in the Gaussian likelihood approach. In con- trast, the SBI approach does not use an explicit covariance matrix for inference and can therefore capture and model these correla- tions implicitly. 3.2 Application to real data In practice, we found the covariance in Equation (15) failed to provide an adequate model of the uncertainties for the real data inversions. Since we did not encounter such issues during the syn- thetic experiments, we hypothesise that this is the result of further, unmodelled uncertainties. For instance, 3-D heterogenieties will likely cause theory errors in the observations that are not captured by our forward model and Earth structure perturbations. In order to stabilise the inversions, we added a small diagonal regularisation term to the covariance matrix: b C = C d + C t (m) + C reg ,(17) where C reg = αI. Each covariance term in Equation (17) is visu- alised in Figure 1. We choose α empirically to be a small fraction of the maximum magnitude of the theory errors: α = ε max C t (m), where ε is chosen empirically. This approach has the advantage of encoding the heuristic that the unmodelled theory errors will be related to the magnitude of the modelled theory errors. This regu- larisation term is a form of covariance shrinkage (Ledoit & Wolf, 2004), a widely used approach for mitigating covariance mismod- elling effects. We note that our covariance b C takes a similar form to Pha . m & Tkal ˇ ci ́ c (2021), where instead a diagonal contribution to the covari- ance was motivated by and interpreted as data noise. We argue that this interpretation is unphysical, in the sense that uncorrelated noise is not a realistic model of the data noise covariance (Duputel et al., 2012; Vasyura-Bathke et al., 2021; Saoulis et al., 2025). Instead, the diagonal term should be understood as a way of regularising against unmodelled theory errors. 4 APPLYING SBI TO THEORY ERRORS 4.1 Neural density estimation A central requirement of SBI is the ability to accurately model target probability densities such as the posterior p(m | D), which may exhibit strong non-Gaussian structure, complex correlations, and multimodality. Neural density estimators (NDEs) address this problem by learning flexible, parameterised distributions that can be trained with standard ML techniques. Normalizing flows provide a particularly effective class of NDEs for SBI (Alsing et al., 2019). The key idea is to represent a complex target density by learning a transformation that maps samples from a simple base distribution to the target. Concretely, let z ∼ p 0 (z) denote a draw from a simple base distribution (typi- cally a standard multivariate Gaussian). A normalizing flow defines an invertible transformation f φ , parametrised by neural networks, such that m = f φ (z; D),(18) where explicit dependence on the observations D allows the trans- formation to model the conditional density q φ (m| D). The result- ing density over m is obtained via the change-of-variables formula, logq φ (m| D) = logp 0 (z) + log det ∂f −1 φ ∂m .(19) with z = f −1 φ (m; D). This formulation yields a simple, compu- tationally tractable expression for the loss in Equation (14), allow- ing the neural network parameters of the transformation f φ to be learned directly by optimisation. In practice, the transformation f φ is constructed as a composi- tion of multiple simple transformations f i φ (see Figure 2c, allowing the overall model to represent highly complex densities while re- maining tractable for both density evaluation and sampling. This hinges on designing f φ with an easily computable Jacobian (the ∂f −1 φ /∂m term in Equation (19)). Different normalizing flow ar- chitectures correspond to different choices of these transformations and their parameterisations, trading off expressivity, computational efficiency, and numerical stability. In this work, we employ masked autoregressive flows (MAFs; Papamakarios et al., 2017) and neural spline flows (Durkan et al., 2019), which are both well-established flow classes that have demonstrated strong performance in SBI and related conditional density estimation tasks. We defer to prior work for a more de- tailed exposition of normalizing flows and their variations (e.g., Papamakarios et al., 2021). 4.2 Compression algorithms In practice, modelling the distribution p(m | D) directly is challenging. This largely results from the high-dimensionality of full-waveform seismic observations, which contain complex latent relationships that are non-trivial to model. Prior work in SBI ad- dresses this issue by compressing the high dimensional observa- tions D into some lower dimensional summary statistics t. Once the compression algorithm is defined, the training datasetD (Equation 12) can be compressed to a set pairs of source parameters m and compressed observations t. SBI then proceeds as usual, modelling the posterior proxy p(m | t). A diagram of this workflow is shown in Figure 2. Note that while this is not the “true” posterior p(m| D), this approach has the advantage of con- structing a posterior solution p(m | t) that is statistically consis- tent with the true source parameters m. This is not the case when significant approximations are made in standard likelihood-based inference techniques. We propose two compression frameworks for applying SBI in the context of theory errors. 4.2.1 Optimal score compression The first relies on a physically-motivated compression algo- rithm to reduce the dimensionality of observations D to a set of 6A. A. Saoulis et al. 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 p 0 (z) using learnable transformations f i φ . informative summary statistics t. The goal of a compression algo- rithm is to preserve as much discriminative power in the summaries t with respect to the source parameters m as possible. One popular and theoretically motivated approach is optimal score compression (Alsing et al., 2017; Alsing & Wandelt, 2018; Alsing et al., 2019). Optimal score compression requires a known, analytic expres- sion for the likelihood p(D | m). For instance, we may re-use our approximate Gaussian likelihood from Equation (8) with some un- known covariance C: p(D| m) =N D g(m, Ω ∗ ), C).(20) Score compression then proceeds by Taylor expanding the likeli- hood p(D | m) about a fiducial set of source parameters m ∗ to first order. The data sensitivities with respect to the source param- eters, G m ∗ = ∇ m ∗ D can then be used to construct an optimal compression algorithm (Alsing et al., 2017). Specifically, Alsing & Wandelt (2018) demonstrate that one such optimal algorithm is t = m ∗ + F −1 ∗ G T m ∗ C −1 (D−μ ∗ ),(21) whereμ ∗ is the fiducial data vectorμ ∗ = g(m ∗ , Ω ∗ ), and the Fisher matrix F ∗ = G T m ∗ C −1 G m ∗ . Equation (21) has a particu- larly simple interpretation: the summary t is the (local) maximum likelihood estimate (MLE) of the source parameters m MLE given an observation D (Tarantola, 2005). For a detailed mathematical derivation we defer to prior work (Alsing et al., 2017; Alsing & Wandelt, 2018), including a previous application of optimal score compression to seismic source inversion (Saoulis et al., 2025). We provide some discussion on our exclusion of higher order terms in Section S1. There are, however, some limitations of optimal score com- pression as applied to seismic source inversions. As highlighted in Saoulis et al. (2025), nonlinearity in the forward model leads to rapid degradation of compression algorithm, and therefore the use- ful information in summaries t. While the forward model is linear with respect to the moment tensor, perturbations in the Earth struc- ture introduce nonlinearities (Pha . m & Tkal ˇ ci ́ c, 2021; Pha . m et al., 2024). We therefore adopt the two-stage approach utilised in Saoulis et al. (2025). First, we utilise the approximate Gaussian-likelihood to find a maximum likelihood estimate of the source parameters m MLE . This is achieved through an iterative least squares algo- rithm: at each iteration, we estimate the local theory covariance b C t (m), and then apply the maximum likelihood estimate in Equa- tion (21). The second step is to then truncate the prior p(m) over a region where the compression is accurate. We follow the proce- dure taken in Saoulis et al. (2025) to use a Gaussian estimate of the posterior to propose a prior that should support the target pos- terior p(m | t). In this approximation, the posterior covariance is given by the inverse Fisher information matrix F −1 (Tarantola, 2005; Alsing et al., 2017). The diagonal of this covariance pro- vides an estimate of the characteristic per-parameter uncertainties σ, which we use to define a weakly informative truncation of the prior. We define a uniform prior centred on m MLE with width Nσ along each parameter direction, where N is a free parameter to be tuned: p Fisher (m;N ) =U (m MLE −N·σ, m MLE +N·σ)∗p(m). (22) SBI then proceeds as usual by sampling from the truncated prior region m∼ p Fisher (m;N ). This approach has several downsides. For one, at all stages it relies on a Gaussian likelihood to compute a local maximum likelihood estimate. While Pha . m & Tkal ˇ ci ́ c (2021); Pha . m et al. (2024) demonstrated that the approximate Gaussian likelihood is sufficiently accurate to recover an improved MLE, it may also in- Theory errors in moment tensor inversion with SBI7 troduce bias. By then truncating the prior around a potentially bi- ased estimate of the MLE, we run the risk of biasing the entire posterior inference process. We probe the degree to which this ef- fect can bias our posteriors in Section 6. Nonetheless, an improved formulation for utilising SBI is desirable. 4.2.2 Machine learning-based compression Optimal score compression is appealing as a generic and sim- ple compression algorithm. However, as its core assumptions break down — local linearity and a well-specified likelihood — its effec- tiveness at preserving information in the compressed representation t becomes restricted (Gerardi et al., 2024; Park et al., 2025; Saoulis et al., 2025). Machine learning-based compression has now received sig- nificant attention across the field of SBI (Charnock et al., 2018; Prelogovi ́ c & Mesinger, 2024; Lanzieri et al., 2025; Lehman et al., 2025). The highly flexible parametrisation of deep neural networks ensures that they can learn globally valid, expressive compressed representations. Unlike in unsupervised compression (e.g. varia- tional autoencoders), the goal in SBI is to preserve as much dis- criminative power in the compressed statistics t with respect to the model parameters m. Jeffrey et al. (2021) demonstrated that the optimal objective function to achieve this amounts to training the compression model and NDE together to predict the posterior density. Specifically, given a compression model F θ and a density estimation model q φ , the networks can be trained jointly end-to-end on the neural poste- rior estimation (NPE) objective function from Equation (14): L =−E (m,D)∼D [logq φ (m| F θ (D))].(23) One can then use the trained model F θ to produce informative, non- linear compressed summaries of the observation for downstream tasks, or use the chained models q φ (m | F θ (D)) for posterior es- timation directly. For the purposes of this work the two networks can be thought of as a single model performing compression and posterior inference jointly, as in Figure 2. The remaining challenge is to design a compression architec- ture, F θ , that is capable of extracting the relevant information from multi-station seismic observations D. 4.2.3 Deep learning architecture for multi-station full-waveform data Deep learning for processing multi-station seismic data has received significant attention. A broad range of model architec- tures have been explored, such as graph neural networks (GNNs McBrearty & Beroza, 2022, 2023; Si et al., 2024a; Hourcade et al., 2025), transformer-based architectures (M ̈ unchmeyer et al., 2021; Song et al., 2025; Huang & Zhang, 2026), Fourier neural opera- tors (Sun et al., 2023), and aggregation across foundation model embeddings (Si et al., 2024b). In this work, we design an archi- tecture inspired by M ̈ unchmeyer et al. (2021), which combines a per-station CNN-based feature extractor with a transformer-based approach to perform multi-station feature aggregation. A simplified illustration of the architecture is presented in Fig- ure 2. Each station’s three-component seismogram is first processed independently by a shared CNN, which extracts station-specific temporal features while enforcing equivariance across the seismic network. To stabilise training and preserve physically meaningful amplitude information, waveforms are normalised per trace and augmented with an explicit log-amplitude feature (M ̈ unchmeyer et al., 2021). The resulting feature sequences are combined with sinusoidal temporal embeddings and station-location embeddings derived from the array geometry, enabling the model to reason jointly about waveform content, timing, and spatial configuration. Aggregation across stations and time is performed using an axial transformer (Ho et al., 2020). Rather than applying full self- attention over all station tokens as in M ̈ unchmeyer et al. (2021), the model factorises attention into separate station-wise and time- wise operations, preserving the temporal structure within the per- station CNN-extracted features. The attention mechanism in the transformer architecture is invariant under permutations of stations, and can be trained with a varying, arbitrary number of stations. It is therefore a highly generalisable feature extractor, that could in principle be applied to varying station configurations. As stated in the section above, this network is trained end-to- end with a NDE head network, which models the posterior distri- bution. This relies on a compressed summary of the seismic ob- servations, t, to model perform density modelling logq φ (m | t). We therefore perform aggregation of the transformer embeddings via learned query tokens, and mean pool them into a single fixed- dimensional embedding. This is passed through a small feedfor- ward neural network to produce the final compressed summary t. The exact implementation can be found in the provided software. 4.3 Incorporating 3-D Earth structure uncertainty Although we consider 1D Earth model uncertainty that affects all stations in the same manner, in many practical applications, it is often necessary to handle 3-D Earth structure uncertainty in mo- ment tensor inversions. To first order, the lateral heterogeneities in- duce station-specific timeshifts in the waveforms (as used for e.g. teleseismic station corrections in the ISC-EHB catalogue; Engdahl et al., 1998). Prior work often pre-processes the observations by aligning them with the best-fitting synthetics, often by maximising per-station correlations. Ideally, though, we should incorporate this 3-D Earth structure uncertainty in moment tensor solutions them- selves. One approach treats station-specific time-shifts as a proxy for all unresolved Earth structure uncertainties and includes them as extra parameters in the inversion (Hu et al., 2023, 2025). These “nuisance” parameters can then be marginalised over for the de- sired focal mechanism (following Equation (2)). However, as noted in Section 2, inference over the joint distribution of source param- eters and nuisance parameters increases the dimensionality of the problem, resulting in significantly increased computational costs for traditional MCMC-based inference procedures (for instance, Hu et al. 2025 runs over 5× 10 6 forward model evaluations per inversion). The approach taken in Hu et al. (2023) is useful in a likelihood-based analysis because attempting to pre-marginalise the likelihood over station-specific time-shifts (as done for 1-D Earth structure uncertainty in Equations 4 - 8) can further degrade the accuracy of the inversion. However, in theory SBI allows us to add these station-specific time shifts as extra nuisance parameters that we simulate in the prior predictive model of Equation (11). These can then be marginalised along with the 1-D Earth structure uncertainty with the flexible empirical density estimation formula- tion. In practice, our optimal score compression formulation in Section 4.2.1 relies on a relatively accurate Gaussian likelihood, which makes this challenging. However, the deep learning-based compression can straightforwardly handle arbitrary station-specific 8A. A. Saoulis et al. time-shifts as an extra source of uncertainty. For the real data exam- ples in Section 7, we therefore modify the prior predictive model for data generation to include time-shift nuisance parameters: m∼ p(m),Ω∼ p(Ω), t∼ p(t), D∼ p(D| m, Ω, t), (24) where each t i is a station-specific time-shift parameter drawn from an assumed prior distribution p(t). These nuisance parameters ac- count for uncertainties in absolute timing and are marginalised over implicitly during training. 5 EVALUATION 5.1 Covariance goodness-of-fit Given the Gaussian approximation adopted in Equation (8), it is important to assess whether the inferred covariance structure provides an adequate statistical description of the theory errors. If the Gaussian assumption holds, the χ 2 statistic χ 2 = (D− g(m, Ω)) T C −1 (D− g(m, Ω))(25) follows a χ 2 distribution with n dof = dim(D) degrees of freedom, or equivalently a reduced statistic χ 2 red = χ 2 /n dof . This offers a simple and computationally inexpensive diagnostic of the approxi- mation quality. In practice, χ 2 values can be computed for an ensemble of synthetic observations from the prior predictive model. The empir- ical distribution of χ 2 values may then be compared against the analytic expectation. Then, any systematic discrepancies between the observed and expected χ 2 red distribution, such as skewed or heavy-tailed behaviour, indicate a breakdown of the Gaussian co- variance approximation (Tarantola, 2005). This comparison can be performed visually using the reduced-χ 2 distribution, or quantita- tively via the Kolmogorov–Smirnov (KS) statistic. The KS statistic measures the maximum absolute difference between the empirical and theoretical cumulative distributions, D KS = max F emp − F χ 2 red .(26) 5.2 Moment tensor inversion quality To assess whether a given moment tensor inversion approach produces posteriors consistent with the true source parameters, we quantify both bias and uncertainty in the inferred solutions. For each focal mechanism parameter, we compute the posterior mean and measure the absolute bias relative to the known true value. In addition, we compute the posterior standard deviation for each pa- rameter, which provides a measure of how tightly each approach constrains each focal mechanism parameter. More discriminatory are calibration tests. These check whether the inferred posteriors are statistically consistent with the true solution over many repeated artificial inversions. This enables a practitioner to identify whether an inversion procedure is system- atically biased or fails to capture the true uncertainty in the poste- rior. Approaches for evaluating posterior calibration have received a lot of attention in the field of SBI (Talts et al., 2018; Hermans et al., 2021; Lueckmann et al., 2021), since calibration is an essen- tial pre-condition for reliable and trust-worthy posterior inference. Estimating empirical coverage in high-dimensional parameter spaces can be computationally demanding, as many existing ap- proaches require explicit density estimation in order to construct credible regions around the true parameters for quantile evaluation (e.g. Hermans et al., 2021). Lemos et al. (2023) introduced an al- ternative method, termed Tests of Accuracy with Random Points (TARP), which enables substantially more efficient coverage esti- mation. Rather than explicitly computing credible regions, TARP evaluates the fraction of posterior samples lying between a ran- domly drawn reference point and the true model parameters. Re- peating this procedure across multiple inversions yields an empir- ical distribution of credibility levels, which can be directly com- pared to the expected uniform distribution under perfect calibra- tion (see Saoulis et al., 2025, for a more detailed explanation, with examples for moment tensor inversion). In this work, we employ TARP to assess whether both Gaussian-likelihood inversions and SBI-based methods produce well-calibrated posterior distributions. 6 SYNTHETIC EXPERIMENTS 6.1 Experimental details 6.1.1 Forward modelling For our synthetic experiments, we use the source–receiver ge- ometry shown in Figure 3. This configuration corresponds to the network geometry of the Long Valley Caldera event studied later in this manuscript (LV2; see Section 7.1). We perform a suite of synthetic tests with moderate magnitude events, M W ∼ 5, where theory errors dominate over data errors. We follow the experimental and processing setup of Pha . m & Tkal ˇ ci ́ c (2021). All Green’s functions are computed using the Com- puter Programs in Seismology code (Herrmann, 2013). Unless oth- erwise stated, we use a source depth of 5 km and all synthetics are band-passed filtered between 20-50 s. All data are sampled at 1 Hz, and we select a time window of 200 s after the origin time for all stations. We use the SoCal 1-D Earth model as a fiducial Earth model Ω ∗ (Dreger & Helmberger, 1990), and follow Pha . m & Tkal ˇ ci ́ c (2021) in parametrising the Earth model prior p(Ω) in terms of a perturbation parameter κ. Starting from Ω ∗ , we generate prior samples by independently applying fractional perturbations to the V p , V s , and layer depths in each layer. The parameter κ speci- fies the mean percentage amplitude of these perturbations, such that larger κ corresponds to greater deviation from the fiducial model. For each experiment, we compute N = 300 sets of Greens functions for different Ω (n) ∼ p(Ω). We may then compute local estimates C t (m) using the Monte Carlo sampling procedure spec- ified in Equation (10). Forallexperiments,weparameterisetheseis- mic source using the symmetric moment tensor m= (M x ,M y ,M z ,M xy ,M xz ,M yz ),andadoptauni- form prior over its six independent components, p(m) = U [−4 × 10 16 , 4 × 10 16 ]. All density modelling and sampling algorithms are performed in the moment tensor component basis. For evaluation and visualisation, we use the Tape & Tape (2015) parametrisation to represent the isotropic and compensated linear vector dipole (CLVD) components δ and γ, along with the moment magnitude measure M w and best fitting double-couple source orientation angles strike, dip and rake. The sensitivity kernels are computed using 5-point derivative stencils. We perform iterative least squares to find an estimate of m MLE for every event for the Gaussian likelihood approach and score compression-based SBI approach. This is achieved by re- peatedly estimating the local sensitivity kernels G(m) and covari- ance C t (m), and applying the MLE estimation formula in Equa- tion (21). We perform 8 iterations, and select the iteration m MLE Theory errors in moment tensor inversion with SBI9 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 Pha . m & Tkal ˇ ci ́ c (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). that minimises the χ 2 (misfit) of the event. In practice, this was generally not the final iteration as we found the optimisation prob- lem was not perfectly convex or stable. 6.1.2 Gaussian likelihood sampling We draw comparisons between the results of moment tensor inversions using the approximate Gaussian likelihood against those using SBI. We compute a local estimate of the Gaussian covari- ance C(m MLE ) as described in Equation (15) (or Equation (17) for real data), and use MCMC to draw samples from the posterior. We use emcee (Foreman-Mackey et al., 2013) to implement a simple, embarrassingly parallel Metropolis-Hastings sampling strategy. We manually optimised the step-size of the burn-in and final sampling stage for each experimental setup to ensure that the posterior so- lutions converged. We parallelise this over 30 cores (Intel Xeon Gold 6230 @2.10 GHz) and run approximately 100, 000 forward model evaluations for each Gaussian likelihood-based inversion. Final posterior samples are computed with a per-walker burn-in fraction of 0.4. 6.1.3 SBI implementation details For the score compression approach, we use the prior trunca- tion procedure in Equation 22 with N = 15, finding this struck a good balance between ensuring well-calibrated but discrimina- tive posteriors for all experiments (see Saoulis et al., 2025, who found similar results). We run 8, 000 forward model evaluations from the prior predictive model Equation (11), and use the result- ing compressed datasetD = m, t to train the NDE to perform posterior estimation (with 10% of data held back as a validation set for model selection). We use a MAF architecture with 5 transfor- mation layers for the NDE, finding it trained quickly on CPUs and was sufficiently flexible to accurately model the target posteriors. Training was performed until a stopping criteria was reached on the validation set. The total duration of MLE estimation of m, dataset generation and training took around 15 minutes on average for each inversion, using parallelisation over 30 cores. For the deep learning-based compression approach, we per- formed a short manual optimisation to improve the performance of the architecture. The per-station 1-D CNN processes 3 component, 200 s duration seismic traces. The network consists of a stack of 6 1-D convolutional layers with uniform kernel sizes, with early strided convolutions used to progressively downsample the tempo- ral dimension while increasing the feature dimensionality. This de- sign yields per-station features with a temporal resolution of around 30 time samples and 128 feature channels per station. We then add sinuisoidal embeddings that encode the station location for each station feature block, and time embeddings encoding the time in- dex of each of the 128-dimensional feature channels. We build the axial transformer with multiple layers consisting of multi-head attention (4 heads) with model dimension of 256. The transformer learns 8 query tokens, which query the processed seis- mic features over 4 attention layers. After the 4 layers, these queries are mean pooled and processed by a 2 layer feedforward network to produce a final summary embedding t with dimension 128. This summary is passed to a rational-quadratic neural spline flow (RQ- NSF) with 5 transformation layers to perform density estimation. We found that this NDE architecture significantly improved per- formance over the MAF used for the score compression approach, likely as a result of facilitating improved gradient flow back through the feature extractor network. The deep learning-based inference model was trained us- ing 100, 000 forward model evaluations from the prior predictive model in Equation 11, with 10% held out as validation data. Train- ing was performed for 300 epochs with a batch size of 128, and the lowest validation loss model was selected. Dataset generation took approximately 10 minutes and training took approximately 12 hours on a single NVIDIA RTX A6000 GPU. Note that, once trained, this model was globally applicable for any moment tensor m, unlike the other two approaches. All models were trained using the Adam optimiser. For the simple MAF NDE for score compression, we used a learning rate of 1× 10 −4 . For the deep learning model, we used a linear learning rate ramp-up to 1× 10 −4 combined with a cosine decay schedule to 1 × 10 −5 , finding that annealing the learning rate toward the 10A. A. Saoulis et al. Figure 4. A comparison between the analytically expected χ 2 red 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 Gaus- sian 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 likeli- hood becomes a very poor approximation. end of training led to slightly improved models. We used a small weight decay of 1× 10 −4 and a dropout rate of 10%, and found that models generally did not overfit during training. All the precise architectural details can also be found in the provided software. 6.2 Probing the validity of the Gaussian assumption We begin by using the χ 2 test from Equation (25) to evaluate the accuracy of a Gaussian approximation of the theory errors, en- coded by the Gaussian covariance C t (m). This provides a straight- forward way of evaluating the accuracy of the Gaussian assump- tion without needing to perform inference. We estimate the theory covariance C t (m ∗ ) for a representative moment tensor m ∗ un- der varying levels of Earth structure uncertainty, κ ∈ [0.1, 1, 3, 5] and compute the χ 2 red distribution for an ensemble of 2000 random events using the LV2 source-receiver geometry in Figure 3a. Figure 4 shows that under even minor Earth structure pertur- bations, on the order of 1% of the velocity model, the Gaussian assumption of a local estimate of C t (m ∗ ) breaks down. This pro- vides evidence that under minor Earth structure uncertainty, Gaus- sian approximations of the theory errors may lead to biased poste- rior inference. We present some examples of the effect of velocity struc- ture perturbations in Figure 5, alongside the per-trace χ 2 red statis- tic. We find that, intuitively, stations at greater epicentral distances are more sensitive to velocity model uncertainties. The resulting time-shifts are not always perfectly modelled by the Gaussian co- variance C t , leading to occasional traces with higher χ 2 red . How- ever, traces with large changes to waveform amplitude (and low time-shift relative to the fiducial observation) appear to cause the most severe modelling issues. One potential explanation is that the covariance estimate C t (m) is dominated by the variance induced by time-shifts. This leads to higher modelled variance in regions of the waveform with higher gradients (i.e. when particle velocity dx/dt is large), and lower variance in regions where the gradient Figure 5. Examples of the per-trace χ 2 red statistic under Earth structure per- turbations for κ = 3%. Higher values of χ 2 red are denoted by a brighter- copper colour, with significant departure from those expected under a Gaus- sian likelihood. Even under modest perturbations in the Earth structure, there are significant non-Gaussian effects in the observations. approaches zero (i.e. the turning points of the waveform). This be- haviour is demonstrated in Figure 1(b). The modelled covariance, which is a first-order approximation that does not model higher or- der statistics, performs poorly on rarer traces with large amplitude differences and low time-shifts from the fiducial. In summary, we find that even under minor velocity model uncertainty, a Gaussian covariance provides a poor statistical fit of the observed variability in the observations. However, the degree to which this matters in a practical sense is best tested through actual moment tensor inversions. 6.3 Synthetic inversions 6.3.1 Case study 1: Gaussian likelihood vs. SBI under high uncertainty We run a large number of synthetic inversions under varying levels of Earth structure uncertainty. To begin with, we explore how each of the three approaches performs under moderate (κ = 5%) uncertainty in the Earth velocity model. In Figs. 6 and 7, we show two examples of moment tensor inversion with each of the three approaches: Gaussian likelihood- based MCMC, SBI with optimal score compression and SBI with a deep-learning based compression model. Figure 6 provides an ex- ample where the Gaussian likelihood approximation is sufficiently accurate to provide a good solution, leading to posterior contours that are consistent with the artificial solution. The Gaussian likeli- hood approach yields posteriors with the lowest uncertainty, which we interpret below. The ML-based compression approach yields significantly tighter posterior contours than the score compression algorithm. The deep learning compression architecture is therefore capable of extracting more information from an observation than the linearised, gradient-based score compression, indicating that there are non-negligible non-linearities in the forward problem. However, Figure 7 presents a synthetic example where the Gaussian likelihood assumption leads to failure for both the Gaus- Theory errors in moment tensor inversion with SBI11 Figure 6. Posterior inference on an artificial example comparing the results using a Gaussian likelihood (red) against the two SBI frameworks: score com- pression 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 γ,δ,M W 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). sian likelihood MCMC approach, and the score compression SBI approach. Both of these approaches yield moment tensor solutions which are inconsistent with the artificial event. In this instance, the iterative least squares procedure for estimating a local theory co- variance C t (m ∗ ) fails to converge since the Gaussian likelihood is not sufficiently accurate. Since both approaches rely on an ac- curate m ∗ (the score compression SBI approach through the prior truncation procedure in Equation (22)), a severely biased estimate of m ∗ causes these approaches to fail. On the other hand, the ML- based compression algorithm has no such limitations, and yields a posterior that is consistent with the artificial event. We systematically evaluate how each moment tensor inversion approach performs. We perform 400 artificial moment tensor inver- sions, with each approach, drawing moment tensors from the same uniform prior. We then perform coverage testing on each of the re- sulting posterior sets using TARP, as described in Section 5. The results of this coverage testing are presented in Figure 8. Figure 8 demonstrates that only the SBI approaches can produce posteriors that are consistent with the true moment tensor solutions, even under high Earth structure uncertainty. The main issue uncov- ered by this test is that the Gaussian likelihood approach yields severely overconfident, and mildly biased, moment tensor posteri- 12A. A. Saoulis et al. 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. ors. On the other hand, the score compression-based SBI approach produces much better calibrated posteriors with only a small degree of bias. We can interpret these results with the help of Figs. 6 and 7. The Gaussian likelihood approach yields much lower uncertainties, but these often fail to cover the true solution, leading to overcon- fidence. On the other hand, the score compression-based SBI ap- proach produces much greater uncertainties (as in Figure 6) which are generally consistent with the true solution. However, both ap- proaches demonstrate inference bias: Figure 7 provides an example of how inaccuracies in the Gaussian likelihood can introduce this bias. The ML-based compression SBI approach exhibits a mild de- gree of under-confidence, producing conservative uncertainty es- timates that tend to overestimate its own uncertainty (by around 10%). This is not ideal since it indicates that the model could produce tighter constraints on the moment tensor: i.e., the deep- learning based compression algorithm is extracting more informa- tion from the observations that are not fully taken advantage of by the NDE head. From a practical perspective, however, this under- confidence is far preferable to overconfidence, in the sense that at Theory errors in moment tensor inversion with SBI13 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 over- confident (∼ 50%) and slightly biased posterior contours, resulting from inaccuracies in the approximate likelihood. The score compression SBI ap- proach 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 conser- vative (∼ 10%) but nonetheless trustworthy and unbiased. The ±[1, 2]σ uncertainty regions are estimated from bootstrapping the coverage results 50 times. least one can be sure the modelled posterior covers the true solu- tion. There is a wide body of literature that has worked towards ensuring SBI approaches yield conservative or calibrated posteri- ors, rather than over-confident posteriors (Delaunoy et al., 2022; Hermans et al., 2021). We also systematically evaluate the modelled uncertainties and biases of each approach over the 400 artificial inversions. The results for each approach are shown in Table 1. The Gaussian like- lihood approach produces low uncertainties, to some degree driven by its over-confidence. We also find that a Gaussian likelihood pro- duces systematic posterior biases, particularly for M W where the magnitude is significantly under-estimated for a large number of events. The biases on the focal mechanism parameters γ and δ are minor, however, and we did not find much evidence of particular mechanisms being over- or under-represented. Between the SBI approaches, we find that the deep-learning based compression produces significantly tighter (around a factor of×2) constraints on the source parameters. This advantage stems from the fact that the deep-learning model is capable of extracting more information from the observations than score compression. In addition, since the deep learning-based compression algorithm does not rely on an inaccurate Gaussian likelihood, it produces pos- teriors with lower systematic bias, particularly for the focal mech- anism parameters γ,δ, and M W . 6.3.2 Case study 2: Gaussian likelihood inversion quality under different experimental configurations We next explore how the Gaussian likelihood assumption fares under varying experimental configurations. We test the following scenarios with input source models randomly sampled from the prior: (i) differing levels of Earth model uncertainty, (i) introducing shorter period data in the filtering window, (i) inversions with a more balanced station configuration. We begin by exploring model calibration under three levels of Earth structure uncertainty, κ ∈ [1, 3, 5]% from scenario (i). Fig. Figure 9 shows the results of 300 artificial inversions under each of these uncertainty scenarios. We find that while even under mild un- certainty (κ = 1%), a Gaussian likelihood yields significantly over- confident posteriors that underestimate the moment tensor compo- nent uncertainty by around 30%. This increases to a 50% underes- timation by κ = 5%. In addition, we find the Gaussian likelihood approach leads to systematic bias in the moment tensor solutions across all levels of Earth structure uncertainty. This is shown by the asymmetry of the calibration curve, reflecting the biases in the focal mechanism uncovered in Section 6.3.1. We next consider shorter period data by filtering between 6- 50s from scenario (i). We expect that this should exacerbate mod- elling errors due to stronger sensitivity to smaller scale structures and, more importantly, lead to more non-Gaussian effects in the theory errors. Figure 9 demonstrates that the Gaussian likelihood assumption is significantly degraded in this context, leading to poorer calibration relative to the (already problematic) baseline 20 - 50 s period band. We find that the Gaussian likelihood assumption now underestimates its moment tensor component uncertainties by around a factor of×2, and there is a persistent and significant bias. We show a typical example of posterior inference in Fig. S1. Finally, we repeat this procedure for scenario (i), performing repeated inversions using a balanced network with good azimuthal coverage around the event, with the exact configuration shown in Fig. S2. Figure 9 finds a similar degree of overconfidence to the baseline case (again around 50% overconfident). This indicates that good station coverage does not mitigate against a poorly specified likelihood in terms of producing well-calibrated moment tensor un- certainties. However, it is important to note that in such scenarios, the moment tensor is recovered with higher precision. Better station coverage is therefore still desirable, since incorrect uncertainty es- timation may not affect the overall interpretation of the focal mech- anism solution. We can nonetheless infer that the limitations with the likelihood-based approach uncovered in this section cannot be attributed to poor station coverage. By comparison, we find that the SBI approaches produce reli- able posteriors under scenarios (i) and (i), with calibration curves for the SBI approaches shown in Figs. S3 & S4. However, including higher frequency information in the waveforms, as done in scenario (i), exacerbates the bias in the score compression SBI approach. This is likely for the same reasons discussed in Section 6.3.1 (i.e., biased MLE estimates of m due to an incorrect Gaussian likeli- hood), though calibration is still much better than using the Gaus- sian likelihood directly. The ML-based compression approach re- mains reliable, unbiased and slightly conservative for most of the scenarios tested, though we find its performance also degrades on the short period data. We provide some comment on this in Fig. S4. We also test to what extent model misspecification affects the ML-based compression approach in Fig. S4, and find that while it 14A. A. Saoulis et al. MethodCompressionMetric γδ M w (10 −2 )strikediprake SBI Optimal score σ1.512.100.962.561.912.94 Bias-0.170.49-0.230.05-0.240.29 Deep learning σ0.710.840.371.250.951.52 Bias-0.060.06-0.01-0.230.06-0.26 Gaussian- σ0.660.870.450.790.711.04 Bias-0.150.15-0.320.07-0.180.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. M w values are reported in units of 10 −2 , all other values in units of degrees. induces predictable miscalibration, the method does not catastroph- ically fail. We present an in-depth analysis showing the posterior accura- cies over the 300 synthetic inversions for several of the scenarios tested in Figs. S5, S6 & S7. We also repeat the posterior width and bias analysis of Table 1 for the κ = 3% uncertainty and shorter period inversions in Tables S1 & S2. 6.3.3 Case study 3: shallow, isotropic events Shallow isotropic sources have proven theoretically challeng- ing to fully resolve using regional surface waves, even when ac- counting for uncertainty in the Earth structure (K ˇ r ́ ı ˇ zov ́ a et al., 2013; Hu et al., 2023; Chiang et al., 2025; Pha . m et al., 2025). This is because fundamentally different source mechanisms, such as isotropic and vertical CLVD, could produce nearly identical re- gional surface waves (Kawakatsu, 1996; Ford et al., 2012; Hu et al., 2023). The ambiguity could result in much larger solution uncer- tainty for non-double-couple source types than for seismic sources at greater depths. Here we test SBI’s ability to retrieve truthworthy uncertainty estimates in the challenging setting. To do so, we study a synthetic scenario using the same station configuration as in the LV2 event (Figure 3) and consider highly isotropic sources with a source depth of 500 m. We set the Earth structure uncertainty to κ = 5%, and test each of the three approaches on 10 randomly sampled isotropic events. The results of this analysis are shown in Figure 10. We focus only on the lune plot in order to evaluate how well each approach resolves the focal mechanism components of each event. We find several instances where the Gaussian likelihood-based approach (and, by extension, the score compression SBI approach) signifi- cantly underestimates the isotropic component of the event. These can be seen most clearly for lunes (b), (e), (h), and (i), where the maximum a posteriori point for both a Gaussian likelihood-based analysis and the score compression approach substantially under- estimates the isotropic component. In contrast, the ML-based com- pression approach to SBI produces posterior contours that more faithfully recover the isotropic component of the event. It is worth noting that the uniform prior of the moment tensor components assumed here leads to very low prior density toward the pure isotropic sources. This is true for the Tape parametrisation as well (Tape & Tape, 2015). This induces a preference for smaller isotropic components, particularly when the observational uncer- tainties are poorly modelled by the likelihood. Prior volume effects may therefore explain the strong isotropic—double-couple trade- off visible for many inversions in Figure 10. Nonetheless, we may reasonably conclude that the more accurate (implicit) likelihood modelling performed by the ML-based approach ensures improved posterior solutions, regardless of the prior chosen. These prelimi- nary results, not yet fully explored in this study, are promising for future real data applications in the context of non-proliferation seis- mology (Ford et al., 2012; Chiang et al., 2018). 6.4 A note on computational costs As described in Section 6.1, the ML-based compression al- gorithm had a significant drawback relative to the score compres- sion approach: for each source-receiver configuration, it cost 12 GPU hours to train the deep learning model, in conjunction with the larger number of forward simulations required for training. However, once trained, the ML-based approach is highly ef- ficient, and can perform posterior inference on hundreds of events within seconds. This advantage was felt directly in the synthetic tests in Section 6, where for each experimental configuration the ML-based compression training and inference time was fixed at 12 hours. However, for the other two approaches the local covariance m MLE estimation and all inference overheads had to be repeated for every event. For the 300 synthetic inversions required to produce reliable coverage statistics, this amounted to 3-4 days of compute using 30 CPUs for each experimental configuration. For each ex- periment this totaled ∼ 25× 10 6 forward model evaluations for the Gaussian likelihood approach and ∼ 3.5× 10 6 for the score compression approach; 1 to 2 orders of magnitude more than the 0.1× 10 6 forward model evaluations for the ML-based approach. It would be quite restrictive to need to re-train a deep-learning compression model for each event under study. Fortunately, while we do not investigate this here, prior studies using similar archi- tectures have demonstrated the ability to generalise across varying source locations and source-receiver geometries (e.g., M ̈ unchmeyer et al. 2021; Hourcade et al. 2025). A single ML model could there- fore be trained to perform inversions in a given region even under varying source-receiver geometries. In principle, this could enable extremely efficient Bayesian source mechanism determination for whole earthquake catalogues, with several order of magnitudes re- duction in computational cost relative to standard methods. 7 APPLICATION TO REAL DATA We repeat the inversion procedures for the two real events in Figure 3. All data are processed as in the synthetic experiments, us- ing a 20 - 50 s period filtering band and a 200 s window following the origin time of the event. The per-station component noise vari- ance is estimated by using the window prior to the event, and then used to rescale the data covariance matrices C d . For the Gaussian likelihood and the score compression SBI Theory errors in moment tensor inversion with SBI15 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. approaches we use prior moment tensor solutions to estimate station-specific time-shifts and align the observations with the syn- thetics. We then performed iterative least squares for a best fitting moment tensor estimate m MLE . Finding that this did not yield per- fectly stable or recoverable best fitting solutions, we also performed a single round of MCMC, which mildly improved (but did not en- sure) the stability of the inversion process. As this did not occur during the synthetic tests, we attribute these issues to further un- modelled theory errors. For the ML-based compression approach, we perform no alignment and implicitly marginalise over station specific time- shifts using the approach detailed in Section 4.3. We use a simple, heuristic prior over time-shifts that absorbs errors in the absolute origin time with a Gaussian distribution, and per-station corrections that are uniformly distributed: t 0 ∼N (0,σ 2 ot ), δ i ∼U (−∆ sc , ∆ sc ), t i = t 0 +δ i . (27) We use p(t;σ ot = 2 s, ∆ sc = 3 s), which is empirically moti- vated by the best fitting time-shifts found in prior regional stud- ies (Pha . m & Tkal ˇ ci ́ c, 2021; Hu et al., 2025). We did not perform any detailed analysis regarding this choice, but found that the ML- based compression solutions were not significantly impacted by small changes to σ ot and ∆ sc . 7.1 Long Valley Caldera M W 4.9, 11/22/1997 We reanalyse a well-studied event at the Long Valley Caldera in California, USA at 17:20:35 on 11/22/1997, referred to as LV2 in Minson & Dreger (2008). Prior studies have indicated a significant non-double-couple component likely related to volcanic processes at the caldera (Dreger et al., 2000; Pha . m & Tkal ˇ ci ́ c, 2021). We use the SoCal 1-D Earth model as in the synthetic examples, and use κ = 5% as the fractional uncertainty in the layer depths and layer velocities V p and V s . We use a theory error regularisation factor of ε = 0.005 in Equation (17). We present the results of each of the three approaches: a Gaussian likelihood based inversion, SBI with score compression, and SBI with ML-based compression in Figure 11. We also show the focal mechanism solution found in Pha . m & Tkal ˇ ci ́ c (2021), who introduced the Gaussian likelihood methodology we reproduce here in Figure 11c. As in the synthetic examples, we find close agreement be- tween the Gaussian likelihood and score compression-based SBI approach for the maximum a posteriori solution of the focal mecha- nism. Again, however, the SBI approach yields substantially higher uncertainties, which our results in Section 6 suggest are a better characterisation of the true uncertainty. On the other hand, the ML- based approach yields slightly different moment tensor solutions, with a fractionally higher isotropic component than the other two approaches. We found this higher isotropic component persisted across ML training runs, as well as when time-shifts were man- ually corrected as in the two other approaches (i.e. by training a model with p(t;σ ot = 0 s, ∆ sc = 0 s) and aligning the observation with the best fitting synthetics). All solutions yield similar trade- offs between component pairs, even when they slightly disagree on the best solution or uncertainty widths (note e.g. a strong positive tradeoff between the magnitude M w and isotropic component γ). We show a random subset of 10,000 posterior predictive checks in Figure 11a. These posterior predictive checks show that both SBI approaches yield plausible fits to the data. For some sta- tions (e.g. KCC), the ML-based compression approach yields better fits to the data, while for others (e.g. ORV) the score compression approach seems to fit better. We verify quantitatively in S3 that the ML-based moment tensor ensemble produces synthetics that are competitive with or improve over the the other two approaches in terms of station-component power content, envelope misfit, and χ 2 fit, which all provide somewhat independent probes of waveform fit. We interpret this as providing evidence that all approaches pro- vide reasonable agreement with the data. 16A. A. Saoulis et al. 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. 7.2 Zagreb, Croatia, M W 5.3 2020/03/22 The 2020 Zagreb earthquake had a significant human toll, causing one fatality, many injuries, and reported economic losses exceeding several billion euros (Novak et al., 2020; Atali ́ c et al., 2021). The focal mechanism of the earthquake has been analysed in several prior works (e.g., Marku ˇ si ́ c et al., 2020; Hu et al., 2025). Hu et al. (2025) used a station-specific time-shift approach to treat un- resolved Earth structure uncertainty, and retrieved a double-couple source-type mechanism largely consistent with the tectonic back- ground in the region. To aid interpretation, we use stations utilised in their analysis, and use the same 2-D composite model (made up of three different 1-D Earth models depending on the station location) used in their work. We present details of the composite model in Fig. S8, and defer to prior work for a more detailed exposition of the tectonic and geophysical background of the region (e.g., Herak et al., 2009; Stip ˇ cevi ́ c et al., 2011, 2020; Hu et al., 2025). We apply κ = 3% perturbations to the velocity models, motivated by the approximate uncertainty levels discovered in Stip ˇ cevi ́ c et al. (2020). These are applied to each of the 1-D Earth models independently. This leads to a problem of greater complexity (and realism) than the LV2 event forward modelling, since structure perturbations are now only az- imuthally (rather than globally) correlated. The results of all three approaches are presented in Figure 12. We find close agreement of a double-couple mechanism with the Gaussian likelihood approach and the score compression SBI ap- proach with Hu et al. (2025), though there are minor differences in fault orientation and magnitude. We found some instability in the exact maximum likelihood solution m MLE even with a large theory regularisation term ε = 0.05, indicating significant unmodelled er- rors. Although broadly compatible with the other two approaches and with the solution of Hu et al. (2025), the ML-based compres- sion approach shows evidence of a minor CLVD component to the event, and prefers different strike and rake orientations (with dis- crepancies between 5−10 ◦ ). The ML-based compression approach also has larger uncertainties than in the previous examples, which we hypothesise results from the more challenging to model Earth structure uncertainties, with independent 1-D Earth models causing uncorrelated errors at different stations. We verified that training with a larger dataset of simulations reduced these uncertainties, so conclude this effect is related to increased inherent difficulties in modelling the composite model uncertainties. The ML-based compression yields a plausible solution — it is consistent with a very low (∼ 2 ◦ ) CLVD component within the 95% confidence interval. It is also possible that nearby faults with differing orientations may have been activated during or con- tributed to the main shock, introducing a small modelled non- double-couple component. More thorough investigation is war- ranted to better understand this fractional CLVD component. While the ML-based approach performs a more comprehensive treatment of the uncertainties (i.e., by marginalising over station-specific time-shifts and not making any assumptions about the form of the errors), it is not possible to rule out that unmodelled effects may in- troduce bias in the resulting solution. As such, we encourage some caution in further interpretation of this result. Again, a quantitative analysis of the posterior predictive checks in Table S4 suggests that all approaches provide compa- rable fits to the data. However, we find that the posterior predictive checks in Figure 12a are in some sense under-dispersive and do not fully cover the observations. This is an indication of model mis- specification; i.e., the theory errors induced by our 1-D Earth struc- ture perturbations underestimate or do not capture the true mod- elling uncertainty. 8 DISCUSSION This work has focused on uncertainties arising from 1-D Earth structure, with azimuthal variations treated in a simplified manner Theory errors in moment tensor inversion with SBI17 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 Pha . m & Tkal ˇ ci ́ c (2021). through station-specific time shifts. An important next step is to rig- orously assess how the different inversion approaches perform on more realistic data. The score-compression SBI approach is closely tied to a Gaussian likelihood, whose behaviour and robustness on real data is better understood. Nonetheless, 3-D Earth structure ef- fects, particularly at lower periods, will likely exacerbate the Gaus- sian likelihood mismodelling issues uncovered in this work. The reliability of the fully data-driven deep learning compres- sion remains unclear. In particular, unmodelled effects such as 3-D Earth structure and site response may lead to more severe out-of- distribution errors, which could underlie the minor moment tensor solution differences uncovered in Section 7 (Van Amersfoort et al., 2020; Schmitt et al., 2023; Pierre et al., 2026). A more comprehen- sive evaluation is therefore warranted, for example by testing the deep learning approach on 3-D synthetics that also include realistic observational contaminants. Future work should attempt to model more realistic Earth structure uncertainties, such as 3-D heterogenieties (as in e.g., Pha . m et al., 2024) and anisotropy. These could be derived from tomographic model ensembles (Chiang et al., 2025), for which SBI is better suited to deal with than likelihood-based analysis (as ex- plained in Section 2). Beyond Earth structure, there are many more modelling phenomena that SBI could be applied to: source time functions (Vall ́ e et al., 2011; St ̈ ahler & Sigloch, 2014), source lo- cation uncertainties (and receiver uncertainties in the case of ocean bottom seismometers; Stachnik et al., 2012; Trabattoni et al., 2020; 18A. A. Saoulis et al. 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). Lindner et al., 2023), composite or finite source models, or poten- tially even errors in the forward modelling itself. All these effects could be treated either as model parameters to be inverted for or “nuisance” parameters to marginalise over. We explored two extremes of compression for SBI: a generic compression algorithm applied on the raw full waveforms, and a fully ML-based data-driven approach. These both had limitations: the former approach was brittle under minor non-linearities, while the latter requires designing and training a tailored deep learn- ing architecture. In practice, the field of SBI has heavily relied on “two-step” compression approaches (Alsing & Wandelt, 2018; Als- ing et al., 2019; Cranmer et al., 2020; Gatti et al., 2024), which first extracts more stable physically motivated summaries from the observations before compressing them further. For full-waveform seismic observations, there are a wide range of options to con- sider: power spectra, time-frequency representations or wavelet co- efficients (Cesca et al., 2006; Fichtner et al., 2008; Vavry ˇ cuk & K ̈ uhn, 2012); peak amplitude, waveform envelope or coda statistics (Nakahara, 2008; Eulenfeld et al., 2022); arrival times, first-motion polarities and S to P amplitude ratios (Hardebeck & Shearer, 2002; Skoumal et al., 2024, see also Song et al., 2025 for a relevant ML- based example). Indeed, our purely data-driven approach could no Theory errors in moment tensor inversion with SBI19 doubt also benefit from some hybrid: for instance, including phys- ically motivated summaries (Makinen et al., 2024; Jeffrey et al., 2025) or explicit simulator feedback (e.g., Graikos et al. 2022; Chung et al. 2023; Holzschuh & Thuerey 2024; see Thuerey et al., 2021 for a review of methods), to incorporate information from the forward model during inference. This work utilised the neural posterior estimation (NPE) framework, which trains a model of the posterior. This allows for rapid posterior sampling, but is somewhat restrictive in that it bakes in the prior used for dataset generation. Neural likelihood estima- tion (NLE) effectively learns probabilistic emulator of the forward model (including all uncertainty effects), and once trained can be used to probe the effects of the prior on moment tensor solutions (Papamakarios et al., 2019). Likelihoods across independent mea- surements can also be multiplied, which could enable the appli- cation of SBI to joint inversions of full-waveform data with first- motion polarities (Hamidbeygi et al., 2023; Chiang et al., 2025) or InSAR data (Delouis et al., 2002; Weston et al., 2011). 9 CONCLUSIONS This manuscript has demonstrated that accurate Bayesian mo- ment tensor inversions under model uncertainty rely on well- characterised likelihoods. Standard Gaussian likelihood approx- imations introduce bias and overconfidence in the solutions. Simulation-based inference (SBI) produces better models of the likelihood using machine learning, enabling accurate and calibrated moment tensor inversions under challenging sources of uncertainty. We showed that Gaussian likelihood approximations break down under very minor velocity model uncertainty. Predictably, as the degree of velocity structure uncertainty increases from 1% to 5% more pathologies are introduced in the Gaussian likelihood- based moment tensor solutions. These problems are exacerbated when treating higher frequency data. We also showed that Gaussian likelihood inversions of shallow isotropic sources may significantly underestimate the isotropic component of the source. We introduced two frameworks for empirical likelihood mod- elling with SBI. One leveraged the generic and lightweight score compression algorithm, which uses a Taylor expansion of the for- ward model to reduce the dimensionality of full-waveform seismic observations. Score compression-based SBI produced much better calibrated moment tensor posterior solutions than a Gaussian like- lihood approach, albeit with very occasional biases. The second framework used a fully data-driven approach, using a deep learning model to compress the observations and model the posterior. The latter framework is significantly more flexible and can straightfor- wardly incorporate more complex sources of uncertainty. Impor- tantly, the flexibility and efficiency of deep learning-based SBI has the potential to enable rapid and robust Bayesian moment tensor inversion at the catalogue level. ACKNOWLEDGMENTS We thank Iva Dasovi ́ c and Kre ˇ simir Kuk at the Croatian Seis- mograph Network, as well as Jinyin Hu, for providing expertise and access to the data for the Zagreb 2020 earthquake. We are grateful to Davide Piras, Alessio Spurio Mancini, and Benjamin Joachimi for helpful discussions in preparation for this work. AAS was sup- ported by the STFC UCL Centre for Doctoral Training in Data In- tensive Science (grant ST/W00674X/1) and by departmental and industry contributions. AAS was also supported by the A. G. Lev- entis Foundation educational grant scheme. T.-S. P. acknowledges financial support from the Australian Research Council through a Discovery Early Career Research Award (DE230100025). DATA AVAILABILITY Python code used to produce the results of this study is avail- able at https://github.com/asaoulis/seismo-sbi. We used several widely available Python packages for processing the data and train- ing the neural networks: obspy (Beyreuther et al., 2010), PyTorch (Paszke et al., 2019), PyTorch Lightning (Falcon, 2019), sbi (Tejero-Cantero et al., 2020), and nflows (Durkan et al., 2020a). References Abercrombie, R. E. & Ekstr ̈ om, G., 2001. Earthquake slip on oceanic transform faults, Nature, 410(6824), 74–77. Akuhara, T., Shinohara, M., Yamada, T., Azuma, R., Hino, R., Obana, K., Takahashi, T., Fujie, G., Kodaira, S., Murai, Y., et al., 2025. Non-double-couple components of the 2024 noto earth- quake aftershocks: influence on focal mechanism estimation, Earth, Planets and Space, 77(1), 145. Alsing, J. & Wandelt, B., 2018. Generalized massive optimal data compression, MNRAS, 476, 60–64, doi: 10.1093/mnrasl/sly029. Alsing, J., Wandelt, B., & Feeney, S., 2017. Massive optimal data compression and density estimation for scalable, likelihood-free inference in cosmology, MNRAS, 000, 1–14. Alsing, J., Charnock, T., Feeney, S., & Wandelt, B., 2019. Fast likelihood-free cosmology with neural density estimators and ac- tive learning, Monthly Notices of the Royal Astronomical Soci- ety, 488(3), 4440–4458. Atali ́ c, J., Uro ˇ s, M., ˇ Savor Novak, M., Dem ˇ si ́ c, M., & Nastev, M., 2021. The mw5. 4 zagreb (croatia) earthquake of march 22, 2020: impacts and response, Bulletin of Earthquake Engineer- ing, 19(9), 3461–3489. Atlas Collaboration, 2025.An implementation of neural simulation-based inference for parameter estimation in atlas, Re- ports on Progress in Physics, 88(6), 067801. Beyreuther, M., Barsch, R., Krischer, L., Megies, T., Behr, Y., & Wassermann, J., 2010. Obspy: A python toolbox for seismology, Seismological Research Letters, 81(3), 530–533. Blom, N., Hardalupas, P.-S., & Rawlinson, N., 2023. Mitigat- ing the effect of errors in source parameters on seismic (wave- form) tomography, Geophysical Journal International, 232(2), 810–828. Cesca, S., Buforn, E., & Dahm, T., 2006. Amplitude spectra mo- ment tensor inversion of shallow earthquakes in spain, Geophys- ical Journal International, 166(2), 839–854. Charnock, T., Lavaux, G., & Wandelt, B. D., 2018.Auto- matic physical inference with information maximizing neural networks, Physical Review D, 97(8), 083004. Chiang, A., Ichinose, G. A., Dreger, D. S., Ford, S. R., Matzel, E. M., Myers, S. C., & Walter, W., 2018. Moment tensor source- type analysis for the democratic people’s republic of korea– declared nuclear explosions (2006–2017) and 3 september 2017 collapse event, Seismological Research Letters, 89(6), 2152– 2165. Chiang, A., Ford, S. R., Pasyanos, M. E., & Simmons, N. A., 20A. A. Saoulis et al. 2025. Bayesian inference for the seismic moment tensor us- ing regional waveforms and teleseismic-p polarities with a data- derived distribution of velocity models and source locations, Bul- letin of the Seismological Society of America. Chung, H., Kim, J., Mccann, M. T., Klasky, M. L., & Ye, J. C., 2023. Diffusion posterior sampling for general noisy inverse problems, in The Eleventh International Conference on Learn- ing Representations. Cranmer, K., Brehmer, J., & Louppe, G., 2020.The fron- tier of simulation-based inference, Proceedings of the National Academy of Sciences of the United States of America, 117, 30055–30062, doi: 10.1073/PNAS.1912789117. Dax, M., Green, S. R., Gair, J., Gupte, N., P ̈ urrer, M., Raymond, V., Wildberger, J., Macke, J. H., Buonanno, A., & Sch ̈ olkopf, B., 2025. Real-time inference for binary neutron star mergers using machine learning, Nature, 639(8053), 49–53. Deistler, M., Goncalves, P. J., & Macke, J. H., 2022. Truncated proposals for scalable and hassle-free simulation-based infer- ence, in Advances in Neural Information Processing Systems, vol. 35, p. 23135–23149, Curran Associates, Inc. Deistler, M., Boelts, J., Steinbach, P., Moss, G., Moreau, T., Gloeckler, M., Rodriguez, P. L. C., Linhart, J., Lappalainen, J. K., Miller, B. K., Goncalves, P. J., Lueckmann, J.-M., Schr ̈ oder, C., & Macke, J. H., 2025. Simulation-based inference: A practical guide, arXiv. Delaunoy, A., Hermans, J., Rozet, F., Wehenkel, A., & Louppe, G., 2022. Towards reliable simulation-based inference with bal- anced neural ratio estimation, Advances in Neural Information Processing Systems, 35, 20025–20037. Delouis, B., Giardini, D., Lundgren, P., & Salichon, J., 2002. Joint inversion of insar, gps, teleseismic, and strong-motion data for the spatial and temporal distribution of earthquake slip: Applica- tion to the 1999 izmit mainshock, Bulletin of the Seismological Society of America, 92(1), 278–299. Dreger, D. S. & Helmberger, D. V., 1990. Broadband modeling of local earthquakes, Bulletin of the Seismological Society of Amer- ica, 80(5), 1162–1179. Dreger, D. S., Tkalcic, H., & Johnston, M., 2000. Dilational pro- cesses accompanying earthquakes in the long valley caldera, Sci- ence, 288(5463), 122–125. Duputel, Z., Rivera, L., Fukahata, Y., & Kanamori, H., 2012. Un- certainty estimations for seismic source inversions, Geophysical Journal International, 190(2), 1243–1256. Duputel, Z., Agram, P. S., Simons, M., Minson, S. E., & Beck, J. L., 2014. Accounting for prediction uncertainty when inferring subsurface fault slip, Geophysical journal international, 197(1), 464–482. Duputel, Z., Jiang, J., Jolivet, R., Simons, M., Rivera, L., Am- puero, J.-P., Riel, B., Owen, S. E., Moore, A. W., Samsonov, S. V., et al., 2015. The iquique earthquake sequence of april 2014: Bayesian modeling accounting for prediction uncertainty, Geophysical Research Letters, 42(19), 7949–7957. Durkan, C., Bekasov, A., Murray, I., & Papamakarios, G., 2019. Neural spline flows, in Advances in Neural Information Process- ing Systems, vol. 32, Curran Associates, Inc. Durkan, C., Bekasov, A., Murray, I., & Papamakarios, G., 2020a. nflows: normalizing flows in PyTorch, doi: 10.5281/zen- odo.4296287. Durkan, C., Murray, I., & Papamakarios, G., 2020b. On con- trastive learning for likelihood-free inference, in Proceedings of the 37th International Conference on Machine Learning, vol. 119 of Proceedings of Machine Learning Research, p. 2771– 2781, PMLR. Dziewonski, A. M. & Woodhouse, J. H., 1983. An experiment in systematic study of global seismicity: Centroid-moment tensor solutions for 201 moderate and large earthquakes of 1981, Jour- nal of Geophysical Research: Solid Earth, 88(B4), 3247–3271. Engdahl, E. R., van der Hilst, R., & Buland, R., 1998. Global teleseismic earthquake relocation with improved travel times and procedures for depth determination, Bulletin of the Seismologi- cal Society of America, 88(3), 722–743. Eulenfeld, T., Dahm, T., Heimann, S., & Wegler, U., 2022. Fast and robust earthquake source spectra and moment magnitudes from envelope inversion, Bulletin of the Seismological Society of America, 112(2), 878–893. Falcon, W., 2019. Pytorch lightning. Ferreira, A., Weston, J., & Funning, G., 2011.Global com- pilation of interferometric synthetic aperture radar earthquake source models: 2. effects of 3-d earth structure, Journal of Geo- physical Research: Solid Earth, 116(B8). Ferreira, A. M. & Woodhouse, J. H., 2006. Long-period seismic source inversions using global tomographic models, Geophysi- cal Journal International, 166(3), 1178–1192. Fichtner, A., Kennett, B. L., Igel, H., & Bunge, H.-P., 2008. Theoretical background for continental-and global-scale full- waveform inversion in the time–frequency domain, Geophysical Journal International, 175(2), 665–685. Ford, S. R., Walter, W. R., & Dreger, D. S., 2012. Event dis- crimination using regional moment tensors with teleseismic-p constraints, Bulletin of the Seismological Society of America, 102(2), 867–872. Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J., 2013. emcee: the mcmc hammer, Publications of the Astronom- ical Society of the Pacific, 125(925), 306. Gatti, M., Jeffrey, N., Whiteway, L., Williamson, J., Jain, B., Ajani, V., Anbajagane, D., Giannini, G., Zhou, C., Porredon, A., et al., 2024. Dark energy survey year 3 results: Simulation-based cosmological inference with wavelet harmonics, scattering trans- forms, and moments of weak lensing mass maps. validation on simulations, Physical Review D, 109(6), 063534. Gerardi, F., Cuceu, A., Joachimi, B., Nadathur, S., & Font-Ribera, A., 2024. Optimal data compression for lyman-α forest cosmol- ogy, Monthly Notices of the Royal Astronomical Society, 528(2), 2667–2678. Gloeckler, M., Deistler, M., Weilbach, C. D., Wood, F., & Macke, J. H., 2024. All-in-one simulation-based inference, in Inter- national Conference on Machine Learning, p. 15735–15766, PMLR. Gonc ̧alves, P. J., Lueckmann, J.-M., Deistler, M., Nonnenmacher, M., ̈ Ocal, K., Bassetto, G., Chintaluri, C., Podlaski, W. F., Had- dad, S. A., Vogels, T. P., et al., 2020. Training deep neural density estimators to identify mechanistic models of neural dynamics, elife, 9, e56261. Graikos, A., Malkin, N., Jojic, N., & Samaras, D., 2022. Diffusion models as plug-and-play priors, Advances in Neural Information Processing Systems, 35, 14715–14728. Greenberg, D. S., Nonnenmacher, M., & Macke, J. H., 2019. Au- tomatic posterior transformation for likelihood-free inference. Hallo, M. & Gallovi ˇ c, F., 2016. Fast and cheap approximation of green function uncertainty for waveform-based earthquake source inversions, Geophysical Journal International, 207(2), 1012–1029. Hamidbeygi, M., Vasyura-Bathke, H., Dettmer, J., Eaton, D. W., & Dosso, S. E., 2023. Bayesian estimation of non-linear centroid Theory errors in moment tensor inversion with SBI21 moment tensors using multiple seismic data sets, Geophysical Journal International, 235(3), 2948–2961. Hardebeck, J. L. & Shearer, P. M., 2002. A new method for de- termining first-motion focal mechanisms, Bulletin of the Seismo- logical Society of America, 92(6), 2264–2276. Hejrani, B., Tkal ˇ ci ́ c, H., & Fichtner, A., 2017. Centroid moment tensor catalogue using a 3-d continental scale earth model: Ap- plication to earthquakes in papua new guinea and the solomon islands, Journal of Geophysical Research: Solid Earth, 122(7), 5517–5543. Herak, D., Herak, M., & Tomljenovi ́ c, B., 2009. Seismicity and earthquake focal mechanisms in north-western croatia, Tectono- physics, 465(1-4), 212–220. Hermans, J., Begy, V., & Louppe, G., 2020.Likelihood-free MCMC with amortized approximate ratio estimators, in Pro- ceedings of the 37th International Conference on Machine Learning, vol. 119 of Proceedings of Machine Learning Re- search, p. 4239–4248, PMLR. Hermans, J., Delaunoy, A., Rozet, F., Wehenkel, A., & Louppe, G., 2021. Averting a crisis in simulation-based inference, stat, 1050, 14. Herrmann, R. B., 2013. Computer programs in seismology: An evolving tool for instruction and research, Seismological Re- search Letters, 84(6), 1081–1088. Hingee, M., Tkal ˇ ci ́ c, H., Fichtner, A., & Sambridge, M., 2011. Seismic moment tensor inversion using a 3-d structural model: applications for the australian region, Geophysical Journal In- ternational, 184(2), 949–964. Hj ̈ orleifsd ́ ottir, V. & Ekstr ̈ om, G., 2010.Effects of three- dimensional earth structure on cmt earthquake parameters, Physics of the Earth and Planetary Interiors, 179(3-4), 178–190. Ho, J., Kalchbrenner, N., Weissenborn, D., & Salimans, T., 2020. Axial attention in multidimensional transformers. Holzschuh, B. & Thuerey, N., 2024.Flow matching for posterior inference with simulator feedback, arXiv preprint arXiv:2410.22573. Hourcade, C., Juhel, K., & Bletery, Q., 2025. Pegsgraph: A graph neural network for fast earthquake characterization based on prompt elastogravity signals, Journal of Geophysical Research: Machine Learning and Computation, 2(1), e2024JH000360. Hu, J., Pha . m, T.-S., & Tkal ˇ ci ́ c, H., 2023. Seismic moment tensor inversion with theory errors from 2-d earth structure: implica- tions for the 2009–2017 dprk nuclear blasts, Geophysical Jour- nal International, 235(3), 2035–2054. Hu, J., Tkalc ˇ ic ́ , H., Pha . m, T.-S., Herak, M., Dasovi ́ c, I., & Br ˇ ci ́ c, M. M., 2025. Bayesian reassessment of seismic moment tensors and their uncertainties in the adriatic sea region, Seismica, 4(2). Huang, X. & Zhang, Y., 2026. Msep-tformer: a multitask source estimation parameter transformer network for earthquake moni- toring, Geophysical Journal International, 244(1), ggaf453. Jeffrey, N., Alsing, J., & Lanusse, F., 2021. Likelihood-free infer- ence with neural compression of des sv weak lensing map statis- tics, Monthly Notices of the Royal Astronomical Society, 501(1), 954–969. Jeffrey, N., Whiteway, L., Gatti, M., Williamson, J., Alsing, J., Porredon, A., Prat, J., Doux, C., Jain, B., Chang, C., et al., 2025. Dark energy survey year 3 results: likelihood-free, simulation- based w cdm inference with neural compression of weak-lensing map statistics, Monthly Notices of the Royal Astronomical Soci- ety, 536(2), 1303–1322. Kawakatsu, H., 1996. Observability of the isotropic component of a moment tensor, Geophysical Journal International, 126(2), 525–544. Kolb, J. & Leki ́ c, V., 2014. Receiver function deconvolution us- ing transdimensional hierarchical bayesian inference, Geophysi- cal Journal International, 197(3), 1719–1735. K ˇ r ́ ı ˇ zov ́ a, D., Zahradn ́ ık, J., & Kiratzi, A., 2013. Resolvability of isotropic component in regional seismic moment tensor inver- sion, Bulletin of the Seismological Society of America, 103(4), 2460–2473. Kubo, H., Naoi, M., & Kano, M., 2024. Recent advances in earth- quake seismology using machine learning, Earth, Planets and Space, 76(1), 36. Lanzieri, D., Zeghal, J., Makinen, T. L., Boucaud, A., Starck, J.- L., & Lanusse, F., 2025. Optimal neural summarization for full- field weak lensing cosmological implicit inference, Astronomy & Astrophysics, 697, A162. Ledoit, O. & Wolf, M., 2004. A well-conditioned estimator for large-dimensional covariance matrices, Journal of multivariate analysis, 88(2), 365–411. Lehman, K., Krippendorf, S., Weller, J., & Dolag, K., 2025. Learning optimal summary statistics of galaxy catalogs with sbi, Journal of Cosmology and Astroparticle Physics, 2025(12), 032. Lemos, P., Coogan, A., Hezaveh, Y., & Perreault-Levasseur, L., 2023. Sampling-based accuracy testing of posterior estimators for general inference, arXiv preprint arXiv:2302.03026. Lindner, M., Rietbrock, A., Bie, L., Goes, S., Collier, J., Rychert, C., Harmon, N., Hicks, S. P., Henstock, T., & Group, V. W., 2023. Bayesian regional moment tensor from ocean bottom seismo- grams recorded in the lesser antilles: Implications for regional stress field, Geophysical Journal International, 233(2), 1036– 1054. Lueckmann, J.-M., Bassetto, G., Karaletsos, T., & Macke, J. H., 2019. Likelihood-free inference with emulator networks, in Sym- posium on Advances in Approximate Bayesian Inference, p. 32– 53, PMLR. Lueckmann, J.-M., Boelts, J., Greenberg, D., Goncalves, P., & Macke, J., 2021. Benchmarking simulation-based inference, in International conference on artificial intelligence and statistics, p. 343–351, PMLR. Makinen, T. L., Sui, C., Wandelt, B. D., Porqueres, N., & Heavens, A., 2024. Hybrid summary statistics, arXiv preprint arXiv:2410.07548. Marku ˇ si ́ c, S., Stanko, D., Korbar, T., Beli ́ c, N., Penava, D., & Ko- rdi ́ c, B., 2020. The zagreb (croatia) m5. 5 earthquake on 22 march 2020, Geosciences, 10(7), 252. McBrearty, I. W. & Beroza, G. C., 2022. Earthquake location and magnitude estimation with graph neural networks, in 2022 IEEE international conference on image processing (ICIP), p. 3858– 3862, IEEE. McBrearty, I. W. & Beroza, G. C., 2023. Earthquake phase asso- ciation with graph neural networks, Bulletin of the Seismological Society of America, 113(2), 524–547. Minson, S., Simons, M., & Beck, J., 2013. Bayesian inversion for finite fault earthquake source models i—theory and algorithm, Geophysical Journal International, 194(3), 1701–1726. Minson, S. E. & Dreger, D. S., 2008.Stable inversions for complete moment tensors, Geophysical Journal International, 174(2), 585–592. Minson, S. E., Simons, M., Beck, J., Ortega, F., Jiang, J., Owen, S., Moore, A., Inbal, A., & Sladen, A., 2014. Bayesian in- version for finite fault earthquake source models–i: the 2011 great tohoku-oki, japan earthquake, Geophysical Journal Inter- national, 198(2), 922–940. 22A. A. Saoulis et al. Mousavi, S. M. & Beroza, G. C., 2022. Deep-learning seismology, Science, 377(6607), eabm4470. Mousavi, S. M., Ellsworth, W. L., Zhu, W., Chuang, L. Y., & Beroza, G. C., 2020. Earthquake transformer—an attentive deep- learning model for simultaneous earthquake detection and phase picking, Nature communications, 11(1), 3952. M ̈ unchmeyer, J., Bindi, D., Leser, U., & Tilmann, F., 2021. The transformer earthquake alerting model: A new versatile approach to earthquake early warning, Geophysical Journal International, 225(1), 646–656. Nakahara, H., 2008. Seismogram envelope inversion for high- frequency seismic energy radiation from moderate-to-large earthquakes, Advances in Geophysics, 50, 401–426. Novak, M. S., Uros, M., Atalic, J., Herak, M., Demsic, M., Ban- icek, M., Lazarevic, D., Bijelic, N., Crnogorac, M., & Todoric, M., 2020. Zagreb earthquake of 22 march 2020-preliminary re- port on seismologic aspects and damage to buildings. Papamakarios, G. & Murray, I., 2016. Fast ε-free inference of simulation models with bayesian conditional density estimation, in Advances in Neural Information Processing Systems, vol. 29, Curran Associates, Inc. Papamakarios, G., Pavlakou, T., & Murray, I., 2017. Masked au- toregressive flow for density estimation, in Advances in Neural Information Processing Systems, vol. 30, Curran Associates, Inc. Papamakarios, G., Sterratt, D., & Murray, I., 2019. Sequential neural likelihood: Fast likelihood-free inference with autoregres- sive flows, in The 22nd international conference on artificial in- telligence and statistics, p. 837–848, PMLR. Papamakarios, G., Nalisnick, E., Rezende, D. J., Mohamed, S., & Lakshminarayanan, B., 2021. Normalizing flows for probabilis- tic modeling and inference, J. Mach. Learn. Res., 22(1). Park, M., Gatti, M., & Jain, B., 2025. Dimensionality reduction techniques for statistical inference in cosmology, Physical Re- view D, 111(6), 063523. Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., et al., 2019. Pytorch: An imperative style, high-performance deep learning library, Advances in neural information processing systems, 32. Petersen, G. M., Cesca, S., Heimann, S., Niemz, P., Dahm, T., K ̈ uhn, D., Kummerow, J., & Plenefisch, T., 2021. Regional cen- troid moment tensor inversion of small to moderate earthquakes in the alps using the dense alparray seismic network: challenges and seismotectonic insights, Solid earth, 12(6), 1233–1257. Pha . m, T.-S. & Tkal ˇ ci ́ c, H., 2021. Toward improving point-source moment-tensor inference by incorporating 1d earth model’s uncertainty: Implications for the long valley caldera earth- quakes, Journal of Geophysical Research: Solid Earth, 126(11), e2021JB022477. Pha . m, T.-S., Tkal ˇ ci ́ c, H., Hu, J., & Kim, S., 2024. Towards a new standard for seismic moment tensor inversion containing 3- d earth structure uncertainty, Geophysical Journal International, 238(3), 1840–1853. Pha . m, T.-S., Tkal ˇ ci ́ c, H., Hu, J., & Wei, Z., 2025. On the global centroid moment tensor achievements and the next generation earthquake catalogs, Physics of the Earth and Planetary Interi- ors, p. 107490. Pierre, S., R ́ egaldo-Saint Blancard, B., Hahn, C., & Eickenberg, M., 2026.Mitigating model misspecification in simulation- based inference for galaxy clustering, Physical Review D, 113(4), 043536. Poppeliers, C. & Preston, L., 2021. The effects of earth model uncertainty on the inversion of seismic data for seismic source functions, Geophysical Journal International, 224(1), 100–120. Poppeliers, C. & Preston, L., 2022. An efficient method to prop- agate model uncertainty when inverting seismic data for time domain seismic moment tensors, Geophysical Journal Interna- tional, 231(2), 1221–1232. Prelogovi ́ c, D. & Mesinger, A., 2024. How informative are sum- maries of the cosmic 21 cm signal?, Astronomy & Astrophysics, 688, A199. R ̈ osler, B., Stein, S., & Spencer, B. D., 2021. Uncertainties in seis- mic moment tensors inferred from differences between global catalogs, Seismological Research Letters, 92(6), 3698–3711. R ̈ osler, B., Spencer, B. D., & Stein, S., 2024a. Which global moment tensor catalog provides the most precise non-double- couple components?, Seismological Research Letters, 95(4), 2444–2451. R ̈ osler, B., Stein, S., Ringler, A., & Vack ́ a ˇ r, J., 2024b. Apparent non-double-couple components as artifacts of moment tensor in- version, Seismica, 3(1). Saoulis, A., Piras, D., Spurio Mancini, A., Joachimi, B., & Fer- reira, A., 2025. Full-waveform earthquake source inversion us- ing simulation-based inference, Geophysical Journal Interna- tional, 241(3), 1740–1761. Sawade, L., Beller, S., Lei, W., & Tromp, J., 2022. Global centroid moment tensor solutions in a heterogeneous earth: The cmt3d catalogue, Geophysical Journal International, 231(3), 1727– 1738. Schmitt, M., B ̈ urkner, P.-C., K ̈ othe, U., & Radev, S. T., 2023. De- tecting model misspecification in amortized bayesian inference with neural networks, in Dagm german conference on pattern recognition, p. 541–557, Springer. Si, X., Wu, X., Li, Z., Wang, S., & Zhu, J., 2024a. An all-in- one seismic phase picking, location, and association network for multi-task multi-station earthquake monitoring, Communica- tions Earth & Environment, 5(1), 22. Si, X., Wu, X., Sheng, H., Zhu, J., & Li, Z., 2024b. Seisclip: A seismology foundation model pre-trained by multimodal data for multipurpose seismic feature extraction, IEEE Transactions on Geoscience and Remote Sensing, 62, 1–13. Simut ̇ e, S., Boehm, C., Krischer, L., Gokhberg, A., Vall ́ e, M., & Fichtner, A., 2023. Bayesian seismic source inversion with a 3-d earth model of the japanese islands, Journal of Geophysical Research: Solid Earth, 128(1), e2022JB024231. Skoumal, R. J., Hardebeck, J. L., & Shearer, P. M., 2024. Skhash: A python package for computing earthquake focal mechanisms, Seismological Research Letters, 95(4), 2519–2526. Sokos, E. & Zahradn ́ ık, J., 2013. Evaluating centroid-moment- tensor uncertainty in the new version of isola software, Seismo- logical Research Letters, 84(4), 656–665. Song, J.-H., Kim, S., Rhie, J., & Park, D., 2022. Moment tensor solutions for earthquakes in the southern korean peninsula us- ing three-dimensional seismic waveform simulations, Frontiers in Earth Science, 10, 945022. Song, X., Men-Andrin, M., Ellsworth, W. L., & Beroza, G. C., 2025. Foconet: transformer-based focal-mechanism determina- tion, Authorea Preprints. Stachnik, J., Sheehan, A. F., Zietlow, D., Yang, Z., Collins, J., & Ferris, A., 2012. Determination of new zealand ocean bottom seismometer orientation via rayleigh-wave polarization, Seismo- logical Research Letters, 83(4), 704–713. St ̈ ahler, S. C. & Sigloch, K., 2014. Fully probabilistic seismic source inversion–part 1: Efficient parameterisation, Solid Earth, 5(2), 1055–1069. Theory errors in moment tensor inversion with SBI23 St ̈ ahler, S. C. & Sigloch, K., 2016. Fully probabilistic seismic source inversion–part 2: Modelling errors and station covari- ances, Solid Earth, 7(6), 1521–1536. Stip ˇ cevi ́ c, J., Tkal ˇ ci ́ c, H., Herak, M., Marku ˇ si ́ c, S., & Herak, D., 2011. Crustal and uppermost mantle structure beneath the exter- nal dinarides, croatia, determined from teleseismic receiver func- tions, Geophysical journal international, 185(3), 1103–1119. Stip ˇ cevi ́ c, J., Herak, M., Molinari, I., Dasovi ́ c, I., Tkal ˇ ci ́ c, H., & Gosar, A., 2020.Crustal thickness beneath the dinarides and surrounding areas from receiver functions, Tectonics, 39(3), e2019TC005872. Stockman, S., Lawson, D. J., & Werner, M. J., 2024. Sb-etas: using simulation based inference for scalable, likelihood-free in- ference for the etas model of earthquake occurrences, Statistics and Computing, 34(5), 174. Sun, H., Ross, Z. E., Zhu, W., & Azizzadenesheli, K., 2023. Phase neural operator for multi-station picking of seismic ar- rivals, Geophysical Research Letters, 50(24), e2023GL106434. Takemura, S., Okuwaki, R., Kubota, T., Shiomi, K., Kimura, T., & Noda, A., 2020. Centroid moment tensor inversions of offshore earthquakes using a three-dimensional velocity structure model: slip distributions on the plate boundary along the nankai trough, Geophysical Journal International, 222(2), 1109–1125. Talts, S., Betancourt, M., Simpson, D., Vehtari, A., & Gelman, A., 2018. Validating bayesian inference algorithms with simulation- based calibration, arXiv preprint arXiv:1804.06788. Tape, W. & Tape, C., 2015. A uniform parametrization of moment tensors, Geophysical Journal International, 202(3), 2074–2081. Tarantola, A., 2005. Inverse problem theory and methods for model parameter estimation, SIAM. Tarantola, A., Valette, B., et al., 1982. Inverse problems= quest for information, Journal of geophysics, 50(1), 159–170. Tejero-Cantero, A., Boelts, J., Deistler, M., Lueckmann, J.-M., Durkan, C., Gonc ̧alves, P., Greenberg, D., & Macke, J., 2020. sbi: A toolkit for simulation-based inference, Journal of Open Source Software, 5(52), 2505. Thuerey, N., Holl, P., Mueller, M., Schnell, P., Trost, F., & Um, K., 2021.Physics-based deep learning, arXiv preprint arXiv:2109.05237. Thurin, J., Modrak, R., Tape, C., McPherson, A., Rodr ́ ıguez- Cardozo, F., Kintner, J., Ding, L., Liu, Q., & Braunmiller, J., 2025. Mtuq: a framework for estimating moment tensors, point forces, and their uncertainties, Geophysical Journal Interna- tional, 241(2), 1373–1390. Trabattoni, A., Barruol, G., Dreo, R., Boudraa, A., & Fontaine, F., 2020. Orienting and locating ocean-bottom seismometers from ship noise analysis, Geophysical Journal International, 220(3), 1774–1790. Vack ́ a ˇ r, J., Burj ́ anek, J., Gallovi ˇ c, F., Zahradn ́ ık, J., & Clinton, J., 2017.Bayesian isola: new tool for automated centroid moment tensor inversion, Geophysical Journal International, 210(2), 693–705. Valentine, A. P. & Trampert, J., 2012. Assessing the uncertainties on seismic source parameters: Towards realistic error estimates for centroid-moment-tensor determinations, Physics of the Earth and Planetary Interiors, 210, 36–49. Vall ́ e, M., Charl ́ ety, J., Ferreira, A. M., Delouis, B., & Vergoz, J., 2011. Scardec: a new technique for the rapid determination of seismic moment magnitude, focal mechanism and source time functions for large earthquakes using body-wave deconvolution, Geophysical Journal International, 184(1), 338–358. Van Amersfoort, J., Smith, L., Teh, Y. W., & Gal, Y., 2020. Un- certainty estimation using a single deep deterministic neural network, in International conference on machine learning, p. 9690–9700, PMLR. Vasyura-Bathke, H., Dettmer, J., Dutta, R., Mai, P. M., & Jonsson, S., 2021. Accounting for theory errors with empirical bayesian noise models in nonlinear centroid moment tensor estimation, Geophysical Journal International, 225(2), 1412–1431. Vavry ˇ cuk, V. & K ̈ uhn, D., 2012. Moment tensor inversion of waveforms: a two-step time-frequency approach, Geophysical Journal International, 190(3), 1761–1776. von Wietersheim-Kramsta, M., Lin, K., Tessore, N., Joachimi, B., Loureiro, A., Reischke, R., & Wright, A. H., 2025. Kids-sbi: Simulation-based inference analysis of kids-1000 cosmic shear, Astronomy & Astrophysics, 694, A223. W ́ eber, Z., 2006. Probabilistic local waveform inversion for mo- ment tensor and hypocentral location, Geophysical Journal In- ternational, 165(2), 607–621. Weston, J., Ferreira, A., & Funning, G., 2011.Global com- pilation of interferometric synthetic aperture radar earthquake source models: 1. comparisons with seismic catalogs, Journal of Geophysical Research: Solid Earth, 116(B8). Yagi, Y. & Fukahata, Y., 2011. Introduction of uncertainty of green’s function into waveform inversion for seismic source pro- cesses, Geophysical Journal International, 186(2), 711–720. Zahradn ́ ık, J. & Sokos, E., 2025. Isola2024: Assessing and un- derstanding uncertainties of full moment tensors, Seismological Research Letters, 96(4), 2647–2659. Zammit-Mangion, A., Sainsbury-Dale, M., & Huser, R., 2025. Neural methods for amortized inference, Annual Review of Statistics and Its Application, 12(1), 311–335. Zhu, W. & Beroza, G. C., 2019. Phasenet: a deep-neural-network- based seismic arrival-time picking method, Geophysical Journal International, 216(1), 261–273. Zhu, W., Mousavi, S. M., & Beroza, G. C., 2019. Seismic signal denoising and decomposition using deep neural networks, IEEE Transactions on Geoscience and Remote Sensing, 57(11), 9476– 9488. Supporting Information The Supporting Information for this study contains additional figures and explanations that provide further details on the method- ology, results, and discussions presented in the main manuscript. Below is a brief description of each figure in the Supporting Infor- mation: Section S1: A discussion of higher order terms in the optimal score compression scheme, which we neglect in our approach. Figure S1: Example of posterior inference and posterior predictive checks for an example shorter period 6− 50 s synthetic inversion. Figure S2: Diagram of the artificial balanced receiver configura- tion for one of the synthetic experiments. Figure S3: Calibration performance for the score compression SBI approaches for the various synthetic experiments addressed in Sec- tion 6. Figure S4: Same as in Fig. S3, but for the ML-based compression SBI approach. Also probes the effect of ML-based compression calibration under model misspecification. Figure S5: An analysis of the full distribution of per-parameter posterior accuracies for the 300 calibration test inversions with κ = 3% for each of the three approaches. Figure S6: Same as Fig. S5 but for κ = 5%. 24A. A. Saoulis et al. Figure S7: Same as Fig. S5 but for κ = 5% and shorter period data between 6− 50 s. Table S1: The mean posterior standard deviation and bias per- parameter for each approach, as in Table 1, for κ = 3%. Table S2: The mean posterior standard deviation and bias per- parameter for each approach, as in Table 1, for κ = 5% and shorter period data between 6− 50 s. Figure S8: The 2-D composite model for the Croatia event adapted from Hu et al. (2025), as described in the main text. Table S3: Results from a quantitative analysis of the posterior pre- dictive checks for each of the three appraoches for the LV2 event. Table S4: Same as Table3 for the Croatia event. The full Supporting Information file is available online, linked with this article. Supporting Information 20 March 2026 S1 OPTIMAL SCORE COMPRESSION In the main text, we quoted the following result for optimal score compression: t = m ∗ + F −1 ∗ G T m ∗ C −1 (D−μ ∗ ),(S1) where the compressed statistic t can be computed from a given observation D using the local sensitivity kernels G m ∗ and covariance estimate C m ∗ (note that F ∗ = G T m ∗ C −1 m ∗ G m ∗ ). However, if the covariance matrix C depends on model parameters m, as it does for theory errors, the full result from Alsing et al. (2017) includes extra terms that rely on the gradients of the covariance ∇ m C: t = m ∗ + F −1 ∗ " G ∗ T C −1 ∗ (D−μ ∗ ) + 1 2 (D−μ ∗ ) T C −1 ∗ (∇C ∗ ) C −1 ∗ (D−μ ∗ )− 1 2 tr C −1 ∗ ∇C ∗ # . (S2) The first two terms correspond to Equation (S1), while the latter two terms depend on non-zero gradients of the covariance∇C ∗ . This correction is quadratic in the data and is therefore a higher- order correction term. In preparation for this work, we implemented and explored the results using the complete result for an optimal statistic t. However, we found that this gave very poor perfor- mance, with the extra two terms completely dominating the MLE estimate t. This was driven by very large gradient terms∇ m C ∗ ; i.e. the theory covariance is very sensitive to perturbations to the moment tensor parameters. We therefore hypothesise that the strong model parameter dependency of the theory covariance C t (m) causes the simple, linearised relationship assumed in the Taylor expansion to break down. 2 A. A. Saoulis et al. For this reason, throughout this work we adopt the simplified optimal score statistic given by Equa- tion (S1), and neglect the covariance-gradient terms in Equation (S2). In practice, we find that this approximation yields significantly more stable behaviour while retaining the desired sensitivity to variations in the model parameters. S2 SYNTHETIC EXPERIMENTS S2.1 Short period inversion example In Figure S1 we present an example of the three inversion approaches for a shorter period filtering band between 6− 50 s, which we analysed in the main text. This illustrates the issues with the likelihood-based approaches to inversion, which can occasionally lead to biased moment tensor solutions and underestimation of uncertainty. The degree of bias shown in Figure S1 is not uncommon for the Gaussian likelihood approach; the z-score reaches z > 3 for several parameters. Shown systematically for the 300 calibration test inversions in Figure S7, this occurs roughly 10− 30% of the time, depending on the parameter. S2.2 Balanced network configuration In the main text, we explored how each inversion approach performed with a more azimuthally balanced network configuration. The network is shown in Figure S2. The results in the main text demonstrated that even with very good station coverage, Gaussian likelihood assumptions still yield uncalibrated, and therefore misleading, posterior uncertainties. S2.3 Further calibration tests of the SBI approaches In Figs. S3 and S4 we show the calibration performance of the two SBI formalisms over various experiments. Figure S3 repeats the varying experimental configurations from the main text using the score compression SBI approach. We find that across all experiments, the SBI approach is much better calibrated than a Gaussian likelihood-based inversion. We find minor bias for κ = 5% as well as for the short period inversions, which we commented on in the main text. Nonetheless, the calibration statistics are acceptable and far improved over the likelihood-based alternative. Supporting Information3 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 estima- tion (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 pa- rameters γ,δ,M W 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). 4 A. A. Saoulis et al. 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 S4(a) shows the performance of the ML-based compression SBI approach for the most challenging examples. For the varying Earth structure uncertainty levels κ ∈ [3, 5]%, we find a minor indication of underconfident (or conservative) posteriors. We discuss this in the main text. However, for the short-period inversions we find an indication of bias. We found that training stability for the short period inversions was relatively poor, and the deep learning model began to overfit quite quickly. We attribute this to the more challenging modelling required to process the short period data, and flag this issue as an avenue for future work. In particular, improvements to model architecture or larger, more varied training datasets may be required to mitigate these issues. Figure S4(b) shows a simple test of the effects of model misspecification on the ML-based compression approach. Specifically, we train the model on a given level of 1-D Earth structure uncertainty κ trained , and apply it to artificial observations produced by a different level of uncer- tainty κ applied . The results are relatively intuitive: models trained on a higher (lower) degree of Earth model uncertainty over(under)-estimate the uncertainty when applied to observations pro- duced with lower (higher) Earth structure uncertainty. The encouraging result of these tests is that applying the model to slightly out-of-distribution data does not lead to model collapse or serious Supporting Information5 Figure S3. The calibration performance of the SBI with optimal score compression under varied exper- imental 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 net- work 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 exper- imental 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. 6 A. A. Saoulis et al. bias. However, as discussed in the main text, real observations will have very different imprints of model misspecification than those explored here, so more comprehensive testing is required in future work. S2.4 Quantitative analysis of moment tensor solutions from each approach Here, we take a more in-depth look at the results of the 300 artificial inversions performed for each experimental configuration. We compute the z-score, which measures how far the true parameter lies from the posterior mean when scaled by the posterior uncertainty. For a given parameter m, the z-score is defined as z = μ m − m true σ m ,(S3) where μ and σ denote the posterior mean and standard deviation, respectively. If the posterior is well calibrated, the z-scores computed across repeated synthetic experiments should be distributed according to a standard normal distribution,N (0, 1). Deviations from this behaviour indicate ei- ther biased posterior estimates or misestimation of posterior uncertainties. While this assumes that the posterior is Gaussian (which is not the case in our examples), the z-score still provides a in- tuitive and interpretable diagnostic for the degree of discrepancy between the inferred moment tensor solutions and the true source. We show the distributions of bias and z-score for each parameter independently for different experimental configurations in Figs. S5, S6, and S7. Figure S5 shows that, for κ = 3% Earth structure uncertainty, the degree of bias in the inferred solutions is similar across all approaches. The ML-based compression approach performs only slightly better in producing precise posterior solutions. However, the Gaussian likelihood approach produces poor uncertainty estimates that lead to over-confidence. The two SBI-based approaches produce better calibrated uncertainties, with closer alignment between the posterior solution z- scores and the expected statistical distribution in Figure S5(c). We see evidence of the mild bias of the score compression SBI approach, and the under-confidence (conservativeness) of the ML- based compression approach, which was discussed in the main text. Supporting Information7 At κ = 5%, Figure S6 shows that the degree of bias in both the Gaussian likelihood approach and the score compression approach increases significantly. While the score compression approach largely accounts for this with larger uncertainties, the Gaussian likelihood approach is very over- confident. On the other hand, the ML-based approach produces posterior solutions that are much more precise (i.e. they have a much narrower distribution of bias). These differences are apparent across all the focal mechanism parameters. We find similar results for the shorter period inversions, with filtering between 6− 50 s, in Figure S7. Again, the ML-based approach is much more precise than the alternatives in recovering the true source parameters. However, as highlighted in the main text, we also find some minor indications of bias and miscalibration: for instance, the CLVD component γ is systematically under-estimated. We also report the mean per-parameter biases and posterior standard deviations of the κ = 3% inversions in Table S1 and for the 6− 50 s filtering band inversions in Table S2 (κ = 5% is given in the main text). For κ = 3%, the Gaussian likelihood and score compression SBI biases are small but significant. However, for the shorter period inversions the Gaussian likelihood approach has major biases. These underscore the findings above and in the main text. In particular, M w is systematically underestimated by over one standard deviation of the posterior (σ = 0.37× 10 −2 , bias = 0.50× 10 −2 ). The positive isotropic component is also systematically under-estimated by around 50% of σ. 8 A. A. Saoulis et al. Figure S5. A detailed look at the posterior solutions for the 300 artificial inversions used to perform cover- age 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. S3 REAL DATA INVERSIONS S3.1 Croatia 2-D composite Earth model Figure S8 shows the 2-D composite Earth model used to analyse the Zagreb, Croatia event described in the main text. It has been adapted from Hu et al. (2025) to include 1-D velocity model perturbations. The simplistic treatment of velocity model perturbations used in this work (with Supporting Information9 Figure S6. A detailed look at the posterior solutions for the 300 artificial inversions used to perform cover- age 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. independent random perturbations applied to layer thickness, V p and V s ) can lead to unphysical effects such as occasional slower layers at greater depths. We therefore smoothed Model 4 and reduced the number of layers relative to the model used in Hu et al. (2025); Figure S8 shows the smoothed Model 4 used in this work. 10 A. A. Saoulis et al. Figure S7. A detailed look at the posterior solutions for the 300 artificial inversions used to perform cov- erage 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. These unrealistic perturbations could be addressed in future work, ideally using data-derived tomographic ensembles (Chiang et al. 2025). Supporting Information11 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). S3.2 Quantitative evaluation of the posterior predictive checks For both the real event inversions in the main manuscript, we perform 10,000 posterior pre- dictive simulations. We use these to quantify how closely each of the modelled posteriors can reproduce the observations. We emphasise that this analysis is not intended as a conventional posterior predictive check in the sense of assessing whether the observed data are typical under the posterior predictive distri- bution. Instead, our goal is to evaluate the capability of each inferred posterior to reproduce the MethodCompressionMetric γδ M w (10 −2 )strikediprake SBI Optimal score σ0.420.770.280.640.490.78 Bias-0.06-0.07-0.02-0.09-0.03-0.02 Deep learning σ0.530.650.260.690.560.87 Bias-0.01-0.04-0.01-0.020.05-0.01 Gaussian- σ0.260.420.170.380.300.47 Bias-0.050.01-0.08-0.07-0.11-0.09 Table S1. Average posterior standard deviations (σ) and biases for each inference method over the 300 arti- ficial 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. M w values are reported in units of 10 −2 , all other values in units of degrees. 12 A. A. Saoulis et al. MethodCompressionMetric γδ M w (10 −2 )strikediprake SBI Optimal score σ1.491.771.123.021.953.04 Bias-0.01-0.570.34-0.43-0.140.16 Deep learning σ0.751.010.401.120.991.50 Bias-0.21-0.02-0.05-0.090.10-0.07 Gaussian- σ0.520.680.370.920.580.90 Bias0.03-0.290.50-0.330.11-0.26 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. salient features of the observations when marginalising over a large number of nuisance param- eters. This is useful since our ML-based compression treats station-specific time-shifts as extra nuisance parameters during the inversion Concretely, we evaluate a set of complementary distance or misfit metrics between each reali- sation and the observed data. For a given metric, we then retain the subset of posterior predictive samples that minimise that metric, and report the mean misfit over this “best-matching” subset. This procedure provides a lower-bound estimate on the discrepancy achievable by each posterior under a given notion of similarity, and allows us to compare different inference strategies in terms of their ability to reproduce the observations. We use the following metrics to serve as somewhat independent checks on waveform similarity: (i) Reduced chi-squared (χ 2 red ). Covariance-weighted waveform misfit (see main text for de- tails). (i) Correlation misfit. One minus the Pearson correlation between observed and synthetic waveforms to quantify phase and shape agreement. (i) Envelope misfit. Mean absolute relative difference between the signal envelopes. (iv) Power misfit. Relative difference in total per-station-component waveform power. (v) Mean squared error (MSE). Unweighted elementwise ℓ 2 misfit between observed and synthetic waveforms. Supporting Information13 Methodχ 2 red Correlation misfitEnvelope misfitPower misfitMSE Score compression SBI0.5630.0520.4210.0013.72e-13 ML compression SBI0.5900.0950.4550.0015.04e-13 Gaussian likelihood0.6060.1100.4190.0017.09e-13 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. Across both Tables S3 and S4, we find the SBI approaches produce moment tensor solutions that are competitive with or improve upon the Gaussian likelihood approach in recovering the observation. In Table S3, the Gaussian likelihood approach produces worse fitting “closest” solu- tions than the SBI approaches across all metrics but the envelope misfit. Interestingly, the score compression SBI approach produces better fitting synthetics than the ML-based approach across all metrics, despite its lower effective model dimensionality (since the ML-based compression approach has greater expressivity through its sampling over station-specific time-shifts). This pro- vides some evidence that for the LV2 event, the score compression approach gives marginally more plausible solutions than the ML-based compression approach. In Table S4, the results are more mixed, with each approach performing better or worse de- pending on the metric used. This indicates that all approaches produce roughly comparable solu- tions in terms of reproducing the observation. We emphasise that these are not traditional posterior predictive checks, which we avoid since we used a different (broader) prior for the ML-based compression approach. In particular, models Methodχ 2 red Correlation misfitEnvelope misfitPower misfitMSE Score compression SBI2.6480.0550.4820.0081.53e-11 ML compression SBI2.8080.0830.4580.0051.90e-11 Gaussian likelihood3.5750.1820.4380.0023.97e-11 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. 14 A. A. Saoulis et al. with broad posterior predictive variance may achieve low best-subset misfits despite assigning low probability to observation-like realisations. For this reason, the results in Tables S3-S4 are best interpreted as a diagnostic of model expressivity and solution ability to reproduce the observations. They are sufficient to show that the ML-based compression SBI solutions yields comparable fits to the observations. Supporting Information15 REFERENCES Alsing, J., Wandelt, B., & Feeney, S., 2017. Massive optimal data compression and density estimation for scalable, likelihood-free inference in cosmology, MNRAS, 000, 1–14. Chiang, A., Ford, S. R., Pasyanos, M. E., & Simmons, N. A., 2025. Bayesian inference for the seismic moment tensor using regional waveforms and teleseismic-p polarities with a data-derived distribution of velocity models and source locations, Bulletin of the Seismological Society of America. Hu, J., Tkalc ˇ ic ́ , H., Pha . m, T.-S., Herak, M., Dasovi ́ c, I., & Br ˇ ci ́ c, M. M., 2025. Bayesian reassessment of seismic moment tensors and their uncertainties in the adriatic sea region, Seismica, 4(2).