MTUQ: a framework for estimating moment tensors, point forces, and their uncertainties

SUMMARY

We introduce MTUQ, an open-source Python package for seismic source estimation and uncertainty quantification, emphasizing flexibility and operational scalability. MTUQ provides MPI-parallelized grid search and global optimization capabilities, compatibility with 1-D and 3-D Green’s function database formats, customizable data processing, C-accelerated waveform and first-motion polarity misfit functions, and utilities for plotting seismic waveforms and visualizing misfit and likelihood surfaces. Applicability to a range of full- and constrained-moment tensor, point force, and centroid inversion problems is possible via a documented application programming interface, accompanied by example scripts and integration tests. We demonstrate the software using three different types of seismic events: (1) a 2009 intraslab earthquake near Anchorage, Alaska; (2) an episode of the 2021 Barry Arm landslide in Alaska; and (3) the 2017 Democratic People’s Republic of Korea underground nuclear test. With these events, we illustrate the well-known complementary character of body waves, surface waves, and polarities for constraining source parameters. We also convey the distinct misfit patterns that arise from each individual data type, the importance of uncertainty quantification for detecting multimodal or otherwise poorly constrained solutions, and the software’s flexible, modular design.

1 INTRODUCTION

Accurate estimation of seismic source parameters is important for understanding regional and global tectonic processes, assessing seismic hazards, and monitoring compliance with nuclear test-ban treaties. Two primary mathematical models used to represent seismic point sources are (1) the moment tensor, a second-order, symmetric tensor |$\mathbf {M} \in \mathbb {R}^{3\times 3}$|⁠, suitable for representing small- to moderate-sized seismic events (Aki & Richards 2002); and (2) the point force |$\mathbf {f} \in \mathbb {R}^3$|⁠, appropriate for mass movements such as landslides and volcanic eruptions (Kanamori & Given 1982). Operationally, hundreds of inversions for such sources are carried out each day (e.g., Hayes et al. 2009; Herrmann et al. 2011; Ekström et al. 2012).

The moment tensor |$\mathbf {M}$| is defined as a symmetric second-order tensor that can be expressed as

where |$\mathbf {U}$| is the orthonormal matrix of eigenvectors (⁠|$\mathbf {UU}^T = \mathbf {I}$|⁠), which define the source orientation, and |$\mathbf {\Lambda }$| is the diagonal matrix of eigenvalues (⁠|$\lambda _1 \ge \lambda _2 \ge \lambda _3$|⁠), which describe the source type and magnitude. This provides a natural framework for source parametrization and visualization.

Accurate characterization of seismic sources remains challenging, despite the relatively low number of parameters required to define moment tensors or point forces. Low signal-to-noise ratios, inaccurate Earth models (e.g. Ford et al. 2009a; Duputel et al. 2014; Phạm & Tkal̆cić 2021; Hu et al. 2023) or tradeoffs between related source types (Ford et al. 2010) give rise to multimodal or otherwise challenging misfit spaces. The misfit space for small or sparsely observed events can exhibit multiple minima, causing conventional local-optimization methods to fail. Even for larger events recorded by many stations, the misfit surfaces over source type or orientation may exhibit complicated patterns that error bars, ellipses or other simple methods fail to capture. Uncertainty analyses established by early studies provided a means to detect and quantify such misfit complexity (e.g. Riedesel & Jordan 1989; Valentine & Trampert 2012; Käufl et al. 2014; Stähler & Sigloch 2014). More recent geometric parametrizations (Tape & Tape 2012, 2015, 2016) have made moment tensor uncertainty quantification easier, facilitating visualization of higher dimensional spaces and reducing distortion in the projections.

We present MTUQ (Moment Tensor Uncertainty Quantification), a Python toolbox for seismic source estimation and uncertainty quantification that uses these new parametrizations to implement moment tensor and force parameter searches. MTUQ combines source parameter exploration with a flexible software design, in which users can adapt example scripts for common inversion tasks or implement fully customized workflows, drawing from an API that supports a variety of research avenues and operational applications.

MTUQ has been used in several studies (Thurin et al. 2022; Delbridge et al. 2023; Kintner et al. 2023; Kumar et al. 2023; Thurin & Tape 2023), yet its design and capabilities have not yet been thoroughly discussed. Below, we first recall the formulation of source estimation as a function minimization problem, with optional ‘cut-and-paste’ time-shift corrections to account for Earth model errors (Zhu & Helmberger 1996). We then describe MTUQ’s structure, features, and usage. Finally, we apply it to three test cases: (1) a moment tensor inversion for an intraslab earthquake; (2) a point force inversion for an episode of the Barry Arm landslide in south-central Alaska; and (3) a moment tensor inversion for the 2017 Democratic People’s Republic of Korea nuclear test.

2 POINT SOURCE ESTIMATION

Point source estimation can be framed as a function minimization problem, in which we determine the seismic source |$\mathbf {s}$|—which is typically either a moment tensor |$\mathbf {M}$| or a force |$\mathbf {f}$|—that minimizes the misfit

between observed data |$u_i^{\rm obs}(\mathbf {x}_r,t)$| and synthetic data |$u_i(\mathbf {x}_r,t;\mathbf {s})$|⁠, with N denoting the number of stations, T the length of data window over time t, and |$u_i(\mathbf {x}_r,t)$| the ith component of the displacement at the rth station location |$\mathbf {x}_r$|⁠. The misfit function here is written using the |$L_2$|-norm of the residuals, and other norms can alternatively be used.

