SH Wave Seismogram Synthesis by the Finite Element Method

A method for the finite element computation of the s)·nthetic seismo­ gran1s of SH wave propagation in a two-dimensional (2-D) medium is pro­ posed. The finite element method has the ad,rantage of specifying the arbi­ trar)' density and velocity field in the medium. To complete the waveform modeling for a n earthquake, the transformation of line source snapshots and seismograms of the 2-D finite element modeling to the point source responses is critically necessary. The accuracy of the finite element model­ ing is verified b)1 comparing it with the generalized ray method and \\1ith a finite difference method which computes the 2-D wave field from a differ­ ent approach. The aim of this study is to demonstrate the successfulness of seismogram synthesis by the finite element method and to propose the potential applications of this method in the Taiwan area. (Ke)· words: SH wave, Finite element method, Seismogram synthesis, 2-D modeling)


INTRODUCTION
Wav'e propagation problems in co1nplicated geological structt1res have generated consid erable interests in rec.ent )'ears. The full \vave equation modeling with numerical methods such as the finite dit ' ferences (FD), the t111ite ele1nent (FE) and the pseudo-spect14 um (PS) methods is a suitable approach. H<)\\lev'er, ffi(1deling wave propagation with three-di1nensi()nal (3-D) struc tures is, in fact, difficult to handle with present-day1 computers. These methods require a tre-1nend()US amount ()fcc)rnputer time and also allo\\1 C)nl)1 a limited 11umber ot· \\1avelengths to propagate. F1· on1 this p<)int <-lf \1iew, 2-D modeling is more practical to use. The '1pplication of� 2-D ffi()deling has been pr()posed in many· problerns t)f· · seisn1olt)gical interest Lamer, 1970� Dravinski, 1983 ). Ho\\'ever, the Green's t � uncti()n of 2-D (line source) is actually differ ent to that of 3-D (poir1t S()Llrce), so applicatior1 is limited to earthquake \vaveform synthesis, such as the gene. ralized 14ay theo�y (GRT) method (Heaton and Helmberger, 1977), WKBJ method (Chapman, 1978) and the frequency-wavenumbe1· (FK) met.hod (Wang and Herrmann, 1980). Recent1y, Vi dale et al. ( 1985) ·proposed a new approach to seismic w·ave propagation mode . ling including laterally \'arying structures \vitl1 the 2-D FD scheme a11d applying it to wa\ief()fin synthesis. 258 TAO, Vol. 7, No. 3, Septem. ber 1996 . . They generated the seismic waves in source area by the GRT and computed the wave propagation by the 2-D · FD method. Finally, they calibrated the propagated energy from 2-D to 3-D by a suitable. transformation discussed below. This approach is success·ful in computing a 3-D wavefield for a 2-D lateral heterogene. ous model and is widely applied to ground motion simulation (Vidale et al., 1985;Ho-Liu and Helmberger, 1989). In this study, a new FE ap proach to model the 2-D Green's functions based on equivalent body force. si s proposed. The couples that we . re required to obtain equivalent forces for a generally oriented displacement discontinuity (Aki and Richards, 1980) were directly implemented. This simplified the nu merical computation for dislocation sources. The obtained Green's functions were easy to modify for the linear moment tensor inversions algorithms (Stump and Johnson, 1977;Kikuchi and Kanamori, 1991;Huang, 1994 ). Furthermore, a numerical technique of computing convo lution on grid points of the 2-D model is proposed to simulate the point source wavefield snapshots, which have not been discussed before. (1) c where G3 is the 3-D Gree . n's function, 8 (t) is the Dirac delta function and V 3 is the 3-D Laplacian ope. rator in Cartesian coordinates. The wave speed c in equation (1) may be an arbitrary function of spatial coordinates, i.e., c(x,y· ,z). For constant c, the solution of equation ( 1) is: The significant properties of G1 are: (1) its amplitude decays (geometric spreading) as where G2 is the 2-D Green's function and V 2 is the 2-D Laplacian operator. Physically, equa tion (3) describes a wavefield emanated from a point source in space with coordinates x and z, which are equivalent to a line source along the y-axis in space with coordinates x, y and z. The same conditions hold as in the 3-D case discussed above. Thus, the 2-D Green's function has the fallowing form: where H(t) is the Heaviside step f ' unction. The t· eatures of G2 are: (1) its an1plitude decays as l /( � r), (2) its propagating vv·avef ' orm is different from the source pulse and includes a low decay tail. The differences betv,1een the 2-D and 3-D Green's functio11s are shown in Figure 1.

