Simultaneous Determination of Earthquake Source Parameters Using Far-Field P Waves : Focal Mechanism , Seismic Moment , Rupture Length and Rupture Velocity

1 Institute of Geophysics, National Central University, Chung-Li, Taiwan, ROC 2 Seismology Center, Central Weather Bureau, Taipei, Taiwan, ROC 3 Center of General Education, Vanung University, Chung-Li, Taiwan, ROC * Corresponding author address: Dr. Tzu-Wei Lin, Seismology Center, Central Weather Bureau, Taipei, Taiwan, ROC; E-mail: daynight@scman.cwb.gov.tw We have developed a new method for the simultaneous determination of the earthquake source parameters: seismic moment, focal mechanism, rupture length, and rupture velocity, using far-field P-waves. The pseudo radiation patterns of P-waves are obtained from each station, accomplished through a grid search of the source time function and simple inversion. The observed waveforms are corrected for instrument response, geometrical spreading and surface effects. These corrections are made for the corresponding ray parameter of each station based on IASP91 velocity model. To obtain the best solution for the focal mechanism, pseudo radiation patterns are compared with the theoretical radiation patterns calculated for the given focal mechanisms obtained through the grid search. The source time functions obtained from all available stations for an earthquake are used to estimate the rupture length and rupture velocity based on rupture directivity. We have applied this method to three moderateto large-sized global earthquakes. The results agree well with those of previous studies, suggesting that this method is applicable for seismic hazard assessment. A moderate Taiwan earthquake is analyzed to demonstrate that the method can indeed be used to rapidly determine the source parameters of earthquakes in Taiwan. (


INTRODUCTION
The source parameters of an earthquake are usually its magnitude, source time function, focal mechanism and seismic moment.Amongst these parameters, the source time function includes information regarding an earthquake's rupture length and rupture velocity, characteristics that can be determined through rupture directivity analysis based on the finiteness source theory (cf.Ben-Menahem 1961;Aki and Richards 2002).Such information is essential for understanding the physical behavior of earthquake rupturing.
Presently, moment tensor inversion (Harvard CMT solutions and NEIC fast moment tensors) is the standard technique used to routinely determine the focal mechanism and seismic moment of earthquakes (Dziewonski et al. 1981;Sipkin 1994).However, pulse width of the source time function at each station needs to be taken as a particular constant in the moment tensor inversion method to solve for the optimal fault plane solution of an earthquake.Traditionally, rupture directivity is mainly based on the variation of the source time function with the station azimuth, which is generally derived from forward waveform modeling, the deconvolution of diffracted P-waves, and the empirical Green's function of surface waves (e.g., Kanamori and Stewart 1976;Chung and Kanamori 1980;Ruff and Kanamori 1983;Ammon et al. 1993).Although the empirical Green's function method can efficiently estimate rupture directivity for large earthquakes, it has strict constraints.The focal mechanism and source depth of the main shock must be analogous to those of the aftershocks (e.g., see Ammon et al. 1993).In addition, the rupture directivity of moderate-sized earthquakes is generally more difficult to estimate by a method such as the empirical Green's function.Moreover, forward waveform modeling is restricted by the correct fault plane solution and source depth, factors which can directly affect the calculation of theoretical far-field seismograms and the estimation of rupture directivity (e.g., Kanamori and Stewart 1976;Chung and Kanamori 1980).
In view of the inconveniences mentioned above, we propose a novel far-field P-waveform inversion procedure that can systemically determine source parameters of an earthquake.This procedure does not resemble the moment tensor inversion or forward waveform modeling; and it can rapidly and efficiently estimate the fault plane solution, seismic moment and rupture directivity for moderate-to large-sized earthquakes.We first determine the variation of the source time function with the station azimuth, and pseudo radiation patterns, i.e., the actual radiation patterns multiplied by the seismic moment, for each given station (see Section 2).The focal mechanism and seismic moment are then calculated from pseudo radiation patterns made on the basis of a grid search for the strike, dip and slip angles of an earthquake fault.Variations in source duration vs. station azimuth can also be used to estimate rupture directivity during earthquakes through several simple fault models with various rupture propagation types (cf.Ben-Menahem 1961;Chung and Kanamori 1980).

METHODOLOGY
Following the notations of Kanamori and Stewart (1976) and Okal (1992), a far-field P-wave displacement can be expressed by: where u t p ( ) is the far-field seismogram; M 0 is the seismic moment; α h , β h , and ρ h are the P-wave velocity, S-wave velocity and density, respectively, in the focal region; g( ) ∆ is the geometrical spreading factor as a function of epicentral distance and focal depth; and r is the radius of the Earth (about 6371 km).R p , R pP , R sP and are the radiation patterns for the P-, pP-, and sP-waves, respectively (cf.Aki and Richards 2002); V pP represents the reflection coefficient of upgoing and downgoing P-waves at the surface; and V sP is the reflection coeffi- cient of upgoing S-and downgoing P-waves.i h and j h represent the take-off angles of the Pand S-waves from the source, and i o is the incident angle for the P-wave at the free surface.

C i p o
( ) is the free surface effect at the receiver as a function of the incident angle; f t ( ) is the source time function, which is specified to be an isosceles triangle with an area of 1.0 in this study; t p , t pP , and t sP are the travel times for the P-, pP-, and sP-waves, respectively, which can be calculated according to a given Earth velocity model.Q t ( ) is the attenuation factor of the Earth based on Azimi's law (cf.Yoshida 1988), and I t ( ) is instrument response.After correcting for instrument response, equation (1) can be rewritten as: where u t c P ( ) is the far-field seismogram corrected for instrument response, and is the combination of several specific terms.
We rearrange equation (2) in matrix form, Ax b = , as follows:  * ( ) where n is the number of total data points in the time series.The vector x denotes pseudo radiation patterns, i.e., the actual radiation patterns ( R p , R pP , R sP ) multiplied by the seismic moment M 0 , which is the objective to be solved.Matrix A gives synthetic waveforms for the P-, pP-, and sP-waves, exclusive of the seismic moment and the radiation patterns.The attenuation factor Q t ( ) is applied by taking t * .=1 0 (the travel time over the attenuation factor) for the P-wave in this study (e.g., Yoshida 1988;Okal 1992).Vector b represents observed waveforms after corrections for instrument response and P ara , which is calculated once a unique corresponding ray parameter is determined.The ray parameter is calculated in this study using the IASP91 velocity model.
For a given station, with a known ray path, we can use equation ( 3) to invert the pseudo radiation patterns by a given source time function with specific source duration.In other words this equation can be used to simultaneously derive information regarding the source, including rupture directivity and focal mechanism.After testing the source time functions with various source durations, we can determine the optimal source duration and corresponding pseudo radiation patterns, when the misfit between the observed and theoretical seismograms reaches a minimum value.The misfit R 1 is defined as: where ( − ⋅ is the difference between the observed and theoretical seismograms at the ith data point.Other notations are the same as in equation ( 3).
In the inversion process, to search for the optimal source duration, the upper bound of the source duration is chosen by an empirical relationship between the seismic moment and the source duration (cf.Furumoto and Nakanishi 1983).The lower bound of the source duration is dependent on the filtering range.Equation ( 3) is used to obtain the optimal source duration and pseudo radiation patterns for all available stations.The pseudo radiation patterns can then be used to derive the fault plane solution; and the variation of the source duration with the station azimuth can be used to estimate the rupture directivity of an event.First, the focal mechanism is decided by comparing the pseudo radiation patterns with the theoretical radiation patterns, for all available stations.The theoretical radiation patterns are calculated on the basis of the formulas of Aki and Richards (2002).We search for values for strike of 0° -360°, dip of 0°-90°a nd rake of -180° -180° at an interval of 2°.The optimal fault plane solution and seismic moment are determined through a linear regression y ax = , where y and x are the pseudo and theoretical radiation patterns of all available stations and a is the slope which implicitly indicates the seismic moment.There are two constraints in the above-mentioned linear regression: one is that the fitting line has to go through the origin, and the other is that a positive slope is required.In this study, only the pseudo radiation pattern of the P-wave is adopted because the estimation of the fault plane solution obtained using the P-wave is more stable than that of combining the P-, pP-and sP-waves.The misfit between the theoretical and pseudo radiation patterns is defined as: , where y j and x j are the pseudo and theoretical radiation patterns of the P-wave at the jth station, respectively, and m is the number of available stations.
In addition, the variation of the source duration vs. the station azimuth can be used to estimate the rupture directivity of an earthquake as (cf.Ben-Menahem 1961; Chung and Kanamori 1980): where T SDT is the source duration observed for a given station, L is the rupture length of the fault, V is a parameter indicating the rupture velocity, C is a parameter indicating the P-wave velocity in the source area, and Φ is a parameter related to Θ, which is the angle measured between the rupture direction and the ray radiating from the initial source point to the station.For an earthquake with unilateral faulting, , where V R is the rupture velocity and V P is the P-wave velocity.In addition, for one event with bilateral for a unilateral rupture and a positive slope ( ) for a bilateral rupture.The values of intercepts for both types of faulting are positive.These conditions are always helpful in analyzing the rupture directivity of an earthquake.Using equation ( 4), we can estimate its intercept and the slope; and the rupture length and rupture velocity can be estimated with the P-wave velocity given in the source area.In addition, rupture directivity analysis can also be used to distinguish faulting type and fault plane utilizing optimal linear fit from equation (4).

DATA AND RESULTS
Broadband teleseismic P-waves for stations with an epicentral distance of 30° -90° are retrieved from the IRIS database.This is to avoid near-field effects and shadow zone interference.After correcting for instrument response, the recordings are transferred into displacement seismograms.In this study, we extracted the P-waveforms, including the P-, pP-, and sP-waves, with time lengths of about 30 -55 sec (i.e., 5 sec before and 25 -50 sec after the P-arrival), for earthquakes with moment magnitudes ( M w ) of 5.8 -7.3.The time length of the P-wave mainly depends on the size of the earthquake, and can be roughly estimated by an empirical formula (see Furumoto and Nakanishi 1983).We employed a band-pass filter to clarify the P-waveforms but not to change the original displacement waveforms too much.
We tested the proposed method, as described in Section 2, by applying it to data for three moderate to large global earthquakes as listed in Table 1.We further apply the technique to a local earthquake in Taiwan (also see Table 1) in order to confirm whether the technique can be used to provide rapid determination of source parameters of local earthquakes for the purposes of hazard assessment.The estimated fault parameters of the four earthquakes are delineated below.
Table 1.Basic information for the four moderate-to large-sized earthquakes analyzed in this study.3.1 Northridge, California, Earthquake of January 17, 1994 ( M w 6.7) The 1994 Northridge earthquake occurred on January 17, 1994 and was located at longitude of 118.54°W, latitude of 34.21°N and a depth of 18.4 km (Table 1).Figure 1 maps the distribution of stations valuable for analyzing the integral rupture characteristics of the earthquake.The P-wave velocity, S-wave velocity and density in the source area are 5.80 km sec 1 − , 3.36 km sec 1 − , and 2.45 g cm 3 − , respectively.The duration of the P-wave displacement used in the analysis is 30 sec (5 sec before and 25 sec after P-arrival), and the waveforms are bandpass filtered with corner frequencies of 0.05 and 0.5 Hz. Figure 2 also shows a comparisons between the synthetic and observed waveforms made through the proposed inversion process.
Next the pseudo radiation patterns and the source duration for each station are derived by equation (3).A comparison of theoretical radiation patterns with pseudo radiation patterns via a grid search of the fault status at 2° intervals shows that both radiation patterns exhibit a good linear relationship (see Fig. 3a).The determined double coupled focal mechanisms for strike, dip, and rake are: (113°, 35°, 102°) and (278°, 56°, 82°), respectively, with a seismic moment of 8 93 10 18 .× Nm ( M w 6.6).These results compare well with the results from the Harvard CMT; Wald et al. (1996); and Thio and Kanamori (1996), as shown in Fig. 3b.From the    distribution of aftershocks (Fig. 3b), which is limited to a region of about 30 × 30 km (Hudnut et al. 1996;Wald et al. 1996), the fault plane can be confirmed to be 113°, 35°, and 102°.We calculated cos Θ from the source to each station, according to the determined fault plane, and then used equation ( 4) to link cos Θ with the source duration.The rupture directivity analysis of the 1994 Northridge earthquake shows that the earthquake was a unilateral faulting event with an optimal rupture angle of 109° measured counterclockwise from the strike of the fault on the fault plane (Fig. 4).The rupture length is then about 17.2 km with a rupture velocity of 2.95 km sec 1 − (about 0.88 β h ) and source duration of about 5.8 sec.These results are in good agreement with the distribution of aftershocks and previous studies based on GPS and seismic data (Hudnut et al. 1996;Wald et al. 1996;Thio and Kanamori 1996).
3.2 Landers, California, Earthquake of June 28, 1992 ( M w 7. 3) The basic source parameters of the 1992 Landers earthquake are listed in Table 1.Recordings from sixteen stations are displayed in Figs. 5 and   P-wave waveforms for the 1992 Landers earthquake corresponding to the optimal focal mechanism derived from this study.
and Harvard CMT (Fig. 7b).Rupture directivity analysis confirms that the 1992 Landers earthquake was a unilateral faulting event with the best fault plane of 341°, 80°, and -178° (see Fig. 8).This is consistent with the distribution of aftershocks as shown in Fig. 7b.The rupture length and rupture velocity are estimated to be about 98 km and 3.02 km sec 1 − , with source duration of 32.4 sec.These values are again comparable to the results of other studies (Wald and Heaton 1994;Dreger 1994;Cohee and Beroza 1994;Unruh et al. 1994;Lazarte et al. 1994).

Northern Iran Earthquake of May 28, 2004 ( M w 6.3)
Seismograms from the 53 stations used are displayed in Figs. 9 and 10.A band-pass filter of between 0.05 and 1.0 Hz is applied to the waveforms after correction for instrument response.The P-wave velocity, S-wave velocity, and density in the source area are 5.80 km sec 1 − , 3.36 km sec 1 − , and 2.45 g cm 3 − , respectively.Observed and synthetic seismograms are compared in Fig. 10.
The relationship between the pseudo and theoretical radiation patterns is demonstrated in Fig. 11a.
The results in Fig. 11b show that the best focal mechanism has a strike of 95°, a dip of 31°,  a rake of 58°, and a seismic moment of 2 51 10 18 .× Nm, similar to the results from several other institutes (Harvard CMT; USGS; ERI).According to the rupture directivity analysis in Fig. 12, the earthquake was a unilateral faulting event, with a rupture length of 6.3 km, a rupture velocity of 2.74 km sec 1 − , and source durations of about 2.3 sec.

3.4
Chia-Yi, Taiwan, Earthquake of October 22 , 1999 ( M w 5.8) The above study focusing on moderate and large global earthquakes shows that reliable results can be obtained by the proposed method.We now apply the technique to the 1999 Chia-Yi, Taiwan, earthquake.The earthquake occurred at a depth of about 17 km, with a Pwave velocity, S-wave velocity, and density of 5.80 km sec 1 − , 3.36 km sec 1 − , and 2.45 g cm 3 − .
Seismic recordings from forty-four stations are used for the analysis (see Figs. 13 and 14).The waveforms are band-pass filtered from 0.05 to 1.0 Hz.The solution of the proposed inversion process (see Fig. 15a) gives the following results for the focal mechanism of the earthquake: strike = 179°, dip = 55°, rake = 80°, and a seismic moment of 1 18 10 18 .× Nm, which are comparable with routine reports (USGS; Havard CMT; CWB; IES BATS) and with a study      made by Chang (2004), as shown in Fig. 15b.According the distribution of the aftershocks, the fault plane would be in strike, dip, and slip: 179°, 55°, and 80°, respectively, to be the rupture plane (see Fig. 15b), which is used to analyze rupture directivity (see Fig. 16).The rupture length, the rupture velocity, and the source duration are then estimated to be about 7.4 km, 1.80 km sec 1 − , and 4.1 sec, with unilateral faulting.

DISCUSSION AND CONCLUSION
It is of great importance when clarifying the rupture characteristics of earthquake faulting to derive fault parameters and rupture directivity.In order to systemically estimate the fault parameters and rupture directivity, we propose an inversion process with P-waveforms.The proposed inversion process is different from previous methods, such as the empirical Green's function method (e.g., Ammon et al. 1993), the deconvolution of diffracted P-waves method (e.g., Ruff and Kanamori 1983), and forward waveform modeling (e.g., Kanamori and Stewart 1976;Chung and Kanamori 1980).We determine the best focal mechanism of an earthquake through a grid search method.This also differs from the moment tensor inversion method (e.g., Dziewonski et al. 1981;Sipkin 1994).The proposed method does not estimate the source depth, as in the study based on USGS and Harvard CMT reports.In the proposed inversion process a simple method based on the global mean velocity model (IASP91 earth model) is adopted to calculate far-field P-wave seismograms.However, if a more complicated crustal model is used it does not significantly change the estimated parameters.
The proposed inversion process has been successfully applied to four earthquakes with M w 5.8 -7.3 (Table 1).The estimated fault parameters and rupture directivity for the four earthquakes are quite consistent with those made by previous studies.All four earthquakes, except for the 1992 Landers earthquake which had a strike-slip mechanism, were found to have thrust mechanisms (Table 2).The rupture directivity analysis only confirms the best fault plane of the 1992 Landers earthquake, the rest can not be clearly done.In other words, the double-couple rupture directivity analysis method is efficient at distinguishing the best fault plane of strike-slip earthquakes.It is difficult to determine the rupture plane of earthquakes with thrust mechanisms due to the similarity of the rupture directivities of the two fault planes.Hence, the fault planes of the other three must be fixed by the distribution of aftershocks.In additional, the four earthquakes are all revealed to have unilateral rupturing; and the direction is measured counterclockwise to the strike of the fault on the fault plane.
The rupture lengths and source durations of the four earthquakes lengthen as the seismic moment (or magnitude) increases.This is consistent with global observations made by Furumoto and Nakanishi (1983), who gave the empirical relationship between the seismic moment, source duration and rupture length of thrust-type earthquakes of M w 7.0 or greater.The rupture length of the 1992 Landers earthquake estimated in this study is 98.0 km and the source duration is 32.4 sec, slightly different from the values of about 84.0 km and 46.8 sec calculated from the empirical relationship by Furumoto and Nakanishi (1983).The other three earthquakes with M w < 7.0 have shorter rupture lengths and source durations than those extrapolated from Furumoto and Nakanishi's empirical formulas.This might imply that for relations between fault parameters that deviate from Furumoto and Nakanishi's observations for earthquakes of M w < 7.0 the physical behavior of moderate-sized earthquakes might be different from that of large-sized ones.It can be seen that the proposed method can systemically and efficiently determine the fault parameters of moderate-sized earthquakes, so that, scaling law for moderate-sized ones can be derived.This is useful in providing more information for understanding the physical behavior of earthquakes.
In this study, we present a promising approach to determine the focal mechanism, seismic moment and rupture length of an earthquake.Future efforts related to this work will focus on developing a semi-automatic data process based on the method to quickly estimate source parameters.Moreover, for seismic hazard assessment, this method can be efficiently applied to regional earthquakes on account of proper local velocity structure model and real-time data transmission.
Acknowledgments Special thanks are due to IRIS for providing high-quality seismic waveform data.The authors also thank Dr. Paul Wessel and Walter H. F. Smith for developing the graphic software GMT.This study was funded by the National Science Council, R.O.C. under Grant Nos.NSC91-2119-M-432-001 and NSC92-2119-M-238-002. Table 2. Major source parameters for the four moderate-to large-sized earthquakes obtained in this study.
Note: *1 indicates the rupture angle is calculated in a counterclockwise direction to the strike of the fault on the fault plane.
Note: *1 represents data from IRIS, *2 stands for data from Harvard CMT, *3 stands for data from IIESC-INSN and *4 indicates data from Central Weather Bureau, Taiwan.

Fig. 1 .
Fig. 1.Distribution of the 34 available stations for the 1994 Northridge earthquake.

Fig. 2 .
Fig. 2. Comparison of the observed (thick lines) and synthetic (thin lines)P-wave waveforms for the 1994 Northridge earthquake corresponding to the optimal focal mechanism derived from this study.

Fig. 3 .
Fig. 3.For the 1994 Northridge earthquake: (a) Correlation between the pseudo and theoretical radiation patterns of P-wave for all analyzed stations.(b) The best fault plane solution obtained by this study in comparison with other results (Harvard CMT and USGS).The shadow circles indicate the distribution of the aftershocks occurring within one year after the main shock.Two profiles at the bottom display the projective cross-sections perpendicular to the strike directions of orthogonal fault planes.

Fig. 4 .
Fig. 4. Estimation of the rupture parameters made by using a least square method and by comparing individual source duration times with cos Θ for all stations.For the 1994 Northridge earthquake: (a) Plane 1-orthogonal fault plane.(b) Plane 2-orthogonal fault plane.
Fig. 5. Distribution of the 16 available stations for the 1992 Lander earthquake.

Fig. 6 .
Fig. 6.Comparison of the observed (thick lines) and synthetic (thin lines)P-wave waveforms for the 1992 Landers earthquake corresponding to the optimal focal mechanism derived from this study.

Fig. 7 .
Fig. 7.For the 1992 Landers earthquake: (a) Correlation between the pseudo and theoretical radiation patterns of P-wave for all analyzed stations.(b) The best fault plane solution obtained by this study in comparison with other results (Harvard CMT and USGS).The shadow circles indicate the distribution of the aftershocks occurring within one year after the main shock.Two profiles at the bottom display the projective cross-sections perpendicular to the strike directions of orthogonal fault planes.

Fig. 8 .Fig. 9 .
Fig. 8. Estimation of the rupture parameters made by using a least square method and by comparing individual source duration times with cos Θ for all stations.For the 1992 Landers earthquake: (a) Plane 1-orthogonal fault plane.(b) Plane 2-orthogonal fault plane.

Fig. 10 .
Fig. 10.Comparison of the observed (thick lines) and synthetic (thin lines) P-wave waveforms for the 2004Iran earthquake corresponding to the optimal focal mechanism derived from this study.

Fig. 11 .
Fig. 11.For the 2004 Iran earthquake: (a) Correlation between the pseudo and theoretical radiation patterns of P-wave for all analyzed stations.(b) The best fault plane solution obtained by this study in comparison with other results (Harvard CMT, USGS and ERI).The shadow circles indicate the distribution of the aftershocks occurring within two months after the main shock.Two profiles at the bottom display the projective cross-sections perpendicular to the strike directions of orthogonal fault planes.

Fig. 12 .Fig. 14 .
Fig. 12. Estimation of the rupture parameters made by using a least square method and by comparing individual source duration times with cos Θ for all stations.For the 2004 Iran earthquake: (a) Plane 1-orthogonal fault plane.(b) Plane 2-orthogonal fault plane.

Fig. 15 .
Fig. 15.For the 1999 Chia-Yi earthquake: (a) Correlation between the pseudo and theoretical radiation patterns of P-wave for all analyzed stations.(b) The best fault plane solution obtained by this study in comparison with other results (Harvard CMT, USGS, CWB and BATS).The shadow circles indicate the distribution of the aftershocks occurring within one month after the main shock.Two profiles at the bottom display the projective cross-sections perpendicular to the strike directions of orthogonal fault planes.

Fig. 16 .
Fig. 16.Estimation of the rupture parameters made by using a least square method and by comparing individual source duration times with cos Θ for all stations.For the 1999 Chia-Yi earthquake: (a) Plane 1-orthogonal fault plane.(b) Plane 2-f orthogonal fault plane.