The displacement response |$\mathbf {u}$| to a moment tensor |$\mathbf {M}$| or a point force |$\mathbf {f}$| at source location |$\mathbf {x}_s$| and origin time |$t_0$| is given by (Shearer 2009)

or

respectively, where |$S(t)$| is the source–time function and G is the impulse response (i.e. Green’s functions) of the Earth model |$\mathbf {V}$|⁠.

In terms of Green’s functions, the misfit function (eq. 2) can be written concisely as

where |$\mathbf {u}^{\rm syn}$| is obtained from the linear combination |$\mathcal {G} \mathbf {s}$|⁠, and the temporal dependence and source time convolution having been made implicit. This concise formulation shows that the forward problem can be carried out using pre-computed Green’s function databases, as the synthetics become just a linear combination of Green’s functions and source terms.

2.1 Uncertainty in the source characterization problem

Uncertainty in source estimation problems arises from data noise and model errors, while the source–receiver geometry influences the ability to resolve the source by affecting data coverage and sensitivity (e.g. azimuthal gap, stations aligned with the nodal plane of the source). While data noise is embedded in the observed data |$\mathbf {u}^{\rm obs}$|⁠, model error arises from limitations of |$\mathbf {V}$| which will not fully capture the Earth’s physical properties, as well as numerical errors in computing synthetics. In practice, Earth model limitations can be partly alleviated through the so-called ‘cut-and-paste’ (CAP) method (e.g. Zhao & Helmberger 1994; Zhu & Helmberger 1996; Zhu & Ben-Zion 2013) at local and regional scales (typically within 0–|$20^\circ$| distance). CAP works by partitioning waveforms into filtered and windowed body and surface wave phases and evaluating the misfit (5) between observations and time-shifted synthetics. Because at long periods, model error manifests primarily through phase misalignment, the CAP method can provide well-behaved long-period measurements and has become widely adopted in regional seismology. A common approach for determining the necessary time-shifts is to maximize the normalized cross-correlation coefficient:

The maximum of |$\chi (t)$| usually occurs when the observed and synthetic waveforms are correctly aligned, but this may not be the case for nodal components or when cycle skipping occurs. MTUQ allows users to evaluate the misfit function (eq. 5) either in its original form or using the CAP method (i.e. with time-shifted synthetics based on eq. 6). Because these misfit function evaluations are the most computationally expensive step in the whole parameter estimation problem, MTUQ provides highly optimized misfit function implementations, as described in Section 3.

2.2 Parameter space exploration strategies

A simple but powerful way to estimate seismic point sources is to evaluate the misfit at regular intervals over the source parameter space. Such searches may include source type, source orientation, magnitude, depth, centroid, and origin time parameters, resulting in grids of up to 10 dimensions. Evaluating the misfit function over the grid yields the minimum misfit solution along with the accompanying local and global structure of the misfit surface. In practice, we use the moment tensor parametrization proposed by Tape & Tape (2015) to construct grids that are either regularly spaced or randomly drawn from a uniform distribution. Unlike the |$M_{ij}$| parametrization, Tape & Tape (2015)’s approach involves bounded ranges for source type and orientation, making it possible to avoid arbitrary lower and upper limits (the source magnitude is still unbounded). Evaluating the misfit function over individual grid points is an embarrassingly parallel problem, making it easy to implement on parallel systems, although limits posed by dimensionality make high-dimensional searches impractical.

For high-dimensional searches, we can turn to more efficient sampling strategies, which tend to prioritize low-misfit over high-misfit subsets of the parameter space. Among strategies that retain a grid structure, Xu et al. (2021) proposed a double-grid approach with coarse-to-fine grid search strategies. Another approach suggested by Lee et al. (2011) takes advantage of a coarse grid search and uses the best solution as the starting point for a Gauss–Newton iterative optimization method.

Another common strategy is to cast the moment tensor inversion problem in the Bayesian framework, using function-informed methods (e.g. Vackář et al. 2017; Pugh & White 2018; Phạm & Tkal̆cić 2021; Hu et al. 2023) such as classical Markov chain Monte Carlo (MCMC) algorithms, or gradient-informed methods (Fichtner & Simute 2018), such as Hamiltonian Monte Carlo (HMC). Given the challenges associated with convergence in Bayesian inversion, grid searches remain valuable as they provide a comprehensive benchmark for evaluating the accuracy and performance of these probabilistic approaches, particularly when exploring high-dimensional misfit surfaces. Grid searches are also the basis of formal uncertainty quantification which require a uniform sampling of the moment tensor space (Tape & Tape 2016). As such, MTUQ’s primary approach for source parameter estimation is based on exhaustive grid searches, offering a systematic and complete exploration of the parameter space.

Beyond grid searches, MTUQ’s modular design accommodates advanced solvers as optional plugins. For example, the Covariance Matrix Adaptation Evolution Strategy (CMA-ES; Hansen & Ostermeier 2001; Hansen et al. 2015), a gradient-free global optimization method with self-tuning that has performed well in blind tests (Akimoto et al. 2010; Glasmachers et al. 2010), has been successfully applied in prior source studies (Thurin et al. 2022; Thurin & Tape 2023). It allows addressing higher dimensional problems, such as joint inversion for source and centroid parameters or multiple subevents characterization. MTUQ’s modular structure has also allowed us to pair it with a Hamiltonian Monte Carlo sampler (Ding et al. 2021) through SEISGen.

In this study, we expand on previous efforts (Thurin et al. 2022; Thurin & Tape 2023) that applied new optimization methods to moment tensor inversion. Below (Section 4.3), we present a combined workflow that uses CMA-ES, which we found to be especially promising, for initial source parameters estimation, followed by grid searches for quality control and uncertainty quantification.

