Moving from 1-D to 3-D velocity model: automated waveform-based earthquake moment tensor inversion in the Los Angeles region

SUMMARY

Earthquake focal mechanisms put primary control on the distribution of ground motion, and also bear on the stress state of the crust. Most routine focal mechanism catalogues still use 1-D velocity models in inversions, which may introduce large uncertainties in regions with strong lateral velocity heterogeneities. In this study, we develop an automated waveform-based inversion approach to determine the moment tensors of small-to-medium-sized earthquakes using 3-D velocity models. We apply our approach in the Los Angeles region to produce a new moment tensor catalogue with a completeness of M_L ≥ 3.5. The inversions using the Southern California Earthquake Center Community Velocity Model (3D CVM-S4.26) significantly reduces the moment tensor uncertainties, mainly owing to the accuracy of the 3-D velocity model in predicting both the phases and the amplitudes of the observed seismograms. By comparing the full moment tensor solutions obtained using 1-D and 3-D velocity models, we show that the percentages of non-double-couple components decrease dramatically with the usage of 3-D velocity model, suggesting that large fractions of non-double-couple components from 1-D inversions are artifacts caused by unmodelled 3-D velocity structures. The new catalogue also features more accurate focal depths and moment magnitudes. Our highly accurate, efficient and automatic inversion approach can be expanded in other regions, and can be easily implemented in near real-time system.

North America, Waveform inversion, Earthquake source observations, Wave propagation

1 INTRODUCTION

A more accurate and comprehensive earthquake moment tensor (or focal mechanism) catalogue is important to constraining the regional stress field, understanding tectonic processes and potentially mitigating seismic hazard. Focal mechanisms are routinely computed for large earthquakes (e.g. M > 5.0) using global seismic data (e.g. the U.S. Geological Survey [USGS], the Global Centroid Moment Tensor [GCMT]) and for smaller earthquakes ( M > 4.0) using regional seismic data (e.g. the Southern California Earthquake Data Center [SCEDC], the Japan Meteorological Agency [JMA]). Usually, these focal mechanisms are determined by inverting either the P-wave first-motion polarities, S/P amplitude ratios, or seismic body and/or surface waveforms, and all of them assume simple 1-D Earth velocity models (Dziewonski et al.1981; Kawakatsu 1995; Pasyanos et al.1996; Zhu & Helmberger 1996; Hardebeck & Shearer 2002; Clinton et al.2006; Ekstrom et al.2012; Yang et al.2012). However, recent tomographic studies demonstrate the existence of strong 3-D velocity heterogeneities in the Earth, ranging from global to regional and local scales (Romanowicz 2003; Tape et al.2009). Thus, the inversions based on simplified 1-D velocity models may lead to biases and large uncertainties in focal mechanism solutions, as the complicated Earth 3-D structure effects cannot be adequately quantified.

Numerous efforts have been proposed to account for the effects of 3-D velocity structure on waveform-based moment tensor inversion, including using relative long-period waveforms (Duputel et al.2012; Ekstrom et al.2012), breaking seismograms into segments and allowing time-shifts during fitting the observations with synthetics (Zhao & Helmberger 1994; Zhu & Helmberger 1996) and applying path-dependent corrections (Tan & Helmberger 2007; Ekstrom et al.2012; Chu et al.2014; Wang et al.2017). However, in regions with highly heterogeneous Earth structures, the aforementioned approaches are not adequate. The 3-D structure effects may become particularly ineluctable in moment tensor inversion for small-to-medium-sized earthquakes, as these earthquakes usually only have high signal-to-noise ratio (SNR) waveforms in relatively short periods (Tan & Helmberger 2007). Alternatively, ‘careful’ data selection—by excluding the ‘complicated’ waveforms—have been commonly used to achieve stable inversion results. However, the data selection procedures are usually highly variable and subjective. In addition, a simple exclusion of ‘complicated’ waveforms will inevitable lower the azimuth coverage and increase the uncertainty. This effect is particularly strong in regions with complex velocity structures, such as the Los Angeles region, where a large portion of seismic stations are located within or near the basin (Fig. 1).

Figure 1.

Seismicity and velocity structure in the study region. (a) Purple dots are seismicity from 1981 to 2018 with magnitude above 2.0 from Hauksson et al. (2012). The triangles are the permanent broadband seismic stations in this area, a large portion of which are located within or near the Los Angeles basin. (b) The latest version of the Southern California Earthquake Center (SCEC) Community Velocity Model (CVM)—CVM-S4.26 (Lee et al.2014b)—along the profile in (a). The standard 1-D Southern California velocity model (SoCal) (Hadley & Kanamori 1977) (red) and the 1-D velocity models (black) extracted from the 3-D CVM-S4.26 model (Lee et al.2014b) are compared to demonstrate the strong lateral heterogeneities in this area. (c) Waveform comparison among real data, the synthetics generated using the 1-D SoCal and the 3-D velocity models, indicating the strong 3-D structure effects. The location of earthquake and locations of stations are given in orange star and triangles in (a).

With the recent progress in seismic tomographic studies using traveltime, ambient noise and full waveform inversion, high-resolution 3-D Earth models become available in various areas. Meanwhile, with the advances in computation power, 3-D numerical simulation of seismic wave propagation has become routine. Thus, it is natural to consider using 3-D velocity models in source inversion (Graves & Wald 2001; Liu et al.2004; Zhao et al.2006; Chao et al.2011; Lee et al.2011; Zhu & Zhou 2016; Hejrani et al.2017). Incorporating realistic 3-D velocity models in source inversion (1) can better separate the source signals from the complex 3-D propagation effects; (2) enable us to use more seismic records for better azimuthal coverage and (3) allow us to extend the waveform inversion to higher frequencies to obtain smaller earthquakes’ focal mechanism.

