Automatic determination of the optimal MT solution with thoroughly scanned

To efficiently achieve stable and reliable results for moment tensor inversion, the first important and necessary step is to avoid manual data selection before the

inversion. In this study, we only consider stations whose average signal-to-noise (S/N) ratios of three-component waveforms are greater than 2.0 in the desired frequency band (0.01~0.09 Hz). We take 150 s time windows before and after the P arrival to estimate the spectra and calculate average point-by-point S/N with 5-point moving time-window smoothing. We exclude stations with epicentral distances smaller than 30 km to reduce the mislocation effect. Since many studies typically use three or more three-component broadband stations for regional moment tensor inversion (Dreger and Helmberger, 1993; Zhu and Helmberger, 1996; Kubo et al., 2002; Pondrelli et al., 2002;

Rueda and Mezcua, 2005; Clinton et al., 2006; Ristau, 2008), we selected a minimum of 3 to a maximum of 7 S/N-qualified stations to start the inversion based on three different criteria: (1) the best azimuthal coverage, (2) the highest S/N, and (3) the shortest epicentral distances. For earthquakes inside a homogeneous network, the most logical and easiest approach for criterion (1) is to divide stations evenly into 8 half quadrants according to azimuth and pick at least one S/N-qualified station in each sector. Unfortunately, this simple method is not appropriate for Taiwan due to restricted azimuthal coverage in most cases since the majority of earthquakes occur near the coast and eastern offshore region. To accommodate this natural condition, we retain the stations at the minimum and maximum azimuth, and use these two azimuths as boundaries to divide the remaining stations into three even sectors. We then pick one

to three stations randomly in each sector so that criterion (1) is fulfilled. For criteria (2) and (3), stations are chosen according solely to their mean S/N and epicentral

distances.

It is normally expected that better azimuthal coverage improves reliability of MT solutions. The problem of selecting stations solely based on azimuthal coverage in the Taiwan area is that the selected stations near the coast or on small islands are

occasionally contaminated by long-period ambient noises even though the average S/N might be slightly better than threshold value. Criteria (2) and (3) not only provide other possible choices of stations, but the consistencies among resulting solutions can further ensure the MT solutions are insignificantly biased by those stations with notable noises near the shore or at far distances. Note that the quality (S/N) of regional waveforms is usually site-dependent and unrelated to epicentral distances; therefore, the three different criteria are complementary.

Once stations are selected, we perform the inversions with different parameter settings. To accommodate the velocity variations in the Taiwan region and to improve the stability of inversion with the longest possible periods, we apply three different frequency bands and four velocity models for all selected stations in the inversions. For each earthquake, three consecutive passbands are applied from the followings 5

standard frequency ranges: 0.01~0.04 Hz (25~100 s), 0.02~0.06 Hz (16~50 s), 0.03~0.08 Hz (12~33 s), 0.04~0.09 Hz (10~25 s), and 0.05~0.15 Hz (7~20 s). The choice depends on the local magnitude reported by CWB. For earthquakes with a local magnitude above 5, within 3.5 to 5, and below 3.5, the lowest corner frequency starts from 0.01, 0.02, and 0.03 Hz, respectively.

For the velocity models, it is unrealistic to create a high-resolution 1D velocity

model which reflects all the path effects for each event-station pair because the seismic waves usually travel through distinct structures in Taiwan (especially along the

east-west direction). Considering the waveforms for MT inversion have been filtered for relatively long periods, we start with a simple three-layer model averaged from a 3D tomography model of Taiwan (Rau and Wu, 1995) as a reference. This model, with Moho depth set at 40 km, has been used for 20 years by the current BATS MT

inversion, and simulations in waveforms are reasonably good (Kao et al., 1998a, 2002;

Kao and Jian, 1999; Kao and Chen, 2000; Liang et al., 2004). We then further consider the effect of regional variations of Moho depths from 20 to 50 km under different geologic units that were revealed by receiver-function and tomography studies (Wang et al., 2010; Kuo-Chen et al., 2012; Huang et al., 2014). Since the dominant periods used are more sensitive to structure deeper than 20 km relative to the shallow crust within 10 km (Dahlen and Tromp, 1998), we create a series of velocity models by changing the Moho at 5 different depths ranging from 25 to 45 km with intervals of 5 km. Each candidate potentially represents the average/pseudo structure under the paths covered by the event. Moho depths between 30 and 45 km are used for inland

earthquakes, while a shallower range from 25 to 40 km is considered for offshore cases.

After carefully examining the inversion results, we confirmed that trying different Moho-depth models can effectively reduce misfit and the compensated linear vector dipole (CLVD) component. Following Dreger (2003), we allow Green’s functions to shift by ±2 seconds for each station to lower the misfit due to the possible mislocation and the imperfectness of the velocity model. For fast moment tensor inversion, the Green’s functions of each 1D model are pre-calculated and stored in the databases with a 1 km grid size in distance and focal depth.