3 SOFTWARE IMPLEMENTATION

MTUQ provides a flexible Python library for seismic source estimation and uncertainty quantification. Along with conventional |$l_1$| and |$l_2$| misfit evaluation, it implements the CAP algorithm of Zhu & Ben-Zion (2013) and its extensions (e.g. Alvizuri & Tape 2016; Silwal & Tape 2016; Alvizuri et al. 2018). Building on ObsPy and Instaseis data structures (Krischer et al. 2015; van Driel et al. 2015), MTUQ offers data misfit and data processing functions, local and remote Green’s function database clients, grid search and global optimization capabilities, and a range of visualization utilities. Customization is possible by either replacing default functions with user-defined ones or through a flexible frontend/backend architecture using object-oriented programming (e.g. subclassing).

MTUQ users define body-wave and surface-wave windows separately, allowing each set of windows to be weighted and normalized on its own. This approach prevents the inversion from being dominated by the large amplitudes of surface waves, while also allowing more control and flexibility on the misfit definition. Distinct pre-processing steps, such as bandpass filtering and time-shift corrections, can also be applied to each wave type. As shown by Zhu & Helmberger (1996), treating body and surface waves individually leverages their respective sensitivities to different parts of the Earth’s structure and accounts for their strengths within differing frequency bands.

Waveform misfit functions are implemented with Cython acceleration (Behnel et al. 2011) and optional MPI parallelization, making them suitable for either consumer-grade or high-performance computing (HPC) systems and adaptable to both research and operational needs. MTUQ includes several implementations of misfit functions, ranging from readable, pure Python versions to optimized C/Python versions. The pure Python versions are useful for verifying the correctness of the more optimized ones. In addition to evaluating waveform misfit, MTUQ can use first-motion polarities comparing observed and predicted first motions to further constrain the source estimation. Misfit objects define various user-supplied options, such as the norm to use (L1, L2 or hybrid L1-L2), allowable time-shifts, and time-shift groupings (e.g. separating Rayleigh and Love waves by components). This flexible structure offers users fine-grained control over how different waveform types are treated in the inversion.

MTUQ supports the SAC data format (Goldstein et al. 2003) and makes extensive use of ObsPy data structures, with most operations on seismic traces handled by ObsPy Stream() and Trace() methods. In addition to providing seismic data, users can optionally provide a data weight file, which controls the relative contribution of individual stations and components to the overall data misfit. This weight file also allows for tuning body wave arrival times, static corrections (per-station correction to account for path effects larger than the allowable time-shift), and other windowing settings. Vertical and radial P waves, as well as vertical, radial, and transverse surface waves, can be adjusted independently. The effect of geometrical spreading can also be compensated via an amplitude scaling relationship based on distance (Zhu & Helmberger 1996), which helps reduce the misfit contribution from amplitude decay that is otherwise unaccounted for. The weight files are designed with data sharing in mind, so that, when used with MTUQ library functions, data processing and misfit evaluation results can be reproduced.

Interfaces are provided for computing synthetics through elastic wave solvers and Green’s function databases (Table 1), including FK (Zhu & Rivera 2002), AxiSEM (Nissen-Meyer et al. 2014), SPECFEM3D_Cartesian (Komatitsch et al. 2004), SPECFEM3D_GLOBE (Komatitsch & Tromp 2002) and strain Green’s functions generated with SEISGen (Ding et al. 2021). If no local Green’s function databases are available, MTUQ can fetch Green’s tensors from the IRIS syngine web service (Krischer et al. 2017), which provides three-component AxiSEM-generated Green’s functions for periods as low as 2 s for various 1-D Earth models.

Finally, utilities are provided for visualizing MTUQ grid search and global optimization results. Currently utilities are based on Matplotlib and optionally on Generic Mapping Tools (Wessel & Smith 1991) and its Python interface PyGMT (Tian et al. 2024), with support for additional graphics backends possible through a frontend/backend architecture. Available functions include those for plotting misfit, likelihood, variance reduction, and tradeoffs (source-type versus magnitude, source-type versus depth) over point source parameters, as well as for visualizing how related station and trace attributes (such as misfit, time-shifts, normalized cross-correlation coefficient, and amplitude ratios) vary geographically. Using parametrizations from Tape & Tape (2012, 2015), we can plot marginal likelihood distributions over the moment tensor’s angular distance (Tape & Tape 2016) and within the so-called v-w space to quantify the probability that the true solution lies in the vicinity of the minimum misfit region.

The general MTUQ workflow is represented in Fig. 1.

4 APPLICATIONS

We apply the MTUQ methods and software to three seismic events, highlighting some important practical issues. These examples include an intraslab earthquake, a landslide, and the underground explosion resulting from the 2017 DPRK nuclear test. For the first event, we compare results obtained by searching over different combinations of source parameters and illustrate the main visualization capabilities. The second event demonstrates the code’s applicability to non-moment tensor cases, particularly in events involving mass movements. Finally, the third example illustrates the flexible programmatic interface in a detailed comparison of source type derived from different seismic wave types and polarity data.

4.1 Case 1—2009 intraslab earthquake near Anchorage, Alaska

We begin with an event included in the MTUQ examples directory: a |$M_{\rm w}\,$|4.6 intraslab earthquake that occurred on 2009 April 7, north of Anchorage (⁠|$61.45^\circ$|N, |$149.74^\circ$|W) at 20:12:55 UTC. This event was chosen for several reasons, including its magnitude (appropriate for a point source approximation), its data coverage, and because it was previously studied in Silwal & Tape (2016), which can serve as a baseline for comparison. Incidentally, this event also shares similarities in terms of location, orientation, and source type with the |$M_{\rm w}\,$|7.1 earthquake that occurred on 2018 November 30 (West et al. 2019).