In this study, we discuss the importance and feasibility of using 3-D velocity model in automated earthquake moment tensor inversion in the Los Angeles region, as the Los Angeles region is prone to earthquakes and is located within a tectonically complex region with strong 3-D velocity heterogeneities in the crust (e.g. the deep sedimentary basins, the Moho lateral variations, Lee et al.2014b; Shaw et al.2015,Fig. 1). A series of 3-D community velocity models (CVMs) have been proposed in the region by synthesizing earthquake/noise tomography studies, seismic reflection and refraction surveys, well logs and geologic studies (Süss & Shaw 2003; Tape et al.2009; Lee et al.2014b; Shaw et al.2015; Small et al.2017). Ground motion simulations demonstrate that the 3-D synthetic seismograms computed using the CVMs fit the observed data reasonably well up to 2 s (Komatitsch et al.2004; Taborda & Bielak 2013; Lee et al.2014a). Former studies had attempted to include the CVMs in source inversion (Liu et al.2004; Zhao et al.2006; Lee et al.2011). For example, Liu et al. (2004) used the spectral element method and the Fréchet derivatives to develop a centroid moment tensor (CMT) inversion procedure for southern California earthquakes. However, limited by the computational cost, Liu et al. (2004) demonstrated with case study for three earthquakes, and did not expand to an automated procedure to produce a catalogue in southern California. Zhao et al. (2006) proposed to use strain Green's tensor (SGT) and source–receiver reciprocity approach in moment tensor inversion to reduce the number of 3-D Green's function calculation. Lee et al. (2011) applied this approach and developed an automated procedure to study M_L > 3.0 earthquakes’ CMT solutions in southern California. However, compared with ∼2345 M_L > 3.0 earthquakes in the region, only 165 earthquakes’ CMT solutions have been resolved in Lee et al. (2011)’s catalogue. This incompleteness is likely due to their simplified data selection procedure and the usage of a constant low-pass filter of 0.2 Hz for all the magnitude earthquakes.

To achieve a more accurate and comprehensive earthquake moment tensor catalogue, we develop a waveform-based approach to automatically determine the small-to-moderate-sized (M_L > 3.5) earthquakes’ moment tensor solutions and focal depth, by combing 3-D velocity models, the SGT method (Zhao et al.2006) and the generalized Cut-and-Paste 3-D (gCAP3D) method (Zhu & Zhou 2016) with automatic data quality control. In the following sections, we first describe our automated inversion approach. We apply our approach in the Los Angeles region and discuss the importance of using 3-D velocity model in source inversion. We also compare the full moment tensor solutions obtained using 1-D and 3-D models to investigate the effect of velocity model on non-double-couple components. We then compare our results with the current catalogues and discuss their similarities and discrepancies. By using 3-D waveform simulations, we demonstrate that our approach can refine the existing moment tensor catalogues in southern California.

2 AUTOMATED EARTHQUAKE MOMENT TENSOR INVERSION

In waveform-based moment tensor inversion, we aim to find the best moment tensor solutions (⁠|${M_{ij}}$|⁠) by minimizing the waveform misfits between observations and synthetics. The synthetic seismograms of a point source can be expressed as

$$\begin{equation*} {u_n}{\rm{\ }}\left( {x,t;{x_0}} \right) = {M_{ij}}{\rm{\ \ }}\partial _j^{\rm{^{\prime}}}{G_{ni}}\left( {x,t;{x_0}} \right) \end{equation*}$$

(Aki & Richards 2002), where |${u_n}$| is the n-component displacement field at location x and |${M_{ij}}$| is the second-order moment tensor at source |${x_0}$|⁠. The Green's tensor (G, or Green's function), describing the Earth's impulse response from |${x_0}$| to x, can be computed for a given Earth velocity model. In this paper, we aim to use 3-D Green's functions, which are significantly more expensive to compute than 1-D Green's functions. In the next subsections, we will discuss the Green's function calculation and our inversion approach.

2.1 Strain Green's tensor and reciprocity

Calculation of 3-D Green's functions requires considerable computation depending on the numerical accuracy for the synthetics, the dimension of study area, and the number of stations. Moreover, 3-D Green's function is a non-linear function of the source location. In the forward simulations (Fig. 2a), the number of simulations will be proportional to the number of sources, making it impractical when the number of earthquake locations to be tested is large (e.g. a large number of earthquakes in a broad region). Alternatively, the strain Green's tensor (SGT) and source–receiver reciprocity approach (Fig. 2b) is more efficient in calculating the 3-D Green's function (Eisner & Clayton 2001; Zhao et al.2006). Taking into the symmetry of the moment tensor and the reciprocity property of the Green's tensor, the displacement at x due to an earthquake at x₀ can be written as:

$$\begin{equation*} {u_n}{\rm{\ }}\left( {x,t;{x_0}} \right) = \frac{1}{2}\left[ {\partial _j^{\rm{^{\prime}}}{G_{in}}\left( {{x_0},t;x} \right) + \partial _i^{\rm{^{\prime}}}{G_{jn}}\left( {{x_0},t;x} \right)} \right]{M_{ij}}, \end{equation*}$$

where the |${H_{jin}}\ ( {{x_0},t;x} ) = \frac{1}{2}\ [ {\partial _j^{\prime}{G_{in}}( {{x_0},t;x} ) + \partial _i^{\prime}{G_{jn}}( {{x_0},t;x} )} ]$| is defined as SGT, which describes the wavefield from the receiver x to the source x₀ (Zhao et al.2006). Thus, the displacement at the receiver can be obtained by computing the wave propagation with the ‘virtual source’ acting at the receiver (Fig. 2b). As a consequence, the number of simulations in 3-D Green's function calculation is proportional to the number of receivers, regardless of the number of potential sources. Thus, the usage of source–receiver reciprocity approach is preferred when the number of source grids (80 001 in this study) is much larger than the number of stations (32 in this study). Especially, the source–receiver reciprocity approach is particularly suitable for (near) real-time moment tensor inversion, as the synthetic seismograms for arbitrary locations/mechanisms of sources within the study area can be quickly simulated (within seconds) with the pre-calculated SGT database. This approach has been discussed in Zhao et al. (2006) and extensively used in source inversion (Zhao et al.2006; Chao et al.2011; Lee et al.2011; Zhu & Zhou 2016).

