A Study on Source Rupture of the 1978 Lan-Hsu , Southeastern Taiwan Earthquake

, Seven teleseismic long period P-wave records from World-Wide Standard Seismo­ graph Network (WWSSN) are used to infer the source rupture process of the July 23, 1978 Lanhsu, eastern Taiwan earthqua k e through the body-wave iterative deconvolution method (Kikuchi and Kanamori, 1982; 1986). Results show that the rupture sequence consists of one main event and one largest subevent near the hypocenter determined by Pezzopane and Wesnousky (1989) and nwnerous smaller ones at a distance from the hypocenter. Comparison of the spatial distribution of mainshock rupture sequence with aftershocks shows that aftershocks are located in the area which did not break before the occurrence of the above-mentioned events. The total seismic moment of the three largest events is 8.76 x 1026 dyne-cm.


INTRODUCTION 1
Taiwan is located at the convergent boundary of the Eurasian and Philippine Sea plates.It is proposed that the Eurasian plate subducts southeastward beneath the Philippine Sea plate at about 21 S N along the Manila Trench, whereas the Philippine Sea plate subducts northeastward beneath the Eurasian plate along the Ryukyu trench at 23.5" N (Tsai, 1986).
For the past 4 million years, the collision between the Luzon island arc and the Chinese continental margin of the northern Eurasian plate has resulted in folded thrust sheets in central Taiwan and volcanic arc in the eastern coast of Taiwan.The Longitudinal Valley fault is assumed to be the suture of the arc-continent collision (Barri er and Angelier, 1986).
From the data of the Taiwan Telemetered Seismographic Network (TTSN), it is located at 22°21.11'N and 121°19.72'E (see Figure la) and its focal depth is 6.1 km.The hypocenter determined by Lee and Tsai is named as LT hypocenter in this work.Lee and Tsai (1981) reported that about one and a half years before the mainshock, microearthquake activity remarkably increased in number and clustered in the space toward the ensuing focal zone.
and some larger aftershocks have strike-slip fault focal mechanisms.Giardini et al. (1985)   determined the source mechanism of this earthquake through the centroid-moment-tensor method.Pezzopane and Wesnousky (1 989) used generalized ray theory to invert the source rupture process.The hypocenter detennined by Pezzopane and Wesnousky is named the PW hypocenter in this work.The focal mechanisms and source parameters reported by these authors are listed in Table 1 and   21-50 L--'--l---L�'--...L--'--'---'-�"'---'--'--'   Sea planes, it is significant to invert the distribution of mechanical prorerties in the source area to infer its relevance to regional tectonics.In this study, a body-wave iterative deconvolution method will be used to invert the source space-time function.