We perform three different inversions using a double-couple grid (DC), a deviatoric grid (DEV), and a full moment tensor grid (FULL), each including a magnitude search. The DC grid corresponds to the centre point of the eigenvalue lune, the DEV grid covers the equator of the lune, and the FULL grid spans the entire lune (Fig. 2c). This progression from DC to DEV to FULL grids reflects an increasing dimensionality of the source parameter space. The centroid depth is fixed at 33 km, based on the Alaska Earthquake Center catalogue. We removed the instrument response and rotated the horizontal components to the radial and transverse directions for data processing. For this and all the following examples, we used ObsPy for data pre-processing and applied causal (one-pass) filters exclusively.

We use a 15 s window for P waves, bandpass filtered between [3, 10] s and a 150 s window for surface waves, bandpass filtered between [16, 40] s. Both windows begin at the predicted (or manually picked) arrival time, determined using Taup (Crotwell et al. 1999) with a specified Earth model (in this case, ak135f Montagner & Kennett 1996). The allowable time-shifts are |$\pm$| 2 s for P waves and |$\pm$| 10 s for surface waves, with Rayleigh and Love waves treated independently because they travel at different group velocities. We search for magnitude from |$M_{\rm w}\,$|4.4 to |$M_{\rm w}\,$|4.7 over regular increments of |$M_{\rm w}\,$|0.1.

In Fig. 2, we show the waveform comparison figure headers output by MTUQ for the three cases (DC, DEV, FULL). The headers list the best-fitting source and the inversion settings, including the solver and the Earth model. Because these settings are the same for all cases, differences between the best-fitting solution are related to the dimensionality of the source parameter space. We see that waveform misfit decreases with increasing grid dimensionality, as expected from the addition of free parameters. These results were obtained using the misfit function:

where |$i = 1, ..., n$| denotes the n source-receiver pairs, |$\mathcal {G}_{s,i}$| represents the i-th Green’s function for the surface waves time-windows, |$\mathcal {G}_{p,i}$| represents the i-th Green’s function for the P-wave time-windows, and |$\mathbf {M}$| is the moment tensor source term. |${{\bf d}}^{\rm obs}_{s,i}$| represents the i-th surface wave data and |${{\bf d}}^{\rm obs}_{p,i}$| is the i-th observed P-wave data. The terms |$\mathbf {w}_{s,i} \in \mathbb {R}^3$| and |$\mathbf {w}_{p,i} \in \mathbb {R}^2$| contain scalar weights applied to the residuals of the surface wave and P-wave windows, respectively.

Fig. 3 illustrates the misfit-parameter tradeoffs as a function of double-couple orientation. For each pair of angles in a misfit map, we consider all possible third angle values and all possible magnitude values, and then we plot the lowest misfit associated with the best-fitting parameters. (Thus, the plots are not cross-sections of the misfit function through parameter space.) The optimal solution (best orientation/angular triplet) is marked by a green circle across the three quadrants and corresponds to the value shown in the waveform plot figure header in Fig. 2.

MTUQ also allows for the simultaneous inversion of both source parameters and location by expanding the grid dimension to include the source’s coordinates in space. In this case, the full source-parameter grid will be evaluated for each source location (latitude, longitude, depth), enabling visualization of tradeoffs between source location and source parameters as shown in Fig. 4 for depth. As searching for depth (or latitude and longitude) increases the dimensionality of the grid, it directly impacts the computational cost of the grid search method. In our experience, querying Green’s functions is typically not the limiting factor for grid search approaches, as the grid size is predefined and manageable. However, this can become a bottleneck for other optimization methods, such as CMA-ES or Hamiltonian Monte Carlo, where a larger number of sampled location points is often required. An example for conducting a centroid search is also bundled with the code but is not presented in this paper. For the sake of simplicity, we keep the catalogue depth of 33 km for the subsequent results with this example event.

When moving to the full moment tensor grid, the source parameter space includes not only orientation parameters, but also source type parameters, which can be visualized using the eigenvalue lune (Tape & Tape 2012). The misfit map over the eigenvalue lune contains insights about the tradeoffs between different source types, which can be useful for source discrimination (e.g. identifying explosions or collapses) or to assess the distance of an inversion outcome to the best double couple solution (Tape & Tape 2019). Similar to the double-couple misfit maps, the lune misfit map is constructed by considering the best misfit value for each location on the lune, taking into account all possible orientations and magnitudes. The L2-norm based misfit map on the eigenvalue lune is displayed in Fig. 5.

MTUQ also supports custom misfit functions |$\mathcal {C}(\mathbf {M})$|⁠, such as joint waveform and first motion polarity misfit:

where |$\mathcal {C}_w(\mathbf {M})$| is the least-squares waveform misfit, normalized by the data norm |$||\mathbf {u}||_2$|⁠, |$\mathcal {C}_p(\mathbf {M})$| is the first-motion polarity misfit (expressed as the discrete number of mismatching polarities) over the number of observed first-motions |$n_p$|⁠, and |$\alpha$| is a weighting factor that controls the relative contribution of waveform and polarity misfit terms, with |$\alpha = 1$| corresponding to purely normalized waveform misfit and |$\alpha = 0$| corresponding to purely polarity misfit. In Fig. 6, we show the results for the DC, DEV, and FULL cases, with |$\alpha = 1$| and |$\alpha = 0.5$|⁠. These results show the benefits of including first-motion polarities as an extra constraint in the inversion and highlight the advantages of customizable misfit functions. The results also differ from those obtained using the absolute misfit in Fig. 2, as a consequence of modifying the misfit function, highlighting the importance of data normalization and weighting.

