Citation: JafarGandomi, A., and Takenaka, H., 2010, Three-Component 1D Viscoelastic FDM for Plane-Wave Oblique Incidence, in Advances in Geosciences, Volume20 (Solid Earth), edited by Kenji Satake, World Scientific Publishing Company, Singapore, 311-322.

THREE-COMPONENT 1D VISCOELASTIC FDM FOR PLANEWAVE OBLIQUE INCIDENCE ARASH JAFARGANDOMI Department of Earth and Planetary Sciences, Kyushu University, 6-10-1 Hakozaki, Higashi-ku, Fukuoka, 812-8581, Japan E-mail: [email protected]

HIROSHI TAKENAKA Department of Earth and Planetary Sciences, Kyushu University, 6-10-1 Hakozaki, Higashi-ku, Fukuoka, 812-8581, Japan E-mail: [email protected]

We propose a very efficient scheme for modeling three-component seismic plane waves in vertically heterogeneous attenuative media using the finite-difference (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. 29, 30] However, for gradient velocity or randomly heterogeneous structures the semianalytical 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

1

2

structure model will inflate 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. 11, 13, 16, 25, 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 plane-wave 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 Takenaka [24] for perfectly elastic media. JafarGandomi and Takenaka [13] 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 constant [17, 4] or any arbitrary frequencydependent [1] quality factor over the seismic frequency range of interest. Emmerich and Korn [9] showed that a generalized Maxwell body can 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

3

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 Kristek [18]). In this study we use relaxation functions of GZBs [4] as (t ) = M R (1 −

1 L

L

(1 − l =1

τ εl −t / τ σ )e ) H (t ) , τ σ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, M R is the relaxed modulus, and H(t) is Heaviside step function. Introducing the memory variables [5, 6] based on a formulation of Robertsson et al. [22], we derive the 3D viscoelastodynamic equation as follows:

ρ∂ t ui = ∂ jσ ij + f i , (2) ∂ tσ ij = (Π − 2M )∂ k u k δ ij + M (∂ i u j + ∂ j ui ) + rij ,

(3) ∂ t rijl = −

1

τ σl

[r + (∆π l ij

l

)

]

− 2∆µ l ∂ k u k δ ij + ∆µ l (∂ i u j + ∂ j ui ) ,

(4) rij =

1 L

L

rijl , l =1

(5) where u j , σ ij , f i and rijl 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; Π = ( λR + 2 µ R ) 1 −

1 L τ Sl τ Pl 1 L (1 − εl ) , M = µ R 1 − (1 − εl ) , L l =1 τσ L l =1 τσ

(6) are the unrelaxed moduli; and

4

∆π l = (λR + 2µ R )

τ εPl τ Sl − 1 , ∆µ = µ R εl − 1 , l τσ τσ l

(7) 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 Takenaka [13]). 3. Viscoelastodynamic equations for plane waves We now consider an incident plane-wave problem in three-dimensional Cartesian coordinates ( x1 , x 2 , x 3 ) 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. Chapman [28]) both with respect to the two horizontal coordinates produces a family of threecomponent 1D viscoelastic wave equations that govern the evolution of the wavefields in the intercept time-slowness (“plane-wave”) domain, as follows: ∂ t u1 = − −

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

p1 (Λ + M ) p12 r − p 2 1+ 11 ∆S ∆P ∆S

∂ t u2 = −

1+

2(Λ + M ) p12 (Λ + M ) p1 p22 r , r12 − 22 ∆P ∆P∆S

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

(Λ + M ) p22 r − p1 1 + 2(Λ + M ) p22 r − (Λ + M ) p2 p12 r , p2 1+ 22 12 11 ∆S ∆P ∆S ∆P ∆P∆S Mp Mp2 p p 1 ∂ t u3 = − 1 ∂ 3u1 − ∂ 3u2 + ∂ 3 33 − 1 r13 − 2 r23 , −

S

S

S

S

S

Λp 2 + 2Mp12 p (Λ + M )(2 Mp12 + Λp 2 ) ∂ σ ∂ t σ 11 = Λ 1 + ∂ 3u 3 − 1 Π + 3 13 ∆P ∆S ∆P −

p2 (Λ + M )( 2Mp12 + Λp 2 ) ∂ σ Λ+ 3 23 ∆S ∆P

+ 1+ +

p12 (Λ + M )( 2Mp12 + Λp 2 ) Π+ ∆S ∆P

r11

