One of the key problems in seismology is the understanding of earthquake source processes. The kinematic source properties of an earthquake include the seismic moment, fault geometry, dislocation across the fault plane and the spatio-temporal distributions of the slip. When an earthquake occurred, besides the origin time and hypocenter location, the fundamental information of focal mechanism describing the fault geometry and dislocation across the fault plane is directly related to the response of the tectonic stress accumulation and release at the fault zone. Other than the focal mechanism, the rupture parameters including the rupture extent, speed, direction, duration and stress drop would control the behavior of the rupture front and reveal the physical properties of the fault itself. The estimated seismic moment and kinematic source rupture characteristics are the basis for further simulations of dynamic rupture scenario of an earthquake.
Understanding of the lithospheric deformation and tectonic stress regimes relies on continuous long-term observations of seismicity accompanied with a complete and reliable focal mechanism catalog. Routine determinations of focal mechanisms of global earthquakes with moment magnitude Mw5.5 or regional earthquakes in seismically vulnerable areas had been developed and performed by many state
agencies or research institutes, such as Global Centroid Moment Tensor (Global CMT) project (Ekström et al., 2012), the United States Geological Survey (Spikin, 1986), GEOFON (Bormann, 2012), and GEOSCOPE (Vallée et al., 2011). In general, the earthquake focal mechanism can be represented by 5 or 6 independent components of the moment tensor (MT) associated with equivalent body-force couples acting on the
point-source location. Then the moment tensor solution of an earthquake can be simply obtained by a linear inversion of the observed waveform data and theoretical Green’s functions. However, the challenges of retrieving stable and reliable moment tensor solutions reside in the intrinsic noises, imperfect coverage of data and inaccuracy of the Green’s functions because of an inappropriate velocity model and mislocated centroid position used. In recent years, owing to the speedy progress of the computational capacities and data transmission of seismic networks, the real-time moment tensor inversion becomes possible to quickly gain the basic earthquake information for characterizing and responding to potential seismic hazards. Lee et al.
(2014) developed a real-time moment tensor monitoring system (RMT) to
automatically detect and simultaneously determine the epicenters and moment tensors of earthquakes in Taiwan, achieved at the expense of enormous computational
resources. On the other hand, studies of seismogenic structures and seismotectonics demand a complete earthquake catalog that contains reliable and well-resolved
moment tensor solutions. It is thus necessary to carry out a comprehensive exploration of several inversion settings, including the data, centroid location, velocity structure, and others that could lead to unstable or ambiguous inversion results but often were neglected in real-time solutions.
In the first part of my PhD thesis, I focus on establishing a new AutoBATS moment tensor inversion algorithm mainly for abundant regional seismicity in Taiwan.
The goal of the study is to report both near real-time and stable MTs by taking into account the uncertainties of the inversion settings mentioned previously. Different from Lee et al. (2014) conducting a grid search of candidate earthquake locations using continuous longer-period broadband seismic data at a few BATS stations, I prefer to
take the best use of high-precision hypocenters to lower the mislocation effects on the MT solutions, which are determined by much denser CWBSN (Central Weather Bureau Seismograph Network). Instead of fixing the inversion settings for all the earthquakes that occurred in tectonically distinct settings, I would like to deliver the final moment tensor solutions together with the validation of the stability or reliability of the solutions by scanning four different Moho-depth models representative of the structures in the Taiwan region, three strategies of data selection, and three filtering frequency bands for all the used data.
In Chapter 2 , I present the new AutoBATS moment tensor inversion approach and its application to regional earthquakes in the Taiwan area. The reliability of the inversion scheme is verified by making the comparison of our newly obtained CMT catalog with the Global CMT solutions and the previous BATS catalog. Benefiting from this fast, reliable inversion scheme, I also update the regional moment tensor catalog for Taiwan earthquakes from 1996 to 2016. A new relation between moment magnitude (MW) and local magnitude (ML) is also examined based on our complete AutoBATS MT catalog. Figure 1.1 illustrates the flow chart of the newly developed AutoBATS moment tensor inversion method which has already been implemented in the real-time moment tensor report system of Data Management Center of IES (DMC-IES).
Figure 1.1 Flowchart of new real-time AutoBATS MT inversion scheme.
In general, a point-source approximation for major and great earthquakes is not sufficient to characterize the earthquake rupture propagating away on the finite-size fault plane. Therefore, the second part of the thesis is devoted to investigate the rupture characteristics of large earthquakes using high-density, high-frequency seismic array data.
Over the past several decades, numerous studies have employed different waveform inversion methods such as the multiple point source moment tensor
inversion (Kao et al., 2000; Chen and Wen, 2015) and finite fault inversion (Antolik et al., 1996; Ji et al., 2002; Ye et al., 2012) to unravel the complex rupture history of many devastating earthquakes. These methods intend to recover the details of the rupture process by modeling a finite earthquake source as multiple subevents or discretized subfaults on the fault plane with variable source parameters. Although the
forward modelling theory that forms the foundation for the source inversion and the numerical inversion algorithm have been developed for years and become routinely applied to global large earthquakes, the finite-source inversion is known to be a notoriously ill-posed and underdetermined geophysical inverse problem (Ide, 2015) with too many unknown, poorly constrained model parameters involved in the
inversion. As a result, it leads to highly non-unique results that strongly depend on the imposed a priori constraints such as the fault geometry, rupture speed, rise time, and etc. In addition, the relative complexity of high-frequency seismograms often affected by strong crustal velocity heterogeneity but not accounted properly in the waveform inversion limits the use of shorter period waveform data in the inversion.
Contributed by the worldwide deployment of dense seismic arrays such as USArray, ORFEUS (European network) and Hi-net (Japanese network), a high-resolution back projection (BP) technique alternative to the finite source waveform inversion have been widely applied to imaging the rupture propagation history of large earthquakes. In contrast to the inversion approach, the BP method requires neither the adjustment of a priori damping and smoothing parameters imposed on the inversion models nor the knowledge of the Green’s functions, fault geometry, and rupture kinematics. It basically retrieves and stack coherent, high-frequency teleseismic seismic wavefield recorded across the array and project them back in time to track the positions of the strongest energy radiators among subfault grids throughout the rupture duration. As the onsets of the wavetrains used in the BP imaging are
aligned properly, the errors in predicted travel times due to the heterogeneous velocity structure between the source and receiver are negligible. Ishii et al. (2005) first
demonstrated the efficacy of the BP method to image a ~1300-km long rupture
propagating northward for the 2004 Sumatra-Andaman megathrust earthquake. Since then, it has become an effective approach to illuminate the rupture process for major earthquakes, such as the 2008 Sichuan earthquake, 2010 Chilean earthquake, and intermediate- and deep-focus earthquakes (Xu et al., 2009; Kiser et al., 2011b; Kiser and Ishii, 2012; Koper et al., 2012).
The classical BP approach, similar to the time-domain beamforming technique, shifts the traveltime moveout from a trial subfault grid relative to the hypocenter and stacks the moveout-corrected waveforms across the array to locate the energy radiation point associated with the maximum amplitude of the stacked waveforms integrated over certain time intervals. However, Ishii et al. (2007) noticed the spatial and temporal artifacts of the ghost stacks from the synthetic BP experiments that could cause the biased estimate of the time and location of the energy radiators. Such the ghost swimming artifact mainly arises from nonstationary seismic wavetrains within which the amplitudes are highly variable such that the weak arrivals in the beamforming image can be obscured by the movement of large-amplitude arrivals (Walker and Shearer, 2009; Xu et al., 2009; Meng et al., 2012a, 2012b). Meng et al. (2012b) developed a more viable BP method for non-stationary seismic waveforms based on the subspace-based method, the Multiple Signal Classification (MUSIC) algorithm. It has been demonstrated that the MUSIC BP technique provides superior resolution of the source images compared to that obtained with the classical time-domain BP method by adopting the “reference window” strategy to mitigate the artificial
swimming effect efficaciously. Despite of the improvement of the BP method, it is still difficult to entirely suppress the smearing effect as illustrated in Figure 1.2b because of the intrinsic trade-off between travel time and distance. Therefore, most of the BP
imaging studies is restricted to 2D applications only, i.e. the candidate radiation sources confined in an assumed horizontal plane positioned at the focal depth of the earthquake.
On May 24, 2013, the largest deep earthquake ever recorded by modern seismic networks struck about 609 km beneath the Sea of Okhotsk. The focal mechanism from the Global CMT solution indicates this deep earthquake ruptured on either a
sub-horizontal or sub-vertical fault plane. The finite-fault inversion reaches the inconclusive fault-plane orientation because the results of assuming either one of the nodal planes yield equal fits to the observed and synthetic waveforms (Ye et al., 2013).
On the other hand, Chen et al. (2014) and Zhan et al. (2014) performed the multiple point source inversion and directivity analysis but obtained opposite conclusions on the actual fault plane. In Chapter 3 , I attempt to develop a three-dimensional (3D) MUSIC BP method (Figure 1.2a) to resolve the ambiguity between the sub-vertical and
sub-horizontal fault plane. To facilitate suppressing the smearing effects in the BP image and retrieve the stable rupture propagating process in 3D, I include the 3D BP images obtained by depth phases, pP, which travel upward to the surface and turning back to the distant stations and have the take-off angle different from the down-going direct P-wave. Since the smearing zone is parallel to the ray path from the source to the receiver array, combining the BP images obtained with different phases would help reduce the smearing area (Figure 1.2b).
Figure 1.2 Illustration of a 3D back projection method and the smearing effects.
(a) Configuration of a 3D source volume discretized into grids of subfaults and the source-to-array geometry. (b) The elongated smearing zones expected to be seen in the BP images obtained with teleseismic P and depth phases.
The combined P- and pP- wave 3D BP image of the 2013 Okhotsk earthquake seen from the North American (NA) seismic network shows a two-stage, anti-parallel rupture. The rupture propagates about 30-40 km northeastward during the first stage; it then jumps to greater depths and propagate at least 80 km in the reverse (S-ward) direction. The two anti-parallel ruptures are separated by a vertical distance of about 10-15 km. The average rupture speeds during the two stages are 3.0-3.3 and 4.25-4.80 km/s, respectively. The initial NE-ward rupture ends near the northern tip of the
subducted Pacific slab that most likely has been heated up by the ambient warm mantle and asthenospheric flow around the slab edge (Peyton et al., 2001; Park et al., 2002;
Davaille and Lees, 2004a). This may cause a relatively slower speed of the first NE-ward propagating rupture as well as prevent it from continuously growing.
In order to understand the resolution and limitations of the 3D MUSIC BP
technique, I perform a suite of synthetic experiments that mimic several possible rupture scenarios for the Okhotsk earthquake, including the rupture propagating vertically up- or down-ward, and horizontally perpendicular or parallel to the source-to-array azimuth. All of these synthetic results show the capability of the 3D P-wave BP MUSIC method to recover the hypothetical subevents. As expected, the combined P- and pP-waves BP images can help stabilize the resolved rupture image and locate the energy radiators with much smaller spatial uncertainties than that obtained with P waves only.
As the BP method mainly uses high-frequency teleseismic wavefield for source imaging, there are few applications to shallow earthquakes with magnitude less than 7 that could also produce destructive hazards. There are on average one to three
earthquakes with magnitude greater than 6 that hit the Taiwan region every year, and some of them have caused severe damages and fatalities such as the 2016 Meinong earthquake and 2018 Hualian earthquake. To better understand the rupture kinematics of these earthquakes and corresponding seismogenic fault zones, I applied the MUSIC BP imaging and directivity analysis to unravel the rupture process of two strong earthquakes, the 2016 Meinong and 2010 Jiashian earthquakes that recently occurred nearby in southern Taiwan with comparable magnitudes and focal mechanisms.
In Chapter 4, I present the first application of the MUSIC BP method to quantify the details of the rupture process for a moderate earthquake, the 2016 Mw 6.4 Meinong earthquake that severely hit southwestern Taiwan, a region with relatively low
seismicity in the past and caused the deaths of over 100 people in the populated Tainan city about 30 km west of the epicenter. I use high-frequency P and depth-phase sP wavetrains from the dense European (EU) and Australian (AU) seismic networks. The
P-wave BP results either from the EU or AU array show a consistent NW-ward rupture extending about 16 km with a single-peak radiation burst and duration of ~7 s. Besides, through forward waveform modeling, the slowness-backazimuth analysis, and
radiation pattern of depth-phase sP/pP amplitude ratios, I confirm the second displacement pulse that begins ~8 s after the onset of the first P-wave pulse mainly consists of the depth-phase pP and sP energy. Therefore, I conduct the same BP imaging using the aligned sP wavetrains. The integrated P- and sP-wave BP images resulting from the EU and AU array data are overall in accordance with each other showing a unilateral subhorizontal rupture propagating northwest. To verify the
kinematic rupture properties derived from the BP imaging, I also employ the directivity analysis assuming a kinematic unilateral rupture model using global broadband
displacement waveforms.
As demonstrated by the study of the Meinong earthquake, the MUSIC BP
technique has shown great success in unraveling the rupture behavior of moderate-size earthquakes generally with relatively short rupture length about tens of kilometers and source duration about a few seconds. Therefore, I attempt to apply the same approach to studying the rupture process of the 2010 Jiashian earthquake which shares several seismological similarities with the 2016 Meinong earthquake, including the nearby epicenter, similar focal mechanism and comparable seismic moment. The results are summarized in Chapter 5.
According to the integrated P- and sP-wave BP images obtained with the EU and AU array data, I found the Jiashian earthquake predominantly ruptured toward
northwest as the Meinong earthquake did. However, unlike the Meinong event with a nearly horizontal rupture, the Jiashian earthquake ruptured obliquely in both the
along-strike and updip direction, which is about 16-20° from the horizontal plane. The estimated rupture length and duration is about 10-11 km and 4-5 s, respectively, which give the average rupture speed similar to that of the Meinong earthquake. Similarly, I conduct the directivity analysis to validate our BP results. I perform a grid search in the space of the rupture length and velocity, and the strike and dip of the rupture
propagation direction and find the optimal solution that yields a global minimum misfit of observed and predicted apparent rupture durations. The F-test shows that the
additional testing variables, the rupture strike and dip, are statistically significant to improve the fit between the observed and predicted source time duration.
Motivated by the similarity of the seismic moment tensor and rupture properties between the Jiashian and Meinong earthquake as well as the spatial distributions of the mainshock and aftershock sequences, I proposed a blind subsurface fault striking NW and dipping to NE (hereafter referred as to the JSMN fault) is the causative fault responsible for the occurrence of the two earthquakes. Besides, I investigate the temporal b-value variations and static stress changes induced by the two events. By compiling all the observations and analyzing results, I propose the Jiashian and Meinong earthquakes are the doublet events which have occurred as a result of failure of the two strong asperities on the JSMN fault. After continuous stress buildup on the fault for at least 45 years mainly exerted by the plate convergence, the accumulated stress exceeding the strength of the smaller asperity on the deeper eastern end of the fault released and initiated the rupture of the Jianshian earthquake. The static Coulomb stress transfer triggered by the slip of the Jianshian event suggests an additional stress increase in the vicinity of the shallower western asperity. With additional tectonic stress loading for six years, the failure of the shallower and stronger asperity took place
and initiated the rupture of the Meinong earthquake.
Chapter 2 A new automatic full-waveform regional moment tensor