METHOD
The inversion technique of body-wave iterative deconvolution method used in this study was developed by Kikuchi andKanamori (1982 and1986) and Kikuchi andFuk1,to (1985 and1987).(1) where J0 (t) is the observed seismogram at the j-th station, the time t is measured from the initial arrival at each station, and Mis the total number of stations.In the next step, the observed waveform in Eq. ( 1) is replaced by the residual waveform: The point source parameters (m2, t2, Xi. y2) are detemtined by Eq. (2) through the same procedure as above.Thus, after N iterations, we obtain N point source parameters (m;, ti, x1, y), i = 1, 2, ... , N. The resultant synthetic wavefoim is given by: N lj(t) = L.m1 W(t-t; ;x;.yJ i=I The moment rate function M0 ( t ) and total seismic moment M0 are given as: and ( 5 ) where u; (t) is the source time function of an individual source with a unit seismic moment.
1brough the above procedure, the temporal-spatial distribution of each point source is iteratively determined.

DATA
The vertical-component records of teleseismic long-period P-waves at seven WWSSN stations have been analyzed.The epicentral distances are ranged from 30° to 90° to avoid the complicated phases from the crust and upper mantle and the diffracted phases from the core-mantle boundary.The related parameters for the seven stations are listed in Table 2.
The seismic signal has been filtered by instrument and aff ected by the local geological structure near station, especially for the later-arr iving phases.Hence, only few motions of the record have been taken into account for the reliability of data.A far field waveform can be represented as:  (1985) might not be an appropriate initial model for the inversion of source ruptures from the short-period data by the Kikuchi and Kanamori's technique.The values of differences between the observed amplitudes and the two calculated ones from the Pezzopane and Wesnousky's as well as the Lee and Tsai's parameters are in the range of 0. 15 to 0.40.
For practical computation, it is important to take an adequate segment of seismograms for each station to cover all of the P, pP and sP phases generated from all of the point sources during the source rupture but it should include the other reflected and refracted phases.Unfortunately, there is not any definite way to determine source rupture duration before the inversion procedure.The only approximate way is to estimate it from earthquake magnitude.We, therefore, proposed that larger earthquakes should have longer duration time.As some previous studies suggested (Kikuchi and Fukao 1985;1987), the main subevents are expected to happen within 20 seconds after the initial rupture for an earth quake with seismic moment in order of 10 26 dyne-cm.After further consideration of the  travel time of the sP phase, a 60-second rupture duration was selected.
The difference ' in the arrival times between the second phase ( pP phase) and the first phase ( P phase) and that between the third phase ( sP phase) and the first phase ( P phase) were measured.The former is• denoted by 11T;,pp and the latter by 11T;; pp.  is not used.The two initial fault-plane models used are determined by Lee and Tsai (1981) and Pezzopane and Wesnousky (1989).Their geometric structures are shown in Figure 4. ,; The area of the fault planes are inferred from the aftershocks.Totally, 18 x 12 and 12 x 12 grid points with an equal spacing of 5 hn are defined for the two models to include all aftershocks under distinctive fault plane solutions.Time history of individual subevent is a .
trapezoid time function with two time constants: rise time and rupture duration, which are determined by Pezzopane and Wesnousky (1989).
For far-field wavefonn synthesis, the geological structure is simplified to a half-space.
Tue rupture velocity is taked as 3.4 hn I sec, which is about 0.9 times of half-space S-wave velocity of the source area.

RESULTS
Tue Kikuchi and Kanamori 's inversion method can be used for both single-station data and multi-station data.In this work, we study the source mechanism as inverted first through the single-station inversion, then through the multi-station inversion.

Single-station inversion
Tue iteration errors for the seven stations through the single-station inversion based on   o-RAR  The subevents of the group at depth are, in general, larger than those in the shallow zone.
The seismic moments for subevents inverted after 20 iterations are listed in Table 3.
The seismic moments obtained from the data of Pacific stations are almost half of those obtained from the data of Eurasian stations, and 3.5 times of that reported by Pezzopane and Wesnousky (1989) and 1.8 times of that by Giardini et al. (1985),

Multi-station inversion
In the single-station inversion, the distinction between the results from the data of Pacific stations and Eurasian stations is obvious.Therefore, we have performed inversions for these two individual data sets and finally a joint inversion for all of these data.Among the data, as mentioned before, the fitting of synthetic seismograms with observed ones at NDI aren't good, and we have therefore excluded the data of this station for the multi station inversion.The iteration errors for the multi-station inversions are shown in Figure  Results show that some large multi-rupture events appear at the boundary grids.Since the area of each grid is 25 km 2 , it is impossible to accumulate enough energy for multi ruptures in such a small area within a short period of time.So, we have extended the fault plane to a larger size, such as 100 x 60 km 2 , then resumed the inversion.The result is shown in Figure 9, in which the prescribed boundary events remarkably move along .thedirection of boundary extension.It indicates an existing error in the inversion of spatial distribution caused by the requirement of matching the synthetic waveform with the observed one.This error will be reduced when the fault plane is extended, but an unreasonable inversion results.So the inversion results for the boundary events are less reliable.
Figure 10 shows the variations of inverted seismic moment with the number of itera tions for Pacific and Eurasian data sets.For Pacific data set, the seismic moment almost monotonically decreases with iteration.A similar variation also exists for Eurasian data set except for the appearance of an abnormal value at the sixth iteration.Before the seventeenth iteration, the values of seismic moment from Pacific data set are smaller to those from Eurasian data set, but vice verse after the seventeenth iteration.However, after the tenth iteration, the values of seismic moment from the two data sets are very similar.

DISCUSSION
The spatial pattern of subevents inverted from the data of Eurasian stations (Figure Sb) is more complicated than that from the data of Pacific stations (Figure Sa).One possible reason for this distinction is due to the different crustal structures of these two regions: a continental type for Eurasian stations and an oceanic type for Pacific stations.In general, the oceanic crust has a thickness of about 5 to 10 km, while the continental crust has a thickness of about 30 to SO km.The travel times of body waves would be different within the two crusts.Besides, the continental crust composition is more complex than the oceanic crust.The former consists of two main layers: the upper layer with granite and the lower layer with basalt; in contrast, the latter is essentially a layer of basalt.The shallow part of the continental crust is quite complicated and may be composed of nonflat structures.However, this shallow part is not thick and does not have a remarkable influence on the intermediate- .J.._ ._.. From Table 3, the resultant seismic moments after 20 iterations are larger than those reported by other authors (e.g.Pezzopane and Wesnousky, 1989;Giardini et al.,• 1985).In our inversion process, all the data recorded in the truncated seismogram are deemed to be generated directly from the source rupture.However, some refracted and reflected signals might exist.These signals would yield small subevents, thus increasing the heterogeneity According to the asperity concept (Kanamori and Stewart,1976), the main :fractured zone during the mainshock is regarded as a strong asperity with a breaking strength higher tll an a critical value.This zone breaks when the accumulated stress is greater than the critical state.Aftershock activity is related to the mainshock rupture through a certain mechanism, for instance stress corrosion as indicated by Yamashita and Knopoff (1987).They proposed two models to interpret aftershock activity.In the first model, aftershocks are assumed to be caused by subsequent slip on a8perities of a fault which are locked during the fracture in the mainshock.In the second one, aftershocks are caused by a catastrophic coalescence of small fractures with the fracture surface of the mainshock.Stress corrosion cracking is assumed to be the cause of aftershock activity for both models.In Figure 11, the first model is applied to interpret aftershock activity in the case where the small subevents inverted from the seismograms are also considered to be a part of the mainshock.On the other hand, the second model works as these small subevents are considered to be after shocks which immediately occurred after the mainshock or 'ghost' events through the inversion process due to the contamination of reflected and/or refracted waves.Mendoza and Hartzell (1988) also mentioned that the spatial distribution of aftershocks reflects either a continuation of slip in the outer regions of the areas of maximum coseismic displacement or the activation of subsidiary faults within the volume surrounding the boundaries of the mainshock rupture.
shown in Figure lb.

Fig. 1 .
Fig. 1.(a) The location of July 23, 1978 Lanhsu earthquake and its surrounding area.The 'solid star' denotes the epicenter.(b) Fault plane solutions of the Lanhsu earthquake: dash lines by Lee and Tsai (1981), solid lines by Pezzopane and Wesnousky (1989), and dotted lines by Giardini et al. (1985).
The source process is described as a sequence of point dislocations, each being characterized by seismic moment (mi), onset time {ti) and coordinates (x;, Y i) on the fault plane, where we take the x-axis along the fault strike and the y-axis_ along the dip.The considered fault plane is divided into nx x n y grids, and adjacent two grid points are sepa rated by L\x in both x-and y-directions.The fault mechanism and source time history are assumed to be identical for all point sources.The synthetic wavelet (i.e. Green function) at the j-th station, generated from .a unit point source (x i, Y i) at time ti, is denoted by "'.; (ti; xi, Yi), and the synthetic record generated from a point source (m;.t1, xi, y;) is represented as mi HJ (t; ; x;, Y;).In the first iteration, the point source (m1, t1, x1, y1) is determined through the minimization of the error: 6) where p, c, S(t) are density, • velocity and far-field source time function respectively.The first-motion amplitude is proportional to ( R ef r0), where Re is the radiation pattern in terms of fault parameters, take-off angle of the seismic ray and azimuth angle of the station (Aki and Richards, 1980).The observed first-motion amplitudes and calculated values from the three fault plane solutions determined by Pezzopane and Wesnousky (1989), Giardini et al. (1985) and Lee and Tsai (1981) are shown in Figure 2. It can be found that the variation of amplitude with distance of the observed data is similar to those of calculated data under the parameters b:y Pezzopane and Wesnousky (1989) and Lee and Tsai (1981) in spite of the match of the values of amplitude, but does differ from that under the parameters by Giardini et al. (1985).The fault plane solution, determined by Giardini et al. from the seismograms of periods from 45 to 135 sec., shows the long-period source mechanism.Thus the fault plane solution determined by Giardini et al.

Fig. 2 .
Fig. 2. Observed and calculated first-motion amplitudes nonnalized with respect to that at station MUN: dotted line for observed data and solid lines for calculated values (square, cross and triangle symbols denote the values based on the focal parameters by Lee and Tsai (1981), Pezzopane and Wesnousky (1989), and Giaridini et al. (1985), respectively).

Fig. 3 .
Fig. 3. Differences of arr ival times of the second and third wavelets with respect to that of the first wavelets observed at seven stations NDI, MUN, RIV, KEV, NUR, AFI and RAR (from left to right).
Before iteration, the initial model must be first constructed.Since the dominant periods of the seismograms used ranged from 15 to 27 sec., the fault-plane solution obtained byGiardini et al. (1985), which is based on the seismograms with periods greater than 40 sec., Fig. 4. Geometry of models: (a) based on Lee and Tsai's source parameters and (b) based on Pezzopane and Wesnousky's source parameters.

Fig. 5 .DISTANCE
Fig. 5. Variations of error with iteration: (a) for seven stations through single-station inversion based on Lee and Tsai's source param eters; (b) for seven stations through single-station inversion based on Pezzopane and Wesnousky's source parameters; (c) for Pacific stations, Eurasian stations and all stations through multi-station inversion based on Pezzopane and Wesnousky's source parameters.

Fig. 8 .
Fig. 8. Temporal-spatial distributions of subevents based on Pezzopane and Wesnousky's model through multi-station inversion: (a) for the data of four Pacific stations; (b) for the data of two Eurasian stations; (c) for six data described by (a) and (b).
5c.The errors for the two individual sets are approximately constant (about 0.15) after 5 iterations, but the error for the joint data decreases with the number of iteration and is about 0.25 after 20 iterations.Inverted results are shown in Figure8.The spatial pattern of subevents from the data of Pacific stations is similar to those from each data set of these stations by single-station inversion.Both results show that a main subevent occurs first after 2.0 sec.followed by five smaller subevents.Two larger subevents are located near the PW hypocenter, and smaller ones are distributed over the fault plane.The spatial pattern in verted from the data set of Eurasian stations shows that three subevents occur at 15.0, 3S.O and 53.0 sec.Comparison of Figure7and FigureSbshows that the source time function in the multi-station inversion is an average superposition of those in single-station inversions of the data of KEV and NUR.This phenomenon is particularly remarkable in the case of joint data inversion (i.e.FigureSc).

Fig. 9 .
Fig. 9. Temporal-spatial distributions of subevents through multi-station inversion for extended source area with dimension of 100 x 60 km2 from Pezzopane and Wesnousky's model: (a) for four Pacific stations and (b) for two Eurasian stations.

Fig. 10 .
Fig. 10.Variations of seismic moment in terms of number of iteration: 'solid line with star' for Pacific data set and 'dotted line with square' for Eurasian data set.

Fig. 11 .
Fig. 11.Spatial distribution of mainshock rupture events (in shaded circles) and aftershocks (in integers to represent their duration magnitudes).

Table 2 .
Station parameters