May 10, 2010 11:58

AOGS - SE

9in x 6in

b951-v20-ch20

Advances in Geosciences Vol. 20: Solid Earth (2008) Ed. Kenji Satake c World Scientific Publishing Company


THREE-COMPONENT 1D VISCOELASTIC FDM FOR PLANE-WAVE OBLIQUE INCIDENCE ARASH JAFARGANDOMI∗ and HIROSHI TAKENAKA† Department of Earth and Planetary Sciences, Kyushu University, 6-10-1 Hakozaki, Higashi-ku, Fukuoka, 812-8581, Japan ∗ [email protected] † [email protected]

We propose a very efficient scheme for modeling three-component seismic plane waves in vertically heterogeneous attenuative media using the finitedifference (FD) method in the time domain. The scheme is able to represent P , SV, and SH waves incoming from different direction. The algorithm can also calculate the plane-wave responses of media for different incident angels. In the algorithm, neglecting lateral heterogeneity, the wave equation is rewritten for plane waves by applying a 2D Radon transform to the 3D general wave equation. QP and QS are incorporated via generalized Zener body rheological models for viscoelasticity to represent anelastic attenuation. Our FD scheme uses a 1D grid, which leads to a significant reduction in computation time and memory requirements as compared to the corresponding 3D or 2D computations. It may be an efficient tool for pre- or post-analysis of local structural effects including anelastic attenuation, before or after the large-scale seismic wave simulation. To demonstrate the ability and efficiency of the scheme, we calculate a synthetic vertical seismic profile (VSP) for a borehole with a highly heterogeneous velocity model and frequency dependent attenuation model in the Kanto area of Japan.

1. Introduction Computation of the responses of horizontally layered media can be calculated by semi-analytical methods such as the propagator matrix method [e.g. Refs. 29 and 30]. However, for gradient velocity or randomly heterogeneous structures the semi-analytical methods face difficulties. In such cases numerical methods such as finite-difference method have to be used. In the semi-analytical methods, increasing numbers of layers and observation points for a specific size of the structure model will inflate

299

May 10, 2010 11:58

300

AOGS - SE

9in x 6in

b951-v20-ch20

A. Jafargandomi and H. Takenaka

the computation time, while this is not the case in the finite-difference method. The finite-difference (FD) method in the time domain is one of the most common techniques used for simulating seismic wave propagation [e.g. Refs. 11, 13, 16, 25 and 26]. This technique is popular because it is simple and easy to program. In a large-scale 3D simulation using a realistic structure model, even only one job of the computation requires very long computation time and huge computer memory, and thus such jobs could not be run so many times. It may often lead to insufficient tuning of the structure model parameters and/or insufficient interpretation of the simulation results. In order to avoid these insufficiencies without increasing number of trial large 3D computations, we need efficient tools which support construction of the structure model, check of the effects of the changes of the model parameters and so on. In this study we propose an efficient tool for pre- or post-analysis of local structural effects including anelastic attenuation, before or after large-scale seismic wave simulations. That is a multi-component 1D viscoelastic finite-difference scheme for planewave simulation, which efficiently calculates seismic responses of vertically heterogeneous attenuative media to an oblique plane-wave incidence. The proposed algorithm of this paper was first developed by Tanaka and Takenaka24 for perfectly elastic media. JafarGandomi and Takenaka13 then extended the algorithm to the viscoelastic media for P-SV plane-wave modeling. In the subsequent sections we describe the viscoelastodynamic equation in the time domain that governs 3D seismic wave propagation with anelastic attenuation, and derive three-component 1D plane-wave equations from the 3D viscoelastodynamic equation. We then describe the FD scheme exploited here and show a numerical example in a homogeneous medium to confirm the accuracy of the algorithm. To demonstrate the efficiency of the algorithm, we also apply the technique to generate a synthetic VSP for a deep borehole in Kanto area of Japan.