Figure 2.

Schematic diagram illustrating the (a) forwarding simulation and (b) source–receiver reciprocity approach. The source–receiver reciprocity states that the locations of source and receiver can be switched and the exactly same seismic response will be observed.

2.2 Automated gCAP3D inversion scheme

Once we have the SGTs, the moment tensor can be obtained by minimizing the misfits between observations and synthetics. In this study, we use gCAP3D (Zhu & Zhou 2016) as our main driver for inversion. Following the principles of the CAP methodology (Zhao & Helmberger 1994; Zhu & Helmberger 1996), gCAP3D breaks the seismograms into P and S/Surface waves, and model them simultaneously by allowing different time-shifts between observations and synthetics. The time-shifts can accommodate inaccuracies in the assumed velocity models and earthquake locations. The amplitude ratios between P and S/Surface waves allow us to better constrain focal depths. We refer the readers to Zhu & Zhou (2016) for details of the inversion method. Here we focus on proposing an automated inversion approach that can be applied in near real-time and produce a routine moment tensor catalogue. We revise the gCAP3D with an iterative scheme with automatic data quality control.

In our inversion, we first cut the seismograms into five segments (P waves on vertical and radial components, and S/Surface waves on vertical, radial and tangential components), and we conduct SNR control in the frequency bands of interest. The lengths of segments and frequency bands depend on the magnitude of the earthquake, which will be discussed in the next section. We then conduct an initial inversion using the observations with high SNR (larger than 2.5 in this study) data. The initial inversion is in general quite stable, given the high quality data from the Southern California Seismic Network (SCSN) and the high performance of the SCEC CVM. But the inversion occasionally falls into faulty solutions, especially for small earthquakes, due to either strong noise or ‘complicated’ waveforms. Thus, we use both the waveform cross-correlation coefficient (CC) and the relative segment misfit (RSM) to quantify the waveform fitting for each segment (Fig. 3b). The CC is defined as:

$$\begin{equation*} CC\ = \frac{{\mathop \smallint \nolimits_{{t_1}}^{{t_2}} u\left( t \right)s\left( {t + \tau } \right){\rm d}t}}{{\sqrt {\mathop \smallint \nolimits_{{t_1}}^{{t_2}} {u^2}\left( t \right){\rm d}t\mathop \smallint \nolimits_{{t_1}}^{{t_2}} {s^2}\left( t \right){\rm d}t} }}, \end{equation*}$$

where |$u( t )$| and |$s( t )$| represent observation and synthetics, respectively, |${t_1}$| and |${t_2}$| define the time window and the |$\tau $| refers to the time-shift. The RSM is defined as the ratio of the segment misfit to the average misfit (aveMISFIT):

$$\begin{equation*} \mathit{ RS{M_i}} = \frac{{\mathit{ MISFI{T_i}}}}{{\mathit{ aveMISFIT}}},\,\,\mathit{ MISFI{T_i}} = {w_i} \times {\left\| {u\left( t \right) - s\left( {t - \tau } \right)} \right\|^2}, \end{equation*}$$

where the |${w_i}$| is the weighting to the segment. We use a threshold CC, which starts from a low value (e.g. 30 per cent) and slowly increases (e.g. 5 per cent in each step) with the number of iteration, to reject the low-CC data. In addition, we automatically reduce the weighting factor of the waveform segments that have large RSM (e.g. larger than 3 times of aveMISFIT) (Fig. 3b). Our automatic reweighting algorithm is summarized in Fig. 3 and Fig. S1. This reweighting step is important, because the basin stations always have large amplitude and the 3-D velocity model may not capture the basin effects fully. Lowering the weighting factor of segments with large RSM can avoid the overfitting of basin stations without sacrificing the fits to the other stations. After removing the low-CC segments and reducing the weighting factor of the large-RSM segments, the source inversion is performed again in the next iteration (Fig. 3b). The iteration procedure will be stopped if the all the waveform CC reach the final threshold (e.g. 70 per cent in this study) or after a certain number of iterations (e.g. 8 iterations). We emphasize that the data selection is governed by the weighting factor of waveform segments. We only need to read the data once and change the weighting factor in each iteration to speed up the inversion without redundant data reads.

Figure 3.

The flow chart of the proposed automated moment tensor inversion scheme. (a) The system consists of model setup and inversion. The model setup is conducted once before earthquakes, while the inversion part will run for individual earthquake. The inversion part consists of modules for data downloading, data pre-processing, Green's functions preparation, inversion (shown in b), and results broadcasting. (b) The iterative inversion scheme with automatic data quality control. The data quality control is governed by the weighting factor of waveform segments (SegW). We use both the waveform cross-correlation coefficient (CC) and the relative segment misfit (RSM) to quantify the waveform fitting for each segment. The weighting factor of each segment changes iteratively until the inversions converge. More details about the automated inversion scheme can be found in section 2.

We also invert for the earthquake centroid depths in this study, which are usually less constrained than the epicentre. The best depth is determined by repeating the inversion at different focal depths, each of which is independent with the automatic data quality control (Fig. 3b). In principle, the best focal depth and the best moment tensor solutions will have the lowest misfit and can fit the most observations (numSeg, the number of segments that have CC larger than the threshold of CC). Thus, we select the best focal depth by minimizing the averaged misfit and maximizing the number of segments that can be fitted:

$$\begin{equation*} \ dept{h_{best}} = \ min\left\{ {\left[ {\frac{\mathit{ {aveMISFIT}}}{{min\left( \mathit{ {aveMISFIT}} \right)}} - 1} \right]\! \times \!\ 100 - numSeg} \right\}. \end{equation*}$$

3 APPLICATION IN THE LOS ANGELES REGION

