# Global centroid moment tensor solutions in a heterogeneous earth: the CMT3D catalogue

Lucas Sawade, Stephen Beller, Wenjie Lei, Jeroen Tromp

► To cite this version:

Lucas Sawade, Stephen Beller, Wenjie Lei, Jeroen Tromp. Global centroid moment tensor solutions in a heterogeneous earth: the CMT3D catalogue. Geophysical Journal International, 2022, 231 (3), pp.1727-1738. 10.1093/gji/ggac280. hal-04920704

HAL Id: hal-04920704

https://hal.science/hal-04920704v1

Submitted on 4 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Global centroid moment tensor solutions in a heterogeneous earth: the CMT3D catalogue

Lucas Sawade<sup>1</sup>,<sup>1</sup> Stephen Beller,<sup>1,2</sup> Wenjie Lei<sup>1,3</sup> and Jeroen Tromp<sup>1</sup>

Accepted 2022 July 21. Received 2022 July 5; in original form 2021 December 30

## SUMMARY

For over 40 yr, the global centroid-moment tensor (GCMT) project has determined location and source parameters for globally recorded earthquakes larger than magnitude 5.0. The GCMT database remains a trusted staple for the geophysical community. Its point-source moment-tensor solutions are the result of inversions that model long-period observed seismic waveforms via normal-mode summation for a 1-D reference earth model, augmented by path corrections to capture 3-D variations in surface wave phase speeds, and to account for crustal structure. While this methodology remains essentially unchanged for the ongoing GCMT catalogue, source inversions based on waveform modelling in low-resolution 3-D earth models have revealed small but persistent biases in the standard modelling approach. Keeping pace with the increased capacity and demands of global tomography requires a revised catalogue of centroid-moment tensors (CMT), automatically and reproducibly computed using Green's functions from a state-of-the-art 3-D earth model. In this paper, we modify the current procedure for the full-waveform inversion of seismic traces for the six moment-tensor parameters, centroid latitude, longitude, depth and centroid time of global earthquakes. We take the GCMT solutions as a point of departure but update them to account for the effects of a heterogeneous earth, using the global 3-D wave speed model GLAD-M25. We generate synthetic seismograms from Green's functions computed by the spectral-element method in the 3-D model, select observed seismic data and remove their instrument response, process synthetic and observed data, select segments of observed and synthetic data based on similarity, and invert for new model parameters of the earthquake's centroid location, time and moment tensor. The events in our new, preliminary database containing 9382 global event solutions, called CMT3D for '3-D centroid-moment tensors', are on average 4 km shallower, about 1 s earlier, about 5 per cent larger in scalar moment, and more double-couple in nature than in the GCMT catalogue. We discuss in detail the geographical and statistical distributions of the updated solutions, and place them in the context of earlier work. We plan to disseminate our CMT3D solutions via the online ShakeMovie platform.

**Key words:** Waveform inversion; Computational seismology; Seismicity and tectonics.

## 1 INTRODUCTION

Imaging Earth's deep interior using global seismic tomography relies on accurate characterizations of earthquake source mechanisms. Beginning with the first sets of centroid-moment tensors (CMTs; Dziewoński et al. 1981) and the first global tomographic waveform inversions (Dziewonski 1984; Woodhouse & Dziewoński 1984), images of Earth's interior have been continuously improving (Romanowicz 2008; French & Romanowicz 2014; Moulik & Ekström 2014; Bozdağ et al. 2016; Fichtner et al. 2018; Lei et al. 2020). This progress owes to advances in computational resources and

infrastructure (Peter et al. 2011; Krischer et al. 2016), numerical algorithms that take advantage of such resources (Komatitsch & Tromp 2002a, b; Leng et al. 2019) and flexible data selection and processing tools (Maggi et al. 2009; Krischer et al. 2015). Yet, the global CMT (GCMT) project continues to be based on normal-mode summation in a 1-D reference earth model with 3-D path corrections for surface wave phase speeds and crustal structure (Dziewoński et al. 1981; Hjörleifsdóttir & Ekström 2010; Ekström et al. 2012).

A comprehensive study of the effects of unmodelled 3-D structure in the GCMT algorithm found biases in the solutions (Hjörleifsdóttir

<footnote n="1"><sup>1</sup>Department of Geosciences, Princeton University, Princeton, NJ 08544, USA. E-mail: lsawade@princeton.edu</footnote>

<footnote n="2"><sup>2</sup>Géoazur, Université Côte d'Azur, 06905 Sophia Antipolis, France</footnote>

<footnote n="3"><sup>3</sup>Lamont-Doherty Earth Observatory of Columbia University, New York, Palisades, NY 10964, USA</footnote>
& Ekström 2010). Both the centroid time and depth are systemati-
cally too large when comparing GCMT results to synthetic forward
modelled data in a 3-D regime. In other words, an inversion for a
new CMT solution using the GCMT solution as a starting model in
a fully 3-D heterogeneous regime is expected to both decrease the
centroid time and the depth by ~1 s and ~5 km, respectively; the
trade-off between the two parameters is negligible. Furthermore,
the algorithm also exhibits a trade-off between scalar moment and
tectonic setting. In regions with thicker crust, the GCMT algorithm
overestimates the scalar moment by up to 30 per cent, whereas in
regions where the crust is thin, for example at mid-ocean ridges, the
scalar moment is underestimated by up to 15 per cent. This differ-
ence is slightly mitigated using crustal corrections, but not resolved
entirely (Hjörleifsdóttir & Ekström 2010).

Despite its biases, the GCMT solution catalogue has become a
staple in the global seismologist's tool box due to its robustness and
quality. Changing the workflow and computational approach is chal-
lenging. Forward modelling of earthquakes in a 3-D earth model
is computationally costly. Inverting for moment tensor, hypocentre
and centroid time requires 10 simulations per earthquake using the
methods proposed by Liu *et al.* (2004). Attempts have been made
to utilize the adjoint method for source inversion, thereby reduc-
ing the number of simulations to two per earthquake per iteration,
but Kim *et al.* (2011) showed that the point source is ill-suited for
back-propagation of the residual if the structure of Earth's interior is
insufficiently known. Thus, no attempts have been made to create a
routine workflow for the inversion of earthquake source parameters
using full 3-D heterogeneous Green's functions because the com-
putational cost has remained too high to justify correcting GCMT
biases due to lower accuracy forward modelling.

Here, we build a database that takes GCMT solutions as a starting
model for the sources but uses a 3-D earth model for the wave propa-
gation. Aiming to exceed the quality of GCMT solutions but keeping
the robustness, we mimic the GCMT workflow as much as possi-
ble, except for using forward modelling algorithms and waveform
data containers that support high-performance-computing (HPC).
A remaining question is which earth model should be used for the
forward computation of the synthetics. While there are a great deal
of good models (many shown in Moulik & Ekström 2014; French
& Romanowicz 2014; Hosseini et al. 2018; Lei et al. 2020), model
GLAD-M25 (Lei et al. 2020) was ultimately chosen for two reasons.
First, we have at our disposal synthetic data for thousands of events
from the GCMT catalogue. Secondly, this means we can conduct
the source inversions in an earth model that is compatible with our
tomographic studies, which is desirable (Valentine & Woodhouse
2010).

We use SPECFEM3D_GLOBE, with its well-documented capabil-
ity of including 3-D crustal structure, topography and bathymetry,
ocean loading, rotation and self-gravitation, while still being more
computationally efficient than other options (Komatitsch & Tromp
2002b; Leng et al. 2019; Lei et al. 2020). As such, our inversions
will have the potential to become a routine algorithm and an ad-
ditional functionality to the ShakeMovie framework (Tromp et al.
2010). Currently, ShakeMovie provides the seismology community
with routinely calculated 3-D synthetic seismograms for every new
GCMT solution, as they are added to the catalogue. Hence, it is nat-
ural for us to provide routine updates to these source mechanisms,
taking them as the starting solutions, but using Green's functions
computed in a 3-D heterogeneous earth model.

We describe the events, stations and forward modelled data used
in the inversion in Sections 2, 3 and 4. Then, in Section 5, we address
data selection and measurement, showing where we encounter the

bias. In Section 6, we delineate how we can reduce this systematic
misfit. Results and discussions are presented in Sections 7 and 8.

## 2 GCMTs