2. Three-Dimensional Viscoelastodynamic Equation Distribution of relaxation mechanisms of rheological models in generalized form can be used to obtain a nearly constant4,17 or any arbitrary frequencydependent1 quality factor over the seismic frequency range of interest. Emmerich and Korn9 showed that a generalized Maxwell body can

May 10, 2010 11:58

AOGS - SE

9in x 6in

b951-v20-ch20

Three-Component 1D Viscoelastic FDM for Plane-Wave Oblique Incidence

301

approximately represent viscoelastic behavior of earth materials for timedomain computation of anelastic seismic waves. Carcione et al.5,6 also used a parallel connection of several standard linear solids (Zener bodies) to introduce a theory of a generalized Zener body (GZB) into a time-domain modeling of seismic waves. The resulting relaxation functions of generalized Maxwell and Zener models are equivalent (e.g. Moczo and Kristek18 ). In this study we use relaxation functions of GZBs4 as
  L  1 τεl −t/τσl H(t), (1) Ψ(t) = MR 1 − 1− l e L τσ l=1

where t is time, τσl and τεl are stress and strain relaxation times of the lth Zener body, L is number of relaxation mechanism, MR is the relaxed modulus, and H(t) is Heaviside step function. Introducing the memory variables5,6 based on a formulation of Robertsson et al.,22 we derive the 3D viscoelastodynamic equation as follows: ρ∂t u˙ i = ∂j σij + fi ,

(2)

∂t σij = (Π − 2M )∂k u˙ k δij + M (∂i u˙ j + ∂j u˙ i ) + rij ,  1  l l rij + (∆π l − 2∆µl )∂k u˙ k δij +∆µl (∂i u˙ j + ∂j u˙ i ) , ∂t rij =− τσl

(3) (4)

L

1 l rij = rij , L

(5)

l=1

l where u˙ j , σij , fi and rij are particle velocities, stress components, body forces and memory variables, respectively; i, j, k = (1, 2, 3); δij is Kronecker’s delta; ∂t and ∂i indicate temporal and spatial derivatives;   L  1 τεP l Π = (λR + 2µR ) 1 − 1− l , L τσ l=1   L  1 τεSl M = µR 1 − 1− l , (6) L τσ l=1

are the unrelaxed moduli; and  Pl  τε l ∆π = (λR + 2µR ) −1 , τσl

 ∆µ = µR l

 τεSl −1 , τσl

(7)

May 10, 2010 11:58

AOGS - SE

9in x 6in

b951-v20-ch20

A. Jafargandomi and H. Takenaka

302

which come from the difference between the relaxed and unrelaxed moduli. The superscripts P and S denote P and S waves, respectively. The relaxed moduli are determined by use of P and S wave velocities at a reference frequency (f0 ) (e.g. JafarGandomi and Takenaka13 ).

3. Viscoelastodynamic Equations for Plane Waves We now consider an incident plane-wave problem in three-dimensional Cartesian coordinates (x1 , x2 , x3 ) with downward positive x3 -direction for a 1D heterogeneous medium where material properties depend only on depth. With this simplification a Radon transform (slant stack) (e.g. Chapman28 ) both with respect to the two horizontal coordinates produces a family of three-component 1D viscoelastic wave equations that govern the evolution of the wavefields in the intercept time-slowness (“plane-wave”) domain, as follows:

 Λp1 1 (Λ + M )p1 p2 (Λ + M )p21 ∂3 u˙ 3 + ∂3 σ23 1+ ∂3 σ13 + ∆P ∆S ∆p ∆P ∆S



 p1 p2 (Λ + M )p21 2(Λ + M )p21 − 1+ r11 − 1+ r12 ∆S ∆P ∆S ∆P

∂t u˙ 1 = −

(Λ + M )p1 p22 r22 , ∆P ∆S

 Λp2 1 (Λ + M )p1 p2 (Λ + M )p22 ∂3 u˙ 3 + ∂3 σ13 ∂t u˙ 2 = − 1+ ∂3 σ23 + ∆P ∆S ∆P ∆P ∆S



 p2 p1 (Λ + M )p22 2(Λ + M )p22 − 1+ r22 − 1+ r12 ∆S ∆P ∆S ∆P −