We apply our automated inversion approach in the Los Angeles region (Fig. 4a). We use a 4th order staggered-grid 3-D finite difference method (FD, Graves 1996) and the SCEC CVM-S4.26a (Lee et al.2014b) to calculate the SGT database. Here, the topography, bathymetry and the oceans are not included in our FD simulation, as the current 3-D FD code (Graves 1996) cannot handle the effects of surface topography and water accurately. Previous studies showed that the effect of topography and water could be significant when topography is high and water is deep (Komatitsch et al.2004; Zhou et al.2016). Therefore, the full complexity of the 3-D model (topography, bathymetry and the oceans) will need be included in the future studies. Our FD simulations are carried out in a volume of |$120\ {\rm {km}}\ \times \ 180\ {\rm {km}}\ \times \ 45\ {\rm {km}}$|⁠, with the implementation of the absorbing boundary condition (Clayton & Engquist 1977) to limit the potential broader effects. With a grid spacing of 100 m (Fig. 4a) and a sampling interval of 0.005 s, our FD simulations result in numerical accuracy up to 0.5 Hz. We then save the SGT database on a mesh of |$2\ {\rm {km}}\ \times \ 2\ {\rm {km}}\ \times \ 1\ {\rm {km}}$| grids (down to 20 km depth by considering the seismogenic depth) as potential earthquake locations (Fig. 4b). Using a modern computer cluster at the High Performance Computing Center at California Institute of Technology, it takes about 17 hr to finish a single station SGT calculation by using 480 CPUs, resulting in 274 GB per station storage on a hard disk. We select 32 high-quality broadband SCSN stations evenly distributed in the study region (Fig. 4a) to construct the whole SGT database, resulting ∼7.8 TB of SGT storage. Table 1 summarizes the parameters related to the FD simulations. Generating the whole SGT database on the dense grid is still a computationally expensive task, but the computational cost is one-time. Once we have the SGT database, the SGT at any source location within the study area can be extracted either from an exact grid match or using an interpolation method. In this study, we search for the nearby 72 grids that bracket the source location to conduct 3-D cubic-spline interpolation (Fig. 4b). In general, it takes about 7 s to obtain a single station SGT. Fig. 4(c) shows the comparison of a SGT calculated from an exact forward simulation with that calculated by interpolating the pre-established SGT database. The two calculations yield exactly the same result, confirming the accuracy and efficiency of the 3-D SGT approach.

Figure 4.

Model setup in this study. (a) Purple dots are earthquakes (M_L> = 3.5) investigated in this study, and the yellow triangles are seismic stations used in inversions. Our Strain Green's Tensor (SGT) simulations are carried out in a volume of |$120{\rm{\ }}km{\rm{\ }} \times {\rm{\ }}180{\rm{\ km\ }} \times \ 45{\rm{\ }}km$|⁠. The grid spacing is 100 m. (b) We mesh the study region into |$2{\rm{\ }}{\rm {km}}{\rm{\ }} \times {\rm{\ }}2{\rm{\ }}{\rm {km}}{\rm{\ }} \times {\rm{\ }}1{\rm{\ }}{\rm {km}}$| grids (down to 20 km depth by considering the seismogenic depth) to save the SGT database. Thus, the SGT at any source location within the study area can be extracted in merely seconds using an interpolation method. We search for the nearby 72 grids that bracket the source location to conduct 3-D cubic-spline interpolation. We use less grids in vertical direction as the vertical direction is highly sampled (1 km spacing) and has less grids (20 grids in total). (c) Comparison of the SGT calculated from a forward simulation (black line) with that calculated by interpolating the grid-based SGT database (red line) confirms the accuracy of the 3-D Green's functions used in our inversions.

Table 1.

Model parameters used in this study.

Model setup
Model	CVM-S4.2.6.M01	SCEC-CVM model
VPmin	1.8	Minimum Vp used in calculation (km s–1)
VSmin	0.25	Minimum Vs used in calculation (km s–1)
DENmin	2.2	Minimum density used in calculation (g cm–3)
Model_Lat	34.0	Latitude of model origin
Model_Lon	−117.9	Longitude of model origin
Model_Rot	0	Rotation of Cartesian model coordinates
Model_Xlen	180	x-axis (positive east) dimension of model space (km)
Model_Ylen	120	y-axis (positive south) dimension of model space (km)
Model_Zlen	45	z-axis (positive down) dimension of model space (km)
Model_H	0.1 km	Spatial grid size (km)
Qsfrac	100	Scaling factor for Qs, at each grid point Qs = Vs*Qsfrac
Qpfrac	50	Scaling factor for Qp, at each grid point Qp = Vp*Qpfrac
Nt	40 000	Number of time steps to compute
Dt	0.005	Time sampling interval (s)
Strain Green's Tensor's Storage setup
SGT_X	2	SGT grid size in x-axis (km)
SGT_Y	2	SGT grid size in y-axis (km)
SGT_Z	1	SGT grid size in z-axis (km)
Simulation Summary
N_CPUs	480	Number of CPUs
Simulation Cost	9.5 hr	per station for three-point body forces
Storage Cost	72 G	per station for three-point body forces

