Large fluctuation of wave amplitude produced by small fluctuation of velocity structure

Physics of the Earth and Planetary Interiors 120 Ž2000. 201–217 www.elsevier.comrlocaterpepi

Large fluctuation of wave amplitude produced by small fluctuation of velocity structure Mitsuyuki Hoshiba) Meteorological Res. Inst., Japan, Nagamine 1-1, Tsukuba, Ibaraki-ken, 3050052, Japan Received 1 January 1998; accepted 1 August 1999

Abstract Seismic velocity structure of the real Earth is described by a combination of deterministic large scale structure and small scale fluctuation. When waves propagate in media having the small scale fluctuation, the fluctuation makes the rays bend, introduces focussing and defocusing, and produces fluctuations of wave amplitude. To estimate the size of wave amplitude fluctuation that is caused by velocity fluctuation, the variances of wave amplitudes are evaluated as a function of propagation distance from numerical simulations of scalar wave propagation in 3D random media having Gaussian or exponential autocorrelation functions. The simulations are performed using a phase screen method, which is a numerical procedure for solving the parabolic wave equation. The following conclusions are obtained from the simulations that use two types of incident wave: stationary monochromatic sine wave and transient Ricker wavelet. For velocity fluctuation as small as several percent, the variance of wave amplitude is larger than the square of its average in some cases of high frequency waves. These results show the importance of small error and fine spatial resolution of velocity structure model for proper modeling of wave amplitude. When the fluctuation of the velocity structure is large, the wave amplitude variance increases with the propagation distance until it reaches some peak value after which it decreases with the increase of propagation distance. The variance for the Ricker wavelet incidence is smaller than that for a monochromatic sine wave. Because of the stochastic dispersion due to the small fluctuation of velocity structure, the maximum amplitude decreases with the propagation distance even if there is neither an intrinsic absorption nor back-scattering for the Ricker wavelet incidence. q 2000 Elsevier Science B.V. All rights reserved. Keywords: Random media; Scattering; Attenuation; Wave propagation

1. Introduction Conventional deterministic models of seismic velocity in the Earth represent large scale structures, and the wave propagation in such structures had

)

Tel.: q81-298-52-9382; fax: q81-298-51-3730. E-mail address: [email protected] ŽM. Hoshiba..

been well studied in seismology. An example is the wave propagation in horizontally stratified media Že.g., Kennett, 1983.. Wave amplitude for high frequencies larger than 1 Hz can sometimes be modeled using the horizontally stratified media; however, seismic velocity structure of the real Earth is described by a combination of the large scale features and small scale spatial fluctuation, as has been observed in deep boreholes Že.g., Wu et al., 1994;

0031-9201r00r$ - see front matter q 2000 Elsevier Science B.V. All rights reserved. PII: S 0 0 3 1 - 9 2 0 1 Ž 9 9 . 0 0 1 6 5 - X

202

M. Hoshibar Physics of the Earth and Planetary Interiors 120 (2000) 201–217

Jones and Holliger, 1997.. When waves propagate through media having small fluctuation in velocity, the spatial fluctuation makes the rays bend and causes focussing if the velocity is smaller than that of the surroundings, and defocusing if it is larger ŽFig. 1.. Focussing makes amplitude large and defocusing makes it small, so the wave amplitudes are influenced by small scale fluctuations that are superimposed on a large scale feature. Small scale fluctuation produces fluctuations of wave amplitude. Coda waves of local small earthquakes are considered to be waves that are scattered by inhomogeneities in the lithosphere ŽAki, 1969; Aki and Chouet, 1975., and the inhomogeneity is considered to be the small scale fluctuation of velocity ŽSato, 1984; Yoshimoto et al., 1997.. Wu and Aki Ž1988. compiled measurements of fractional velocity fluctuation, and their figure shows that the RMS fractional velocity fluctuation lies between 0.01 and 0.1 for the lithosphere to the upper mantle. Many models have been proposed to relate the scattering characteristics and the coda wave amplitude Že.g., Sato, 1977; Hoshiba, 1991, 1995; Zeng et al., 1991., and a number of authors have investigated the attenuation

of wave amplitude due to the scattering Že.g., Sato, 1984, 1989; Fehler et al., 1992; Mayeda, 1992; Hoshiba, 1993; Yoshimoto et al., 1993.. They have focused only on the average amplitude of coda waves, not the fluctuation of wave amplitude Žor the variance of the distribution of the amplitude.. Fig. 2 shows an example of attenuation of observed wave amplitude for 2–4 Hz after corrections of geometrical spreading, source power and site amplification factor. Though the average of the wave amplitude decreases with distance, large variances Žthat is, fluctuation. of wave amplitude are observed. As mentioned above, the spatial fluctuation of the wave amplitude should be expected if inhomogeneity exists in the lithosphere. This paper investigates the relationship between the small scale fluctuation of velocity structure and the variance of the distribution of the wave amplitude propagating through the fluctuating media. The effects of inhomogeneous media on wave propagation are characterized by Born’s or Rytov’s theory Že.g., Ishimaru, 1978. if the scattering is weak. If scattering is sufficiently strong, asymptotic theory, e.g., strong fluctuation theory ŽIshimaru,

Fig. 1. Schematic illustration of the propagation of a wave through a fluctuating medium. The inhomogeneity makes rays bend, introduces focussing and defocusing, and produces fluctuation of travel time and wave amplitude.

M. Hoshibar Physics of the Earth and Planetary Interiors 120 (2000) 201–217