−

M p1 M p2 1 p1 p2 ∂3 u˙ 1 − ∂3 u˙ 2 + ∂3 σ33 − r13 − r23 , ∆S ∆S ∆S ∆S ∆S  

 Λp2 + 2M p21 p1 (Λ + M )(2M p21 + Λp2 ) = Λ 1+ ∂3 u˙ 3 − Π+ ∆P ∆S ∆P

 p2 (Λ + M )(2M p21 + Λp2 ) × ∂3 σ13 − Λ+ ∂3 σ23 ∆S ∆P

∂t u˙ 3 = − ∂t σ11

(Λ + M )p2 p21 r11 , ∆P ∆S

May 10, 2010 11:58

AOGS - SE

9in x 6in

b951-v20-ch20

Three-Component 1D Viscoelastic FDM for Plane-Wave Oblique Incidence

∂t σ22

∂t σ33



 p21 (Λ + M )(2M p21 + Λp2 ) + 1+ Π+ r11 ∆S ∆P

 p2 (Λ + M )(2M p21 + Λp2 ) + 2 Λ+ r22 ∆S ∆P   2p1 p2 2M p21 + Λp2 + (Λ + M ) 1 + r12 , ∆S ∆P  

 Λp2 + 2M p22 p2 (Λ + M )(2M p22 + Λp2 ) = Λ 1+ ∂3 u˙ 3 − Π+ ∆P ∆S ∆P

 2 2 p1 (Λ + M )(2M p2 + Λp ) × ∂3 σ23 − Λ+ ∂3 σ13 ∆S ∆P

 p22 (Λ + M )(2M p22 + Λp2 ) + 1+ Π+ r22 ∆S ∆P

 p21 (Λ + M )(2M p22 + Λp2 ) + Λ+ r11 ∆S ∆P   2p1 p2 2M p22 + Λp2 r12 , + (Λ + M ) 1 + ∆S ∆P   Λ 2 p2 Λp1 Λp2 = Π+ ∂3 σ13 − ∂3 σ23 + r33 ∂3 u˙ 3 − ∆P ∆p ∆P Λp21 2Λp1 p2 Λp22 r11 + r12 + r13 , ∆P ∆P ∆P   M 2 p1 p2 M p21 M p1 = ∂3 u˙ 2 + M 1 + ∂3 σ33 ∂3 u˙ 1 − ∆S ∆S ∆S   M p21 M p1 p2 + 1+ r23 , r13 + ∆S ∆S

 2(Λ + M )p21 2ΛM p1 p2 M p2 1+ ∂3 σ13 = ∂3 u˙ 3 − ∆P ∆S ∆P

 M p1 2(Λ + M )p22 − 1+ ∂3 σ23 ∆S ∆P 



  M 2 2(Λ + M )p21 2(Λ + M )p22 2 + 1+ p 1+ + p1 1 + r12 ∆S 2 ∆P ∆P



 M p1 p2 M p1 p2 2(Λ + M )p21 2(Λ + M )p22 + 1+ r11 + 1+ r22 , ∆S ∆P ∆S ∆P +

∂t σ13

∂t σ12

303

May 10, 2010 11:58

9in x 6in

b951-v20-ch20

A. Jafargandomi and H. Takenaka

304

∂t σ23

AOGS - SE

  M 2 p1 p2 M p22 M p2 = ∂3 u˙ 1 + M 1 + ∂3 σ33 ∂3 u˙ 2 − ∆S ∆S ∆S   M p22 M p1 p2 + 1+ r13 , r23+ + ∆S ∆S 1 τσl  l  × r11 + (∆π l − 2∆µl )(−p1 ∂t u˙ 1 − p2 ∂t u˙ 2 + ∂3 u˙ 3 ) − 2p1 ∆µl ∂t u˙ 1

l =− ∂t r11

l =− ∂t r22