We focus on earthquakes with local magnitude larger than 3.5 in the SCEDC catalogue in this area from 2000 to 2018, resulting in 78 earthquakes to be studied (Fig. 4a). For each earthquake, we use the Auto-gCAP3D approach discussed in the previous section to obtain the source parameters (Fig. 3). We retrieve the three-component seismic waveform data from the SCEDC Seismogram Transfer Program (STP). The source information, including the origin time, location, and magnitude, is adopted from the double-difference relocated catalogue by Hauksson et al. (2012). The P and S arrivals are obtained from the SCEDC STP when available (Hutton et al.2010), or we estimate the theoretical arrival times using the hypocentre and station locations with the 1-D velocity model of SoCal (Hadley & Kanamori 1977). We then conduct a series of data pre-processing, including removing instrument response to ground velocity, resampling, marking the P/S arrivals, and rotation (Fig. 3a). Meantime, the corresponding SGT at each station is extracted from the pre-calculated SGT database using the 3-D cubic-spline interpolation method (Fig. 4c). In our inversion, parameters such as frequency bands and time window of P and S/Surface waves are chosen based on the magnitude of the target earthquake (Fig. 5). These values are determined empirically, but they have proven to be robust and appropriate for earthquakes in the Los Angeles region. The earthquake source duration used in inversion is based on the scalar moment of the event (Ekstrom et al.2012), as our results are less sensitive to source durations due to the usage of relative long-period waveforms. We conduct grid search (5° searching interval) to find optimal focal parameters (strike, dip and rake). Typically, the solutions converge after eight iterations. We also run our inversion using different focal depths. The best focal depth is obtained by minimizing the averaged misfit and maximizing the number of segments. In the last step, we fix the focal depth at preferred values, and use a finer grid (2° searching interval) to output the final moment tensor solutions and the associated uncertainties, which are estimated using a bootstrapping method (Efron & Tibshirani 1991; Zhan et al.2012). Fig. 6 shows an example inversion result of the 2018/08/29 02:23 La Verne M4.4 earthquake, which occurred near the northern edge of the Los Angeles basin with most of the azimuth coverage is provided by stations in the basins to the southwest. By using 3-D Green's functions, we are able to find a well-constrained focal mechanism solution that fits most of the waveforms well, with averaged CC higher than 85 per cent. The bootstrapping results show that the uncertainties of strike, dip, and rake are all within 3°. The preferred centroid depth of 5 km not only provides the best RSM data fitting and waveform correlations, but also allow us to include the maximum number of waveforms in inversion. For all the earthquakes studied, the obtained moment tensor solutions can be obtained from the Supporting Information.

Figure 5.

Empirical inversion parameters used in this study. The band-pass filter used in inversion is described by the high-pass (f1) and low-pass (f2) cut-offs. We use different frequency bands for different magnitude earthquakes because smaller earthquakes usually only have high signal-to-noise ratio waveform in short periods. The time window is defined as −5/−10 s before and Tp/Ts s after the on-set of P/S arrivals. The length of time window varies with magnitude to accompany the different frequency bands of the waveforms used in inversion. These values are proposed based on detailed case studies for varies magnitude earthquakes in the Los Angeles region, but they have proven to be robust and appropriate for our automated inversion in this region.

Figure 6.

Source parameters inversion results for the 2018/08/29 02:33 La Verne M4.4 earthquake. (a) The earthquake location, seismic stations and inversion parameters used in this study. The stations are colored by the averaged waveform cross-correlation coefficients (CC), which serves as proxies for the quality of inversion. The upper right beachball shows the final fault plane solutions, with the uncertainties (95 per cent confidence level) estimated by a bootstrapping method. (b) The averaged CC, the number of segments with CC larger than the threshold of CC, and the waveform misfit as a function of focal depth and focal mechanism, indicating a well-resolved focal depth of 5 km. (c) Waveform fitting at different stations. The black and red lines indicate data and synthetics, respectively. The grey and blue lines represent the waveforms discarded in our automatic data selection. The station names are indicated on the left of waveform pairs along with the distance in km and azimuth in degree. The waveform CC (lower) and time-shits (upper) are shown around the waveform pairs.

4 EFFECT OF VELOCITY MODELS IN INVERSION

To understand and demonstrate the effect of velocity models on moment tensor inversion, we compare the inversion results obtained using 1-D and 3-D velocity models.

We use a moderate earthquake (the 2014/03/29 21:32 La Habra M_w 4.1) occurred close to the central of study area as an example (Fig. 7a). We invert for the moment tensor solutions and focal depths using either the 1-D Southern California Model (1D-SoCal) (Hadley & Kanamori 1977), or the 1-D regional basin model extracted from the Crust 1.0 (1D-Basin) (Laske et al.2013), or the 3-D CVM-S4.26a model (Lee et al.2014b). We keep the inversion parameters and automatic data selection criteria same in all three cases for fair comparisons. As expected, using the 3-D velocity model produces the best waveform fits and can include the most seismic observations in inversions (quantified by the number of segments that have CC larger than the threshold of CC; Fig. 7a). The inversion using the 1D-SoCal model yields a strike-slip mechanism, with large uncertainties shown by the bootstrapping results (Fig. 7b). Even though the inversion using the 1D-SoCal model converges to a minimum focal depth of 4.5 km, the number of segments that can be fitted (red line in Fig. 7b) is very small, indicating the current inversion is unstable. The inversion using the 1D-Basin model results a mixture of strike-slip and thrust mechanisms, with smaller scattering in the bootstrapping solutions (Fig. 7b). Centroid depth is estimated to be 8 km, although 5.5 km provides similar waveform fitting with lower CCs and more discarded data (Fig. 7b). The inversion utilizing the 3-D model suggests that this earthquake is almost a pure thrust event. The uncertainties in focal mechanism also decrease dramatically using the 3-D velocity model (see the corresponding histogram distributions for the strike, dip and rake in Fig. S2), owing to the fact that we can better fit almost all the observations (Fig. 8). In addition, the constraint on depths to be ∼6 km is strong with a sharp misfit curve and clear maxima in both average CC and number of waveforms used.

Figure 7.

Comparison of moment tensor inversion results using 1-D and 3-D velocity models. Panel (a) shows the seismic stations included in the inversion for the 2014/03/29 21:32 La Habra M4.1 earthquake, with at least one of the waveform segments satisfying our data selection criteria. The background color shows the basin thickness in this region (Lee et al.2014b). In general, the 1-D models can only fit relatively simple body and/or surface signals for bedrock stations, not basin sites. (b) The averaged cross-correlation coefficient (CC), the number of segments that can be fitted (those waveforms have CC larger than the threshold of CC), and the waveform misfit as a function of focal depth and focal mechanism. The misfit function curve is more convergent and less fluctuant in the case of using the 3-D velocity model. In addition, the usage of the 3-D velocity model produces clear maxima in both average CC and number of waveforms used.

Figure 8.

Comparison of waveform fits using 1-D and 3-D velocity models. The obtained focal mechanism solutions are shown as beachballs, where the black lines are optimal results, and the grey lines are uncertainties estimated by a bootstrapping method (95 per cent confidence level). Histograms of the cross-correlation coefficients of waveform fits (including those discarded by our automatic data selection) are shown on the top. We display all the waveform fits to demonstrate that the improvement of waveforms fits is universal for all the stations.