p22 (Λ + M )(2 Mp12 + Λp 2 ) r + 2 p1 p2 (Λ + M ) 1 + 2Mp12 + Λp 2 r , Λ+ 22 12 ∆S ∆P ∆S ∆P

5

p Λp 2 + 2Mp22 (Λ + M )( 2Mp22 + Λp 2 ) ∂ σ ∂ 3u3 − 2 Π + 3 23 ∆P ∆S ∆P

∂ tσ 22 = Λ 1 + −

p1 (Λ + M )(2 Mp22 + Λp 2 ) ∂ σ + Λ+ 3 13 ∆S ∆P

+

p12 (Λ + M )(2Mp22 + Λp 2 ) r + 2 p1 p2 (Λ + M ) 1 + 2Mp22 + Λp 2 r , Λ+ 11 12 ∆S ∆P ∆S ∆P

∂ tσ 33 = Π +

1+

p22 (Λ + M )(2Mp22 + Λp 2 ) Π+ ∆S ∆P

Λp 2 Λ2 p 2 Λp Λp Λp 2 2Λp1 p2 ∂ 3u3 − 1 ∂ 3σ 13 − 2 ∂ 3σ 23 + r33 + 1 r11 + r12 + 2 r13 , ∆P ∆p ∆P ∆P ∆P ∆P

∂ tσ 13 =

M 2 p1 p2 Mp12 Mp1 Mp12 Mp1 p2 ∂ 3u2 + M 1 + ∂ 3u1 − ∂ 3σ 33 + 1 + r13 + r23 , ∆S ∆S ∆S ∆S ∆S

∂ tσ 12 =

2ΛMp1 p2 Mp2 2(Λ + M ) p12 Mp1 2(Λ + M ) p22 ∂ 3u3 − 1+ ∂ 3σ 13 − 1+ ∂ 3σ 23 ∆P ∆S ∆P ∆S ∆P M 2 2(Λ + M ) p12 2(Λ + M ) p22 p2 1 + + p12 1 + ∆S ∆P ∆P

+ 1+

Mp1 p2 2(Λ + M ) p12 Mp1 p2 2(Λ + M ) p22 1+ r11 + 1+ r22 , ∆S ∆P ∆S ∆P

+ ∂ tσ 23 =

r12

M 2 p1 p2 Mp22 Mp2 Mp22 Mp1 p2 ∂ 3u1 + M 1 + ∂ 3u2 − ∂3σ 33 + 1 + r23+ + r13 , ∆S ∆S ∆S ∆S ∆S

∂ t r11l = −

1

τ σl

∂ t r22l = −

1

τ σl

[r + (∆π l 11

[r + (∆π l 22

)

− 2∆µ l (− p1∂ t u1 − p 2 ∂ t u2 + ∂ 3 u3 ) − 2 p1 ∆µ l ∂ t u1

l

l

)

]

− 2∆µ l (− p1∂ t u1 − p2 ∂ t u 2 + ∂ 3 u3 ) − 2 p2 ∆µ l ∂ t u 2

]

1

[r l + (∆π l − 2∆µ l )(− p1∂ t u1 − p2 ∂ t u 2 + ∂ 3u3 ) + 2∆µ l ∂ 3u3 ] τ σl 33 1 ∂ t r12l = − l {r12l + ∆µ l [− p1∂ t u 2 − p2 ∂ t u1 ]}, τσ

∂ t r33l = −

∂ t r23l = − ∂ t r13l = −

1

τ σl 1

τ σl

{r

l 23

+ ∆µ l [− p 2 ∂ t u3 + ∂ 3u 2 ] ,

}

{r

l 13

+ ∆µ l [− p1∂ t u3 + ∂ 3 u1 ] ,

}

(8) where p1 and p 2 are the horizontal slownesses in the x1- and x2-directions,

r22

6

respectively, p 2 = p12 + p 22 , ∆ P = ρ − Πp 2 , ∆ S = ρ − Mp 2 , and Λ = Π − 2 M . 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 Takenaka [24] 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 elastic [25, 26, 16, 11] and viscoelastic media [e.g. 22]. Here we use a 1D FD staggered grid scheme with second-order accuracy in time and fourth-order accuracy in space. 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.

7

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

The discretized equations can be derived for three-component, two-component (P-SV) or one-component (SH) plane-wave equations. The P-SV and SH equations can be easily derived by setting p 2 = 0 in the three-component planewave equations (Eq. 8). The discretized equations for P-SV plane waves are described in JafarGnadomi and Takenaka [13] and 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. Robertsson [21] 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 Takenaka [24] 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].