1 l [r + (∆π l − 2∆µl )(−p1 ∂t u˙ 1 − p2 ∂t u˙ 2 + ∂3 u˙ 3 ) − 2p2 ∆µl ∂t u˙ 2 ] τσl 22

l ∂t r33 =−

1 l [r + (∆π l − 2∆µl )(−p1 ∂t u˙ 1 − p2 ∂t u˙ 2 + ∂3 u˙ 3 ) + 2∆µl ∂3 u˙ 3 ] τσl 33

l ∂t r12 =−

1 l {r + ∆µl [−p1 ∂t u˙ 2 − p2 ∂t u˙ 1 ]}, τσl 12

l =− ∂t r23

 1  l r23 + ∆µl [−p2 ∂t u˙ 3 + ∂3 u˙ 2 ] , l τσ

l ∂t r13 =−

1 l {r + ∆µl [−p1 ∂t u˙ 3 + ∂3 u˙ 1 ]}, τσl 13

(8)

where p1 and p2 are the horizontal slownesses in the x1 - and x2 -directions, respectively, p2 = p21 + p22 , ∆P = ρ − Πp2 , ∆S = ρ − M p2 , and Λ = Π − 2M . Derivation process of viscoelastodynamic equations for plane waves is shown in details for P-SV case in JafarGandomi and Takenaka,13 which is similar to that for the three-component case of Eq. (8). Similar plane-wave approaches have also been done by Blanch et al.3 for viscoacoustic media and by Tanaka and Takenaka24 for perfectly elastic media.

4. FDM Computations To solve a set of the 1D equations derived in Sec. 3, we employ a finitedifference time-domain technique with staggered grids in time and space. For numerical solution of the wave equations several staggered grid finitedifference schemes have been proposed for 2D and 3D elastic25,26,16,11 and viscoelastic media [e.g. Ref. 22]. Here we use a 1D FD staggered grid scheme with second-order accuracy in time and fourth-order accuracy in space.

May 10, 2010 11:58

AOGS - SE

9in x 6in

b951-v20-ch20

Three-Component 1D Viscoelastic FDM for Plane-Wave Oblique Incidence

305

Fig. 1. Arrangement of quantities on staggered grid scheme: (a) spatial arrangement Type 1, (b) spatial arrangement Type 2, and (c) temporal arrangement.

4.1. Finite-difference scheme Figure 1 illustrates the position of the wavefield variables on the 1D staggered grids. We propose two types of spatial grids (Type 1 and Type 2) by putting different velocity components on the free surface, for the cases that the free surface condition need be considered. In Type 1 (Fig. 1a), horizontal components of velocity and in Type 2 (Fig. 1b) vertical component are located on the free surface. Depending on the purpose one of the two grid types may be chosen. For example, in exploration seismology the vertical component might be more important and its exact position on the free surface might then be desirable, while in strong-motion seismology more attention is paid to the horizontal components at the free surface. Using the FD scheme the plane-wave equations (Eq. (8)) are discretized. The discretized equations can be derived for three-component, twocomponent (P-SV) or one-component (SH) plane-wave equations. The PSV and SH equations can be easily derived by setting p2 = 0 in the three-component plane-wave equations (Eq. (8)). The discretized equations for P-SV plane waves are described in JafarGnadomi and Takenaka13 and

May 10, 2010 11:58

306

AOGS - SE

9in x 6in

b951-v20-ch20

A. Jafargandomi and H. Takenaka

the peer-reviewed FORTRAN code related to their article is available at http://software.seg.org/2007/0005.