For a more realistic simulation of synthetic waveforms, we assign 0.5, 1.0, and 2.0 seconds as duration of source time function for earthquakes with magnitude smaller than 4, between 4 and 6, and above 6, respectively. Since events associated with volcanic activities have been reported in the Okinawa Trough in northeast Taiwan (Lin et al., 2007), our newly developed MT inversion method includes three different isotropic conditions: free-constrained, limited (≤10%), and zero isotropic components.

For the free-constrained or zero isotropic component, we simply invert either 6 or 5 independent elements of the moment tensor. For limited isotropic constrained inversion, we compose a weighting vector inside the kernel matrix, in which the sum of three isotropic elements equals 0, as shown in following

form: [𝑑

0] = (𝑤 𝑤 𝑤 0 0 0^𝐺 )(𝑚11 𝑚₂₂ 𝑚₃₃ ⋯)^𝑇. The inversion is performed iteratively with an initial weighting factor (w) 1, and then adjusted by increments of 1 until the isotropic component of the moment tensor is reduced to 10% or lower.

The AutoBATS MT inversion starts with auto picking of the first P arrivals for each station (Allen, 1978). The system then performs all inversion processes

concurrently with different combinations of settings/parameters mentioned above using the GNU parallel tool (Tange, 2011). Among those inversion results, the one with the minimum misfit and applicable non-double-couple (non-DC) component is reported as the final/preferred MT solution (see Table 2.1 for classifications). Generally, we expect both the CLVD and isotropic components of the tensor to be low for small- to

moderate-size tectonic earthquakes associated with shear faulting and relatively

uncomplicated ruptures. According to Kagan (2002), an additional non-DC component can be introduced due to the data contamination of background noises, an over

simplified velocity model, or the mislocation of earthquakes; therefore, we set up rigorous criteria for non-DC. Only solutions with an isotropic component  20%, CLVD  30%, and non-DC40% are accepted for reporting. For each resultant MT solution, we assign a quality factor (QF) according to the waveform misfit and non-DC component (Table 1). Although the MT solutions with limited non-DC components were reported, the AutoBATS parameter-scanning algorithm is capable of searching for the non-tectonic seismic events if the non-DC criterion is removed.

In this study, we applied the procedure of AutoBATS MT inversion to determine a total of 3058 solutions for earthquakes with a ML greater than 3.5 between May 1996 and April 2016 in the Taiwan region (Figure 2.1). Figure 2.2 and Figure 2.3 demonstrate how well the focal mechanisms are constrained by the AutoBATS

inversion for a shallow earthquake (M_L=4.9) inland and another intermediate-depth event (ML =6.7) offshore, respectively. In each event, the best solution of focal mechanism and the corresponding stations used are highlighted (Figure 2.2a, and Figure 2.3a). Inversion results from different combined settings (such as Moho depth of models, frequency bands, and centroid depths) are directly compared in terms of their misfits and non-DC components. For the 2015/05/25 inland earthquake, the distribution of misfit and non-DC contours shows that filtering at a lower frequency band (0.02~0.06 Hz) and a model with a deeper Moho of 40 km yields the smallest misfit for this normal faulting event with a centroid depth of 10 km (Figure 2.2b). Most inverted tensors are close to pure normal-faulting as indicated by the triangle diagram of focal mechanism (Figure 2.2d). The triangle diagram, developed by Frohlich (2001), is an effective way to demonstrate the fault-type partitioning and the frequency

distribution of the focal mechanisms from the inversion results of all tested parameters.

Figure 2.2e further gives the distribution of the Kagan angles (Kagan, 2003), K-angle hereafter, the smallest 3D rotation angle measured between the best focal mechanism and other solutions. It can be shown that the differences between solutions are more sensitive to the focal depths than the frequency band or velocity model for this earthquake. The overall small misfit is also evident by the excellent agreements between the observed and synthetic waveforms calculated from the final solutions (Figure 2.2c). The uncertainties of moment magnitude, focal depth, and CLVD component are assessed through the standard deviation (S.D.) from all scanned resultant inversions. This example reveals that the non-DC component of the low misfit (<0.4) solutions can vary significantly with Moho depth and frequency band (Figure 2.2b). The variability of solutions also confirms the necessity of scanning various settings while performing a regional moment tensor inversion.

Table 2.1 Quality classifications of the misfit and non-DC component in the released AutoBATS MT catalog

Misfit Category Non-DC Comp. (%) Class

<0.3 A <10 1

0.3~0.5 B 10~20 2

0.5~0.7 C 20~30 3

>0.7* D >30* 4

* Program rejects solutions with non-DC component larger than 40% or misfit greater than 0.75.

Figure 2.2 MT inversion results for the earthquake on 2015/05/25 (ML=4.9) in southern Taiwan.

(a) Map showing earthquake epicenter (star) and distribution of 3 station groups (dist: shortest distance, az: best azimuthal coverage, SN: highest S/N) considered in the inversion scanning scheme. The stations used for the final/best solution are highlighted with bolded outlines. The final MT solution is shown in the upper left corner. The original time, epicenter, focal depth and local magnitude reported by CWB are listed at the top. (b) Contour images of misfit (Top) and non-DC (Bottom) for four varied Moho-depth velocity models, three different frequency bands, and 25 scanned focal depths. The crosses mark the best suited combination of parameters used to produce the final MT solution. (c) The fitness between