Fig. 2. Example of observed wave amplitude, A, for 2–4 Hz versus hypocentral distance, r, after corrections of geometrical spreading, source power and site amplification factor using coda normalization method ŽAki, 1980.. These data are collected at central Kushu, Japan.

1978; Shapiro and Kneib, 1993. can be applied. Theoretical works are summarized in, e.g., Chernov Ž1960., Tatarskii Ž1961., Tatarskii et al. Ž1993., Ishimaru Ž1978., and Uscinski Ž1980.. However, many important phenomena in wave propagation lie in the region between weak and strong scattering. The theory to interpolate between these two extreme cases is so difficult that simulations are currently the most powerful and flexible tools to provide an accurate information about wave propagation in inhomogeneous media. The fluctuation of the wave amplitude is an important subject in optics, acoustic, oceanography and meteorology, and it is characterized by a parameter called scintillation index which is defined as the variance of squared amplitude normalized by the square of its average. Martine and Flatte Ž1988, 1990. estimated the scintillation index of the incident monochromatic sine waves for some random media using the phase screen method, which is based on the parabolic approximation. In seismology, the maximum amplitude has often been used to represent the strength of ground motion, and a pulse wave is often used as a model of source time function. In the present paper, we estimate the variance of the maximum amplitude for an incident transient wave ŽRicker wavelet. and compare it with that of monochromatic sine waves.

2. Method of simulation Following Martine and Flatte Ž1988., we simulate the wave propagation in 3D random inhomogeneous

203

media using the phase screen technique, which is based on the parabolic approximation, to estimate the variance of maximum amplitude. In this section, the method is described. Intrinsic absorption of wave energy is not considered in this simulation. Let us imagine waves propagating through a locally isotropic medium ŽAki and Richards, 1980. as schematically shown in Fig.1. For x ) 0, the inhomogeneities are characterized by a correlation length a. Neglecting the conversion scattering between the P and S waves, the S wave propagation is approximated by the scalar wave propagation, 1

Dy

V Ž x.

2

E t2 P u Ž x,t . s 0,

Ž 1.

where D is the Laplacian and V Žx. is the S wave velocity. We introduce the fractional fluctuation of the reciprocal of the velocity square as, 2

2 n Ž x. s

1rV Ž x . y 1rV02

Ž 2.

1rV02

where V0 s ² V Žx.: and n is a real number. Angle brackets mean statistical average. Then, wave Eq. Ž1. is written as:

Dy

1 V02

 1 q 2 n Ž x . 4 E t2 P u Ž x,t . s 0.

Ž 3.

We assume here that the backscattering is negligible and, therefore, restrict propagation to the forward Žqx . direction. The wave field u is written as a superposition of plane waves U propagating in the x direction, u Ž x ,w;t . s

1

`

H d v U Ž x ,w; v . e 2p y`

iŽ k xyw t .

Ž 4.

where w is a displacement vector in the transverse Ž y–z . plane, v is the angular frequency, and k ' vrV0 . Using Laplacian in the transverse plane, D T s E y2 q Ez2 , the wave equation in the frequency domain is described as 2 ik Ex q D T q 2 k 2 n U q Ex2 U s 0.

Ž 5.

For ka 4 1, U should be a slowly varying function of x and varies only over the distance of the scale

204

M. Hoshibar Physics of the Earth and Planetary Interiors 120 (2000) 201–217

size a of the inhomogeneities, we note that < kE x U < 4 < E x2 U <, and therefore neglect the last term, 2 ik Ex q D T q 2 k 2 n U s 0.

Ž 6.

The approximation is called parabolic approximation Že.g., Ishimaru, 1978., and is known to be applicable to the forward propagation even if there is some large angle scattering. Born theory predicts that forward scattering is dominant when ka 4 1 ŽSato and Fehler, 1998.. To solve Eq. Ž6., we use the phase screen method ŽFisk and McCartor, 1991; Fisk et al., 1992; Wu, 1994; Liu and Wu, 1994; Wild and Hudson; 1998.. In this method, the 3D inhomogeneous medium is represented by a sequence of 2D screens perpendicular to the major direction of wave propagation. For homogeneous media Ž n s 0., Eq. Ž6. is expressed, 2 ik Ex q k 2 P TT  U Ž x ,w; v . 4 s 0

Ž 7.

where TT is the 2D Fourier transform of Ž y,z . into k s Ž k y , k z .. Here TUŽ x q D x, w; v .4 is approximately described using TUŽ x, w; v .4

Thus, the procedure to propagate waves in inhomogeneous media over D x is approximately expressed from Eqs. Ž8. and Ž10., U Ž x q D x ,w; v . f Ty1 exp Ž yiAD x . P TT  exp Ž i u Ž x , y, z . . T =U Ž x ,w; v . 4 .

Ž 11 .

The process is repeated for all screens until the observing plane is reached. Martine Ž1992. described the error and limits of application for the phase screen method. In Appendix A, waveforms calculated using the phase screen method are compared with the theoretical waveforms.

3. Random media The procedure to generate inhomogeneous random media will be explained. We write Eq. Ž2. approximately for the case of small fluctuation, 2

2 nŽ x . s

1rV Ž x . y 1rV02 1rV02

f2

V0 y V Ž x . V0

Ž 12 .

TT  U Ž x q D x ,w; v . 4 f exp Ž yiAD x . P TT  U Ž x ,w; v . 4 ,

Ž 8.