A CMT solution is a representation of earthquake point sources that
contains the components of the symmetric moment tensor **M** =
[M<sub>rr</sub>, M<sub>θθ</sub>, M<sub>φφ</sub>, M<sub>rθ</sub>, M<sub>rφ</sub>, M<sub>θφ</sub>], a spatial location **r** = [z, θ, φ]
(deep, colatitude and longitude), a centroid time t<sub>CMT</sub> and a half-
duration *h*<sub>dur</sub> describing a Gaussian moment-rate function or a
source–time function describing the time-dependent energy release
during the event. The word ‘centroid’ is used because a CMT solu-
tion captures the location of the geometrical centre of the seismic
energy radiation.

Moment tensors can be used in a variety of geophysical research.
At local to regional scales, moment tensors have the ability to de-
scribe complex fault geometry (e.g. Vavryčuk *et al.* 2021), and
they are integral to understand the regional stress field by revealing
the fault process of complex earthquake sequences (Wang & Zhan
2020, e.g. Ridgecrest sequence 2019). At global scale, the most
notable use of CMTs is in global waveform tomography.

Forty years ago, the GCMT project began to invert for the first
GCMTs at Harvard University, and was led by Adam Dziewoński
for 20 yr, before being taken over by Göran Ekström. The main re-
search product of this ‘Harvard CMT Project’ is the ‘Harvard CMT
Catalog’. It was transitioned to the Lamont-Doherty Earth Obser-
vatory at Columbia University and continued as the ‘the Global
CMT Project’ and the ‘Global CMT Catalog’. The main goal of the
project is to take care of, improve and grow the existing CMT cat-
alogue. These CMT solutions are distributed through a simple, yet
efficient, web interface (www.globalcmt.org, Ekström *et al.* 2012).

For waveform tomography, it is essential to have an accurate
representation of forward modelled waveforms, which, requires an
accurate representation of the source (Dziewoński *et al.* 1981). A
CMT is by definition the ‘best-possible’ point source of an earth-
quake (Dziewonski *et al.* 1981). Since the first development of the
GCMT workflow, GCMTs have been commonly used in tomogra-
phy workflows (Dziewoński *et al.* 1981; Li & Romanowicz 1996)
and in forward modelling studies of Earth structure at regional and
global scales (Komatitsch & Tromp 2002b; Beller *et al.* 2018; Hos-
seini *et al.* 2020).

Although GCMTs have been used in a variety of global tomo-
graphic studies (French & Romanowicz 2014; Moulik & Ekström
2014; Bozdağ *et al.* 2016; Lei *et al.* 2020)—some of which also
updated the GCMT solutions—the traditional GCMT has certain
shortcomings that can be mitigated by the use of 3-D Green’s func-
tions. First, the current GCMT approach does not account for the
elastic structure of the Earth at the source or the receiver. At the
receiver, this primarily leads to an error in the predicted wave am-
plitude (Abercrombie & Ekström 2001). Secondly, the GCMT pro-
cedure does not account for the elastic focusing and 3-D variations
in the global attenuation structure, also leading to amplitude errors
(see, e.g. Dalton & Ekström 2006). Third, off-great-circle propaga-
tion is not accounted for, potentially leading to arrival angle and/or
take-off angle deviations, and hence orientation of the moment ten-
sor (Foster *et al.* 2014). Finally, the effects of rotation, leading to
Love–Rayleigh wave coupling, are not accommodated.

For the purpose of this paper, we select a subset of 9382 events
from the GCMT catalogue from the year 2000 up until February
2021 with moment magnitudes band-limited to the 5.7 ≤ M<sub>w</sub> ≤ 7.5
range. The geographical distribution of events is shown in Fig. 1(a).
<figure><img src="image_3.png" /><figcaption>**Figure 1.** (a) Selection of 9382 events from the GCMT catalogue which are considered in this article. The events are chosen based on moment magnitude $M_w$, where $5.7 \le M_w \le 7.5$. (b) Geographical distribution of stations and networks used for inversion. (c) P-wave speed structure at 100 km depth according to tomographic model GLAD-M25, which is used in the forward modelling.</figcaption></figure>

The subselection largely removed intracontinental and ridge events, which are mostly of lower magnitude. We chose the events based on requirements for global tomography and the largest magnitude at which an event can be viewed as a point source for global forward modelling using the spectral-element-method.

## 3 DATA

Our inversion aims to reduce the waveform difference between observed and synthetic three-component displacement seismograms. Observed data are downloaded from the Incorporated Research Institutions for Seismology (IRIS) Data Management Center (DMC; Trabant *et al.* 2012). Following Ekström *et al.* (2012), seismograms are downloaded from several seismographic networks, which provide ample global coverage, namely, the Caribbean USGS Network (CU), GEOSCOPE (G), GEOFON (GE), the New China Digital Seismograph Network (IC), the Global Seismic Network (II, IU) and the Mediterranean Very Broadband Seismographic Network (MN). For each station, a window of 11 500 s from the broadband channel (LH\*) is downloaded, starting 300 s prior to the centroid moment time provided by the GCMT starting solution. Three types of seismic waves are used in the inversion, which are defined by both their wave characteristics and period band, namely, (1) body waves ($T > 40$ s), (2) intermediate-period surface waves ($T > 50$ s) and (3) mantle waves (long-period, surface waves that sense the deep-mantle; $T > 125$ s), each chosen for their own benefits. Body waves carry information about the focal sphere by sampling it uniformly in three dimensions. Intermediate-period surface waves are essential for the inversion of smaller earthquakes due to their good signal-to-noise ratio. This type of wave is not considered for events deeper than 300 km, since intermediate surface waves generated by such events are too small. Finally, mantle waves are very-long period in nature and therefore insensitive to the spatial extent of the earthquake; for example a $M_w$ 7.0 earthquake has a surface rupture length of approximately 60 km (Wells & Coppersmith Kevin 1994), which is small compared to the wavelengths of the corresponding mantle waves. Hence, mantle waves are of great value in the inversion of earthquakes with larger magnitudes. Furthermore, we include records of mantle waves up to 4.5 hr because minor- and major-arc surface waves sample the earthquake in two distinct directions, adding further constraints to source inversions (Ekström *et al.* 2012). After a preliminary check of record quality and conversion to displacement by station response removal data are processed as follows.

We largely follow the processing choices made in the GCMT Project to determine the effects that a 3-D earth model has on the

inversion of a CMT (Ekström *et al.* 2012). Each wave type is processed separately. First, data are cut to the desired length in time. Next, data are detrended to eliminate station drift and demeaned. Data are tapered to prevent any ringing in subsequent Fourier analysis, and then filtered, interpolated and rotated to vertical, radial and transverse components with respect to the geometry between event and stations. Tables containing the exact parameters used by the processing are provided in Supplement S.4.2. Each wave type is processed to different record lengths, but interpolated to the same time-steps from the start of the event at a 0.5 Hz sampling rate.

## 4 SYNTHETICS

We calculate synthetic seismograms using model GLAD-M25 (Lei *et al.* 2020), which was constrained by data from 1480 events $5.5 \le M_w \le 7.2$, and the global 3-D Spectral-Element-Method (SEM) software package SPECFEM3D_GLOBE. The routine computation of heterogeneous 3-D synthetics is possible as a result of the highly parallelizable nature of the SEM (Komatitsch & Tromp 2002b) and its symbiosis with Graphics Processing Units (GPU), which have become essential for fast, parallel, seismic forward simulations (Komatitsch *et al.* 2009, 2010; Komatitsch 2011).

We are using a downsampled version of GLAD-M25. Since the minimum period of interest here is 40 s, the spectral-element mesh does not need to resolve nearly as much detail as the mesh used in Lei *et al.* (2020). To save computational resources and simulation time, the cubed-sphere mesh was downsampled from $6 \times 256 \times 256$ (number of chunks × number of spectral elements along the surface of one side of a chunk squared) spectral elements at the surface to $6 \times 128 \times 128$ surface elements. The resampling allows for reduction of required computational resources from 20 to 2 GPU-h for the computation of a 3 hr simulation; for example a 3 hr simulation on a 6 GPU node now takes 20 min. The downsampling utilizes the Lagrange basis functions of the original mesh to find the values of the new, lower-resolution mesh. Low-pass filtered to 40 s, the waveforms of the old and new mesh are indistinguishable to numerical precision, as expected. Except for instrument response removal, synthetic seismograms are processed identically to observed data. After both forward and observed data are in place, we select data windows based on waveform similarity, as we discuss next.

## 5 DATA SELECTION