4.2. Boundary conditions In our scheme the domain boundaries are only at the top and bottom. The lack of artificial side reflections is another advantage of this approach. The free surface is implemented by applying the so-called image methods,16, 21 which extend the grid half a sample above the free surface. Robertsson21 summarized the possibilities for treating particle velocity and stress components above the free surface, which are formally required by the FD scheme as image methods 1, 2, and 3. Tanaka and Takenaka24 showed that using particle velocity as an even function and traction as an odd function with respect to the free surface (image method 2) gives a better result than setting only the tractions as an odd function with respect to the free surface and the particle velocity to zero above the free surface (image method 3). In our scheme we use image method 2. To implement a non-reflecting boundary, we apply the absorbing-zone technique of Cerjan et al.7 and the first-order one-way wave-equation techniques of Clayton and Engquist.8 To choose an appropriate time increment (∆t) we use the relation ∆t < 0.606∆x3 /vmax , where vmax is the highest P -wave apparent velocity in the model. In the case of oblique incidence [vmax = 1/(minimum vertical slowness)]. 4.3. Incidence of initial wave The grid of the proposed algorithm is divided into two main parts: an initial wave zone and a computational domain. For the incident wave, we consider a perfect elastic homogeneous zone. This zone can be set on top of or under the computational domain. There should be no velocity contrast between the initial wave zone and its adjacent layer of the computational domain to prevent artificial reflections and conversions. However, the attenuation contrast may cause a negligibly weak reflection in the case of much highly attenuative models. The thickness of the initial wave zone should be larger than the total width of the incident plane wave and the absorbing zone. The absorbing zone occupies 50 grid points in the examples shown in the subsequent section.

May 10, 2010 11:58

AOGS - SE

9in x 6in

b951-v20-ch20

Three-Component 1D Viscoelastic FDM for Plane-Wave Oblique Incidence

307

5. Numerical Examples We examine the performance of the proposed algorithm with two numerical examples. The first example demonstrates the accuracy of our FD scheme by comparing results from FD computations and an analytical solution in the time domain. For the analytical solution of the plane wave, we use the solution obtained by Blanch et al.2 They incorporate constant Q as a function of frequency analytically via Futterman’s10 expression. In the solution of Blanch et al.,2 the incident angle is not considered; so in the case of oblique incidence in a homogeneous medium, we use apparent velocity as the propagation speed. In the second example, the ability of the algorithm for calculating planewave responses of realistic models for vertically heterogeneous attenuative media is demonstrated by calculating a synthetic vertical seismic profile (VSP) for a borehole in the Kanto area of Japan.

5.1. Accuracy test To check the accuracy of our FD calculations, the result from a numerical simulation is compared with the corresponding analytical solution. The structure model is a homogeneous attenuative medium. The S- and P wave velocities are 1500 m/s and 2000 m/s, respectively, at the reference 3 frequency of f0 = 2 (Hz) and density of 2,000 kg/m . The attenuation model is a constant QP = QS = 50, which is modeled by five Zener bodies for each of QP and QS . In this example, the bottom boundary of the model is set to be nonreflecting. A Ricker wavelet with a central frequency of 2 Hz and is set at a depth of 3,000 m in the initial wave zone as the incident P wave with incident angle of 20◦ and azimuth of 45◦ , and the receiver is placed at a depth of 13,000 m. The computational domain starts from a depth of 3,500 m. The time increment and spatial grid interval are 5 ms and 20 m, respectively. The total number of time steps and grid points are 1,500 and 1,000, respectively. The exact location of the vertical component is half a grid step below the horizontal components, so the vertical propagation distances for horizontal and vertical components are 12,500 m and 12,510 m, respectively. Figure 2 shows the three-components waveforms calculated by the proposed FD scheme, compared with the analytic solution. For better demonstration of viscoelastic effects, the calculated waveforms with the FD scheme in the

May 10, 2010 11:58

AOGS - SE

308

9in x 6in

b951-v20-ch20

A. Jafargandomi and H. Takenaka

Fig. 2. Three components of calculated responses to an incident P wave with incident angle of 20◦ and azimuth of 45◦ .

elastic case are also shown. Clearly there is quite a good agreement between the viscoelastic FD solution and the analytic solution both in amplitude and phase. 5.2. Synthetic VSP The second example includes a simulation of a plane wave propagating through a sedimentary basin in Kanto area via a synthetic VSP. The velocity model is a 2,781-m-deep sedimentary system including sediments and basement in the Kanto area, Japan. The velocity and density structure models (Fig. 3) are models used by JafarGandomi and Takenaka,13 which were constructed from well-log data from the Fuchu borehole (Suzuki et al.23 ) and the velocity model estimated by Yamamizu et al.27 Recently