Point Source Expressions for the 2-D Green's Function
Fol1ov\ling the deri\'atic)n ot· the GRT and being solved by· the Cagniard-de Hoop tech nique (Hel111berger, 1983 )., the displacement potential ot, the point SC)urce, </J1), in a wholespace can be expressed as: \V here ,,. and z are the r�ldi al and \'ertical coordinates, res pee ti ''el y, and .
where p is the ra)' par,tmeter, and 1] == (1 I (�2p 2). Im represents the image part of a com plex variable. Hc)we . ver., ff()l11 the \Vork of Gi1bert and Knopof ' f ( 1961 ), the solution for the line S(1urce, </> 1_, excitation can be expressed as: · 260 TAO, Vol. 7, No. 3, September 1996 Thus, the point source solution (</JP) represented from the line source solution J( t) in equation (6) as J=('1P) <PL can be determined. This approach can be used to simulate the wavefor1r1 of an earthquake using the 2-D numerical Green1s functions.

Finite Element Formulation
The 2-D wave equation, equation (3), can be formulated as the follo\\t·ing finite element equation for shear \\t·ave propagation [Fl is the vector of the external loading fo14ces at source points. To solve the temporal derivativ'es of e . quation (9) on each grid point, in this study, the second derivati · ves are introduced by the use of explicit finite difference approximation (Wylie, 1975), as follows: where the definitions of A and D are the same . as those in e-quation (9). Vis the vector of the global v·elocity, and � t is the time interval used in numerical integration. Thus, the displace ment Don time t + '1.t is computed from the previous computed values of D(t), V(. t -1/2 lit) and A(t). The finite element code development and mesh design are described in more detail in Huang (1989) and Huang and Yeh ( . 1994).

Source Representation for Shear Dislocation Fault
Following the definition of" Helmberger and Hark14ider (1978)� the tangential displace ment U5H can be expressed as: where TSS and TDS are. the vertical tangential strike-slip and the dip-slip displacement Green's functions, respectively. Both A4 and A5 can be expressed as: where () is the strike from the end of the fault, A is the rake angle, and 8 is the dip angle. To introduce the shear dislocation source to numerical mode-ling, Vidale et al. (1985) used the GRT to generate the displacements in a small ring around the source point. Those prescribed displacements drove the FD grids to simulate the Green's functions. In this stud)1, the concept of equivalent body forces (Aki and Richards, 1980) was employed. The direct implementa tion of the couple of body forces on the source grids made the modeling of the explosion and • shear dislocation sources successful, while simplifying the numerical computations.

Comparison With Other 1"1ethods
To test the accuracy L)f. this method, the results are compared \Vith those from other meth ods. Herein, a model \\1ith ()ne lc1yer 0\1er a halt ' -space . that was well tested by the GRT (Apsel and Luco, 1983) is used t � or compa1·ison. The TSS Gre . en's functions ot ' the shear dislocation fault are C()mpared in Figure 3 and the TDS in Figure 4. Following the 1·esults of Vidale et al.

DISCUSSI.ON
The FE approach vvitl1 the 1 i ne t() IJ<)int SC)Ltrce t1-c_1nsf<-)rm hc1s been sho\\'n tc) st1ccesst · 11lly 1nc)del the SH wa\1efo1·ms. A r1u1neric�al test ()f tl1i s rnethod [lg,1i 11st. a11 a!l(,l])1tical 111ettl()d c111d the FDM shows a gol)d agree1nent. The 111c1j()f ctdva11tage of tl1is 1n. ethod is tl1ctt the model can have �trt)it1·ary v·elocity an<l der1sity \'a1· iatio11s, i11clt1di11g (li11pi11g [lJ1d ct1rve l1()t1ndaries, veloc.�ity and density gradients, and (. l lo\V-\ . rel<-1cit. y zc)ne. The 111est1 design t' c)r the FE nt1n1e1·ical g1·ids ca11 he i1-1·egt1la1 · . The FE i11odel c<111 et·t·ect.ivelv be dcf . i11e. d <:ts t'i11e g1·ids i11 a lo�' \!e]c)citv Z()ne ()f ,., '-' ., shall()\\' laye1·s with high \)elocity g1·adients� ctn(l theref' ()fe he defined as coarse grids ii1 11igh velocity ar·eas. f-Io\vev'e1·� SL1ch ch(1ngeci grid sizes ai-e n<.)t al lt)wed by the FD 1nethod. Emp]C))'- • ing this method to propagate earthquake energy, the path effects can be computed for the complex earth model in different scales. The application of this appr· oach may be contributed to strong ground motion modeling for shallow structure amplification, determining source parameters in regional distance using broad band seismograms and stud)1ing slab structure at teleseismic distance. The extension to model the P-SV case is straightforward and is easily modified based on this study. The major disadvantages of this approach are that the computa tions take too much time to propagate earthquake energy numerically and the nature of the high frequency results is limited. Usually, the grid size constraints the highest frequency which is allo\\1ed. Furthermore, the model must be two-dimensional. In principle, the third dimen sion, which is perpendicular to the model, should be considered as an infinite extension. The potential applications of this method to earthquake study in the Taiwan area extend into se. veral branches. First, newly deployed strong motion instruments in the wide area of Tai"W·an provide an excellent data set to study strong motion amplification in shallow basins as 8(Jr-S/1(Jlth HiJ.c111g 265 well as the top()grapl1ic e . ffects in the Central Mc)untain Range. Second, modeling seismo grams from the Tai\v,111 east-west cross section may provide important int .

LINE SOURCE
POINT SOURCE l I I I�.
. . • Fil'?· 6 · . Snapshot of .. the line source response (2-D) and its point source correction (3-D) in \�lhC)le space. The input S()Urce . is a band limited 8 -function as shown in Figure 2 with the explosive type source mechanism. T\\lO hori zontal C<lmp<)nents def. ine the model and a vertical component represents the \:vavefield. amplitude. For the symmetry of the �lavefields, the ampli tudes of' half model <)TI tl1e . t' rc)nt side are . muted.

CONCLUSIONS
A new method to nume1·ically synthesize seismogra1ns was proposed and we1l tested. This method is based l)n the equiv·alence body force concept, in V\1hich only suitable body forces are loaded on the source points of the finite element mesh, and no interface betV\1een source area and numerical mode. I is needed. Furthermore, a simple method to calibrate the line source snapshots to point sou1·ce was developed. These Green's functions are easily substi tuted for the standard Green's functions of the linear moment tensor in\1ersion algorithms, but include 2-D structures. Potential applications of this method for seismological research in Taiwan have also been proposed.