8

To choose an appropriate time increment ( ∆t ) we use the relation ∆t < 0.606∆x3 / vmax , where v max is the highest P-wave apparent velocity in the

model. In the case of oblique incidence [ v max =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. 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’s [10] 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 frequency of f0 = 2 (Hz) and density of 2000 kg/m 3 . The attenuation model is a constant

9

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 3000 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 13000 m. The computational domain starts from a depth of 3500 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 1500 and 1000, 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 12500 m and 12510 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 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.

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

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 2781-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 Kinoshita [14] 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 ) = Q 0 f n : for the Fuchu borehole Q0 = 31 ± 7 , n = 0.76 ± 0.41 (0.4
10

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.

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

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 ([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 Type1 (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,

11

Figure 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°.

major reflections and conversions can be seen between a depth of about 2000 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. Table 1. Values of stress and strain relaxation times used for realization of Q model of Fuchu borehole. l τ Pl τl τ Sl 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

12

6. Conclusion We have proposed an efficient algorithm to calculate the plane-wave response of a vertically heterogeneous attenuative medium. The plane-wave 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 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 Kristek [18] 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, 817–824 (2004). 2. J. O. Blanch, J. O. A. Robertsson, and W. W. Symes, Geophysics 60, 176– 184 (1995). 3. J. O. Blanch, W. W. Symes, and R. J. Versteeg, Society of Exploration Geophysicists Open File Publications No. 4, 13–44 (1998). 4. J. M. Carcione, Elsevier, New York, 390 pp. (2001). 5. J. M. Carcione, D. Kosloff, and R. Kosloff, Geophys. J. Roy. Astron. Soc. 93, 393–407 (1988a). 6. J. M. Carcione, D. Kosloff, and R. Kosloff, Geophys. J. Roy. Astron. Soc. 95, 597–611 (1988b). 7. C. Cerjan, D. Kosloff, R. Kosloff, and M. Reshef, Geophysics 50, 705–708 (1985). 8. R. W. Clayton and B. Engquist, Bull. Seism. Soc. Am. 67, 1529–1540 (1977). 9. H. Emmerich and M. Korn, Geophysics 52, 1252–1264 (1987). 10. W. Futterman, J. Geophys. Res. 67, 5279–5291 (1962). 11. R. W. Graves, Bull. Seism. Soc. Am. 86, 1091–1106(1996).

13

12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30.

R. W. Graves and S. M. Day, Bull. Seism. Soc. Am. 93, 283–300 (2003). A. JafarGandomi and H. Takenaka, Geophysics, 72, H43–H53 (2007). S. Kinoshita, Bull. Seism. Soc. Am. 98, 463–468 (2008). J. Kristek and P. Moczo, Bull. Seism. Soc. Am. 93, 2273–2280 (2003). A. Levander, Geophysics 53, 1425–1436 (1988). H. P. Liu, D. L. Anderson, and H. Kanamori, Geophys. J. Roy. Astron. Soc. 47, 41–58 (1976). P. Moczo and J. Kristek, Geophys. Res. Lett. 32, 1306, doi: 10.1029/2004GL021598 (2005). P. Moczo, J. Kristek, V. Vavycuk, R. J. Archuleta, and L. Halda, Bull. Seism. Soc. Am. 92, 3042–3066 (2002). K. B. Olsen, S. M. Day, and C. R. Bradley, Bull. Seism. Soc. Am. 93, 627– 638 (2003). J. O. A. Robertsson, Geophysics 61, 1921-1934 (1996). J. O. A. Robertsson, J. O. Blanch, and W. W. Symes, Geophysics 59, 1444–1456 (1994). H. Suzuki, R. Ikeda, T. Mikoshiba, S. Kinoshita, H. Sato, and H. Takahashi, Review of Research for Disaster Prevention 65, 1–162 (1981). H. Tanaka and H. Takenaka, Zisin 57, 343–354 (2005). J. Virieux, Geophysics 49, 1933–1957 (1984). J. Virieux, Geophysics 51, 889–901 (1986). F. Yamamizu, H. Takahashi, N. Goto, and Y. Ohta, Zisin 34, 465–479 (1981). C. H. Chapman, Geophys. J. Roy. Astr. Soc. 66, 445–453 (1981). W. T. Thomson, J. Appl. Phys. 21, 89–93 (1950). N. A. Haskell, Bull. Seism. Soc. Am. 43, 17–34 (1953).