The improvements in the focal mechanism solutions by using the 3-D velocity model come from the significantly better waveform fits (Fig. 8). In the frequency ranges of interest (0.08–0.35 Hz for P waves and 0.07–0.23 Hz for S/Surface waves in Fig. 8, or in Fig. 5 more generally), the 1D-SoCal model fits well the relatively simple body and/or surface waves for bedrock-site stations (Fig. 7a). The 1D-Basin model produces slightly better waveform fits (can fit more stations and have slightly higher CCs) than the 1D-SoCal model (Fig. 8), as the earthquake and the majority of stations are located within or near the edge of the Los Angeles basin. However, neither 1-D models capture the complex and strongly amplified waveforms within the basin (Fig. 8). In comparison, the usage of the 3-D velocity model produces substantially better fits to observed waveforms than those obtained using 1-D velocity models. We emphasize that even for the stations that are discarded in the source inversion (CC < 70 per cent), the synthetic waveforms predicted using 3-D velocity model show better fits to the real data (Fig. 8). More quantitatively, the histogram of CCs between observations and synthetics obtained by the 3-D inversion are systematically concentrated at higher values than that from the 1-D inversions (Fig. 8), indicating the importance of incorporating 3-D velocity model in source inversion, especially in regions with complicated 3-D structures.

5 RESULTS AND DISCUSSION

We obtained moment tensor solutions for 75 earthquakes out of the initial set of 78 events (Fig. 9). The remaining three are early aftershocks or triggered events with noisy long-period seismic signals due to the coda of the main events. Overall, moment tensor solutions can be well determined in our automated inversion approach for almost all the earthquakes with M_L > 3.5. These earthquakes are dominated by thrust and strike-slip mechanism (Fig. 9), consisting with the current crustal deformation field in the Southern California region (Yang & Hauksson 2011; Rollins et al.2018).

Figure 9.

Comparison of focal mechanism solutions from various catalogues. (a) Red beachballs show the focal mechanisms obtained using our automated inversion approach; black beachballs represent the focal mechanisms from the YHS catalogue, and the blue beachballs are focal mechanisms from the SCEDC-CMT catalogue. The characters (B, C or D) marked the upper right corner of the YHS solutions indicate the quality of inversion, and the unmarked focal mechanism solutions are results with quality A (nodal plane uncertainty equal or less than 25°) (Yang et al.2012). The yellow color highlights the earthquakes with significant discrepancies among various catalogues. (b) Histogram of the focal mechanism rotation angle between our results and the YHS catalogue.

5.1 Compare with the YHS and the SCEDC-CMT catalogue