4.2 Case 2—2021 Barry arm landslide

The following example involves a surface mass movement that occurred on the active Barry Arm landslide in south-central Alaska (⁠|$61.24^\circ$|N |$147.96^\circ$|W) on 2021 August 9, at 07:45:50 UTC, featured as Barry Arm 3 in Karasözen & West (2024). This event represents one of several detected movements at Barry Arm, highlighting the site’s ongoing instability. The ability to determine force characteristics from seismic data is therefore a valuable tool for monitoring this and similar unstable slopes, particularly in remote locations where direct observations may be challenging. The point force approximation is generally favoured over a moment tensor for representing a mass moving along an inclined surface. This model was introduced by Kanamori & Given (1982) to explain the surface wave radiation patterns recorded from the Mount St. Helens flank collapse (and subsequent explosion) during the 1980 volcanic eruption.

In this example, we consider data from 35 three-component broad-band seismic stations, bandpass filtered between [16, 40] s, allowing a |$\pm 35$| s time-shift to account for imprecision in the onset time. The source depth is fixed at 0 km depth within the ak135f 1-D Earth model (Montagner & Kennett 1996). Due to the deficiency of short-period energy in the available data, the inversion was performed with surface waves only. Allowing a time-shift greater than a quarter of the dominant period of our waveforms can potentially introduce polarity-reversal biases (Braunmiller et al. 1994), which is present in the misfit map (Fig. 7b). We consider force amplitude in the range of |$10^{[9, 12]}$| Newtons, with the exponent varying by 0.1 increments.

Fig. 7(a) shows a subset of waveforms for the best-fitting point force, which points predominantly downward (⁠|$\theta = 165^\circ$|⁠) with a bearing of |$\phi = 278^\circ$| (measured counterclockwise from east). The downward direction might be explained by the steep slope of the Barry Arm landslide scarp. The best-fitting force is represented by a black diamond at the piercing point on the unit sphere. The associated misfit map is shown in Fig. 7(b). The ambiguity between upward and downward force arises from the polarity-reversal-bias, whereby a sign-flip of the force is accompanied by waveforms that are time-shifted by more than a quarter of a period. Reducing the allowable time-shifts could help mitigate this effect, though it may come at the cost of flexibility in accommodating errors in the Earth model or onset timing, which can be challenging to evaluate in the case of a long duration, long-period signal.

4.3 Case 3—2017 DPRK nuclear test

On 2017 September 3 the Democratic People’s Republic of Korea (DPRK) conducted an underground nuclear test at the Punggye-ri Nuclear Test Site. This test marked the largest (⁠|$M_{\rm w}\,$|5.24; Wang et al. 2018) in a series of six nuclear tests conducted over the previous years (Liu et al. 2018). A second, lower magnitude event was detected 8.5 min after the main shock, likely due to a collapse at the test site (Alvizuri & Tape 2018; Chiang et al. 2018), and background seismicity was observed in the weeks following the test (Tian et al. 2018).

The 2017 event has been extensively studied, as detailed in Table 2. Research efforts have focused on yield estimation (Lay 2022; Stroujkova 2018), event location (Myers et al. 2018; Pasyanos & Myers 2018), and characterization of the source parameters and post-test collapse mechanism using advanced seismic analysis (Yao et al. 2018; Walter et al. 2018; Wang et al. 2018; Liu et al. 2018; Chiang et al. 2018; Alvizuri & Tape 2018; Hong et al. 2019; Hu et al. 2023) and interferometric SAR imaging (Chi-Durán et al. 2021, 2024).

We revisit this event with a larger seismic data set and a computationally efficient workflow based on the CMA-ES algorithm in addition to the standard grid search method.

4.3.1 Data acquisition and processing

The data used for this event include regional body waves, surface waves, polarities, and teleseismic polarities, which have all been included in previous studies, though not yet together (Table 2). While prior studies have used similar data types in various combinations, our workflow separates body- and surface-wave windows, weights each data type independently, and integrates polarities from both regional and teleseismic stations in a single inversion. By analysing each data type separately, we investigate the source constraints each one provides and its respective contributions to the overall misfit surface.

For the body wave and surface wave analyses, we used waveforms from 31 regional stations within |$20^{\circ }$| of the epicentre (⁠|$41.3^{\circ }, 129.078^{\circ }$|⁠) using only publicly available networks. The three-component regional data are bandpass filtered between 4–15 s for the body waves and 15–50 s for the surface waves. We used 25 s P-wave time windows, with the start times determined from the theoretical P waves arrival times using the TauP method (Crotwell et al. 1999) in the ak135f model (Montagner & Kennett 1996). Surface wave time windows were set to 250 s with a start time controlled by a constant linear moveout of 3.8 km s⁻¹.

For the first-motion analysis, we added data from 153 teleseismic stations between |$20^{\circ }$| and |$100^{\circ }$|⁠. We picked 10 first motions from the regional data set and 46 polarities from the teleseismic data set. All picked polarities were upward on the vertical raw (velocity) waveforms.

To efficiently generate synthetic waveforms, we worked with pre-computed Green’s functions simulated using the 1-D model MDJ2 (Ford et al. 2009b), with a 1 s dominant period, |$20^{\circ }$| lateral extent and 700 km depth extent. We used the spectral-element solver AxiSEM (Nissen-Meyer et al. 2007, 2014) for the Green’s function simulations and Instaseis (van Driel et al. 2015) for subsequent synthetic seismogram generation.