where A s Ž k y2 q k z2 .r2 k, and D x is the interval between screens. For a plane wave propagating normal to the x–y plane, that is k y s k z s 0, UŽ x q D x,w; v . equals to UŽ x,w; v . independently of v or D x. This means that the numerical dispersion Žor grid dispersion. due to large grid size is not expected in the phase screen method, which is different from the finite difference method. Phase shift is approximately given locally by neglecting the term of D T in Eq. Ž6., 2 ik Ex q 2 k 2 n P U Ž x ,w; v . s 0.

Ž 9.

Solving Eq. Ž9., we obtain U Ž x q D x ,w; v . f exp i u Ž x , y, z . U Ž x ,w; v . , Ž 10 .

The statistical average of nŽ x . is zero; ² nŽ x .: s 0. Let fluctuation n be statistically homogeneous and isotropic, and have Gaussian or exponential random statistics. Then the randomness is characterized by the autocorrelation function ² nŽ x . nŽ xX .:, which is a function of < x I xX < ² n Ž x . n Ž xX . : s n20 exp Ž y< x y xX < 2ra2 . for Gaussian autocorrelation function ² n Ž x . n Ž xX . : s n20 exp Ž y< x y xX
Ž 13 .

where n 20 ' ² n

¦ž

f where u Ž x, y,z . s knD x.

.

2

Ž x . :s

1 4

V0 y V Ž x . V0

¦ž

2

1rV Ž x . y 1rV02 1rV02

2

/;

2

/;

,

Ž 14 .

M. Hoshibar Physics of the Earth and Planetary Interiors 120 (2000) 201–217

n 0 is the approximation of RMS velocity fluctuation. Fourier transform of the autocorrelation function is the power spectrum SŽ l ., T Ž ² n Ž x . n Ž xX . : . s < N Ž l . < 2r Ž Nx D x P Ny D y P Nz D z . ' S Ž l . ,

Ž 15 . where T is the 3D Fourier transform of x s Ž x, y, z . into l s Ž l x ,l y ,l z ., N Ž l . s T Ž nŽ x .., D x, D y and D z are the sampling interval to make the medium discrete along x-, y- and z-axes, respectively, and Nx , Ny and Nz are the total number of samples in these orthogonal directions. SŽ l . is

Ž 16X . where C is determined from the relation that SX Ž l . also satisfies with Eq. Ž15., that is, T Ž ² n Ž x . n Ž x X . : . s SX Ž l . . Ž 15X . Then,

3

3

2

P SX Ž l . d l y S Ž 0 . r Ž Nx D x P Ny D y P Nz D z . .

H

for exponential function,

Ž 16 .

where l s < l <. Thus, the random fluctuation nŽ x . is given by

Recalling that 1rŽ2 p . 3 P

2 0

H SŽ l .d l s n , Ž 21 .

 S Ž l . P Ž Nx D x P Ny D y P Nz D z . 4

Pexp  i w Ž l . 4 ,

1r2

Ž 17 .

where w Ž l . is a uniform random variable within w0,2p . that satisfies w Žyl . s yw Ž l .. Once n 0 , a and V0 are given, random media are created according to Eq. Ž17.. However, we have a problem for the practical use of Eq. Ž17.. Because N Ž l . is the Fourier transform of nŽ x ., N Ž0. is expressed as N Ž 0 . s n Ž x . d x s ² n Ž x . : Ž Nx D x P Ny D y P Nz D z . .

H

Ž 18 . Because ² nŽ x .: is zero by definition, N Ž0. is zero. On the other hand, from Eq. Ž15. N Ž 0 . s  S Ž 0 . P Ž Nx D x P Ny D y P Nz D z . 4

Ž 20 .

C s n 20r  n 20 y S Ž 0 . r Ž Nx D x P Ny D y P Nz D z . 4 .

n Ž x . s Ty1 < N Ž l . < P exp  i w Ž l . 4 sT

2

SX Ž l . s C P 8p n 20 a3r Ž 1 q a2 k 2 . for l / 0 SX Ž l . s 0 for l s 0 for exponential function

s 1r Ž 2p . P HSX Ž l . d lsC P 1r Ž 2p .

for Gaussian function

y1

Here SŽ0. is not zero as required by Eq. Ž16., so N Ž0. is not zero, which is contrary to Eq. Ž18.. To avoid this problem, we introduce SX Ž l . instead of SŽ l ., SX Ž l . s C P p 3r2 n 20 a 3 exp Ž ya2 l 2r4 . for l / 0 SX Ž l . s 0 for l s 0 for Gaussian function

2 n20 s ² n Ž x . :s² n Ž x . n Ž x . :

S Ž l . s p 3r2 n 20 a3 exp Ž ya2 l 2r4 .

S Ž l . s 8p n 20 a 3r Ž 1 q a 2 l 2 .

205

1r2

X

The above formula Ž16 . satisfies ² nŽ x .: s 0 and ² nŽ x . 2 : s n 20 . SX Ž l . is used instead of SŽ l . for the practical use of Eq. Ž17.. In this paper, simulations are performed for four different a, that is a s 2, 4, 8 and 16 km. We choose Nx s 256 and Ny s Nz s 128, and D x s D y s D z s 0.25 km for a s 2 km, and D x s D y s D z s 0.5 km for a s 4, 8 and 16 km. 4. Simulation results A monochromatic sine wave and a Ricker wavelet of unit amplitude are used for the source time function; they are defined as shown in Fig. 3, R s Ž t . s sin Ž 2p ft . for monochromatic sine wave 2 R r Ž t . s 1 y 2 p f c Ž t y 1rfc .