The moment tensor catalogue in the Southern California region is routinely maintained by the Southern California Earthquake Data Center (http://service.scedc.caltech.edu/eq-catalogs/CMTsearch.php, SCEDC-CMT), by inverting the long-period surface waves (Clinton et al.2006). Alternatively, the SCEDC also release updates to the YHS catalogue (Yang et al.2012), which is computed from high-frequency P-wave polarities and S/P amplitude ratios using the HASH method (Hardebeck & Shearer 2002, 2003). Owing to the usage of high-frequency signals, the YHS catalogue has the capability to resolve focal mechanisms for earthquakes down to M_L 1.0, representing the most comprehensive focal mechanism catalogue in the Southern California region.

We compare the focal mechanism solutions from the routine catalogues with ours to evaluate the quality of our inversions. Fig. 9 summaries the comparison among our Auto-gCAP3D catalogue, the YHS, and the SCEDC-CMT catalogues, for the 75 events we have studied. Generally, the focal mechanism parameters among the three catalogues match each other. To quantify the difference between our solutions and those from the YHS catalogue, we calculate the focal mechanism rotation angle (Kagan 1991), which is defined as the minimum rotation that is needed to transfer from one focal mechanism to the other. Statistics show that the mean rotation angle is around 20°, smaller than the formal uncertainty of 25° for quality A events in the YHS catalogue (Yang et al.2012), indicating good consistency between our catalogue and the YHS catalogue.

However, we found that 6 out of the 75 events have rotation angle larger than 50° (highlighted in yellow in Fig. 9b). For example, our inversion of the 2007/09/02 17:29 Trabuco Canyon M4.7 event suggests a thrust faulting mechanism, while the YHS catalogue characterizes this earthquake as a strike-slip faulting mechanism. Similar differences in the estimated focal mechanism solutions also appeared for several other smaller earthquakes (Fig. 9a). The YHS catalogue classifies the inversion results as quality A/B/C/D based on nodal plane uncertainties (Yang et al.2012). Of the 6 events, most of them have quality-A solutions in the YHS catalogue (Fig. 9a), representing the focal mechanism with a mean nodal plane uncertainty equal to or less than 25°. Therefore, the rotation angles larger than 50° are statistically significant. To understand such differences in solutions without knowing the true focal mechanism, we analyse the waveform fits between observed and synthetic seismograms, which are calculated by fixing the 3-D velocity model and adopting focal mechanism solutions from the different catalogues. As the accuracy of the 3-D velocity model has been evaluated in previous studies (Lee et al.2014a; Taborda et al.2016), the better focal mechanism should provide better waveform fits. As shown in Fig. 10(a), 3-D synthetic waveforms generated using our preferred thrust-faulting mechanism does provide better fits to the observed waveforms than the strike-slip-faulting mechanism. The systematic improvement of waveform fits for all the earthquakes shown in Fig. 10(b) demonstrates that the solutions obtained by our technique are better centroid moment tensor solutions.

Figure 10.

(a) Comparison of waveform fits using different source parameters adopted from different catalogues. The waveforms are displayed in the same manner of Fig. 6. Here, we display all the waveform fits to demonstrate that the improvements of waveform fits are not biased by data selection. (b) Histograms of the waveform-fitting cross-correlation coefficients (CC) obtained by using source parameters from our study (red) and the YHS catalogue (black), indicating a systematic improvement of waveform fits for all the earthquakes. The grey lines in earthquake focal mechanisms indicate the uncertainties estimated by a bootstrapping method (95 per cent confidence level), suggesting that our focal mechanism solutions are well constrained.

The YHS catalogue is computed using high-frequency P-wave polarities and S/P amplitude ratios using the HASH method, by assuming simple 1-D Earth velocity models. The HASH method is easily prone to errors caused by the determination of first-motion polarities (Hardebeck & Shearer 2002). In addition, the extremely complicated Earth 3-D structures, which could modify the estimation of source depths and take-off angles, may also contribute to the errors in focal mechanism solutions (Takemura et al.2016). Finally, the difference in focal mechanism solutions between the YHS catalogue and our results could represent the difference between the initial and total rupture process, as the high-frequency P-wave polarities describe the earthquake initiation and the centroid moment tensor solutions describe the main-rupture process. A more detailed study is required to understand the differences between the first motion and centroid moment tensor solutions, which will be left for further investigations.

5.2 Earthquake focal depth

We also compare the earthquake focal depth obtained by our study with previous studies. The usage of 3-D velocity model and the amplitude ratio between P and S/Surface waves in inversion give us strong constrains on earthquake focal depths (Zhu & Zhou 2016). Comparison of the focal depths determined using our algorithm and those determined by double difference relocation (Hauksson et al.2012), the YHS, and the SCEDC-CMT are shown in Figs 11(a)–(c). Noted that the depths in the YHS catalogue mainly inherit from Hauksson et al. (2012)’s double-difference catalogue. The focal depths obtained by our inversion correlate well with those determined by double-difference relocation, except that our results in general yield slightly shallower depths (∼1–2 km); this is most likely due to the usage of different velocity models in different studies. However, there are large discrepancies between the focal depths obtained by our study and the SCEDC-CMT. The reason mainly comes from the fact that the SCEDC-CMT typically places earthquakes at fix depths, as the SCEDC-CMT inversion is not sensitive to the focal depth (Clinton et al.2006). The comparisons shown in Figs 11(a)–(c) demonstrate that the focal depths obtained by our technique are in good agreement with those determined by double-difference method and are an improvement over those derived by the SCEDC-CMT.

Figure 11.

(a-c) Comparison of earthquake focal depths determined using our method (horizontal axes) and those obtained by double-difference relocation (Hauksson et al.2012), the YHS catalogue (Yang et al.2012), and the SCEDC-CMT catalogue. (d-e) The relationship between M_L and M_w. Dashed grey lines indicate the M_w = M_L, blue lines are M_w–M_L relationship from Clinton et al. (2006), and the red lines are M_w–M_L relationship from linear regression in this study. The M_L appears to systematically overestimate the M_w in the Los Angeles region. (f) Comparison of M_w determined using our method and that obtained by SCEDC-CMT. The M_w determined from our automated inversion generally correlate very well with the M_w determined by SCEDC-CMT. (g-h) The relationship between the reduced-M_L (M_Lr) and M_w. The M_Lr is a simple linear adjustment used in SCEDC to bring the M_L closer to M_w.

5.3 The relationship between M_w and M_L

In this section, we compare the scaling relationship between the moment magnitude (M_w) and local magnitude (M_L) for the 75 earthquakes we have studied. The M_L appears to systematically overestimate the M_w in the Los Angeles region (Fig. 11d), which is consistent with the results inferred from the SCEDC-CMT catalogue (Fig. 11e). We perform linear regression for the M_w versus M_L, which is presented by the relation of |${M_{\rm w}} = \ 0.8560\ \times {M_{\rm L}} + 0.4092$|⁠. Our new M_w–M_L relation is ∼0.3 magnitude unit smaller than the regression relation obtained using the whole southern California earthquakes (Clinton et al.2006,Figs 11d and e), potentially indicating regional variation for the M_w–M_L relationship in southern California.

As there is a systematic discrepancy between the M_w and M_L (Figs 11d and e), the SCEDC introduced the M_L–reduced magnitude (M_Lr) since 2016. The M_Lr is a simple linear adjustment (⁠|${M_{\rm{Lr}}} = \ 0.8530\ \times {M_{\rm L}} + 0.40125$|⁠) to bring the M_L closer to M_w, which is useful to understand the magnitude of earthquakes if the M_w cannot be obtained (http://scedc.caltech.edu/eq-catalogs/change-history.html). The comparison between the M_w and M_Lr shows that the M_w determined from our automated inversion generally agrees with the M_Lr, even though there is slightly larger scatters for the small events (Figs 11g and h). We also compare the M_w determined from our study to that derived from the SCEDC-CMT. Differences between the M_w derived from our study and SCEDC are subtle, yet present, especially for small earthquakes (Fig. 11f). The difference in M_w between two catalogues could be partially explained by the difference in earthquake focal depth, as there is a strong trade-off in moment magnitude and focal depth in long-period waveform inversion.

5.4 The effect of velocity models on non-double-couple components

In previous sections, we limited our moment tensor inversion for pure double-couple (DC) solutions, as expected for shear faulting in shallow tectonic earthquakes. However, a wide variety of processes (e.g. volumetric changes in source area, fault geometry complexity) can cause earthquake mechanisms to have significant non-double-couple components (non-DC) (Julian et al.1998; Ben-Zion & Ampuero 2009; Ross et al.2015; Shi et al.2018). On the other hand, the observed non-DC components in moment tensor solutions can be the artifacts in source inversion due to the unmodeled 3-D or anisotropic velocity structure (Vavryčuk 2004; Šílený 2009; Li et al.2018). In this aspect, it is important to quantify the accuracy and reliability of the non-DC components. In this section, we compare the full moment tensor solutions obtained using 1-D and 3-D models to investigate the effect of velocity model on non-DC components.

We adopt the source decomposition of Zhu & Ben-Zion (2013) and include compensated linear vector dipoles (CLVD) and/or isotropic (ISO) components in our inversion. We use the ratio of scalar moment of each component to quantify the strength of CLVD and ISO components, with the maximum percentage of CLVD and ISO being 25 and 100 per cent, respectively (Zhu & Ben-Zion 2013). We invert for the pure deviatoric moment tensor solutions (DC + CLVD, Fig. 12a) and full moment tensor solutions (DC + CLVD + ISO) (Fig. 12b) using either the 1D-SoCal (Hadley & Kanamori 1977), or the 1D-Basin (Laske et al.2013), or the 3-D CVM-S4.26a model (Lee et al.2014b), by keeping the inversion parameters same in all cases for fair comparisons. Results show that a significant number of earthquakes have strong non-DC components in the inversions using 1-D velocity models (Fig. 12), while the inversions using the 3-D velocity model show that most events have nearly DC mechanisms with limited CLVD and ISO components (Figs 12 and S3). The percentage of non-DC components decrease dramatically with the usage of the 3-D velocity model, suggesting that the large percentage of non-DC components in the 1-D inversions is a result of poorly constrained Earth structure. The comparison of the full moment tensor solutions derived using 1-D and 3-D velocity models reinforce the importance of incorporating 3-D velocity model in inversion to understand the faulting behaviour. We would expect a further decrease in the percentage of non-DC components in source inversion with improvements made to the 3-D velocity model.

Figure 12.

Comparison of (a) deviatoric moment tensor solutions and (b) full moment tensor solutions obtained using different velocity models. The percentage of compensated-linear-vector-dipoles (CLVD) and isotropic (ISO) components is defined using the scalar moment of each component to the total moment, with the maximum percentage of CLVD and ISO being 25 and 100 per cent, respectively (Zhu & Ben-Zion 2013). The percentage of non-double-couple components (CLVD and/or ISO) decrease dramatically with the usage of the 3-D velocity model, suggesting that the large percentage of non-double-couple components in the 1-D inversions mainly comes from the unmodeled 3-D velocity structure. For comparison, we also show the percentage of CLVD components for each earthquake in the SCEDC-CMT catalogue.

Meanwhile, we found relatively large percentage of CLVD and ISO components for several events (deviatoric catalogue and full moment tensor catalogue can be obtained in Supplementary material), corresponding to up to ∼8 per cent and ∼26 per cent of the total moment releases, respectively (Fig. 12b), which are difficult to explain by the error of the 3-D velocity model. The large non-DC components may suggest complex natures of faulting, for example multifaulting occurred close in space and time, fluid injection, rock damage-related-radiation. Further investigations are required to understand the large non-DC components for shallow tectonic earthquakes in the Los Angeles region.

6. CONCLUSION

In this study, we developed a highly automated and efficient procedure to determine the moment tensor solutions for small-to-medium-sized earthquakes using 3-D velocity models. We applied our approach in the Los Angeles region to demonstrate the importance and feasibility of using 3-D velocity model in automated moment tensor inversion. We generated a new moment tensor catalogue in the Los Angeles region with the completeness of M_L ≥ 3.5. By comparing our catalogue with the current catalogues (the SCEDC-CMT and the YHS catalogues), our results show that incorporating 3-D velocity model can refine the existing moment tensor catalogues in this region, resulting in more accurate focal mechanism solutions, focal depth, and moment magnitude. In addition, by comparing the full moment tensor solutions obtained using 1-D and 3-D velocity models, our studies show that the percentage of non-DC components (CLVD and/or ISO) decrease dramatically with the usage of 3-D velocity model, indicating that the large percentage of non-DC components in the 1-D inversions mainly comes from the unmodelled 3-D velocity structure. Furthermore, our proposed catalogue can be used in full waveform inversion to improve the resolution of 3-D Earth models. Our proposed catalogue can also be used in studying the temporal and spatial variations of stress conditions, the depth distribution of seismicity, the distribution of strong ground motion, as well as contributing to mitigating the seismic hazard and risk in the area. Our automatic inversion approach can be expanded in the Southern California region and can be easily implemented in near real-time system in the near future.

SUPPORTING INFORMATION

Figure S1. Automatic reweighting algorithm. We use a threshold CC (Tcc), which starts from a low value and slowly increase with the number of iteration, to reject the low-CC data (Case 3). In addition, we automatically reduce the weighting factor of the waveform segments that have large RSM (Case 2). Note that the CC, the Threshold CC, and the RSM are dynamic parameters, where their values change with the number of iteration.

Figure S2. The bootstrapping results in Fig. 8(c) in the main text can be grouped into two classes, according to the strike of the fault planes. Although there are two types of focal mechanism solutions, our ensemble solutions are dominated by Class 1 solutions (∼99 per cent), indicating that the fault plane solutions are well constrained.

Figure S3. Comparison among the double couple solutions, the deviatoric moment tensor solutions, and the full moment tensor solutions obtained using 3-D velocity model. Results show that most events have nearly double-couple mechanisms with limited non-double-couple components, and there is little difference between the pure double-couple mechanisms and the results from full moment tensor inversions.

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.

ACKNOWLEDGEMENTS

This research is supported by SCEC 2019 award 19011 and Air Force research grant 18C0058. We have used waveforms and parametric data from the Caltech/USGS Southern California Seismic Network (SCSN) (doi: 10.7914/SN/CI) stored at the Southern California Earthquake Data Center (SCEDC, doi:10.7909/C3WD3xH1). The generalized Cut-and-Paste 3-D (gCAP3D) code is downloaded from http://www.eas.slu.edu/People/LZhu/home.html. We thank Robert Graves for providing the 3-D finite-difference code and useful discussions, and Voon Hui Lai for her help in setting up the initial 3-D finite-difference simulations. We thank Egill Hauksson and Jennifer Andrews for useful discussions about the SCEDC catalogue. We are also grateful to Editor Sidao Ni and two anonymous reviewers for suggestions that greatly improved the manuscript. Sac 2000 and GMT (Wessel et al.2013) were used for basic data processing and figure development.