To correct phase misalignment biases (i.e. cycle skipping), we employed a two-step process. First, we performed a preliminary inversion that allowed for large time-shifts across all stations to identify systematic delays. Time-shift corrections were then derived by ensuring geographic consistency in the distribution of time-shifts and by visually inspecting waveform fits for evidence of phase misalignment. This step confirmed that eastward oceanic paths toward Japan required up to 30 s of static corrections (Alvizuri & Tape 2018), while the continental paths—including to station MDJ—did not require any static correction. The custom static corrections were applied to selected stations to limit allowable time-shifts during the final inversion to |$\pm 2$| s for P waves and |$\pm 10$| s for surface waves. The combination of static corrections and shorter allowable time-shifts reduces computational cost and reduces the emergence of local minima having a source mechanism flipped in sign from the true solution.

4.3.2 Method

We use CMA-ES, which was previously used to invert multiple subsources of the 2022 Hunga Tonga-Hunga Ha’apai volcanic eruption (Thurin et al. 2022; Thurin & Tape 2023). Here, we present a simpler analysis using only a single source, which demonstrates how initial CMA-ES model space exploration can be used with subsequent grid searches as part of an efficient combined workflow.

In practice, we find that it is cost-effective to explore a large range of depths and magnitudes using CMA-ES, from which a smaller range can be selected for subsequent grid searches. In addition to being computationally efficient, CMA-ES is well-suited for problems where model parameters (e.g. magnitude, source type, source orientation, and location) have different scales or units. Unlike grid search methods, it operates over continuous parameter spaces rather than a pre-defined discrete set.

Our CMA-ES implementation follows Hansen (2016), with the addition of boundary-handling constraint methods from Biedrzycki (2020). Specifically, we implement a periodic boundary for the strike angle and point force bearing (which varies from |$0^{\circ }$| to |$360^{\circ }$| and can wrap around these limits) and a randomized redraw for other out-of-bound samples. A short review of the CMA-ES algorithm can be found in appendix A of Thurin & Tape (2023). We use a population of 64 samples in each CMA-ES iteration, and run the algorithm for 180 iterations (fixed), for a total of 11 520 function evaluations. We then conduct a fine-scale grid search, which included a total of |$5.977 \times 10^7$| grid points, using the magnitude estimate from the CMA-ES inversion. We fix the source centroid at |$41.30^{\circ }$|N, |$129.08^{\circ }$|E, and 1 km depth.

4.3.3 Results for a moment tensor

The misfit map obtained from 180 iterations of CMA-ES is displayed in Fig. 8, which reveals how CMA-ES focuses preferentially on regions with low misfit. Once CMA-ES identifies a descent direction, most subsequent samples are clustered in the upper portion of the eigenvalue lune, consistent with an isotropic (explosive) source.

The source that minimizes the misfit for our combined data set is an explosion-like moment tensor, with |$M_{\rm w}\,$|5.44, determined by optimizing over a continuous magnitude range of |$M_{\rm w}\,$|[5.1–5.9]. This result then informs the subsequent grid search, narrowing the magnitude range to |$M_{\rm w}\,$|[5.42, 5.44, 5.46]. This focused approach prevents searching over an unnecessary large parameter space. In this way, the CMA-ES results serve as a prior constraint for the grid search. We perform the two-part inversions (CMA-ES and grid-search) using regional body waves (BW), regional surface waves (SW), regional first-motion polarities (RP), and teleseismic first-motion polarities (TP) with equal weights (w = 0.25) to ensure a balanced contribution to the normalized joint misfit.

Fig. 9 shows the best-fitting moment tensor and a selection of waveform fits from the fine grid search. Since each misfit is evaluated independently, MTUQ allows us to combine them in a post-processing stage, enabling the computation of misfit maps for multiple scenarios once the grid search is complete, as shown in Fig. 10. From left to right, the figure shows the normalized misfit for [BW] alone, [SW] alone, the combination of [BW, SW], [BW, SW, RP], and finally, the combination of [BW, SW, RP, TP]. Underneath each misfit map, we plotted the confidence curve (Tape & Tape 2016) for a fixed magnitude (from the optimal value of |$M_{\rm w}\,$|5.44) and a uniform random grid. Specifically, the probability density function is expressed as:

where |$\mathbf {M}_0$| is the best-fitting moment tensor, and |$\mathcal {M}$| represents the sample space.

The parameter |$\sigma$| represents the estimated standard deviation of the data errors, which is estimated differently for waveform and polarity data. For waveform misfits, it is determined by evaluating the standard deviation of the residual of the best-fitting solution. For polarity data, |$\sigma$| is estimated via a Beta-binomial update that incorporates the frequency of observed polarity mismatches. We begin with an uninformative |$Beta(1, 1)$| prior on the mismatch probability, which is then updated to a posterior distribution based on the number of observed mismatches. This procedure regularizes the estimate of |$\sigma$| by incorporating prior uncertainty and thereby prevents |$\sigma$| from collapsing to zero. When combining misfits from different data groups, we normalize each misfit by dividing it by the norm of its corresponding observed data vector, then take the arithmetic mean of these normalized misfits. The total variance is obtained by summing the individual variances, each scaled by its corresponding data norm.