½

5

=exp y p f c Ž t y 1rfc .

½

.

Ž 19 .

for Ricker wavelet

2

5 Ž 22 .

M. Hoshibar Physics of the Earth and Planetary Interiors 120 (2000) 201–217

206

for sine wave m2r

Ž x. s

s

¦Ž I

max y

² Imax : T . 2

where f and f c are the frequency and predominant frequency, respectively. The Ricker wavelet can be considered to vanish outside the interval w0, 2rfc x. The spectra of the two source functions are an impulsive function, that is dŽ f ., for the sine wave, and a truncated function over the interval w0, 3 f c x for the Ricker wavelet, respectively. A plane wave incident upon at x s 0 normal to the y–z plane propagates into the random media defined in the previous section. Fluctuation of wave amplitude caused by the randomness of velocity structure is observed. We consider intensity, I Ž x,w . '
Ž x. s

² Ž I y ² I : T . 2 :T ² I :2T

s

² I 2 : T y ² I :2T ² I :2T

T

2 2 : ² Imax T y ² Imax : T

² Imax :2T

for Ricker wavelet

Fig. 3. Sine wave, R sŽ t ., and Ricker wavelet, R r Ž t ., which are used as the source time function. It is shown in Ža. as a function of time; f s1 Hz for sine wave and f c s1 Hz for Ricker wavelet. The spectra are shown in Žb.. The spectrum is a spike-like function Ždelta function. for sine wave, and limited by 0 and 3 f c for Ricker wavelet.

;

² Imax :2T

Ž 23 .

In the above formula, m 2s is the well known scintillation index in optics, and m s and m r correspond to the standard deviation divided by the average. The initial conditions are I Ž0,w . s Imax Ž0,w . s 1, and m2s Ž0. s m2r Ž0. s 0. Fig. 4 shows ² I : T , m2s , ² Imax : T and m2r depending on x for f s f c s 1 Hz for the random media with n 0 s 0.02, a s 4 km and V0 s 3 kmrs for the Gaussian autocorrelation function. The number of media, which satisfy the three stochastic parameters, is deterministically infinite. In this paper, simulations are performed using five different media which have common stochastic parameters. For sine wave, ² I : T does not depend on propagation distance x, which means energy is conserved because intrinsic absorption is not considered and backscattering is neglected here. The method used in this paper ŽEq. 11. only adjusts the phase so no energy is lost from forward scattered wave field. For the Ricker wavelet, however, ² Imax : T decreases with increasing x, which means that the maximum amplitude decreases with distance even if there is no intrinsic absorption nor backscattering. These phenomena are explained by the stochastic dispersion of body waves in random media ŽMclaughlin and Anderson, 1987; Muller et ¨ al., 1992; Shapiro et al., 1996.. Waves in random media prefer fast paths and, therefore, the apparent velocity of wave propagation is larger than the average velocity, which is called velocity shift by Muller ¨ et al. Ž1992. and Roth et al. Ž1993.. The velocity shifts are frequency dependent because the effects of diffraction vary with different frequency and, therefore, the ray paths and travel times are different between high and low frequency waves. The sine wave has only one frequency but the Ricker wavelet has frequencies in the range w0, 3 f c x. The maximum amplitude of the Ricker wavelet is caused by the constructive interference of many frequencies within w0, 3 f c x. The constructive interference becomes in-

M. Hoshibar Physics of the Earth and Planetary Interiors 120 (2000) 201–217

207

Fig. 4. Intensity and its variance as a function of propagation distance, x. Incident wave has frequency, f s 1 Hz for sine wave, and fc s 1 Hz for Ricker wavelet. Random media have Gaussian autocorrelation function of correlation length, a s 4 km, average velocity, ² n : s 3 kmrs, and fluctuation, n 0 s 0.02. Simulations are performed using five media each having the same random media description. Above: average intensity ² I : T for sine wave, and average maximum intensity ² Imax : T for the Ricker wavelet, where ²: T means a statistical average on the y–z plane. Bottom: m2s and m2r which are the variances of I or Imax normalized by the square of ² I : T or ² Imax : T . Predictions of m2s from the weak scattering theory and strong scattering theory are shown in the figure.

complete due to the propagation in random media because of the frequency dependent dispersion. The temporal width of the Ricker wave is 2rfc at x s 0; however, the width becomes larger with increasing x, that is, the pulse broadening occurs in random media ŽWilliamson, 1975; Sato, 1989.. Maximum amplitude must decrease to maintain energy conservation because of the pulse broadening. The theory of weak scattering gives a relation for m2s using Rytov’s theory, m2s f 2'p n 20 k 2 ax 1 y

ž

1 D

/

tany1 D ,

Ž 24 .

ŽChernov, 1960; Aki, 1973; Ishimaru, 1978., where D s 4 xrŽ ka2 .. Strong scattering theory predicts that m2s approaches to 1 as x approaches infinity. For the case of n 0 s 0.02 ŽFig. 4., m2s and m2r increase monotonously with increasing x without reaching a