Fig. 3. Velocity and density models (left and middle panels) and attenuation model (right panel) of the Fuchu borehole site, Japan.

May 10, 2010 11:58

AOGS - SE

9in x 6in

b951-v20-ch20

Three-Component 1D Viscoelastic FDM for Plane-Wave Oblique Incidence

309

Table 1. Values of stress and strain relaxation times used for realization of Q model of Fuchu borehole. l

τσl

τεSl

τεP l

1 2 3 4 5

0.39978 0.10042 0.02522 0.00634 0.00159

0.70598 0.10195 0.02522 0.00634 0.00159

1.20474 0.10350 0.02523 0.00634 0.00159

Kinoshita14 measured the QP and QS values for some deep boreholes in Kanto basin including Fuchu borehole. His Q estimation is in the form of Q(f ) = Q0 f n : for the Fuchu borehole Q0 = 31 ± 7, n = 0.76 ± 0.41 (0.4 < f < 2 (Hz)) for S wave. He also showed that in the Kanto basin QP /QS ≈ 1.5. For calculating the synthetic VSP, we here use his Q model for S wave and QP = 1.5QS for P wave. The Q model is realized by five relaxation mechanisms, i.e. L = 5 for each of P and S waves. The corresponding stress and strain relaxation times are listed in Table 1. For this FD simulation, the spatial grid spacing and time increment are 10 m and 1 ms, respectively. Since the original velocity and density data were digitized with a sampling interval of 1 m (Ref. 23), an averaging technique proposed by Moczo et al.19 is used to incorporate the subgrid structure of velocities and density into the FD grid. This technique helps to achieve accurate arrival times. In this computation we use spatial grid Type 1 (see Fig. 1). A Ricker wavelet with a central frequency of 2 Hz is used as an incident P wave. P - and S-wave velocities of the model correspond to those at the reference frequency of 2 Hz. The incident angle of incoming wave is 15◦ in the bottom half-space from the vertical, and the azimuth is 45◦ . Receivers are located at 50 levels at intervals of 50 m intervals in the borehole. Figure 4(a) and (b) show the three components of synthetic VSP to the oblique incident wave for elastic and viscoelastic cases, respectively. In Fig. 4, major reflections and conversions can be seen between a depth of about 2,000 m and the free surface. In the viscoelastic case (Fig. 4b) the strong attenuation decreases the amplitudes of reflected waves much more for longer propagation. Figure 4 depicts that the proposed algorithm can clearly calculate the reflected and converted phases on the three components.

May 10, 2010 11:58

310

AOGS - SE

9in x 6in

b951-v20-ch20

A. Jafargandomi and H. Takenaka

Fig. 4. Three-component synthetic VSP for the Fuchu borehole site. (a) Elastic and (b) viscoelastic responses to an oblique incident P wave with incident angle of 15◦ and azimuth of 45◦ .

6. Conclusion We have proposed an efficient algorithm to calculate the plane-wave response of a vertically heterogeneous attenuative medium. The planewave equations were derived by applying 2D Radon transform to the 3D viscoelastodynamic equation for 1D structure. The finite-difference time-domain staggered-grid scheme was employed to solve the plane-wave equations numerically. It may be an efficient tool for pre- and post-analysis of local structural effects including anelastic attenuation, before or after large-scale seismic wave simulations. Our scheme is based on the GZB

May 10, 2010 11:58

AOGS - SE

9in x 6in

b951-v20-ch20

Three-Component 1D Viscoelastic FDM for Plane-Wave Oblique Incidence

311

rheological model for incorporating arbitrary QP and QS , while some FDM code for 3D strong-motion simulations employ the GMB model [e.g.12,15,18 ]. However, since the two rheological models are equivalent and there is a relation between them as Moczo and Kristek18 derived, our scheme may also be useful for pre- and post-analysis of 3D simulation based on the GMB model.