As we add data constraints, the area of low misfit reduces, thereby lowering the uncertainty of the results (Fig. 10). This impacts the shape of the confidence curves, which evolve towards a |$\Gamma$|-like shape having a larger area under the curve and, hence, a larger confidence parameter |$\mathcal {P}_{AV}$|⁠. Although polarities alone may not be sufficient to determine a unique best-fitting solution, they significantly contribute when combined with waveform data (Ford et al. 2012; Chiang et al. 2014; Alvizuri et al. 2018; Chi-Durán et al. 2021, 2024). Moreover, including both regional and teleseismic polarities appears beneficial, as the rightmost misfit map (Fig. 10) has the narrowest and tightest low misfit region of all maps (even though the confidence parameter only increases marginally when including teleseismic polarities). The complete set of figures for the 5 inversions can be found in supplementary Figs S14 to S25.

It is helpful to visualize variations in station attributes on a map, illustrating how azimuths, distance or site effects might affect the inversion outcome. The time-shift maps in Fig. 11 highlight the need for time-shift corrections, as all the paths crossing the Sea of Japan require up to |$-30$| s time-shifts in the MDJ2 model, which is consistent with observations in Dreger et al. (2021). Notably, Rayleigh wave time-shifts demonstrate strong azimuthal coherence across the Sea of Japan, whereas Love wave time-shifts show more variability. This variation in Love wave time-shifts might be attributed to the generally lower signal-to-noise ratio in the horizontal components of the waveforms or to the higher sensitivity to shallow crustal structure of Love waves.

4.3.4 Results for a point force

The previous results demonstrate that an all-positive-eigenvalue moment tensor can fit the considered seismic data sets. Motivated by the ambiguity between shallow explosive-like moment tensors and a vertical downward force applied to the Earth’s surface in volcanic explosions (Kanamori & Given 1982), we perform a point force inversion of the DPRK nuclear test. Using the same data sets as in Section 4.3.3, we perform a grid search over all force orientations, with a force amplitude in the range of |$10^{[10, 14]}$| Newtons, with the exponent varying by 0.1 increments, and again at a fixed depth of 1 km and assuming the MDJ2 1-D Earth model. Since the point force grid is much smaller than the full moment tensor grid, using CMA-ES here is less advantageous, and we can directly search with a fine amplitude increment instead.

The waveform fits obtained from the joint inversion of [BW, SW, RP, TP] and associated misfit maps are shown in the Supplementary Figs S29 and S30. The source that best fits the data set is a downward force, with waveform fits that are qualitatively comparable to the moment tensor solution (albeit a higher absolute waveform misfit).

This result suggests that these two models—moment tensor and point force—can both explain the seismic recordings from anthropogenic explosions and from volcanic explosions, as documented by Thurin & Tape (2023) for the Hunga Tonga-Hunga Ha’apai volcanic eruption. The canonical models (moment tensor, force) in their respective fields contribute to this ambiguity. The main feature that distinguishes these two models are the near-field, high-frequency waveforms, aese two models are the near-fieldnd the Love wave radiation patterns: two-lobed for a point force (Kanamori et al. 1984) and four-lobed for a moment tensor. In practice, resolving this ambiguity is challenging due to limited data availability and to the weak amplitudes typically associated with shallow explosive sources. Note that the two conceptual models have been combined into a compound source model in Hu et al. (2024), where a vertical force and an isotropic moment tensor source are shown to explain the main seismic source of the 2022 Hunga Tonga volcanic sequence, highlighting the interest in discriminating source models in complex geophysical settings.

5 DISCUSSION

We presented three distinct applications: an earthquake, a landslide event, and an underground nuclear explosion. The first two cases draw directly from the examples packaged with the software itself, while the third case illustrates a more involved workflow, combining multiple software elements in a new way to understand the contributions from different data types to overall seismic source constraints.

Layered (1-D) models often produce satisfactory results within the context of the CAP method (particularly if using long-period data), but their limitations become evident when large time-shifts are required to align observed and synthetic waveforms at specific stations or azimuthal ranges, as seen for stations across the Sea of Japan in the DPRK test (Fig. 11). These large time-shifts highlight the inadequacy of using a single 1-D model in capturing all path-dependent effects. Although manual adjustments can mitigate these issues, they underscore the need for more sophisticated modelling approaches. MTUQ is compatible with 3-D Green’s function databases derived from numerical solvers such as SPECFEM3D (or SEISGen for Strain Green’s functions), offering a pathway to improve forward modelling accuracy in regions with complex structures, where the importance of accounting for 3-D effects has been well-demonstrated (Ramos-Martínez & McMechan 2001; Liu et al. 2004; Covellone & Savage 2012; Hejrani et al. 2017; Takemura et al. 2020; Sawade et al. 2022)

Previous studies and ongoing efforts have integrated 3-D Green’s functions with MTUQ using its existing interfaces, both for source-side moment tensor inversions (Rodriguez Cardozo et al. 2022; McPherson et al. 2023) and for receiver-side moment tensor inversions (Ding et al. 2021), which requires saving the time-dependent volumetric wavefield within a target source region. These studies have leveraged 3-D Green’s functions to validate expected improvements, including better waveform fits, reduced time-shifts, and moment tensor solutions that align more closely with double-couple mechanisms for tectonic earthquakes.

However, one challenge with many existing 3-D Earth models is that they rely predominantly on long-period seismic data, and higher resolution models need to be developed to improve waveform fitting in the 15–40 s period range relevant for most regional moment tensor applications. Another challenge is the relatively limited understanding of how uncertainties in 3-D Earth models propagate into moment tensor solutions. While model uncertainties have been studied in 1-D contexts (Yagi & Fukahata 2011; Duputel et al. 2014; Phạm & Tkal̆cić 2021), their treatment in 3-D inversion frameworks is under active study (Pham et al. 2024). Future work in this area is critical to ensure that the improved accuracy provided by 3-D models is accompanied by robust uncertainty quantification. A secondary advantage of incorporating 3-D Earth models is the potential to reduce tradeoffs between the centroid and time-shifts. Working directly with 3-D models should improve joint inversions of centroid and source mechanism parameters (Takemura et al. 2020; Ding et al. 2021).