peak. Because of stochastic dispersion, m2r is smaller than m2s . Fig. 5 is the same as Fig. 4, but for n 0 s 0.05. ² I : T does not depend on x again, and ² Imax : T decreases with x; the rate of decrease is larger than that for n 0 s 0.02. ² Imax : T decreases to be about 0.5 at x s 60 km, which corresponds to apparent Qy1 of 6 = 10y3 . For n 0 s 0.05 ŽFig. 5., m 2s and m 2r increase with x for small x, reaching peaks at around x s 20–40 km, and then decreasing. At x s 20–40 km, m 2r is larger than 1 which means quite large variance and that the standard deviation is larger than the average. At small propagation distance, the fluctuation of wave amplitude increases due to repeated focussing or defocusing. Once wave energy is concentrated, to a certain extent, to certain rays, it is difficult to concentrate more energy into the same rays. The probability of repeated focussing or defocusing more

208

M. Hoshibar Physics of the Earth and Planetary Interiors 120 (2000) 201–217

Fig. 5. Same as Fig. 4, but for n 0 s 0.05.

than several times along a common ray path is very small, therefore, m2s and m2r decrease with distance after reaching peaks. In Fig. 6, the time dependent intensityŽs wamplitudex 2 . along the y-axis is shown for the case of the incident Ricker wavelet. The spatial variation in the waveform shape is small and we can recognize the original form of Ricker wavelet at x s 5.5 km where m2r is small, but the spatial waveform variation is significantly large and we cannot recognize the original waveform any more at x s 28 km where m 2r is larger than 1. The distribution of Imax on the y–z planes at several values of x is shown in Fig. 7 for n 0 s 0.05. The left side uses linear–linear scale and the right side uses log–linear scale. At x s 0.5 km, the waveform is almost the same for all y, z and the fluctuation of the maximum amplitude is little, so that the histogram is almost a spike. At x s 30.5 km where m2r is larger than 1, the range of variation of Imax is larger than one order, that is, the largest Imax is larger than 10 times of the smallest Imax . The inten-

sity distribution broadens for x F 30.5 km as x increases, but the broadening reverses for x ) 30.5 km in Fig. 7, which corresponds to the phenomenon that m 2r reaches a peak at around x s 30 km in Fig.5. The comparison of both scales leads to conclusion that the distribution of Imax is approximately a log–normal distribution rather than a normal distribution. Fig. 8 is the same as Fig. 7, but for the maximum amplitude Žthat is, square root of Imax . instead of Imax . The maximum amplitude is also approximately a log–normal distribution rather than a normal distribution. Mclaughlin and Anderson Ž1987. also found a log–normal distribution in amplitude for random media using a Gaussian beam technique. The left of Fig. 9 shows averaged m2r vs. dimensionless travel distance Ž k c x . for several values of k c a calculated using the results for a s 2, 4, 8, and 16 km, where k c s 2p fcrV0 . The right picture of this figure shows the results when the randomness has exponential autocorrelation function. For n 0 s 0.05, m2r increases with increasing k c x and then

M. Hoshibar Physics of the Earth and Planetary Interiors 120 (2000) 201–217

Fig. 6. An example of time-dependent intensity w s Žamplitude. 2 x for Ricker wavelet source Ž f c s1 Hz. at propagation distance, x s 5.5 Žleft. and 28 km Žright.. The random media have Gaussian autocorrelation function with as 4 km, ² n : s 3 kmrs and n 0 s 0.05.

209

path effect due to the small fluctuation of medium because several percent fluctuation produces large fluctuation of wave amplitude whose distribution is a log–normal. The large scatter of wave amplitude means that a few data of wave amplitude have no meaning for the estimation of Q, or other stochastic parameters representing media; the average of a sufficient large number of data is necessary as suggested by Jannaud et al. Ž1991.. An accurate velocity structure model is important for the discussion of amplitude on a certain observed seismogram, especially for the higher frequency waves. If we use inappropriate models and then compare a synthesized waveform with an observed one, the synthesis must be different from real observation even when the source model is precise. We need an appropriate velocity structure model. How fine and precise the model required for the discus-

reaching a peak, then it decreases for both cases of the Gaussian and exponential autocorrelation function.

5. Discussion As shown in Fig. 2, for the estimates of the Q value the wave amplitudes are measured as a function of propagation distance after corrections of geometrical spreading, source power and site amplification factor. Large scatter of the wave amplitude is often observed Že.g., Hoshiba, 1993; Yoshimoto et al., 1993; Harmsen, 1997., and attributed to source complexity Žsource mechanism, directivity of rupture, or asperity. ŽLiu and Helmberger, 1985. or directivity of site amplification factors. However, Figs. 7 and 8 show that the distribution of observed wave amplitudes looks like a normal distribution on a logarithmic plot, i.e., log–normal distribution. Amplitudes observed in seismic arrays have approximately log–normal distributions Že.g., Chang and von Seggern, 1980.. The large scatter observed in the wave amplitude might be easily explained by the

Fig. 7. Distribution of Ima x at x s 0.5, 8.0, 30.5, and 128 km. Incident Ricker wavelet has frequency f c s1 Hz, and the random media have Gaussian autocorrelation function with as 4 km, ² n : s 3 kmrs and n 0 s 0.05. Left and right sides show the same distribution, but linear–linear plot for left, and log–linear plot for right.

210

M. Hoshibar Physics of the Earth and Planetary Interiors 120 (2000) 201–217

Fig. 8. Same as Fig. 7, but for the distribution of maximum amplitude, which is the square root of Ima x .

sion of wave amplitude be? Two examples will now be considered.