GLAD-M25 is a modern model of the Earth, there are parts of observed traces that are still too complicated to be explained by a global tomographic model and/or have such low signal-to-noise
ratio that a comparison in waveform is only a comparison of noise (Maggi et al. 2009; Lei et al. 2020). Hence, windows have to be selected to assure that we capture how a point-source affects waveforms and avoid explaining noise or waveform anomalies introduced by heterogeneities not present in GLAD-M25. Observed data are windowed according to their similarity with synthetic data for the three different wave types. That is, a set of misfit criteria are met for segments (windows) of the data and synthetic trace-pairs, where data and synthetics are ‘close enough’. We use the window selector FLEXWIN’s Python implementation Pyflex (Maggi et al. 2009; Krischer & Casarotti 2015) to find such windows. Tables containing the exact parameters used for the window selection are presented in Supplement S.5. Pyflex uses two main misfit measurements to select windows, namely, the cross-correlation traveltime delay and the amplitude misfit, among many other measurements, which are beyond the scope of this paper. The cross-correlation traveltime misfit

$$ \Delta t = \underset{\tau}{\operatorname{argmax}} \int_{t_0}^{t_1} s(t + \tau) d(t) dt, \qquad (1) $$

is a measure of delay between two waveforms, s and d, where a positive time shift indicates that the data are positively shifted in time compared to the synthetics. The amplitude measurement

$$ d \ln A = \frac{1}{2} \ln \int_{t_0}^{t_1} d^2(t) dt / \int_{t_0 + \Delta t}^{t_1 + \Delta t} s^2(t) dt, \qquad (2) $$

where $t_0$ and $t_1$ denote the start and end time of the windows, and where $\Delta t$ indicates that we are measuring the power ratio between the observed data and the cross-correlation time-shift corrected synthetics. A positive dln A indicates that the power of the observed data is larger than the power of the synthetic data. In the inversion (Section 6) we use the waveform misfit, which ideally reduces the misfit in both of these categories. In Section 1, we introduced the idea of systematic bias in the misfit between the data and the forward modelled data, which is left to be explained by better source parameters. Fig. 2 shows the fit in traveltime and amplitude between the data and the forward modelled events in the subset of the GCMT catalogue used in this study. The fit is good, but there is a systematic negative time-shift in both body and surface waves, and a positive time shift for mantle waves. The amplitude misfit in the windows is positively shifted in the mean sense, indicating that the GCMT scalar moment is on average lower than the observed scalar moment. A source inversion using 3-D Green’s functions as described above should be able to centre the histograms.

# 6 INVERSION

The model parameters we aim to optimize are defined by the model vector $\mathbf{m} = [M_{rr}, M_{\theta\theta}, M_{\phi\phi}, M_{r\theta}, M_{r\phi}, M_{\theta\phi}, \theta, \phi, z, t_{CMT}]$. In this study, as in the standard GCMT inversion algorithm, we omit the inversion for half-duration, $h_{dur}$, and simply relate it to the scalar moment as done in the GCMT Project by the empirical relationship $h_{dur} = 2.4 \times 10^{-6} M_0^{1/3}$ in s, where $M_0$ is in N·m (Dahlen & Tromp 1998, Section 5.4.5).

## 6.1 Objective function

The objective function consists of a sum of waveform misfits, asymmetries of which are balanced by weights. We first describe the waveform misfit, assemble the objective function including the weights and explain the construction of the balancing weights.

The misfit function consists of normalized squared waveform differences in targeted windows, namely,

$$ \frac{1}{N} W_{w}(t) [d(t) - s(t; \mathbf{m})]^2, \qquad (3) $$

where d denotes an observed seismogram, s a corresponding synthetic seismogram for the current model $\mathbf{m}$, $W_w$, $w = 1, \dots, W$, a time window and N a normalization factor for each trace defined by

$$ N = \sum_{w=1}^{W} \int_{0}^{T} W_{w}(t) d^{2}(t) dt, \qquad (4) $$

where T denotes the length of the observed seismogram. This normalization factor puts all seismograms on an equal footing in terms of their contribution to the overall misfit. Normalized waveform misfits (3) from $k = 1, \dots, K$ wave types on $c = 1, \dots, C_k$ seismographic components recorded at $r = 1, \dots, R_{kc}$ seismographic stations in $w = 1, \dots, W_{kcr}$ time windows may be combined with weights to form the overall objective function:

$$ C(\mathbf{m}) = \frac{1}{2} \sum_{k=1}^{K} \omega_k \sum_{c=1}^{C_k} \omega_c \sum_{r=1}^{R_{kc}} \omega_{kcr} \frac{1}{N_{kcr}} \times \sum_{w=0}^{W_{kcr}} \int_{0}^{T_{kcr}} W_{kcrw}(t) [d_{kcr}(t) - s_{kcr}(t, \mathbf{m})]^2 dt. \qquad (5) $$

Our selection of categorical weights $\omega_k$, seismographic component weights $\omega_c$ and combined geographical and azimuthal receiver weights $\omega_{kcr}$ may be summarized as follows.

### 6.1.1 Categorical weights

Following (Ekström et al. 2012), we assign categorical weights, $\omega_k$, to the waveforms depending on earthquake magnitude. This weighting controls the contributions of the different wave types to the overall waveform misfit and is shown in Fig. S15(a). Here, $k = 1, 2, 3$ label body, surface, and mantle wave categories, respectively. Above magnitude 7.5, body and surface waves are not included at all due to the spatio-temporal effects that large earthquakes have on body and surface waves; this is irrelevant for the present study. From magnitudes 7.5 to 6.5, the body and intermediate surface wave weights are linearly increased until fully weighted ($\omega_{1,2} = 1.0$) at magnitude 6.5. Body and surface waves are then fully weighted ($\omega_{1,2} = 1.0$) between magnitudes 6.5 and 5.0, while the mantle wave weight is linearly decreased from magnitudes 6.5 to 5.5. Mantle waves are not included ($\omega_3 = 0.0$) below a magnitude of 5.5 because they are overwhelmed by noise. The categorical weights are normalized after all weights, such that their average is 1.0.

### 6.1.2 Component weights

We assign a component weight, $\omega_c$, to each seismographic component (R-T-Z). While we implemented the infrastructure to assign component-based weights. For the present study, these weights are set to (1,1,1) for Z-R-T.

### 6.1.3 Geographical weights

The waveform misfits are weighted depending on the geographical station distribution. This is necessary because the Northern Hemisphere is densely populated with seismic stations, whereas the
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Time and amplitude misfits between seismic waveforms of various wave types (body waves, surface waves and long-period ‘mantle’ waves) observed on three components (vertical, radial, and transverse from top row to bottom row) and their corresponding synthetics computed based on spectral-element simulations of wave propagation using GCMT earthquake mechanisms (Ekström et al. 2012) in 3-D earth model GLAD-M25 (Lei et al. 2020). Misfit measurements are cross-correlation time-shifts (first three columns) and amplitude ratios (last three columns) over relevant time windows in appropriate frequency bands. Labels to the left of the histograms indicate the wave types as explained in Section 3 for each row, and the labels along the horizontal axes indicate the components. The histograms are normalized to density, and the top-right, outside label shows the number of measurement windows *N* used to compute the histograms. The GCMT label indicates mean and standard deviation for each histogram.</figcaption></figure>

Southern Hemisphere remains poorly covered (see Fig. 1b). Geographical weights, $\omega_{geo}$, are computed for each wave type and component separately following Ruan et al. (2019):

$$[\omega_i^{geo}]^{-1}(\Delta_0) = \sum_{j=1}^{N} \exp\left[-\left(\frac{\Delta_{ij}}{\Delta_0}\right)^2\right], \quad (6)$$

where $\Delta_{ij}$ is the epicentral distance between stations *i* and *j* (symmetric with 0 diagonal), and $\Delta_0$ is a reference geographical distance that is chosen based on the minimum and maximum number of the vector $\omega_i$, called the condition number in Ruan et al. (2019).

### 6.1.4 Azimuthal weights

Next, the azimuthal station density is computed, so that misfits are balanced with respect to source azimuth, in an attempt to balance sampling of the source radiation pattern. Weights for the receivers inside a specific azimuthal bin *b* are defined as

$$\omega_b^{azi} = F \left[ \frac{1}{N_b} \right]^p, \quad (7)$$

where $N_b$ is the number of receivers of in an azimuthal bin, and *F* is a normalization factor that returns the mean of all weights to 1. The exponent *p* is a matter of preference, but *p* = 0.5 has proven to be a suitable choice. This weight is applied to all receivers within bin *b*,

$$\omega_r^{azi} = \omega_b^{azi} \quad \text{if receiver } r \text{ in bin } b. \quad (8)$$