Acknowledgments We wish to thank editor Prof. K. Satake and reviewer Prof. Y. Wang and an anonymous reviewer for valuable comments. This study was partially supported by JSPS.KAKENHI (19540448) and JSPS.KAKENHI (19201034).

References 1. S. Asvadurov, L. Knizhnerman and J. Pabon, Geophysics 69 (2004) 817–824. 2. J. O. Blanch, J. O. A. Robertsson and W. W. Symes, Geophysics 60 (1995) 176–184. 3. J. O. Blanch, W. W. Symes and R. J. Versteeg, Society of Exploration Geophysicists Open File Publications 4 (1998) 13–44. 4. J. M. Carcione, Elsevier, New York (2001) 390. 5. J. M. Carcione, D. Kosloff and R. Kosloff, Geophys. J. Roy. Astron. Soc. 93 (1988a) 393–407. 6. J. M. Carcione, D. Kosloff and R. Kosloff, Geophys. J. Roy. Astron. Soc. 95 (1988b) 597–611. 7. C. Cerjan, D. Kosloff, R. Kosloff and M. Reshef, Geophysics 50 (1985) 705– 708. 8. R. W. Clayton and B. Engquist, Bull. Seism. Soc. Am. 67 (1977) 1529–1540. 9. H. Emmerich and M. Korn, Geophysics 52 (1987) 1252–1264. 10. W. Futterman, J. Geophys. Res. 67 (1962) 5279–5291. 11. R. W. Graves, Bull. Seism. Soc. Am. 86 (1996) 1091–1106. 12. R. W. Graves and S. M. Day, Bull. Seism. Soc. Am. 93 (2003) 283–300. 13. A. JafarGandomi and H. Takenaka, Geophysics 72 (2007) H43–H53. 14. S. Kinoshita, Bull. Seism. Soc. Am. 98 (2008) 463–468. 15. J. Kristek and P. Moczo, Bull. Seism. Soc. Am. 93 (2003) 2273–2280. 16. A. Levander, Geophysics 53 (1988) 1425–1436. 17. H. P. Liu, D. L. Anderson and H. Kanamori, Geophys. J. Roy. Astron. Soc. 47 (1976) 41–58. 18. P. Moczo and J. Kristek, Geophys. Res. Lett. 32 (2005) 1306, doi: 10.1029/ 2004GL021598 (2005). 19. P. Moczo, J. Kristek, V. Vavycuk, R. J. Archuleta and L. Halda, Bull. Seism. Soc. Am. 92 (2002) 3042–3066.

May 10, 2010 11:58

312

AOGS - SE

9in x 6in

b951-v20-ch20

A. Jafargandomi and H. Takenaka

20. K. B. Olsen, S. M. Day and C. R. Bradley, Bull. Seism. Soc. Am. 93 (2003) 627–638. 21. J. O. A. Robertsson, Geophysics 61 (1996) 1921–1934. 22. J. O. A. Robertsson, J. O. Blanch and W. W. Symes, Geophysics 59 (1994) 1444–1456. 23. H. Suzuki, R. Ikeda, T. Mikoshiba, S. Kinoshita, H. Sato, and H. Takahashi, Review of Research for Disaster Prevention 65 (1981) 1–162. 24. H. Tanaka and H. Takenaka, Zisin 57 (2005) 343–354. 25. J. Virieux, Geophysics 49 (1984) 1933–1957. 26. J. Virieux, Geophysics 51 (1986) 889–901. 27. F. Yamamizu, H. Takahashi, N. Goto and Y. Ohta, Zisin 34 (1981) 465–479. 28. C. H. Chapman, Geophys. J. Roy. Astr. Soc. 66 (1981) 445–453. 29. W. T. Thomson, J. Appl. Phys. 21 (1950) 89–93. 30. N. A. Haskell, Bull. Seism. Soc. Am. 43 (1953) 17–34.