Let us assume that we have a 3D velocity structure model whose resolution is quite small but which has p% error, and suppose that the background velocity is nearly 3 kmrs everywhere. In other words, variation in the velocity structure of less than p% is unmodeled. Even if we simulate the wave propagation and estimate the maximum amplitude or intensity Imax using the velocity model, the small velocity fluctuation introduces errors in the estimates of the wave amplitude. Suppose that we hope to predict the maximum intensity Žthat is square of maximum amplitude. for f c s 1 Hz at x s 30 km Žthat is, k c x f 63. with an error less than 70% which corresponds to m 2r being 0.5 f 0.7 = 0.7 on the assumption that the fluctuation has Gaussian autocorrelation function of a s 4 km Ž k c a f 8.4.. In this case, p must be less than 2% ŽFig. 4, Fig. 9.; if p s 5%, the m 2r exceeds the criterion ŽFig. 5, Fig. 9.. We need a velocity model with less than 2% error for a discussion of wave amplitude. Let us assume that we have a 3D velocity structure model whose mesh size is L and suppose that the background velocity is nearly 3 kmrs everywhere. This means that we do not know a velocity fluctuation of correlation length less than 2 L Ž a F 2 L.. Even if we simulate using the model of mesh

Fig. 9. Averaged m2r over five realizations of random media for each case of a s 2, 4, 8, and 16 km. Instead of x and a, dimensionless k c x and k c a are used. Left of this figure shows m 2r for random media with Gaussian autocorrelation function, and the right shows that for the exponential autocorrelation function.

M. Hoshibar Physics of the Earth and Planetary Interiors 120 (2000) 201–217

211

Fig. 10. Maximum m2r in the range a F 2 L, for random media of n 0 s 0.02 and 0.05 having Gaussian or exponential autocorrelation function, where L corresponds to the mesh size of the velocity structure model.

size L, the small fluctuation introduces errors. Variance of the error corresponds to maximum m 2r for a F 2 L, which is hereafter represented as M 2 , if the spectrum of the velocity fluctuation of length larger than 2 L is ‘white’. Fig. 10 shows relation between L and M 2 . This figure suggests how small L should be to estimate the intensity Imax accurately at a certain observation point. For example, when we hope to estimate Imax with an error less than 80% Žcorresponds to M 2 s 0.8 = 0.8 f 0.6. for f c s 1 Hz at x s 60 km Ž k c x s 125. in a medium whose fluctuation has an exponential autocorrelation function, L is not important if n 0 s 2%, but L must be less than 1 km Žcorresponds to k c L s 2.1. if n 0 s 5%. That is, we need a velocity structure model of mesh size less than 1 km, if the model has 5% error.

6. Conclusion The following conclusions are reached in this paper. Ž1. Even for fluctuations of velocity structure as small as several percent, the standard deviation of wave amplitude is larger than its average for high frequency waves. The relation of the velocity struc-

ture fluctuation with the amplitude fluctuation is evaluated, which suggests how small the error and how fine the resolution are necessary for velocity structure model for the discussion of wave amplitude. Ž2. Because of stochastic dispersion due to the randomness of velocity structure, the maximum amplitude decreases with distance for transient pulse wave even if there is neither intrinsic absorption nor back-scattering. The variance of wave amplitude increases with propagation distance, and then reaches a peak, and after that it decreases when the fluctuation of velocity structure is large.

Acknowledgements The author thanks the two reviewers whose comments were useful for revising the manuscript. A part of this study was done when the author was a visiting scientist of Ludwig -Maximilans-Universitat ¨ Munchen by the fund of the Science and Technology ¨ Agency, Japan. The author thanks F. Scherbaum who accepted him to the university and gave useful comments, M. Fehler for his careful reading of the manuscript and improving the usage of English, and

212

M. Hoshibar Physics of the Earth and Planetary Interiors 120 (2000) 201–217

Fig. 11. Snap shots of the waveforms at z s 0 obtained by the simulation using the parabolic approximation and phase screen method compared with the theoretical solution. Time interval of the snap shots is 1 s, and n 1 s 3.0 kmrs, n 2 s 2.4 kmrs, and f c s 1 Hz.

M. Hoshibar Physics of the Earth and Planetary Interiors 120 (2000) 201–217

Fig. 12. Same as Fig. 11, but for n 1 s 3.0 kmrs, n 2 s 2.4 kmrs, and f c s 2 Hz.

213

214

M. Hoshibar Physics of the Earth and Planetary Interiors 120 (2000) 201–217

Fig. 13. Same as Fig. 11, but for n 1 s 3.0 kmrs, n 2 s 3.6 kmrs, and f c s 1 Hz.

M. Hoshibar Physics of the Earth and Planetary Interiors 120 (2000) 201–217

Fig. 14. Same as Fig. 11, but for n 1 s 3.0 kmrs, n 2 s 3.6 kmrs, and f c s 2 Hz.

215

M. Hoshibar Physics of the Earth and Planetary Interiors 120 (2000) 201–217

216

S. Shapiro, K. Holliger, Y. Zeng, T. Iwata and H. Sato who gave useful suggestion.