Bins with large numbers of stations receive smaller weights <1, whereas bins with a smaller station count receive larger weights >1.

### 6.1.5 Combined receiver weights

Finally, geographical and azimuthal weights are combined such that their mean is 1:

$$\omega_r = \frac{\omega_r^{azi} \omega_r^{geo}}{\frac{1}{R} \sum_{r'=1}^{R} \omega_{r'}^{azi} \omega_{r'}^{geo}}. \quad (9)$$

If several azimuthal bins have similar weights in the geographical weighting, azimuthal weighting corrects this slightly so that weighting with respect to the source is more balanced.

At the final stage, categorical weights are normalized based on the combined weights, such that the average of the categorical weights approaches 1. These final weights are shown in Fig. 3. Note how remote stations receive higher weights and clustered stations lower weights in comparison.

## 6.2 Optimization

The misfit function (5) is minimized with respect to the model parameters **m**. We compute the Fréchet derivatives of the synthetic seismograms with respect to model parameters using the methods discussed in Supplement S.2. The Fréchet derivatives with respect to the moment tensor elements are linear, whereas the derivatives with respect to the location are non-linear and are best approximated in the basis of the forward modelling method, that is the Gauss–Lobatto–Legendre basis. The full derivation of gradients with respect to the source location is presented in Supplement S.2.2.

Given the Fréchet derivatives, both the gradient and the Gauss–Newton Hessian matrix of the cost function can be computed, and we can calculate a Newton step to update our model parameters. Note that this is only feasible due to the relatively small parameter space, as each Fréchet derivative requires one additional simulation. For larger problems, such as seismic tomography, adjoint methods are a better choice (Tromp et al. 2004).

## 6.3 Iterative inversion for depth and centroid time

The inversion for the depth of globally distributed earthquakes is close to, but not entirely, linear. This, aside from the fact that moment tensors are not well-defined at structural boundaries (Woodhouse 1981; Dahlen & Tromp 1998), is a reason why the GCMT project chooses to fix the depth of shallow earthquakes (Ekström et al. 2012). We evaluate the necessity of further constraining the
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Component-wise weight distribution for Event C201712131803A. The event with a moment magnitude of 6.5 is located at a transform-boundary between the African Plate and Antarctica beneath the Southern Atlantic Ocean, as indicated by a black star in the maps. The rows are identified at the very left for each component (Z-R-T). Each column is titled by a corresponding wave type and wave weight. Each map is accompanied by factor R, the ratio of maximum and minimum weight, and S, the sum of the weights, for the wave type/component combination. R_T is the same ratio but for all weights, and S_T the sum of all weights, which naturally is one.</figcaption></figure>

depth by iteration. We focus on depth because it is the least well-constrained parameter (Kim et al. 2011), other than dip-slip components (Kanamori & Given 1982). Furthermore, the original GCMT solution is inverted using a spherically symmetric earth model. When this GCMT solution is used in SPECFEM3D_GLOBE, including topography, the source is placed in the mesh using the depth below the surface topography. Events located in regions with high elevation or deep bathymetry may be significantly displaced. Therefore, to check whether a correction in depth is necessary, an iterative approach to depth optimization is analysed. In this approach, we use a regularized Gauss-Newton algorithm to constrain the depth and centroid time of each event to see whether an extra Newton step significantly increases waveform fit. We are using the same objective function as in eq. (6), and compute Fréchet derivatives using the same methods.

## 6.4 Correction for the amplitude discrepancy

A disadvantage of using a waveform misfit function for source inversions is that discrepancies in phase due to the wave speed model, which the source inversion does not account for, result in a loss of amplitude in the synthetics, and hence the scalar moment (Liu et al. 2004). Liu et al. (2004) introduced the idea of using a balanced cost function, involving a combination of waveforms and envelopes, as a solution to this problem. We found, however, that while this approach does decrease amplitude misfit, it can never fully account for the misfit in amplitude, unless only the envelope misfit is used—in which case phase information is fully lost. Therefore, after the waveform inversion, we choose to scale the scalar moment using

amplitude information of the longest period waveforms to correct for the loss of power in the signal. We use the correlation between cross-correlation-time-shifted synthetics and observed data on windows

$$S = \int_{t_0}^{t_1} s(t + \Delta t) d(t) dt / \int_{t_0 + \Delta t}^{t_1 + \Delta t} s^2(t) dt, \quad (10)$$

where integral annotations are the same as in eq. (2), to correct synthetics by the factor S (see, e.g. Ekström & Busby 2008). Note that we are using mantle waves because the longest-period waves are best at constraining the power (Hjörleifsdóttir & Ekström 2010; Ekström et al. 2012). If histograms are centred on zero, the power is matched. Note that a sacrifice is being made to provide a more accurate scalar moment than waveform optimization can provide. Since the correction does not require the waveform misfit to remain the same or decrease, a slight increase in waveform misfit is observed as compared to the pure waveform misfit optimization. It is important, however, that the overall, final waveform misfit remains smaller than the one computed using the starting model.

## 7 RESULTS: A NEW CATALOGUE OF CMT3D

In the following few paragraphs, we present statistics of the observed changes between the original GCMT catalogue and our inverted CMTs to confront our goal of correcting the sources so that the histograms shown in Fig. 2 are centred. In the future, we will refer to the presented catalogue as the CMT3D catalogue. However,
for the purpose of this paper, we will make the following distinc-
tions to be clear as to what is being compared: CMT3D refers to the
single iteration moment tensor inversion, CMT3D+ refers to the
CMT solution after a subsequent iterative inversion for depth, and
the modifier ‘X’ denotes the correction for the power on waveforms.
So, ‘CMT3D+X’ is the final and published catalogue, referred to as
the CMT3D catalogue in the future. An overview of the catalogue is
shown in Fig. 4. Fig. 5 compares the data fit pre- and post-inversion.

The left three columns of Fig. 5 (shift in centroid time) narrowed
in their time shift distribution, and most, but not all, component-
wavetype combinations moved towards zero. The Z and R com-
ponents of the mantle waves show both an increase in time-shift,
mostly after the iterative depth and centroid time inversion. Overall,
the centroid time is shifted by -1.05 s, with a standard deviation of
1.17 s (see Fig. 4d). In other words, the events are corrected to be
earlier in time than their starting model counterparts. The centroid
time shift is mostly apparent for events with deep crustal roots.
Fig. S8 shows that shallow (<70 km) events within Tibet move to
even shallower depths. Similarly, deep (>300 km) events in South
America are shallower than their GCMT counterpart. Events in the
Northwest to Northeast Pacific remain largely unchanged. Positive
changes in centroid time are only observed in the Hellenic subduc-
tion zone, the Sumatran subduction zone, and events surrounding
the Scotia Plate.

Fig. 5 shows that the power discrepancies are overall corrected to
have a zero-mean on all component-wave type combinations, to be
expected due to the amplitude corrections. Note that as a part of the
amplitude correction procedure, the waveform misfit is increased,
as shown in Fig. S5. However, the overall misfit before inversion
and of our final catalogue remains decreased (see Fig. S6), which
is the desirable result. The scalar moment of the events increases
in the mean sense by 5.61 per cent, but shows a large spread of
~17.84 per cent in standard deviation as seen in Fig. 4(e). The
spatial distribution of changes in scalar moment from the GCMT
catalogue to the CMT3D+X catalogue shows that shallow events
have the largest corrections (Fig. S7). This is an expected result that
can be attributed to the added crustal structure (Hjörleifsdóttir &
Ekström 2010). Events in regions with continental crust are nega-
tively corrected, while events underneath an ocean layer are largely
positively corrected. In the depth range from 23 to 35 km, this effect
is less pronounced and the histogram of scalar moment changes is
well-centred around zero. Deeper events are, by and large, increased
in scalar moment.

Fig. 4(f) shows an overall shallowing of the events by ~3.94 ± 4.41 km. Depth corrections appear to increase to some extent with event depth; see Fig. 4c). With increasing event depth, negative depth corrections increase as well until ~300 km. Deep events below 300 km are shallower overall, but receive a smaller correction than intermediate events; intermediate (70–300 km) events receive the largest corrections. The spatial distribution of depth corrections is shown in Fig. S9 of the supplement. The largest corrections are observed for the deepest subduction zone events, followed by shallow earthquakes beneath continents. All shallow events in subduction zones around the Pacific receive very small corrections until a depth of 35 km, where the corrections increase. It is important to emphasize that all depth ranges show overall negative corrections, with the smallest corrections in the 8–11 km range. Positive changes are rare.

Fig. S10 shows the spatial distribution of horizontal relocations.
Very shallow events around the Pacific tend towards the centre of
the Pacific except for events in South America. For intermediate
shallow events, the horizontal relocations become more haphazard,

but they reveal more systematic changes for the deeper shallow
events in the range 35–70 km. The most striking result here is the
directional change of events in South America. Up until 80 km
depth, the dominant change in direction is eastward, away from
the Pacific, but below 80 km the dominant correction of horizontal
location is westward.

The distribution of source-type (DC versus ± CLVD) shows more
DC components (ε = 0); the standard deviation changes negligibly.
The spatial distribution of the change in source type does not follow
a clear pattern. It should be noted, however, that the changes in the
source type are largest at the surface, where some moment tensors
elements are not well resolved (Kanamori & Given 1982). The
changes in source type become increasingly smaller as a function
of depth. Deep events show almost no change in source type at all.
This trend becomes evident in the histograms in Fig. S12 of the
supplement, where the largest spread is shown for shallower events
and then becomes narrower with depth.

We computed the correlation of spatially binned parameter
changes. Changes in event depth and centroid time are not show-
ing any systematic correlated behaviour, except intermediate depth
subduction zone events; the majority of intermediate depth events
are positively correlated. That is, an increase in depth will increase
centroid time-shift and vice versa (see Fig. S13). Changes in depth
and scalar moment are correlated close the surface, where a de-
crease in depth shows a decrease in scalar moment and vice versa
(see Fig. S14).

8 DISCUSSION

In the following paragraphs, we discuss some corrections seen in
Section 7.

On average, our new CMT3D catalogue tends to correct the
mean delay of simulated synthetic waveforms as compared to their
observed counterpart. Overall, we observe a mean centroid time-
shift of about 1 s (earlier) in accordance with those postulated
by Hjörleifsdóttir & Ekström (2010) and similar to the trend ob-
served in other studies (Lei et al. 2020). Mantle waves, in particular
Rayleigh waves, are often not sufficiently well represented using
SPECFEM3D_GLOBE. This leads to a much smaller number of win-
dows as compared to the other two categories. Consequently, mantle
waves have lesser priority for misfit correction, and we observe an
increase in time misfit for the Z and R components of mantle waves.

The decrease in synthetic waveform amplitude as a result of least
squares waveform inversion has been observed in multiple studies
(Liu et al. 2004; French & Romanowicz 2014; Lei et al. 2020),
and can be attributed to phase incoherence that is not correctable
by pure source inversion. There are different ways of alleviating
(e.g. envelope misfit, Liu et al. 2004), or correcting, these discrep-
ancies (e.g. amplitude, Lei et al. 2020). In this study, we decided
to apply a scaling using the correlation of synthetic and observed
seismograms (Ekström et al. 2006). In the future, one hope is that
there will be eventual convergence in waveform misfit between
both source and model through their joint improvement by sequen-
tial or simultaneous (e.g. Valentine & Woodhouse 2010) inversions.
In such a case, the loss of power in the signal can be minimized
by corrections of phase incoherence, and our correction for scalar
moment should not be required. In any case, for now, this ampli-
tude correction still provides a fairly good assessment of the energy
radiated by an earthquake. Using this information, Hjörleifsdóttir
& Ekström (2010) find that in areas with a thick crust (50 km
as compared to PREM’s 24 km crustal thickness), the GCMT’s
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Statistical overview of the inverted events. (a) Distribution of inverted earthquakes across the globe. (b) Change in source-type distribution. (c) Depth as a function of change in centroid depth $\Delta z$. Depth binned statistics are shown for this plot. The mean and bin width are indicated by the bar with dot, and the left and right bars denote the standard deviation to that. (d) Distribution of changes in centroid time $\Delta t_c$. (e) Distribution of changes in scalar moment 100 · Δ(ln $M_0$). (f) Distribution of changes in centroid depth.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Agreement between seismic waveforms of various wave types (body waves, surface waves, and long-period 'mantle' waves) observed on different components (vertical, radial and transverse from top row to bottom row) and their corresponding synthetics computed based on spectral-element simulation of wave propagation in the 3-D earth model GLAD-M25 (Lei *et al.* 2020) using [GCMT, grey patch; same as line in Fig. 2] earthquake mechanisms (Ekström *et al.* 2012) and the new inverted [CMT3D+X, coloured line] earthquake mechanisms. Measures of agreement are cross-correlation time-shifts (first three columns) and amplitude ratios (last three columns) over relevant time-windows in appropriate frequency bands. Labels on the left of the histograms indicate the wave type, as discussed in Section 3, for each row, and the labels along the horizontal axes indicate the component. The top-right label shows both mean and standard deviation for each histogram.</figcaption></figure>

scalar moment is overestimated, whereas in regions with thin crust (10 km), the scalar moment is underestimated. This leads to two effects that should be observable: (1) the correction of the scalar moment should be positive for regions with thin crust and vice versa and (2) if the distribution of events is leaning towards events with

thinner crust, an overall positive trend in the scalar moment correction should be observed and vice versa. In fact, the geospatial distribution of our results in Fig. S7 clearly confirms the trend that shallow, intracontinental earthquakes receive an overall negative correction, while events below an oceanic crust receive on overall
positive correction, highlighting the influence of 3-D crustal varia-
tions on CMT inversions. Furthermore, most events are located in
the shallow part of subduction zones and ocean ridges. Although
rapid variations in structural parameters do occur in these regions,
the majority of events happen where the crust is rather thin. Hence,
postulates (1) and (2) are confirmed: the spatial distribution tends
to events with an oceanic plate atop (thinner crust), and the overall
correction in scalar moment is, in fact, positive. Previous results
additionally suggested that the scalar moment of mid-ocean ridge
events is underestimated, and the scalar moment has to be corrected
(Hjörleifsdóttir & Ekström 2010), a suggestion also supported by
our study.

The average depth correction of –4 km that we observe is also to be expected. For instance, Hjörleifsdóttir & Ekström (2010) predicted that with the use of a fully 3-D heterogeneous earth model and solver, the depth correction would move events to shallower depth by approximately 5.5 km due to the bias introduced using the traditional GCMT approach. Additionally, in regions where the crust is substantially thicker than in the global PREM 1-D model, the error in the depth estimate is substantially larger by about ≳10 km (Hjörleifsdóttir & Ekström 2010), a statement that is underlined by our study. Indeed, some of the largest depth corrections in our catalogue are obtained for places where the crust is anomalously thick as compared to PREM. For example, the largest negative depth corrections occur for very deep subduction zone events and are most likely an effect of the updated 3-D structure in the upper mantle by GLAD-M25 as compared to the 1-D PREM model. Another interesting feature is the change in depth of strike-slip events close to mid-ocean-ridges. The GCMT algorithm generally overestimates the depth for such events (Hjörleifsdóttir & Ekström 2010). As a result, events at these locations experience a more pronounced shallowing in our inversion, and it attests to the strong influence of shallow 3-D crustal variations in the CMT estimation.

The majority of events in this catalogue already show significant
changes in depth after a single iteration for moment tensor and
hypocentre. However, the damping required to stabilize the inver-
sion for the hypocentre limits how much an event can move in depth
during a single iteration. Hence, the iterative depth and centroid time
inversion requires ~2–3 iterations on average, with a maximum of
six iterations (see Fig. S2). The effect of the iterative procedure
on depth is fairly drastic since the single iteration inversion moves
the events by 1.63 km in depth and 0.11 s in centroid time, while
the iterative procedure will move events by ~2.33 km and ~0.85 s
(see Fig. S4f). We conclude that an iterative depth inversion after a
single iteration for the moment tensor and the hypocentre is indeed
necessary.

The systematic correlation of changes in event depth and scalar
moment for very shallow events is expected, but non-existent
for deeper events. The locally, positive correlation of changes in
depth and centroid time for intermediate depth subduction zone
events is expected for travel times to stations close in epicen-
tral distance, but not for distant stations. Thus, we attribute this
overall positive correlation to updated upper mantle structure in
GLAD-M25.

Interestingly, Morales-Yáñez et al. (2020) shows similar trends
in errors in moment tensors and hypocentres as Hjörleifsdóttir
& Ekström (2010) using a different method also based on
long-period waveforms. The common trends observed by the
W-phase inversion (Kanamori & Rivera 2008; Morales-Yáñez
et al. 2020) and our approach suggest a true bias in source
mechanism and in location introduced by the traditional GCMT
approach.

9 CONCLUSION

We presented a new database of CMTs, CMT3D, inverted using
Green's functions from the fully heterogeneous 3-D earth model
GLAD-M25, using the GCMTs as starting solutions. Based on the
sequence of a single iteration inversions for the moment tensor
and hypocentre, an iterative inversion for depth and centroid time,
and a scalar moment correction, we computed a total of 9382 new
CMTs. Our results largely confirm previous predictions made by
Hjörleifsdóttir & Ekström (2010). The events are generally shal-
lower, but receive more substantial corrections where the crust is
significantly thicker than in the PREM model. Earthquakes that are
located below the ocean only experience slight changes in depth,
albeit negative as well. We attribute these changes to the incorpo-
ration of topography and crustal structure in our Green's function
calculations. A significant shallowing of events in the lower parts of
the upper mantle is attributed to a more accurate upper mantle struc-
ture. The centroid time is observed to be slightly earlier than in the
GCMT catalogue. The systematic bias in time-shifts is addressed
by the use of heterogeneous structure. The geographical correction
pattern is less pronounced for the centroid time than for depth, but
events in the Northern Pacific and South East Asia require smaller
corrections than events elsewhere. The scalar moment is increased
on average as a result of two steps. First, in a pure waveform inver-
sion, the scalar moment decreases by fitting waveforms containing
phase differences that cannot be explained by a point source. Cor-
rection of this amplitude loss results in an overall increase in the
scalar moment. Secondly, the correction of the scalar moment ex-
hibits a strong geographical pattern, where events below or within
continental crust receive negative corrections, and events below or
within oceanic crust receive positive corrections. The changes in
moment magnitude decrease with depth, but the mean correction
remains positive. Very shallow events that decrease in depth tend to
locally decrease in scalar moment; a trend that is not observed for
deeper events. Intermediate depth events show an expected, locally
positive correlation between changes in depth and centroid time.

ACKNOWLEDGMENTS

We thank Frederik J. Simons for extensive feedback on the
manuscript, and Göran Ekström and Qinya Liu for their comments
and suggestions for both methods and manuscript. This research
used resources of the Oak Ridge Leadership Computing Facility,
which is a DOE Office of Science User Facility supported under con-
tract DE-AC05-00OR22725. Additional computational resources
were provided by the Princeton Institute for Computational Science
& Engineering (PICSciE). The open source spectral-element soft-
ware package SPECFEM3D_GLOBE and the seismic measurement
software package FLEXWIN used for this paper are freely available
via the Computational Infrastructure for Geodynamics (CIG; geo-
dynamics.org). This research was supported by NSF grant 1644826.

DATA AVAILABILITY

The event catalogue is made available on the ShakeMovie website
and through a dedicated DOI (Sawade et al. 2021a). We understand
that the resources used for the source inversion are not available to
everyone to reproduce the result of this study. However, we want
to emphasize that the software used for computation is created to
make the study reproducible. The software for this publication is
available through the dedicated DOI (Sawade et al. 2021b).
The raw observed seismograms can be downloaded from IRIS
webservices (Trabant *et al.* 2012) using the provided parameters
(Sawade *et al.* 2021b). The computed synthetic seismograms are
not made available through a similar source due to their storage
volume. On demand, the processed observed and seismic data can
be made available, and for specific events, we can recompute the
seismograms if desired.

Using GLAD-M25 (Lei *et al.* 2020), and specfem3D_globe (Komatitsch & Tromp 2002b), the seismograms can also be reproduced by independent users, using the source models presented in this study. Note that the computation of the synthetic seismograms using the events from the catalogue presented in this study requires an earth model that is fully heterogeneous and incorporates topography, Bathymetry, rotation and ellipticity. Synthetic seismograms are expected to be worse if our source models are used in combination with 1-D forward solvers or solvers that do not consider the above-mentioned physics.

REFERENCES

Abercrombie, R.E. & Ekström, G., 2001. Earthquake slip on oceanic transform faults, *Nature*, **410**, 74–77.

Beller, S., Monteiller, V., Operto, S., Nolet, G., Paul, A. & Zhao, L., 2018. Lithospheric architecture of the South-Western Alps revealed by multi-parameter teleseismic full-waveform inversion, *Geophys. J. Int.*, **212**(2), 1369–1388.

Bozdağ, E., Peter, D., Lefebvre, M., Komatitsch, D., Tromp, J., Hill, J., Podhorszki, N. & Pugmire, D., 2016. Global adjoint tomography: first-generation model, *Geophys. J. Int.*, **207**(3), 1739–1766.

Dahlen, F.A. & Tromp, J., 1998. *Theoretical Global Seismology*, 1st edn, Princeton Univ. Press.

Dalton, C.A. & Ekström, G., 2006. Global models of surface wave attenuation, *J geophys. Res.*, **111**(5), 1–19.

Dziewonski, A.M., 1984. Mapping the lower mantle: determination of lateral heterogeneity in P velocity up to degree and order 6, *J geophys. Res.*, **89**(B7), 5929–5952.

Dziewoński, A.M., Chou, T.A. & Woodhouse, J.H., 1981. Determination of earthquake source parameters from waveform data for studies of global and regional seismicity, *J geophys. Res.*, **86**(B4), 2825–2852.

Dziewonski, A.M., Chou, T.A. & Woodhouse, J.H., 1981. Determination of earthquake source parameters from waveform data for studies of global and regional seismicity, *J geophys. Res.*, **86**(B4), 2825–2852.

Ekström, G. & Busby, R.W., 2008. Measurements of seismometer orientation at USArray transportable array and backbone stations, *Seismol. Res. Lett.*, **79**(4), 554–561.

Ekström, G., Dalton, C.A. & Nettles, M., 2006. Long-period instrument gain at global seismicity stations, *Seismol. Res. Lett.*, **77**(1), 12–22.

Ekström, G., Nettles, M. & Dziewoński, A.M., 2012. The global CMT project 2004-2010: centroid-moment tensors for 13,017 earthquakes, *Phys. Earth planet. Inter.*, **200–201**, 1–9.

Fichtner, A. et al., 2018. The collaborative seismic Earth model: generation 1, *Geophys. Res. Lett.*, **45**(9), 4007–4016.

Foster, A., Ekström, G. & Hjörleifsdóttir, V., 2014. Arrival-angle anomalies across the USArray transportable array, *Earth planet. Sci. Lett.*, **402**(C), 58–68.

French, S.W. & Romanowicz, B.A., 2014. Whole-mantle radially anisotropic shear velocity structure from spectral-element waveform tomography, *Geophys. J. Int.*, **199**(3), 1303–1327.

Frohlich, C. & Davis, S.D., 1999. How well constrained are well-constrained T, B, and P axes in moment tensor catalogs? *J geophys. Res.*, **104**(B3), 4901–4910.

Hansen, P.C., 2000. The L-curve and its use in the numerical treatment of inverse problems, in *Computational Inverse Problems in Electrocardiology*, ed. Johnston, P., pp. 119–142, WIT Press.

Hjörleifsdóttir, V. & Ekström, G., 2010. Effects of three-dimensional Earth structure on CMT earthquake parameters, *Phys. Earth planet. Inter.*, **179**(3-4), 178–190.

Hosseini, K., Matthews, K.J., Sigloch, K., Shephard, G.E., Domeier, M. & Tsekhmistrenko, M., 2018. SubMachine: web-based tools for exploring seismic tomography and other models of Earth's deep interior, *Geochem. Geophys. Geosyst.*, **19**(5), 1464–1483.

Hosseini, K., Sigloch, K., Tsekhmistrenko, M., Zaheri, A., Nissen-Meyer, T. & Igel, H., 2020. Global mantle structure from multifrequency tomography using P, PP and P-diffracted waves, *Geophys. J. Int.*, **220**(1), 96–141.

Kanamori, H. & Given, J.W., 1982. Analysis of long-period seismic waves excited by the May 18, 1980, eruption of Mount St. Helens—a terrestrial monopole? *J. geophys. Res.*, **87**(B7), 5422–5432.

Kanamori, H. & Rivera, L., 2008. Source inversion of W phase: speeding up seismic tsunami warning, *Geophys. J. Int.*, **175**, 222–238.

Kim, Y., Liu, Q. & Tromp, J., 2011. Adjoint centroid-moment tensor inversions, *Geophys. J. Int.*, **186**(1), 264–278.

Komatitsch, D., 2011. Fluid-solid coupling on a cluster of GPU graphics cards for seismic wave propagation, *Comptes Rendus - Mecanique*, **339**(2-3), 125–135.

Komatitsch, D. & Tromp, J., 1999. Introduction to the spectral element method for three-dimensional seismic wave propagation, *Geophys. J. Int.*, **139**(3), 806–822.

Komatitsch, D. & Tromp, J., 2002a. Spectral-element simulations of global seismic wave propagation – I. Validation, *Geophys. J. Int.*, **149**, 390–412.

Komatitsch, D. & Tromp, J., 2002b. Spectral-element simulations of global seismic wave propagation – II. Three-dimensional models, oceans, rotation and self-gravitation, *Geophys. J. Int.*, **150**, 303–318.

Komatitsch, D., Michéa, D. & Erlebacher, G., 2009. Porting a high-order finite-element earthquake modeling application to NVIDIA graphics cards using CUDA, *J Parallel Distr Com.*, **69**(5), 451–460.

Komatitsch, D., Erlebacher, G., Göddeke, D. & Michéa, D., 2010. High-order finite-element seismic wave propagation modeling with MPI on a large GPU cluster, *J Comput. Phys.*, **229**(20), 7692–7714.

Krischer, L. & Casarotti, E., 2015. *pyflex: 0.1.4*.

Krischer, L., Megies, T., Barsch, R., Beyreuther, M., Lecocq, T., Caudron, C. & Wassermann, J., 2015. ObsPy: a bridge for seismology into the scientific python ecosystem, *Comput. Sci. Discov.*, **8**(1), doi:10.1088/1749-4699/8/1/014003.

Krischer, L. et al., 2016. An Adaptable Seismic Data Format, *Geophys. J. Int.*, **207**(2), 1003–1011.

Lei, W. et al., 2020. Global adjoint tomography—model GLAD-M25, *Geophys. J. Int.*, **223**(1), 1–21.

Leng, K., Nissen-Meyer, T., van Driel, M., Hosseini, K. & Al-Attar, D., 2019. AxiSEM3D: broad-band seismic wavefields in 3-D global earth models with undulating discontinuities, *Geophys. J. Int.*, **217**(3), 2125–2146.

Li, J., Zheng, Y., Thomsen, L., Lapen, T.J. & Fang, X., 2018. Deep earthquakes in subducting slabs hosted in highly anisotropic rock fabric, *Nat. Geosci.*, **11**(9), 696–700.

Li, X.-D. & Romanowicz, B., 1996. Global mantle shear velocity model developed using nonlinear asymptotic coupling theory, *J geophys. Res.*, **101**(B10), 22245–22272.

Liu, Q., Plet, J., Komatitsch, D. & Tromp, J., 2004. Spectral-element moment tensor inversions for earthquakes in southern California, *Bull. seism. Soc. Am.*, **94**(5), 1748–1761.

Maggi, A., Tape, C., Chen, M., Chao, D. & Tromp, J., 2009. An automated time-window selection algorithm for seismic tomography, *Geophys. J. Int.*, **178**(1), 257–281.

Menke, W. & Russell, J.B., 2020. Non-double-couple components of the moment tensor in a transversely isotropic medium, *Bull. seism. Soc. Am.*, **110**(3), 1125–1133.

Morales-Yáñez, C., Duputel, Z. & Rivera, L., 2020. Impact of 3-D Earth structure on W-phase CMT parameters, *Geophys. J. Int.*, **223**(2), 1432–1445.

Moulik, P. & Ekström, G., 2014. An anisotropic shear velocity model of the Earth's mantle using normal modes, body waves, surface waves and long-period waveforms, *Geophys. J. Int.*, **199**(3), 1713–1738.
Nocedal, J. & Wright, S.J., 2006. *Numerical Optimization*, Springer Science & Business Media.

Peter, D. et al., 2011. Forward and adjoint simulations of seismic wave propagation on fully unstructured hexahedral meshes, *Geophys. J. Int.*, **186**(2), 721–739.

Romanowicz, B., 2008. Using seismic waves to image Earth's structure, *Nature*, **451**, 266–268.

Ruan, Y., Lei, W., Modrak, R., Örsvuran, R., Bozdağ, E. & Tromp, J., 2019. Balancing unevenly distributed data in seismic tomography: a global adjoint tomography example, *Geophys. J. Int.*, **219**(2), 1225–1236.

Sawade, L., Beller, S., Lei, W. & Tromp, J., 2021a. CMT3D - Catalog. DataSpace. https://doi.org/10.34770/gp8e-sx34.

Sawade, L., Beller, S., Lei, W. & Tromp, J., 2021b. GCMT3D - Software. DataSpace. https://doi.org/10.34770/yctp-3c03.

Tape, W. & Tape, C., 2012. A geometric setting for moment tensors, *Geophys. J. Int.*, **190**(1), 476–498.

Trabant, C., Hutko, A.R., Bahavar, M., Karstens, R., Ahern, T. & Aster, R., 2012. Data Products at the IRIS DMC: stepping stones for research and other applications, *Seismol. Res. Lett.*, **83**(5), 846–854.

Tromp, J., Tape, C. & Liu, Q., 2004. Seismic tomography, adjoint methods, time reversal and banana-doughnut kernels, *Geophys. J. Int.*, **160**(1), 195–216.

Tromp, J. et al., 2010. Near real-time simulations of global CMT earthquakes, *Geophys. J. Int.*, **183**(1), 381–389.

Valentine, A.P. & Woodhouse, J.H., 2010. Reducing errors in seismic tomography: combined inversion for sources and structure, *Geophys. J. Int.*, **180**(2), 847–857.

Vavryčuk, V., Adamová, P., Doubravová, J. & Horálek, J., 2021. Moment tensor catalogue of microearthquakes in West Bohemia from 2008 to 2018, *Earth Syst. Sci. Data Discuss.*, **2021**, 1–24.

Wang, X. & Zhan, Z., 2020. Seismotectonics and fault geometries of the 2019 Ridgecrest sequence: insight from aftershock moment tensor catalog using 3-D Green's functions, *J. geophys. Res.*, **125**(5), 1–15.

Wells, D.L. & Coppersmith Kevin, J., 1994. New empirical relationship between magnitude, rupture length, rupture width, rupture area, and surface displacement, *Bull. seism. Soc. Am.*, **84**(4), 974–1002.

Woodhouse, J.H., 1981. The excitation of long period seismic waves by a source spanning a structural discontinuity, *Geophys. Res. Lett.*, **8**(11), 1129–1131.

Woodhouse, J.H. & Dziewoński, A.M., 1984. Mapping the upper mantle: three-dimensional modeling of Earth structure by inversion of seismic waveforms, *J. geophys. Res.*, **89**(B7), 5953–5986.

## SUPPORTING INFORMATION

Supplementary data are available at GJI online.

**Figure S1.** Distribution of DC versus non-DC components as a function of depth. GCMT denotes the distribution prior to inversion and CMT3D denotes the distribution post-inversion. Two sets of binned statistics are shown in the figure: The bar with dot is showing the bin width as well as the mean source type in the bin; the left and right bars each denote one standard deviation; in grey, the source type statistics prior to inversion and in black the source type statistics after inversion.

**Figure S2.** Histogram of the numbers of iterations required for convergence to a given depth. Note that the scale of the y-axis is logarithmic.

**Figure S3.** Single Iteration waveform inversion for moment tensor and hypocentre. Statistical overview of the catalogue. (a) Distribution of inverted earthquakes across the globe. (b) Change in source-type distribution. (c) Depth as a function of change in centroid depth $\Delta z$. Depth binned statistics are shown for this plot. The mean and bin width are indicated by the bar with dot, and the left and right bars denote the standard deviation to that. (d) Distribution of changes in centroid time $\Delta t_c$. (e) Distribution of changes in scalar moment $100 \cdot \Delta(\ln M_0)$. (f) Distribution of changes in centroid depth.

**Figure S4.** Iterative inversion for depth and centroid time. Statistical overview of the catalogue of the subsequent iterative inversion for depth and centroid time. (a) Distribution of inverted earthquakes across the globe. (b) Change in source-type distribution. (c) Depth as a function of change in centroid depth $\Delta z$. Depth binned statistics are shown for this plot. The mean and bin width are indicated by the bar with dot, and the left and right bars denote the standard deviation to that. (d) Distribution of changes in centroid time $\Delta t_c$. (e) Distribution of changes in scalar moment $100 \cdot \Delta(\ln M_0)$. (f) Distribution of changes in centroid depth.

**Figure S5.** Energy correction before and after. Waveform agreement between seismic waveforms of various wave types (body waves, surface waves, and long-period 'mantle' waves) observed on different components (vertical, radial and transverse from top to bottom row) and their corresponding synthetics computed by spectral-element simulation of wave propagation using the waveform inverted CMT3D+ earthquake mechanisms in the 3-D earth model GLAD-M25 (Lei et al. 2020) are compared with their corresponding corrected counterpart CMT3D+X earthquake mechanism. The misfit is measured as normalized L2-norm misfit. Labels on the left of the histograms indicate the wave type for as mentioned in Section 3 for each row and the labels at on the horizontal axes indicate the component. The top-right label shows both mean, and standard deviation, for each histogram.

**Figure S6.** Final Catalogue. Waveform agreement between seismic waveforms of various wave types (body-wave, surface-wave and long-period 'mantle'-wave) observed on different components (vertical, radial and transverse from top row to bottom row) and their corresponding synthetics computed by spectral-element simulation of wave propagation using GCMT earthquake mechanisms (Ekström et al. 2012) in the 3-D earth model GLAD-M25 (Lei et al. 2020) are compared with their corresponding post-inversion counterpart CMT3D+X. The misfit is measured as normalized L2-norm misfit. Labels on the left of the histograms indicate the wave type for as mentioned in Section 3 for each row and the labels at on the horizontal axes indicate the component. The top-right label shows both mean, and standard deviation, for each histogram.

**Figure S7.** Changes in the sources' scalar moment as a function of depth. We collected the scalar moment changes for nine different depth ranges and plotted their 5×5 degree binned mean value. The title of each subplot indicates the depth range. Each map is accompanied by a histogram of changes that is plotted underneath the map. The histogram's units are in percentages of increase or decrease, and in the top left-hand corner the total number of events in the depth range is shown. A vertical line at 0 in the histogram serves as a guide. While the x-axis limits of the accompanying histograms change for each subplot, the colour bar is the same for all subplots.

**Figure S8.** Changes in the sources' centroid time as a function of depth. We collected the changes in centroid time for nine different depth ranges and plotted their 5×5 degree binned mean value. The title of each subplot indicates the depth range. Each map is accompanied by a histogram of changes that is plotted underneath the map. The histogram's units are in seconds, and in the top left corner the total number of events in the depth range is shown. A vertical line at 0 in the histogram serves as a guide. While the x-axis limits of the accompanying histograms change for each subplot, the colour bar is the same for all subplots.

**Figure S9.** Changes in the source depth as a function of depth. We collected the depth changes for nine different depth ranges and plotted their 5×5 degree binned mean value. The title of each subplot indicates the depth range. Each map is accompanied by a histogram of changes that is plotted underneath the map. The histogram's units
are km, and in the top left corner the total number of evens in the depth range is shown. A vertical line at 0 in the histogram serves as a guide. While the x-axis limits of the accompanying histograms change for each subplot, the colour bar is the same for all subplots.

**Figure S10.** Values of the horizontal relocation as a function of depth. We collected the horizontal relocation for nine different depth ranges and plotted their 5×5 degree binned mean value using the vector mean. The title of each subplot indicates the depth range. Each map is accompanied by a histogram of magnitudes of the relocation underneath the map. The histogram's units are kilometres, and in the top left corner the total number of events in the depth range is shown. A vertical line at 0 in the histogram serves as a guide. While the x-axis limits of the accompanying histograms change for each subplot, the colour bar is the same for all subplots.

**Figure S11.** Changes in the source type in terms of the angles $\omega$ (Frohlich & Davis 1999; Hjörleifsdóttir & Ekström 2010; Tape & Tape 2012) and $\gamma$ (Tape & Tape 2012). (a) Histogram of the change in the moment tensor in terms of the angle $\omega$ between the GCMT tensor and the CMT3D+X tensor; x-axis in degrees and y-axis in counts. The numbers in the upper part of the plot indicate the number of events in each segment divided by the dashed lines. (b) Histograms of angle $\gamma$ before (GCMT) and after inversion (CMT3D+X). While the main text (Fig. 4) indicates the $\epsilon$ parameter due to prominence in literature and usage in the GCMT project, the $\gamma$ angle is a more natural choice. Both indicate the change in zero-trace source type.

**Figure S12.** Changes in the source type as a function of depth. We collected the angle $\omega$ (see Fig. S.11) for nine different depth ranges and plotted their 5×5 degree binned mean value using the vector mean. The title of each subplot indicates the depth range. Each map is accompanied by a histogram of angles that is plotted underneath the map. The x-axis limits of the accompanying histograms are the same for all subplots.

**Figure S13.** Spatially binned correlation values between the changes in depth and time as a function of depth. We collected the changes in depth and centroid time for nine different depth ranges in 5×5 degree bins, and we used the binned values to compute the correlation coefficient. Each bin contains at least three events. The title of each subplot indicates the depth range. Each map is accompanied by a histogram of correlation coefficients that is plotted underneath the map. The x-axis limits of the accompanying histograms are the same for all subplots.

**Figure S14.** Spatially binned correlation values between the changes in depth and the scalar moment as a function of depth. We collected the changes in depth and centroid time for nine different depth ranges in 5×5 degree bins, and we used the binned values to compute the correlation coefficient. Each bin contains at least three events. The title of each subplot indicates the depth range. Each map is accompanied by a histogram of correlation coefficients that is plotted underneath the map. The x-axis limits of the accompanying histograms are the same for all subplots.

**Figure S15.** Weighting and filtering setup depending on magnitude and wave type. (a) Weight contribution depending for each wave type. (b) Filtering in period using SAC cosine tapers. Each magnitude and wave type has a unique cosine taper defined by four corners.

**Figure S16.** Typical event-station plot. The inversion results for an event at a single station. The map in the top-left of the figure shows station (upside-down triangle)—event (star) geometry, and an arrow indicates the relocation with a legend in the very top-left corner of

the figure. Below the map we show before (black) and after (blue) moment tensor with a table below indicating the values for the initial moment tensor in black and the values for new moment tensor; values indicating changes with respect to the initial moment tensor are given the d suffix. Mw is the moment magnitude, S, D and R are strike, dip, and rake, respectively, with the two possible solutions. T, B and P are tension, null, and pressure axes and their plunge and azimuth as indicated by the column identifier. EPS is epsilon, which indicates how double-couple/CLVD the moment tensor is. T, LAT, LON and Z are the time-shift with respect to preliminary determined epicentre, latitude, longitude and depth, respectively. Note that the location parameters are given in terms of the change with respect to the initial solution. The first line of the head indicates network (IU) and station (BFO) as well as geographical location. The second line of the header indicates event ID (C051200C) origin time, and location. The third line of the header shows station-event geometry. EPI is the epicentral distance, AZI is the azimuth from event to station and BAZ is backazimuth (station to event). In the top-right of the figure we have a legend indicating which line in the plots is which. Below the legend, we channel information for the respective traces to the left. Here, BHZ-00 indicates channel and location, W indicates the weight used in the inversion for this particular trace. The value below indicates the maximum displacement used for normalization of the traces. M, C and S, are the normalized least-squares misfit, maximum normalized cross-correlation and time-corrected correlation values, respectively, between the observed and GCMT (red) and CMT3D (blue) synthetics. This is repeated for all three channels (Z, R and T). The traces are plotted in bold in their windows and with some transparency outside of the windows. The windows are green if the waveform misfit has improved, and red if the misfit increased. Below the windows we plot measurements between observed and synthetics for GCMT (red) and CMT3D (blue). The values in bold indicate the preferred value.

**Figure S17.** Selection of 70 randomly chosen events and their respective beachball representations, before (GCMT, red) and after (CMT3D+X, blue) inversion. Each subplot contains the Event ID and location in the label above the beachballs and the corresponding changes below the beachballs. '<' denotes the angle $\omega$ (see Fig. S11). d1nM0 is a fractional value.

**Table S1.** Table showing the processing parameters for body waves.

**Table S2.** Table showing the processing parameters for intermediate surface waves.

**Table S3.** Table showing the processing parameters for the mantle waves.

**Table S4.** Windowing parameters for body waves prior to surface waves. Selection mode is set to body waves. For detailed explanations for each model parameter, see Maggi *et al.* (2009).

**Table S5.** Windowing parameters for body waves after to surface waves. Selection mode is set to 'phase_list:PcP,ScS,PcS,ScP'. For detailed explanations for each model parameter, refers to Maggi *et al.* (2009).

**Table S6.** Windowing parameters for intermediate-surface waves. Selection mode is set to 'surface waves'.

**Table S7.** Windowing parameters for mantle waves. Selection mode is set to 'surface waves'.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.