The Kinematic Source Process of the MW 5.9, 1999 Chia-Yi Taiwan Earthquake from Teleseismic Data Using the Hybrid Blind Deconvolution Method

The kinematics rupture process of the Chia-Yi earthquake occurred on October 22, 1999 (Mw 5.9) in Taiwan is investigated. The hybrid blind deconvolution technique is applied to the teleseismic waves to invert source parameters, including slip amplitudes, rise times and rupture velocities on the fault. From the directivity effect of the fault, the actual fault plane is determined as east dipping. According to the derived ASFT, the total duration of the rupture process is 6.5 sec. A good correlation notices that the larger slip amplitude corresponds to the longer rise time of the subfault. By using the inverted source parameters and combining the stochastic method for finite fault, the strong ground motions of 12 stations at epicentral distances ranging from 3.28 to 29 km are simulated. The results show that the agreements between simulated and recorded spectra are quite satisfying. It means that the inverted source parameters are reliable and the stations where located at near source distance are dominated by source effects. The inhomogeneous distribution of slip and the variable corner frequency could play important roles in the simulation process. Although the source effects are dominant, there are some significant discrepancies existing at stations, implying the site effects are influential.


Introduction
October 1999, at 02:19:0.4 GMT, a strong earthquake of M L =6.4 occurred in the Chia-Yi city near the south of the Chelungpu fault in Taiwan. The event is the larger aftershock of the Chi-Chi earthquake (Mw= 7.6). The Chi-Chi earthquake which ruptured along the Chelungpu fault was the largest earthquake to strike Taiwan in the twentieth century. The Chia-Yi earthquake occurred on a basement-involved reserve fault under western plain. Some parts of the basement faults lie directly beneath cities. In the fact-finding, the event caused at least 12 buildings to collapse in Chia-Yi city [1]. It offers a valuable example of the potential destructiveness of moderate-magnitude earthquake when it occurred in the proximity of densely populated areas.
There exists a good azimuthal covering distribution of the stations in far field used in this study, with regard to the epicenter of the examined event, displayed in Figure 1. The focal mechanism of this event is related thrust fault and the fault plane solution is (strike, dip, rake)=( 46 , 52, 125) by Harvard CMT shown in Figure 2. Since 1991, the Central Weather Bureau (CWB) in Taiwan has constructed an island-wide network composed of more 1218 The Kinematic Source Process of the MW 5.9, 1999 Chia-Yi Taiwan Earthquake from Teleseismic Data Using the Hybrid Blind Deconvolution Method than 700 free-field strong-motion stations [2]. Several stations in the vicinity of the epicenter well recorded and generated seismograms of the event. Therefore, we could compare the recorded with the simulated seismograms to examine the contribution of the seismic source effect to the strong-motions near the epicenter. In general speaking, observed seismic data are generally considered the result of a source time function convoluting with the Green's function, but it is puzzled for seismologists to estimate the reliable Green's function and the source time function. Directly estimating both functions from observed seismic data without any prior information about them is referred to as blind deconvolution. In essence, seismic deconvolution is, indeed, a blind process in that it is difficult to quantitatively measure the source wavelet generated by an earthquake [3]. There are some blind deconvolution algorithms proposed to solve the problems in geophysics. For example, Velis and Ulrych [4] used the fourth-order cumulant function on seismic deconvolution. An iterative algorithm with the Gaussian mixture model was developed by Santamaria et al. [5] to identify the reflectivity sequence on seismic deconvolution. A two-channel blind deconvolution method developed by Zerva et al. [6] was designed to the evaluation of site response. Pflug [7] postulated a relationship between the signal passband and the trispectral domain when the fourth-order cumulant function is applied to seismic deconvolution problem. The generalized blind deconvolution technique offered by Liao and Huang [8] is used to eliminate the ground roll effect effectively and identify the location of reflective signal of a seismic data. Liao and Huang [9] employed the blind deconvolution by Santamaria et al. [5] in conjunction with the water level algorithm to find out the directivity from the apparent source time functions and judge the actual fault plane of the earthquake. Furthermore, Liao and Huang [10] developed a hybrid blind deconvolution (HBD) and combined the GA algorithm to invert the source process of the Alaska earthquake (2002) occurred at Denali fault. In addition, Liao et al. [11] employed the HBD method to investigates the kinematic rupture process of the M L 7.3 Chi-Chi, Taiwan, earthquake on September 21, and discovered the important relationships among the parameters including slip, rise time and rupture velocity distributed on the fault.
In this paper, we use the hybrid blind deconvolution (Liao,et al. [11]) combined with the method introduced by Baumont [12] to investigate the source process of the Chia-Yi earthquake. In order to ascertain the validity of inverted source parameters and observe the source effect affecting on the strong motions, we synthesize the strong motions by using the stochastic method on finite faults. In this method, the source is represented by a rectangular fault plane which is subdivided into a certain number of subfaults, each of which is modeled as a Brune point source. Simulations are created for an observation point by summing the subfaults time series. The stochastic method has been widely applied to different regions to predict the strong ground motions effectively ( [13]- [18]).