The wave field inside is generally represented as, `

Ý Bn P Ž 2 n q 1. i n Pn Ž cos u . jn Ž k 2 r . ,

fi s C Ž f .

ns0

Ž A-3.

Appendix A The comparison between the simulation using parabolic approximation and a theoretical solution will be described in this appendix. Let us consider an example model containing a spherical velocity anomaly inclusion in an otherwise homogeneous structure. The homogeneous structure has velocity, y 1 s 3 kmrs, and the spherical inclusion has y 2 s 2.4 or 3.6 kmrs with radius r 0 s 2 km whose center is Ž x 0 , y 0 , z 0 . s Ž4 km, 4 km, 0.. We consider here an incident plane Ricker wavelet of f c s 1.0 or 2 Hz propagating into the qx direction, so that the size of the anomaly corresponds to k c a f 4.2 and 8.4. Here 0.25 km is used for D x for f c s 1 Hz, and 0.125 km for f c s 2 Hz. The theoretical solution is given as follows. The wave field outside the spherical inclusion is expressed as the summation of the incident wave and scattered wave: fo s f inc q fscat . The incident wave of frequency, f, and amplitude, C Ž f ., propagating in the qx direction is expressed as,

f inc s C Ž f . e i k 1Ž xyx 0 .

where k 2 s vrn 2 . A n and Bn are given from the boundary conditions,

fo s f i at r s r0 and fXo s fXi at r s r0 ,

Ž A-4.

where dash ŽX . means derivation for r, An s

X X jn Ž k 2 r 0 . k 1 jn Ž k 1 r 0 . y jn Ž k 1 r 0 . k 2 jn Ž k 2 r 0 . X X Ž1. h n Ž k 1 r 0 . k 2 jn Ž k 2 r 0 . y jn Ž k 2 r 0 . k 1 hŽ1. n Ž k1 r0 . X

Bn sy

X Ž1. jn Ž k 1 r 0 . k 1 hŽ1. n Ž k 1 r 0 . y h n Ž k 1 r 0 . k 1 jn Ž k 1 r 0 . X Ž1.X hŽ1. n Ž k 1 r 0 . k 2 jn Ž k 2 r 0 . y jn Ž k 2 r 0 . k 1 h n Ž k 1 r 0 .

Ž A-5.

After inverse Fourier transform we obtain the theoretical solution. Comparisons of the simulation using phase screen method with the theoretical solution are given in Figs. 11–14, where snap shots are shown with time interval of 1 s. Results of the simulation are shown with solid lines, and those of the theoretical solution are broken lines. Though they are different for the portions of wide-angle scattering or near strong velocity contrast, they are generally in good agreement.

References

`

sCŽ f .

Ý Ž 2 n q 1. i n Pn Ž cos u . jn Ž k 1 r . , ns0