Another key aspect of MTUQ is its compatibility with multiple inversion strategies, allowing users to go beyond exhaustive grid searches for parameter estimation. While grid searches provide a complete exploration of the parameter space, they become computationally expensive as the dimensionality of the problem increases. For example, performing a joint inversion of source parameters and centroid locations dramatically expands the search space. Alternative optimization strategies like the CMA-ES or HMC offer significant computational advantages: CMA-ES can be used as an efficient alternative to grid search when parallel computational resources are not available. These methods, already integrated into MTUQ via compatible plug-ins, are well-suited for exploring misfit spaces with a few dozen of parameters. Pairing these global optimization approaches with refined grid searches, as illustrated in Case 3 (Section 4.3), enables users to focus computational efforts on the most promising regions of the parameter space. These efficient techniques provide opportunities to explore novel applications, such as multi-event inversions or real-time monitoring.

Another promising avenue for future development lies in expanding MTUQ’s capabilities to handle compound sources (combinations of moment tensor and single force). This extension would potentially allow the framework to discriminate ‘exotic’ event types—such as shallow volcanic explosions, landslides or exotic anthropogenic sources—from regular earthquakes. Such a compound source model has been used to invert a potential model for the 2022 eruption of the Hunga Tonga volcano (Hu et al. 2024) and would enhance MTUQ’s ability to deal with data from complex geophysical processes. These developments are important for improving source discrimination in challenging scenarios where mixed source mechanisms contribute to the observed seismicity, such as volcanic eruptions with mass movements or shallow explosions producing secondary collapses.

6 CONCLUSION

We have reviewed the design and capabilities of MTUQ, an open-source software package for seismic source estimation and uncertainty quantification. Through three test cases, we illustrated how the code can be applied to moment tensor and point source estimation problems. Other recent efforts have illustrated additional capabilities and use cases, such as for earthquake seismology (Mahanta et al. 2024), planetary seismology (Maguire et al. 2023), and volcano seismology (Kintner et al. 2022; Thurin et al. 2022; Thurin & Tape 2023).

We aim to continue adding features, such as better ways of accounting for data and model covariance, support for compound sources (combinations of force and moment tensors), and the ability to model multiple moment tensors from different locations with slight time-shifts to capture multi-event or subevent scenarios. We also hope to continue offering virtual and hybrid workshops to attract and sustain MTUQ users. Previous efforts include the 2022 MTUQ virtual workshop and the tutorials at the 2024 SCOPED workshop (Denolle et al. 2025).

In summary, MTUQ provides a flexible source inversion toolbox aimed at both research and operational applications. For reference, the data sets and scripts used for Cases 1 and 2 above are bundled directly with the package.

ACKNOWLEDGMENTS

We thank Editor Andrew Valentine and the reviewers—Arthur Rodgers and an anonymous reviewer—for their helpful comments, feedback, and suggestions. This project was supported by Air Force Research Laboratory contract HQ0034-20-F-0284 and the TREMOR project by contract FA9453-24-9-0001, with additional support from the National Nuclear Security Administration, Defense Nuclear Nonproliferation Research and Development, and the UAF Geophysical Detection of Nuclear Proliferation (GDNP) University Affiliated Research Center (UARC). MTUQ is featured in the SCOPED project, which is supported by NSF Award OAC 2104052. The manuscript was approved for release with LA-UR 24-30535. We thank Ezgi Karasözen and Michael West for recommending the Barry Arm landslide event in Section 4.2 and for curating the data set. We thank Tom VanDeMark, Aileen Fisher, Robert Walker, and Mike Cleveland for useful discussions that helped improve the usability of the package. We acknowledge the University of Alaska Fairbanks (UAF) Research Computing Systems for an allocation on the `chinook' high-performance computing cluster.

DATA AVAILABILITY

• Scripts for reproducing the grid-search results and figures from Cases 1 to 3, along with the best-fitting moment tensor parameters, weight files, first-motion polarity picks, and data for Case 3, are available in the Zenodo collection at (10.5281/zenodo.13868450).

• The github organization at https://github.com/mtuqorg contains the moment tensor code, MTUQ, including bundled data for Case 1 (Silwal & Tape 2016) and Case 2 (Karasözen & West 2024). MTUQorg contains other modules for MTUQ, including SeisHMC and SEISGen, for performing receiver-side source estimation using 3-D Strain Green’s functions (Ding et al. 2021). The module mtuq_cmaes features the CMAES approach for sampling model parameter space (Section 4.3.2). The package Pysep was used for fetching observed waveforms and for pre-processing; it is available at https://github.com/adjtomo/pysep.

• Seismic waveform data used in this study are from the AK (10.7914/SN/AK), AT (10.7914/SN/AT), AV (10.7914/SN/AV), G (10.18715/GEOSCOPE.G), IC (10.7914/SN/IC), II (10.7914/SN/II), IU (10.7914/SN/IU), JP, KG,KS, TW (10.7914/SN/TW), and YV (10.7914/SN/YV_2006). These DOIs were last accessed from https://www.fdsn.org/networks on March 2025. Waveforms were obtained from the EarthScope Data Management Center.

• For the DPRK inversion (Section 4.3), we did not use the DB network stations featured in Wang et al. (2018), because we were unable to obtain the correct metadata needed to remove the instrument response.

REFERENCES