Hybrid Blind Deconvolution
In seismic exploration, the observed seismic data  [9] is applied to calculate the source time function and Green's function spontaneously from the seismic data directly. There are two major processes contained in the hybrid blind deconvolution. The first process is to estimate the Green's function by the method offered by reference [5]. The blind deconvolution algorithm from reference [5] can provide an inverse filter } { i f , such that the convolution of } { i f with the seismogram } { i z removes the source wavelet, which makes it possible to suppress the source effect and eliminate the diffraction effect. The second way is applying the estimated Green's function derived from reference [5] to the EGF and STF boxes developed by Bertero et al. [19] to get more stable solutions of source time function and Green's function by a series of iteration. In the EGF box, the projected Landweber method in reference [19] with the causality constraint is used to adjust the estimated Green's function. On the other hand, the projected Landweber method with constrains of positivity and of boundedness of the support is used in the SFT box to modulate the estimated source time function and fit more of the physical properties. The projected Landweber method is an iterative process for approximating the solution of constrained least-squares problems. It is written as the iterative form as follow.
where n f is n-th iteration of source time function, G is Green's function, u represents the seismogram, and τ is the relaxation parameter which satisfy the following conditions.

Source Inversion
The apparent source time functions, which calculated by the hybrid blind deconvolution method retrieved from different stations, describe the moment release watched just at different distances and orientations to the earthquake [20]. From above viewpoint, we applied the method in Baumont and Courboulex [21] to invert the kinematic rupture process of the earthquake. The ASTF is modeled as a time series of Gaussians by Baumont and Courboulex [21], where k p T is the travel time delay between the sub-fault k, and the nucleation point to the station, µ is rigidity, dx and dy are sub-fault dimension. Based on the IASPEI 91 model (Kennett and Engdahl [22]) the apparent velocities were calculated and then k p T were calculated by the apparent velocities. In order to optimize the non-linear inversion of ASFTs, we employ GA algorithm to invert the distributions of ( k A , k T , k σ ) on the fault.

Stochastic Method
The radiated seismic energy in the form of elastic wave has general spectral characteristics at a specific site. Considering the shear wave contribution to the strong motion, the acceleration spectrum ) ( f A of shear wave at a distance r from a fault with seismic moment 0 M is defined as [17,23] where θφ R is the radiation pattern, FS the free surface amplification effects, PRTITN is a constant, ρ and β are the density and the shear wave velocity. The anelastic attenuation Q value is represented by a mean frequency-dependent quality factor for the surrounding area. The source spectrum ) , ( The c f and max f in equation 8 are the corner frequency and the high cut-off frequency respectively. In this paper, we use the relationship between rise time and corner frequency [25] to calculate the c f .
It is a brand-new attempt to figure out the c f by considering the c f as variable on the fault plane. There is a different formula to calculate c f , however, the method considered c f as a constant value on the fault plane [26].
is a high frequency cut-off filter and it has the form of the fourth-order Butterworth filter. It is expressed as the following For the generation of a random time series with the characteristics of a typical strong motion, the duration of the ground motion d T is defined T is the rise time and p T is the path duration depending on the epicentral distance empirical relation [27]. An envelope function ) (t w of accelergram is given for a more realistic shape of the acceleration time series. The envelope function of Saragoni and Hart [28] is used and defined as follow η ξ and are both of constants and ) (t H is the unit-step function.

Apparent Source Time Function
In order to determine the direction of the rupture, we choose 10 teleseismic stations (Figure 1) uniformly distributed around the epicenter of the Chia-Yi earthquake.

1220
The Kinematic Source Process of the MW 5.9, 1999 Chia-Yi Taiwan Earthquake from Teleseismic Data Using the Hybrid Blind Deconvolution Method The focal mechanism of the Chia-Yi earthquake is shown in Figure 2. Figure 3 shows the broadband P-wave displacement recordings of the vertical components of the stations. The apparent source time functions estimated by using the hybrid blind deconvolution method (real line) are shown in Figure 4. Taking account into the physical characteristics of an STF, an STF is the function that describes the time history of the seismic moment released on the fault. In other words, an observed STF or an STF retrieved from a specific station only describes the moment release watched just at that distance and orientation to the earthquake [20]. Based on the above viewpoint, we can derive the rupture process and rupture direction from the ASTFs. In the Figure 4, it can be seen that the shapes and amplitudes of the apparent source time functions are variable with the different azimuth angles. The distinctive differences among the shapes and amplitudes of the ASTFs allow us to take conclusions with respect to the emergence of directivity effects. The ASTFs obtained at different stations are composed of at least two major sub-events; it is reasonable to conclude, therefore, that this earthquake was a complicated event. It is trivial that the simple Brune model in seismology is not suitable to explain the complex rupture behaviors of this Chia-Yi earthquake.    For a source viewed from a large distance, the ASTF duration τ changes with the azimuth and can be expressed as (Nakanishi [29]): where 0 τ is the non-azimuthal source duration; k is a constant; φ is the azimuth of the station viewed from the epicenter; and 0 φ is the direction of the propagation of the rupture. The relationships between pulse widths durations of the ASTFs and the azimuths of the stations can be investigated by employing the least-squares method. The parameters ( 0 τ , 0 φ ) is equal to (6.5, 324.3) and the results are shown in Figure 5. If the rupture velocity were assumed 3km/s, the dimension of the 1222 The Kinematic Source Process of the MW 5.9, 1999 Chia-Yi Taiwan Earthquake from Teleseismic Data Using the Hybrid Blind Deconvolution Method rupture plane can be estimated approximately 18km at least. It indicates that the azimuth of the rupture is about 324 o , which suggests that the rupture propagated in the northwest direction. The direction of slip vector of plane 1, one of the focal planes of the earthquake ( Figure 2) agrees with the retrieved result from our study. It follows that plane 1 can justifiably be recognized as the real fault plane.
The east dipping slip model is though as the preferred model. It is well consistent with the result of strong motion inversion in reference [1].

Source Model
With the ASTFs of 10 stations, we can invert the distributions of ( k A , k T , k σ ) along the fault by GA algorithm to optimize the non-linear inversion of the ASFTs. We fix the dimensions of the fault at about 18km × 16km with the epicenter located at the depth of 16km. The fault is divided into equal area 2km × 2 km sub-faults. Figure 6 shows the distribution of the slip amplitudes on the fault. The slip amplitudes range from 0.004m to 2.2m, and their average is about 0.38m. In general, the slip appears to be concentrated in two patches. One of the patches located at about 2 km of the right lower part from the nucleation point, the other located the left upper corner on the fault. Figure 7 displays the distribution of the rupture velocities on the fault. It is worth noting that the region with the higher value of 3.5 km/s is located at about 8-12 km in strike direction of the fault. In this region, it had higher rupture velocities and lower slip amplitudes; therefore, it means that the fault in this region could be relatively weak. It may be attributed to the weakness of low-density rocks in the zone (Ozacar et al., [30]) or the fluids in this region might have such reduced friction that tectonic stress was produced by the relatively thin, brittle part of the crust (Fisher et al. [31]). The distribution of rise time on the fault, as determined in the present study, is shown in Figure 8. The maximum rise time on the individual sub-fault is about 1.1 sec, which is about one-sixth of the total rupture duration (6.5 sec) of the earthquake rupturing process. Comparing distributions of the slip and rise time, we find that the larger slip amplitude corresponds to the larger rise time on each subfault. Figure 9 shows the relationship of the linear regressions between slip amplitudes and rise time. This relationship is quite similar to the result of Chi-Chi earthquake [20]. The comparison of the retrieved ASTFs calculated by the hybrid blind deconvolution (solid line) and the synthetic ASTFs (dashed line) using the inverted kinematic source model is displayed in Figure 2. The correlation coefficients between the synthetic and calculated ASFTs are high and at approximately 0.8~0.98.

Stochastic Simulation
For identifying the validation of the inverted source model, we use the stochastic method to simulate the strong ground motions near the epicenter. There are 12 vertical component accelerograms of different stations simulated. The distribution of 12 stations is shown in Figure 10. The 12 stations at epicentral distances range from 3.28 to 29 km. In above simulation, we assume that the stations are all installed on the rock to detect the source effect on the ground motions near the epicenter. The complexity of the faulting process largely dominates the character of the signal at near source distances, when the site effects can be considered as weak or negligible (Emolo and Zollo [32]). The modeling parameters for the application of the simulated method are listed in Table 1.
The material properties are described by density, ρ , and shear-wave velocity, β , which are given the values 2.8g/cm 3 and 3.7km/sec, respectively. We applied a geometric spreading operator R 1 , and the anelastic attenuation was represented as . The kappa operator, κ , was given the value 0.006 sec for very hard-rock site (Boore and Joyner [33]) in the present study. The goal of this study is to investigate the source effects on the frequency and time domain of the seismic signals; therefore, we assume the path and the local site effects appear as the simplest forms. In Figure 11~13 we present the results of the stochastic simulations and their comparisons with the observed strong ground-motion recordings. Peak ground acceleration are well reproduced, expect the stations CHY030 and CHY034. CHY030 is under-estimated about a factor of 0.75 and CHY030 is over-estimated about a factor about 1.45. The amplitudes of the observed spectra of the stations (CHY009, CHY034, CHY038, CHY039, CHY095) which the epicenters are less than 13km are generally matched well with the observed, especially the station CHY009. It demonstrates that the inverted source model is quite reliable in this study. The phenomenon may be due to the fact that the seismic signal of the station near the epicenter is so largely affected by the source that the other effects of path and site are not crucial. Nevertheless, large discrepancies exist in the lower frequency at some stations (CHY028, CHY029, CHY032, CHY050 and CHY081). The epicentral distances of these stations are ranging from 15km to 28.24km, depending on the source-station distance. The reasons for this misfit may be attributed to inadequacies in the representation of the path or the local site.

Conclusions
We apply the hybrid blind deocnvolution method to derive the temporal-spatial rupture process of the M L 6.4 earthquake which occurred on the Chia-Yi, Taiwan, on October 22, 1999. In this study, the total rupture duration of the main shock is about 6.5 sec. Based on the relationship between the pulse widths of the ASFTs and the azimuths of the stations, we can distinguish the ruptured plane (SE-NW direction) from the focal mechanism. The slip distribution on the fault plane concentrates on two patches. From the higher rupture velocity with the lower slip amplitude, the relative weak region is found on 8-12 km in strike direction of the fault. For the most part, the larger slip amplitudes correspond to larger rise times on the subfaults.
The stochastic method is employed to confirm the reliability of the inverted source model. The result shows that the comparison between simulated and observed of the signal in the station near the epicenter is quite matched. 1226 The Kinematic Source Process of the MW 5.9, 1999 Chia-Yi Taiwan Earthquake from Teleseismic Data Using the Hybrid Blind Deconvolution Method It reveals the precision of the source model and the source effect dominates the seismic signal near the epicenter. The effects of path and site are weighting in the signals of the stations with the epicentral distances increasing. The next steps are applied these simulations to the earthquake disaster preventions and the field in earthquake engineering. There, the results in the paper are playing a crucial role in actual applications.