observations (solid black lines) and synthetic waveforms (blue dash lines) corresponding to the best solution. For each station, the name of stations, azimuth angles, epicentral distances, and the average misfits of three component waveforms are indicated on the left. The number beside each seismogram is the individual misfit. (d) The distribution of appearance (in percentage) on the focal mechanism triangular diagram from all inversion results with different inversion settings. The symbols are colored with the appearance frequency. (e) Contour images of the K-angle measured from the 3D rotation angle between the final focal mechanism and all other MT solutions. The explanations of axes and symbols are the same as in Figure 2.2b.

Figure 2.3 MT inversion results for the intermediate-depth earthquake on 2014/12/10 (ML=6.7) beneath the Ryukyu subduction zone.

The layout is the same as in Figure 2.2

Figure 2.3 shows the result of an intermediate depth earthquake occurring in offshore Northeast Taiwan on 2014/12/10. This earthquake was excluded in the previous BATS MT catalog because the original depth reported by CWB exceeds the BATS MT calculation range. Our new AutoBATS inversion can provide more

constraints to the subduction zone earthquakes. Unlike the shallow event shown in Figure 2, the misfit, non-DC component, and K-angle of this mid-focus offshore earthquake are relatively insensitive to velocity model, frequency bands, and focal depth, resulting in highly coherent focal mechanisms and less varied misfits (Figure 2.3). Because this earthquake is large in magnitude (Mw 6.0), filtering at a lower frequency band produces a smaller misfit as expected (Figure 2.3d). Although the

azimuthal coverage of stations is restricted in one quadrant, our MT solution and centroid depth estimation are both highly consistent with the solution from GCMT (Ekström et al., 2012). For intermediate depth earthquakes, the centroid depths are sometimes less well resolvable within our ±12 km scanning depth range because the misfit values are comparably small (Figure 2.3b). In addition to a reliable initial location from the CWB catalog, it is always recommended to consider other

independent constraints, such as the depth phases pP and sP observed at teleseismic distance, to further confirm the actual depth of the event.

The AutoBATS MT inversion scheme shows superior capability for determining focal mechanisms. This study successfully resolved 3058 MTs (about 56%) out of 5500 CWB-reported earthquakes with acceptable data quality. Moreover, among these MT solutions, 87% have misfit smaller than 0.7 (above Category C in quality) and 99% have non-DC below 30% (above Class 4). Some resolved earthquakes are as small as MW 3.0 (Figure 2.4a). With the new AutoBATS catalog, we can further investigate whether the mean misfit, value of the CLVD component, and consistency of MT solutions (as described by the K-angle) have any relation with earthquake magnitudes. Figure 2.4b shows that the average misfit decreases from 0.7 to 0.4 as the moment magnitude increases from 3.5 to 5.6. This can be simply explained by the higher data quality against background noise for larger events. For earthquakes with M_W ≥5.6, the average misfit increases again. This phenomenon may relate to the complexity of earthquake rupture, which often involves source directivity and multiple fault segments with different slip vectors, particularly in large events. Waveforms excited by complex source ruptures cannot be well simulated by a single focal mechanism, and so larger misfits may exist. In this case, the results of MT inversion

represent the overall effect of averaged rupture for the earthquakes. Like the features observed for misfit, the mean K-angle shows almost exactly the same trend; the values first decrease and then increase with magnitude (Figure 2.4b). The similarity between mean misfit and K-angle suggests that the stability and consistency of MT inversion results are mainly affected by the waveform quality for small-medium earthquakes and rupture complexity for larger earthquakes.

Contrarily, the mean CLVD appear to be a constant, independent of the magnitude of events (see Figure 2.4c). This phenomenon agrees with Kagan (2002), who

concluded that CLVDs are usually caused by artificial effects, except for

well-examined events. As already demonstrated in Figure 2.2b, the CLVD component of a shallow earthquake is easily affected by different inversion settings. For Taiwan earthquakes, we believe that the strong heterogeneity in tectonic structures and the long-period ambient noises for some stations are mostly responsible for the observed high CLVDs. The lower limit of the average K-angle or CLVD (defined by mean value minus a S.D.) for each Mw bin represents those earthquakes having overall consistent MTs or small CLVDs among all scanning settings. On the other hand, the upper bound (mean value plus a S.D.) reveals that the MT solutions are affected tremendously by different inversion settings. Again, the deviations (error bars) of the average K-angle and CLVD highlight the importance and necessity of probing the best solutions through comprehensive scanning on different inversion settings or parameters.

Figure 2.4 Number of events, Average K-angle and CLVD component with respect to Mw.

(a) Histogram showing the numbers of earthquakes with respect to MW (bin size = 0.2). (b) The relation of average misfit and K-angle with respect to MW based on all AutoBATS MT solutions. (c) The relation of average CLVD with respect to MW. Both error bars in (b) and (c) correspond to a standard deviation.