Ž A-1. where k 1 s vry 1 , v s 2p f, r s wŽ x y x 0 . 2 q Ž y y y 0 . 2 q Ž z y z 0 . 2 x1r2 , cos u s Ž x y x 0 .rr, Pn is Legendre function and jn is the first kind spherical Bessel function. Here the factor eyi vt is omitted in Eq. ŽA-1.. C Ž f . is obtained from the Fourier transform of the Ricker wavelet ŽEq. 22.. The scattered wave outgoing from the spherical inclusion is generally expressed as, `

fscat s C Ž f .

Ý

A n P Ž 2 n q 1 . i n Pn Ž cos u . hŽ1. n Ž k1 r . ,

ns0

Ž A-2. where hŽ1. n is a spherical Hankel function.

Aki, K., 1969. Analysis of the seismic coda of local earthquakes as scattered waves. J. Geophys. Res. 74, 615–631. Aki, K., 1973. Scattering of P waves under the Montana Lasa. J. Geophys. Res. 78, 1334–1346. Aki, K., 1980. Attenuation of shear waves in the lithosphere for frequency from 0.05 to 25 Hz. Phys. Earth Planet. Int. 21, 50–60. Aki, K., Chouet, B., 1975. Origin of coda waves: Source, attenuation and scattering effects. J. Geophys. Res. 80, 3322–3342. Aki, K., Richards, P., 1980. Quantitative Seismology. Theory and Method Vols. 1–2 W.H. Freeman, 932 pp. Chang, A.C., von Seggern, D.H., 1980. A study of amplitude anomaly and mb bias at LASA subarrays. J. Geophys. Res. 85, 4811–4828. Chernov, L.A., 1960. Wave Propagation in a Random Media. McGraw-Hill, 168 pp. Fehler, M., Hoshiba, M., Sato, H., Obara, K., 1992. Separation of scattering and intrinsic attenuation for the Kanto-Tokai region, Japan using measurements of S-wave energy vs. hypocentral distance. Geophys. J. Int. 108, 787–800. Fisk, M.D., Charrette, E.E., McCartor, G.D., 1992. A comparison

M. Hoshibar Physics of the Earth and Planetary Interiors 120 (2000) 201–217 of phase screen and finite difference calculations for elastic waves in random media. J. Geophys. Res. 97, 12409–12423. Fisk, M.D., McCartor, G.D., 1991. The phase screen method for vector elastic waves. J. Geophys. Res. 96, 5985–6010. Harmsen, S.C., 1997. Estimating the diminution of shear-wave amplitude with distance: application to the Los Angeles, CA, urban area. Bull. Seis. Soc. Am. 87, 888–903. Hoshiba, M., 1991. Simulation of multiple scattered coda wave excitation based on the energy conservation law. Phys. Earth Planet. Inter. 67, 123–136. Hoshiba, M., 1993. Separation of scattering attenuation and intrinsic absorption in Japan using the multiple lapse time window analysis of full seismogram envelope. J. Geophys. Res. 98, 15809–15824. Hoshiba, M., 1995. Estimation of nonisotropic scattering in western Japan using coda wave envelope: application of a multiple scattering model. J. Geophys. Res. 100, 645–657. Ishimaru, A., 1978. Wave Propagation and Scattering in Random Media Vols. 1 and 2 Academic Press, San Diego, CA, 572 pp. Jannaud, L.R., Adler, P.M., Jacquin, C.G., 1991. Frequency dependence of the Q factor in random media. J. Geophys. Res. 96, 18233–18243. Jones, A.G., Holliger, K., 1997. Spectral analyses of the KTB sonic and density logs using robust nonparametric methods. J. Geophys. Res. 102, 18391–18403. Kennett, B.L.N., 1983. Seismic Wave Propagation in Stratified Media. Cambridge Univ. Press, Cambridge, 339 pp. Liu, H.L., Helmberger, D.V., 1985. The 23:19 aftershock of the 15 October 1979 Imperial valley earthquake: more evidence for an asperity. Bull. Seis. Soc. Am. 75, 689–708. Liu, Y.B., Wu, R.S., 1994. A comparison between phase screen, finite difference, and eigenfunction expansion calculations for scalar waves in inhomogeneous media. Bull. Seim. Soc. Am. 84, 1154–1168. Martine, J., 1992. Simulation of wave propagation in random media: theory and applications. In: Tatarskii, V.I., Ishimaru, A., Zavorotny, V.U. ŽEds.., Wave Propagation in Random Media ŽScintillation.. SPIE Press and Institute of Physics Publishing, pp. 463–486. Martine, J.M., Flatte, S.M., 1988. Intensity image and statistics from numerical simulation of plane wave propagation in 3-D random media. Appl. Opt. 27, 2111–2126. Martine, J.M., Flatte, S.M., 1990. Simulation of point source scintillation through three-dimensional random media. J. Opt. Soc. Am. 7, 838–847. Mayeda, K., Koyanagi, S., Hoshiba, M., Aki, K., Zeng, Y., 1992. A comparative study of scattering, intrinsic and coda Qy1 for Hawaii, Long Valley and central California between 1.5 and 15.0 Hz. J. Geophys. Res. 97, 6643–6659. Mclaughlin, K.L., Anderson, L.M., 1987. Stochastic dispersion of short-period P-waves due to scattering and multipathing. Geophys. J. R. Astron. Soc. 89, 933–963. Muller, G., Roth, M., Korn, M., 1992. Seismic-wave travel time ¨ in random media. Geophys. J. Int. 110, 29–41.

217

Roth, M., Muller, G., Korn, M., 1993. Velocity shift in random ¨ media. Geophys. J. Int. 115, 552–563. Sato, H., 1977. Energy propagation including scattering effects. Single isotropic scattering approximation. J. Phys. Earth 25, 27–41. Sato, H., 1984. Attenuation and envelope formulation of three component seismograms of small local earthquakes in randomly inhomogeneous lithosphere. J. Geophys. Res. 89, 1221–1241. Sato, H., 1989. Broadening of seismogram envelopes in the randomly inhomogeneous lithosphere based on the parabolic approximation: Southeastern Honshu. Jpn. J. Geophys. Res. 94, 17735–17747. Sato, H., Fehler, M.C., 1998. Seismic wave propagation and scattering in the heterogeneous Earth. Springer, New York, 308 pp. Shapiro, S.A., Kneib, G., 1993. Seismic attenuation by scattering: Theory and numerical results. Geophys. J. Int. 114, 373–391. Shapiro, S.A., Shwarz, R., Gold, N., 1996. The effect of random isotropic inhomogeneity on the phase velocity of seismic waves. Geophys. J. Int. 127, 783–794. Tatarskii, V.I., 1961. Wave Propagation in a Turbulent Medium. McGraw-Hill Book, 285 pp. Tatarskii, V.I., Ishimaru, A., Zavorotny, V.U. ŽEds.., Wave Propagation in Random Media ŽScintillation.. SPIE Press and Institute of Physics Publishing, 487 pp. Uscinski, B.J., 1980. Intensity fluctuation spectra in a multiply scattering media. Phys. Earth Planet. Inter. 21, 134–147. Wild, A.G., Hudson, J.A., 1998. A geometrical approach to the elastic complex-screen. J. Geophys. Res. 103, 707–725. Williamson, I.P., 1975. The broadening of pulses due to multipath propagation of radiation. Proc. R. Soc. Lond. A. 342, 131–147. Wu, R.S., 1994. Wide-angle elastic one-way propagation in heterogeneous media and an elastic wave complex-screen method. J. Geophys. Res. 99, 751–766. Wu, R.S., Aki, K., 1988. Introduction: seismic wave scattering in three-dimensionally heterogeneous Earth. Pure Appl. Geophys. 128, 1–6. Wu, R.S., Xu, Z., Li, X.P., 1994. Heterogeneity spectrum and scale-anisotropy in the upper crust revealed by the German continental deep-drilling ŽKTB. holes. Geophys. Res. Let 21, 911–914. Yoshimoto, K., Sato, H., Ohtake, M., 1993. Frequency-dependent attenuation of P and S waves in the Kanto area Japan, based on the coda normalization method. Geophys. J. Int. 114, 165–174. Yoshimoto, K., Sato, H., Ohtake, M., 1997. Short-wavelength crustal heterogeneities in the Nikko area, central Japan, revealed from the three-component seismogram envelope analysis. Phys. Earth Planet. Inter. 104, 63–73. Zeng, Y., Su, F., Aki, K., 1991. Scattering wave energy propagation in a random isotropic scattering medium: I. Theory. J. Geophys. Res. 96, 607–619.