Near-source Ground Motion of a Propagating Rupture Fault From the Finite Element Modeling

A finite element algorithm was developed to simulate wave propagation and near-source ground motion of propagating rupture fault. The modified 'split-node' technique of Melosh and Raefsky was employed to simulate the equivalent body forces system for a double couple without moment. Synthetic seismograms were computed for a suite of dip-slip faultings with various slip directions, and rupture velocities to illustrate the ground motion and coseismic deformation variations. Results from this study showed that, for a thrust fa ult with infinite rupture velocity or down going propagating ruptures, the hanging wall has a stronger ground motion than the foot wall; in contrast, for up going or bi-lateral ruptures, the hanging wall has a weaker shacking than foot wall. For the propagating rupture faults, the 'starting' and 'stopping' phases were numerically reproduced using the finite element method. In addition, observation at the surface, the duration pulse of S-waves were significantly changed between both sides of the fault for the rupture directivity. (


INTRODUCTION
The study of earthquake sources is one of the prime concerns for seismologist.The bulk of knowledge about earthquake sources comes from research through practical observation and theoretical investigations.Theoretical study of the earthquake source rupture process was previously published by Haskell (1969) to compute the ground motion of finite ten sile fault.Boore and Zoback (1974) computed the near-source motion of the inplane finite rupture, while Madariaga (1976) computed the dynamic.rupture of circular fault using the finite difference method.Previous theoretical near-source ground motion modeling has typ ically been restricted to full-space (Luco and Sotiropoulas, 1985) and constant propagation rupture (Bouchon, 1980) models.The application of these methods did, indeed, success fully interpret records from the 1966 Parkfield earthquake (Aki, 1968) and the 1979 Imperial Valley earthquake (Archuleta, 1984).In addition, many studies related to source rupture were proposed (Goldstein and Archuleta, 1991;Chiu, 1991), but their results showed that the source rupture process and fault geometry are more complex than those from teleseismic observation.They demonstrated the fact that the detailed source rupture processes may be retrieved from the near-source strong motion simulation (Hartzell and Heaton, 1983;Wald et aL, 1991).However, the resolutions from those snidies were limited due to combined effects from the laterally varying velocity structure, topographic effects, attenuation and anisotropy.Other than this, detailed knowledge about near-source ground motion behavior is still not well defi ne.Recently, the installation of digital strong motion and broad-band instruments have acquired a great deal of high quality seismic data (Abrahamson et al., 1987; CDMG,  1989; Kanamori et al., 1991; Helmberger et aL, 1993), thereby creating new opportunities to improve the knowledge of source rupture processes.Nevertheless, to extract the source information from near-source observation� accurate wave-form modeling is necessary.
The fi nite element method, which is well known as a robust numerical method for synthetic seismogram simulation, is one of the modeling techniques for the study of the near-source ground motion behavior because of it's natural ability to simulate both complex elastic models and complex sources.Many complex structure can be easily modeled using the finite element method without difficult modifications to the algorithm (Zienkiewicz, 1977).Moreover, the method can be easily employed to approximate the topography.In this study, the authors presented the finite element method to model the ground motion from the fault rupture.Also, the 'split-node' technique (Melosh and Raefsky, 1981) was employed to simulate a complex dislocation rupture source.Although it is not presented in this study, however, this approach lent itself as a possible tool for the complete wave form modeling of complex source geometry and slip distribution in a laterally varying velocity structure.In this study, the near-source ground motion computed from the numerical modeling was analyzed.The results of this study showed that rupture velocity changes have •strong effect on the near field ground motions.

Finite Element Formulation
In the finite element method, wave motion due to external loading force can be described as a set of simultaneous, second-order differential equations results: ( 1 ) where [M] is the global mass matrix, [C] is the damping matrix, and [K] is the global stiffness matrix.[A], [VJ and [D] are the vectors of the global acceleration, velocity and displacement of the model at grid points.[F] is the vector of the external loading forces at where 6.t is the time interval used in numerical integration.kesults from the finite element computation are the full solution of propagating waves; thus, all transmitted and reflected waves, converted waves, surface waves, diffraction, head waves and multiple scattered waves are included.
2.2 Fault Dislocation Formulation Melosh and Raefsky (1981) proved that the static rupture fault may be implemented as an equivalent double-couple force distribution without moment.As a result, they developed a split-node technique to represent the fault dislocation (equation ( 7) in Melosh and Raefsky, 1981): ( 3 ) In this study, the 'split-node' technique was extended to simulate the dynamic dislocation rupture processes via solution of the elastic wave equation.In this case, the static forces [ F] used in equation (3) were replaced by the dynamic forces [F(t)] which were represented as ramp functions (Figure 1) to simulate the rupture propagation by suitable time delays (Yeh and Huang, 1986;Huang, 1989).Similarly, the static deformation [ D ]   cases of this study.

Numerical Implementation
The 2-dimensional, 4-node quadrilateral element was used in this study to modeling wave propagation.For handling the fault plane which slips across the elements, the local assemblage technique to redistribute the equivalent forces to grids was used (Zienkiewicz, 1977).To keep numerical dispersion down to an acceptable level, the grid interval which was less than one over ten wavelengths of the dominant frequency of the S wave was selected.
A damping matrix [ C] was used with a damping coefficient as defined by Archuleta and Frazier (1978) in the 0.1 to 0.2 range.The increasing of the damping coefficient will reduce the ringing with calculated signals, however, it also reduces the signal amplitudes.There is not a suitable theory to select the damping coefficient, the authors suggested a trial and error test in each case.In order to suppress numerical reflections from the artificial boundaries of the fi nite element model, the boundary absorbing method proposed by Smith (1974) was used.This method involved the superposition of independently calculated wave equation solutions with Neumann and Dirichlet conditions.Although the method proposed by Smith (1974) can not cancel the multiple ref l ections between free surface and artificial boundaries, however, primary and comer reflections may completely cancel using this method (see Figure 1 of Huang and Yeh, 1991).A similar procedure for source rupture processes was proposed by Benz and Smith (1987) to model the source rupture.For the purpose of comparing the accuracy of both algorithms, a case proposed by Benz and Smith (1987) was reproduced and is shown in Figure 2. Nearly the same amplitudes and wave forms are present in both profiles.The difference in artificial boundary refl ection, noted as BR in Figure 2, may be due to different numerical approximations.Otherwise, the code in this study was checked against with the analytic solution (Huang, 1992) and the pseudo-spectral method (Huang and Yeh, 1991).This study demonstrated two numerical experiments to simulate near-source ground motion.First, we computed the near-source ground motion features for the model size near 14 times that of the fault dimension (Experiment 1).In this case, the fault was simulated by 12 grid points.The second model on which the rupture fault was discreted by 90 grid points.In this experiment, the dimension of the model is about 2 times larger than the fault dimension.The detailed properties of the rupture directivity and rupture complexity were examined,

3.1Experiment 1
The goal of this experiment was to analyze the amplitude characteristics of ground motion from finite source with a propagating rupture process.For this purpose.a model of 45° dip-slip thrust fault on half-space was used and was shown in Figure 3.The ground motion recorded on the free surface (Figure 4) and wavefield snapshots (Figure 5) were examined.The defi nition of the positive direction of seismograms in Figure 4   A model with the same fault geometry as in Figure 3 but as a bi-lateral propagating rupture is proposed to demonstrate the amplitude characteristics of ground motion.A source time function (Figure 8) was chosen to compute the synthetic seismograms to demonstrate the difference between the bi-lateral and the infinite rupture sources (Figure 4).In this case, the rupture was initiated at the center of fault plane then propagated to two ends.Although the fault type (thrust fault) is the same as.the infinite rupture case, however, the free surface ground motion of the propagating rupture has a complex time history than a infinite rupture velocity case (Figure 4) as shown in Figure 9.In the bi-lateral rupture case, the shacking on the foot wall have larger values than that on the hanging wall.The frequency content of the horizontal acceleration seismograms on foot wall is higher than that on the hanging wall.The detailed seismograms recorded at location 0.57 km for different rupture process are shown in Figure 10.The seismogram of the bi-lateral propagating fault on this point shows a complex wave form but a small peak amplitude than that of the infinite rupture velocity case.

Experiment 2
For a detailed study of the wave form features near and above the rupture fault, a thrust slip model of a large fault dimension, a model as shown in Figure 11 was proposed.Figure 12 shows the acceleration seismograms of a dip-slip faulting with constant rupture velocity.This case presents a upward propagating rupture from the lower end tip (point A) to the upper end tip (point B) with a rupture velocity (Vr) equal to 0.8 S-wave velocity (Vs).The solid lines on Figure 12  However, phases Pl and Sl were clearly separated from P2, S2 at RP.For the pure dip-slip faulting of this case, the vertical acceleration on point LP was strongly enhanced.Although the peak amplitude ratio of LP/RP on a vertical displacement is 3, but, the ratio for vertical acceleration is up to 14.However, the amplitude enhanced on the horizontal component is not as clear as on the vertical one with horizontal displacement ratio as 2.5 and horizontal acceleration ratio as 5.1.
Figure 14 shows the seismic profiles of a pure thrust fault with downward propagated rupture which has the same fault geometry as in Figure 11 but the rupture initiated from point B with downward propagation (Vr= 0.8 Vs) to point A. In this case, the directivity is not as clear as the upward propagation (Figure 12).From the acceleration profiles, the arriving phases P2 and S2 can be clearly separated from .Pl and Sl on both side of fault.However, larger amplitudes and shorter duration were found on the hanging wall of the fault.
Selected seismograms at point MP of Figure 11 for both rupture faulting types are shown in Figure 15.The ground motion amplitude of upward propagated faulting is larger than that from the downward propagated faulting.The total duration (from Pl to $2) of upward propagation faulting is shorter than that of the downward propagation faulting.In    both cases, the arrival phases of P and S waves radiated within the fault segment (between point B and A) cannot be identified from those seismograms.The surface ground motions represented the low frequency features during the rupture front passing.
To demonstrate the difference in the recorded ground motions by changing the rupture velocity for the same fault geometry and propagating direction, in this study, we• present two cases.Both cases are the upward propagating ruptures from point A of Figure 11 to point B. The first one, which is noted as •v-f• case in Figure 16 and 17, increases the rupture velocities from 0.8 to 1.12 S wave velocity, uniformly.The other one, which is noted as •v-• case in Figure 16 and 17, decreases rupture velocity from 0.8 to 0.46 S wave velocity, uniformly.The 'Vr' case, which has the constant rupture velocity (0.8 S wave velocity) was discussed in Figure 12 and 13.We selected two observation point at L1 and Rl (Figure 11).

305
Ll indicates the point of 2.5 km to the left of LP.Rl indicates the• point of 2.5 km to the right of RP. Figure 16 and 17 shows the seismograms recorded at LI and Rl for the •v + ', •v-• and 'Yr' cases.Figure 16 shaw that in the beginning of the seismograms are similar among the three cases where as the later-phases have significantly changed.The ground motion was strongly amplified for the increasing rupture velocity case and deamplified for the decreasing rupture velocity case.The duration of seismograms for the v + case is shorter than vr or v-case.For the v + case, the P wave may be completely separated from the S wave on the acceleration seismograms, because P2 recorded early than Sl.Although i. t is not a point source, the increasing rupture velocity case shows a very simple acceleration seismograms.
The ground motion on the hanging wall, for example point R 1 (Figure 17), shows a smaller amplitude value, for all of the attribute, than that on L1 (Figure 16) for the same case.At Rl, the peak amplitudes of seismograms are similar for the rupture velocity change.The ;teceleration seismograms show clear P-arrival at R 1 (Figure 17) than that at point L1 (Figure 16).

DISCUSSION
The near source ground acceleration from the point dislocation source was found as a wave let type motion (Johnson, 1974).In this study, it was found that the ground motion from the finite instantaneous rupture fault was still similar to that from the point source as shown in Figure 4 and Figure 10 (b).However, the wave form was relativity complex for the bi-lateral rupture faulting, as shown in Figure 10 (a), and the unilateral rupture faulting as shown in Figure 13, 15, 16 and 17.Actually, the rupture effect in the source spectra has become known as directivity (Aki and Richards, 1980;Bullen and Bolt, 1985;Douglas et aL, 1988).Observations showed that the amplitudes and durauon of the ground motion were different for the station on different azimuths of rupture sources (Hartzell, 1989;Ammon et aL, 1992).However, most theoretical discussions for the directivity were only for frequency domain features (Haskell, 1964;Hasegawa, 1974) and the fault types were limited to the static state or uniform and/or constant velocity rupture faulting (Boore et al., 1971;Luco and Sotiropoulos, 1985).Using the method proposed in this study, the directivity phenom ena could be reproduced in time domain.Hence, some limitations of source approximation were released using numerical modeling.For example, the directivity effects from complex fault ruptures such as slip directions, rupture velocities and complex slip offset could then be studied.The major objective• of this study was to develop a modeling technique to the oretically interpret suitable wave form features observed in the recorded seismograms� In this study, the characteristics of the near-source ground motion resulting from the simple two-dimensional dislocation models were examined.The wave fonn complexity has been effectively analyzed from the forward modeling using the finite element method.Many phases can be traced from the seismic profiles and snapshots via numerical computation.The application of the computer program which was developed by this study was successful in the case the May 20, 1985 Hualien earthquake.Wu (1990) and Tai (1990) employed this method to model the strong motion records of mainshock and aftershocks at station Hualien.
The comparison of the calculated and observed motions showed that the characteristics of the fault rupture such as slip amount, direction and rupture propagation direction, velocity and rise time were examined from the wave fonn mod eling (Wu, 1990;Tai, 1990).
In the last few years, although a number of near-source earthquake ground motion records have been obtained, the recorded data are not still sufficient to obtain the relative ground motion in _ the near-source region.Actually, it is hard to analyze the characteristics of complex fault ruptures such as the change of slip velocity, amplitude and rupture directions from a few recorded seismograms.However, as the quality and quantity of seismic recordings increase, seismologists will be able to improve their ability to identify and interpret complex wave forms as functions of structure and source properties from recorded data.The theoret ical development and numerical simulation of the near-source ground motion properties are necessary for better understanding of the characteristics of the near-source earthquake ground motion.
Finally, the authors wish to point out that although there are many potential applications, the split node technique, which is employed in this study, is still limited by its infinite fault width assumption.Thus, only the dip-slip fault can be simulated using the 2-dimensional elastic wave equation and most records observed are not dip-slip fault.The application of this method to the modeling of the ground motion is therefore limited.Some extensions are necessary.

CONCLUSIONS
A method to compute the near field ground motion of a moving dislocation rupture was presented.This method has the ability to model complex rupture directivity.The near source ground motion features from the complex rupture model such as directivity, starting phase, stopping phase and rupture passing phase were successfully simulated.Results showed the near-source ground motion was strongly affected by propagating rupture process.In short, with variable rupture velocities.Identification symbols remain the same as in Figure 16.
for a systematic study of the source rupture process, the finite element method is a useful tool.

Fig. 1 .
Fig. 1.The ramp type displacement source time function which is used in all cases of this study.

Fig. 2 .
Fig. 2. Comparison of seismic profiles computed by Benz and Smith (1987) (a) and this study (b).The dashed lines noted as BR are the arrival times of the boundary reflection phases, Pl and SI are arrival times of the 'starting phases' of P and S waves, while P2 and S2 are th� arri val times of the 'stopping phases'.
is up for the vertical component and right for the horizontal component.The definition remains the same in the other figures of this study.In Figure 4, the vertical acceleration profile has simple wave forms.except for the small range near the fault tip on which the surface waves were interfered by the propagating waves at free surface.The ground motion on the hanging wall have stronger shacking and more complex wave forms than that on the foot wall.Most of the energy which propagated to the left was downward and away from the free surface.Major contribution of energy are from scattering from fault tip.The integrated vertical velocity and displacement seismograms of Figure 4 �ave similar characteristics as shown in Figure 6 except lower frequencies.The vertical displacement on the foot wall was subsidence while on the hanging wall it was lift.For the displacement ramp type source time function used in this study to modeling the dislocation rupture, the permanent deformation which show the final offset of the fault rupture can be measured in each end of the seismograms of the displacement profile (Figure 6(b)).Along the profile on the free surface, the permanent deformation is shown in Figure (7).

Fig. 3 .
Fig. 3.The half-space model used in Experiment 1 of this study.The depth of the upper fault tip is 120 m: the dip angle is 45°.Only the central part, thus AA'. with a range of 1.2 km was selected in the future display.

Fig. 4 .
Fig. 4. Acceleration profi les of the ground surface of Experiment 1.The fault geometry is shown in Figure 3.

Fig. 5 .
Fig. 5. Acceleration snapshots of the model of Figure 3 with an infinite rupture velocity.The elastic wavefields are expressed in a vector format.directivity corresponding to large amplitudes and more coherent signal can be clearly seen on the seismograms at the left hand side of the fault tip as shown in Figure 12.The free surface ground motions at locations which are indicated as LP, MP and RP on Figure 11 are shown in Figure 13.From the acceleration seismograms, phases P2, Sl and S2 arrived at LP nearly the same time and produced a single pulse with a peak amplitude larger than RP and MP.However, phases Pl and Sl were clearly separated from P2, S2 at RP.For the pure dip-slip faulting of this case, the vertical acceleration on point LP was strongly enhanced.Although the peak amplitude ratio of LP/RP on a vertical displacement is 3, but, the ratio for vertical acceleration is up to 14.However, the amplitude enhanced on the horizontal component is not as clear as on the vertical one with horizontal displacement ratio as 2.5 and horizontal acceleration ratio as 5.1.Figure14shows the seismic profiles of a pure thrust fault with downward propagated rupture which has the same fault geometry as in Figure11but the rupture initiated from point B with downward propagation (Vr= 0.8 Vs) to point A. In this case, the directivity is not as clear as the upward propagation (Figure12).From the acceleration profiles, the arriving phases P2 and S2 can be clearly separated from .Pl and Sl on both side of fault.However, larger amplitudes and shorter duration were found on the hanging wall of the fault.Selected seismograms at point MP of Figure11for both rupture faulting types are shown in Figure15.The ground motion amplitude of upward propagated faulting is larger than that from the downward propagated faulting.The total duration (from Pl to $2) of upward propagation faulting is shorter than that of the downward propagation faulting.In The vertical seismic profiles of (a).velocity and (b).displac ement of the same model of Figure3.The profile observe d range was the same as in Figure4.

Fig. 7 .
Fig. 7.The displac ement profile of the fi nal offset of the model of Figure 3 with infinite rupture velocit y.HD represe nted the horizon tal displac ement, while VD the vertica l.Only on the range of 0.3 km on Figure 3, thus BB', which centere d on the upper fault tip was selecte d in this profile .
Fig. 12. Horizontal and vertical component synthetic acceleration seismograms for the fault model of Experiment 2 (Figure 11) with a uniform source rupture velocity (0.8 S-wave velocity) from the lower end (point A) upward prop agating to the upper end (point B).P1, S1 are the P-wave and S-waves which radiated from point A; P2, S2 are the P and S-waves from point B.

Fig. 13 .
Fig. 13.Horizontal and vertical component synthetic acceleration, velocity and displacement seismograms on points LP, MP and RP of Figure 11.Iden tification symbols remain the same as in Figure 12.The number at the right end of each seismogram represents the peak value.

Fig. 14 .
Fig. 14.Horizontal and vertical componen_ t synthetic acceleration seismograrns for a thrust fault with a uniform rupture velocity (0.8 S-wave velocity) from• point B of Figure 11 downward propagating to point A. Pi.S1 are the P-wave and S-wave which radiated-from point.B; P2, S2 are the P and S waves from point A. Identification symbols remain the same as in Figure 12.
Comparison of synthetic acceleration, velocity and displacement seismo grams from the fault model (Figure 11) on point MP with a uniform upward (UW) and downward (DW) propagating rupture.Both rupture models had the same rupture velocity (0.8 S-wave velocity).Identifica tion symbols remain the same as in Figure 12.