Nun1erical Modeling for Acoustic Scattering of 3-D Spherical Wavefronts: lntplications on Near Source Basin Amplification

To understand the complex damage pattern produced by the interaction between a 3-D sedimentary basin and 3-D spherical wavefronts, scattering of seismic waves by 3-D models due to a local point source is investigated by using the Pseudo-Spectrum method. The 3-D wavefields at different time steps are evaluated n111nerically for a corner diffraction and a sediment-filled basin model. 3-D wavefronts of both models are investigated from snapshots over free surface and vertical cross-sections. Numerical results show that model corners generate strong out-of-plane scattering energy which causes strong seismic energy focusing and defocusing in some specific locations. For an incident wave from a point source below the basin, the sediment-filled basin traps wave energy which propagates inward and focuses near the basin center. This ene rg)r extends to the basin b·ottom with decreasing amplitude. Besides, part of the incident energy is blocked by this basin, resulting in low amplitudes at the surrounding rock sites. This blocking effect is not predicted by the plane wave incidence. Comparisons are made with the results from 2-D models, and they show that the 3-D wavefronts from a point source or from in-plane scattering can be approximated by 2-D models; however, wavefronts from out-of-plane scattering cannot be reproduced.


INTRODUCTION
The recent occurrence of the very damaging Northridge earthquake (M=6.6, January, 17, 1994) in Calif omia, with hypocenter 14 km under the San Fernando Valley sedimentary basin, underscore the importance of understanding the interaction of 3-D wavefronts with a 3-D basin. The complexity of the damage pattern in the San Fernando Valley may be, at least in part, a result of the interaction of wavefronts with the basin boundary, which undoubtedly produces localized wavefield focusing, resulting in excessively strong shacking amplification. Prior to the Northridge earthquake, it had been found that data from events nearby .the Taiwan SMART-1 array could not be explained by plane wave incidence onto the base of the Langyang sedimentary basin (Abrahamson et al., 1987). Wen et al. (1994) have also found that the strong ground motion response of the 3-D Taipei basin differs by different source approaches.
. The phenomenon of scattering of seismic waves in heterogeneous medium has attracted the attention of seismologists and engineers for many years. Because of the complexity of the problem, exact analytic solutions are obtained that are restricted to basins of simple geometries (Trifunac, 1971;Wang and Trifunac, 1974). More general cases must be obtained by numerical methods. In the past two decades, different numerical methods for studying scattering effects have been developed. These include the 2-D Aki-Larner method (Aki and Lamer, 1970) and its 3-D extension (Ohori· et al., 1990), the finite difference method (FDM) (Boore, 1972), the finite element method (FEM) (Smith, 1975), the glorified optics method (Hong and Helmberger, 1978), the boundary integral equation method (Dravinski, 1983;Khair et al., 1989;Mossessian and Dravinski, 1990) and the Gaussian bean method (Cerveny et al., 1982). Some hybrid methods which combine the merits of two numerical methods have also been developed ( Van den Berg, 1988). Due to the limitations of computational capacity and storage memory of the computer, these methods are usually applied to 2-D cases. Generally, plane waves with different incident angles are employed (Aki and Lamer, 1970;Dravinski, 1983;Van den Berg, 1988). Only a few cases employ a point source or a finite fault rupture (Johnson, 1984;Benz and Smith, 1988;Dong and �cMechan, 1991). The rarely developed 3-D extensions of these methods are used to test the accuracy of their methods and only a few computed samples of typical 3-D models are found (Khair et al., 1989;Mossessian and Dravinski, 1990;Ohori et al., 1990). Recently, the finite difference method for 3-� models was applied to more realistic earth models (Frankel and Vidale, 1992;Frankel, 1993;· Yomogida and Etgen, 1993). These authors have concentrated on the frequency dependent site responses and time series features on the basin surface. Due to heavy computation, they have generally neglected to analyze the spatial, instantaneous seismic energy distribution inside the basin.
In this paper, the out-of plane scattering of acoustic waves from a local point source in a heterogeneous medium is considered. The present authors have developed a program which bases on the pseudo-spectrum method (PSM) to study acoustic wave propagation in a 3-D medium. The PSM was first developed by computational physicists in the 1970s (Gazdag, 1973;Orszag, 1972). The method is based on the Fourier method to calculate the spatial differentiation of a wave equation. The PSM was applied to explosion seismology in the early 1980's (Gazdag, 1981;Baysal, 1982, Kosloff et al., 1984;Reshef and Kosloff, 1985). It is found that a resolution of ten grid points per smallest period is required for the finite element method or the finite difference method in order to adequately model the wave phenomena, whereas only two grid points are required using the PSM (Fornberg, 1987). The advantages of saving storage memory and getting high frequency resolution are the features of this method. The .method, however, is more efficient when used to calculate 3-D problems than 2-D problems (Fomberg, 1987). Supercomputers (the CRAY Y-MP and CRAY 2) were used to execute the numerical code here, its accuracy having been confi11ned.

253
Compa ri sons with 2-D models are also provided. Nu merical results presented here show · the effects of the point source excited 3-D wavefi eld which are not reprodu ced by a plane wav e inpu t moti on or 2-D models. Thu s, to predi ct the ground motion inside a 3-D ba sin excited by an earthquake ri ght below the basin, the consideration of a point source scattering fi eld is necessary.

For111 ulation
The basic equation governing acoustic wave propaga ti on in 3-D medium is as follows: where t denotes ti me, p( x, y, z) denotes the density, K ( x, y, z) is bulk modules, P( x, y, z, t) denotes the pressure and S(x, y, z, t) denotes the source term (Kosloff and Baysal, 1982). If the density p along the x and y axis is constant, the spatial derivative can be simplifi ed and equation (1) where S' ( x, y, z, t) is the source term ex pressing the di -v ergence of the pressure, and C(x, y, z) is the velocity of the medium. Equation (1) is more general tha n equation (2) bu t involves more nu m erical calcula tion for wave propa ga ti on. To redu ce the numerical computati on, for some cases of constant densities along the x and y axes, equation (2) wa s employed here to compute wave propaga ti on.

Numerical Solution Technique
In order to numerically solve equation (1) or (2), a di screti zati on in spa ce and time is perf ormed. In thi s study, the time marching based on the central difference approxima ti on is •• ca lcu lated. From the centra l difference defi ni ti on (Wylie, 1975), the physical quantity P(t) is ex pressed as:

••
where double dots ex press a double di fferentiation with respect to ti me. The P( t) deduced from equation (2) can be shown to be: a2 P a2 P a2 P · aP 1 aP · {}t 2 x y z p z (4 ) Hence, the wave field at ti me t+ � t is given from the values at time t and t-� t as follows: After the spatial Fourier transformation is applied, the spatial derivation can be evaluated by multiplication in the wave number domain. Thus, is the Fourier transform of P( x ). N is the number of data points, and k is the wavenumber.

+ p -oz p oz
where P( x x , Yy, z) is the 2-D Fourier transform of the P( x, y, z) with respect to x and y.

Boundary Conditions
To avoid the wrap-around at the artificial boundaries, some absorbing boundary con ditions based on the gradual reduction of wave amplitudes in the vicinity of the boundary (Cerjan et al., 1985;Sochacki et al., 1987) have been proposed. Cerjan et al. (1985) suggested a Gaussian type damper for 2-D acoustic wave cases. Sochacki et al. (1987) pro posed 5 different types of dampers (linear, exponent, cubic, exponential and Gaussian) with different dampers for different applications. These suggestions are based on 2-D numerical testing. For effi cient absorption of the numerical reflection from the artificial boundaries of I the 3-D acoustic wave equation, the present authors have rigorously checked the 5 types of dampers with different absorption coefficients. No damper for 3-D problems have been discussed, before. In this study, a damper for 3-D cases have been rigorously tested for its numerical behavior.
In this study, the free boundary is approximated by including a wide zone with zero velocity above the upper surface .of the model to simulate the traction-free condition and to avoid the wraparound propagation waves from the bottom of the model (Reshef et al., 1988 a, b).

Source Time Function
A time dependent pressure change at a single node is used as a point source. The source time function (Figure 1) is taken as a first directive of a Gaussian function as: (Equation 9). Herein, a=3500 and to=O.l sec are used.

( 9 )
where to is the central time of a wavelet, and a is a coefficient governing the time interval from the negative to positive peak of the function. The Fourier transfo11n spectrum of equation (9) is: and is shown in Figure 2. This choice is made as a band-limited spectrum in the frequency domain; however, it i· s suitable for numerical computation .

Coding and Testing
The 3-D acoustic wave propagation program is developed under the UNIX based op erating system. Te sting is made on a DEC-3100 workstation, an ETA-10, an IBM-3090 and CRAY series supercomputers, and all have produced stable results. The perfo1·mance of the PSM is dependent on the speed of the fast Fourier transf or111 (FFf) which is very suitable for parallel and vectorization procedures (Temperton, 1983(Temperton, , 1985. However, the present generation of array and vector processors are quite adequate for computing the FFI'. Usually, general scientific computers support high perfo1·1nance 1-D and/or multiple FFI' sys tem libraries. With these library routines being used, the computation speed of the code in this study is greatly improved. Different FFI' codes from the scientific libraries in different mechanisms are compared. The multiple FFI, routine of the CRAY Y-MP has the best per formance in executing the code. To check the numerical accuracy of the PSM, the calculated results are compared with other methods (Huang and Yeh, 1991;Huang, 1992). One of them is shown in Figure 3. The comparison of the results of the PSM and the Cagniard-deHoop method (Daudt et al. , 1989) illustrates that there are virtually no differences between the resulting waveforms.

Corner Diffraction Case
The proposed model (see Figure 4) includes two media with a low velocity and low density (1.5 km/sec and 2.2 g/cm 3 , respectively) cubic comer inside an otherwise homoge neous medium with a uniforrn velocity of 2.6 km/sec and a density of 2.4 g/cm 3 . The model size, which is represented by grid numbers, is set at 128 x 128 x 128 with space intervals (�x, 6y and �z) of 200 m. The low velocity comer has dimensions of 1.3 km both on the x and y axes and 1.1 km on the z axis. Over the free surface (z= 1.66 km), a fictitious layer has been placed with a velocity of 0.01 km/sec on top of the model which numerically sim ulates a free boundary. An explosive type point source is located at 0.1 km above the comer with coordinates of x=I.3 km, y=l.3 km and z=l.2 km. The source time function is shown as equation (9) with a=3;00 and t0=0.05 sec. The computed 3-D wavefield snapshots at time intervals of 0.40 and 0.44 seconds after the explosion are shown in Figure 4, where the snapshots of the 3-D wavefields as equal amplitude surface have been displayed as shown. In order to show the reflection inside the direct wavefront, the wavefields with a y coordinate of less than 0.8 km are made transparent in Figure 4. Some propagation phases are clearly seen in the snapshots. Two selected vertical wavefield cross-sections at time 0.44 sec are shown in Figure 5. In the cross-section of y=0.8 km ( Figure 5 (a)), which includes the low velocity comer and does not include the source, the phases defined in Figure 4 are clearly seen. It is also evident that, due to the source location above the comer, the transmitted waves from the upper comer surface (phase Ez in Figure 4(b)) are stronger than those from the vertical comer surface (phase Ex in Figure 4(b)). Both Ez and Ex have crossed each other causing a transmitted curved wavefront (phase F in Figure 5 (a)) behind. In the vertical cross-section of y= 1.6 km ( Figure 5 (b) ), which is in the homogeneous medium behind the low velocity comer, the diffraction from the upper surface of the low velocity comer and its free surface reflection (phases C and D) are smaller in amplitude than the direct and reflection waves (phases A and B). The strength of phase D decreases in amplitude farther away from  257 the diffraction comer. The horizontal wavefield cross-sections of z=0.8 km and z=l .6 km are shown in Figure 6. In the cross-section of Figure 6 (a), the major energy transmitted into the low velocity region is from phase Ez. The two phases (Ex and Ey in Figure 6 (a)) which enter this cross-section horizontally from both vertical comer surfaces are identified. The same as in the vertical cross-section, both Ex and Ey in Figure 6(a). have crossed each other causing a curved wavefront behind (phase F). The velocity ratio between the two media measured from the 2-D cross-section (assumed to be the source located in this cross-section) is nearly 0.4 which is different from that for the real velocity ratio of both media (0.576). The difference comes from the source location of the 3-D model far away from this cross-section The apparent velocity ratio on both media is less than the real velocity ratio in this vertical profile. The stratigraphic wavefronts of phases Ez, Ex and Ey are exam ined from the individual 2-D cross-sections of the 3-D wavefields. Ez has the shape of a quarter cone with its axis in the z direction, while both Ex and Ey have the shapes of a quarter cone with their axes in the x and y directions, respectively. The wavefront of phase F is found to be a quarter spherical shape. The reflection waves (phases C and D in Figure 4) are found on the cross-section above the low velocity comer (Figure 6(b)). The vertical axis is defined as the depth from free surface. The definitions of the symbols are the same as those in Fig. 4(b). The amplitude of the spatial wavefield of each panel is normalized individually and the same as the other cross-sections in this study.

259
The 3  . The solid lines ex press the ou tline of the basin. The wavefields with a coordinate y less than 0.8 km are transparent. R: the artificial boundary reflection pha se.

261
and a density of 2.3 g/ cm 3 wi th a good contra st to the bedrock of velocity (2.6 km/ sec) and densi ty (2.4 g/ cm 3 ). The ellipsoid has a 3:2: 1 ratio of long (x), inte11 nediate (y) and short (z) axes. The source used in thi s study is a point source just below the basin (a s shown in Fi gu re 7 ) to s i mu l ate a near-bas i n earth qu ake. The strong artificial boundary reftection is found in the propaga ting wavefield as shown in Fi gure 8. It wa s found tha t the Gaussian damper best redu ced the artificial reflections. The Gaussian damper is ex pressed in the fo1m of exp(-a d2 ), where a is the absorpti on coefficient, and d is the di stance (expressed in gri d points) from the inner boundary of the absorbing layer. In thi s study, ma ny absoiption pa ram eters for study, many absorption parameters for the 3-D acoustic wave equation have been tested, and it has been found that an absorption coefficient of a=0.025 coupled with an absorbing distance of 17 grid points had reduced wraparounds to an acceptable level as shown in Figure 8(a). Comparative results without absorption boundaries are shown in Figure S (b). No optimum values of absorption parameters are obtained in the testing procedure here. Larger values of d coupling with smaller values of a of the absorption parameters produce better absorption of artificial reflection. For computation efficiency, a special test is required for the particular problem of interest.
The 3-D wavefront snapshots of the experimental model at later time steps than in Figure 8 are shown in Figure 9. The wavefront of the direct wave (phase A), the basin bottom reflection (phase B) and the transmitted wave and its reflected phase from the free surface (phase C) are observed. The special wavefront (phase D) of the acoustic wave trapped on the basin is clearly seen in Figure 9(b). The plane snapshots of the free surface wavefield at different time steps are plotted in Figure 10. Two interesting features are observed in these snapshots. First, from the blocking of this low velocity basin, the incident wavefront is disconnected on both sides of the basin boundary, and the incident wave decreases its amplitude as it approaches the basin boundary. Tracing the wavefront in time, there are two regions of the basin boundary that give lower peak amplitude than the incident wave ( Figure  I 0). Second, the basin traps energy that propagates inward and focuses energy near the center of the basin of this model. Two vertical cross-sections of 3-D snapshots are shown in Figure 11. The trapped wavefield is clearly seen in the cross-sections and on different time steps. Phase (D), as trapped energy in the basin, extends to the basin bottom with decreasing amplitude from top to bottom. As phase (D) hits the basin boundary, a strong reflection is generated that keeps most of the wave energy within the basin, much the same as a wave excited in an ellipsoidal swimming pool. Reflections from all parts of the basin boundary approach each other f or1ning regions of wave focusing in a manner dependent on the nature of the basin boundary. Two subsurface horizontal snapshots with different elevations are shown in Figure 12. The wavefield inside the basin includes the horizontal propagation trapped waves (Figure 12(a)), while outside (Figure 12(b)) the basin, no trapped waves are found.

Cross-section Related to 2-D Cases
The 2-D models, with an infinite extension of the third dimension of the 3-D experi mental models, are selected to compute the wavefield for comparison with the 3-D model. Two cross-sections related to Figure 5 are chosen, and the computed wavefi eld of the comer diffraction model (not shown) is very similar to its related 3-D cross-section ( Figure 5(a)); however, the diffraction phases C and D of Figure 5  :. · wavefi eld of Figure 13 (a) is very different from that of the related 3-D cross-section ( Figure  6 (a)). The amplitude of the wavefield in the low velocity comer is stronger than that in the other region; besides, phases Ez, Ex and Ey are not found. Phases Ex and E11 can be reproduced by shifting the source slightly away from the comer as shown in Figure 13 (b ). Phase Ez cannot be reproduced because it propagates into this cross-section· vertically, and the out of plane wavefield cannot be simulated by a 2-D model. The same reason explains , why phases B, C and D of Figure 6 (b) cannot be simulated by a full-space 2-D model. o.   Figures 11 ( c) and ( d), the 3-D snapshot cross-sections are located at x= 1.6 km, and the source is far away from the plane (Figure 7). Due to the lateral heterogeneity in the x direction, even though the shapes of the basin for both cross-sections are similar, the spatial distribution of these wavefronts are very different. This shows the limitation of using 2-D models for a 3-D interpretation. Besides, no corresponding 2-D models can be used to approximate the wavefield snapshots for the free surface or for any horizontal cross-sections produced by the 3-D model (Figure 1 2).

DISCUSSION
In thi s study, through 3-D numerical computati ons, it may be understood that excessively stro ng ground moti ons can be generated by the localized wavefield focusing effect as a result of the interacti on of a 3-D basi n with incoming 3-D spherical wavefronts. Until recently, there were only a few cases for 3-D basin resp onse comp utati ons which mainly involved plane wave incidence. Fi nite fault and point source effects are less often di scussed. Generally, di scussion on site eff ects is sep arated from that on source eff ects, and it usually neglects the wave scatteri ng eff ects induced by di stance variations from source to basir.. However, when the source is close to the basi n, it has been found that, even in such simp le models as those used in thi s study, the wavefields are strongly di sturbed by the configurations of the 3-D direct application of this method at this time. With the increased availability and speed of supercomputers, solutions to this computational problem will be overcome in the near future . •

S. CONCLUSIONS
The numerical findings discussed above lead to the following conclusions: (1) Consideration of the out-of-plane scattering of spherical wavefronts and calibration for apparent velocities is necessary when 2-D results are being applied to infer 3-D cases.
(2) The wavefront-blocking effect produced by the sediment-filled basin gives rise to low shacking amplitudes which cannot be reproduced by a model with a plane wave inci dence.
(3) A basin trapped wave energy is found to propagate inwardly and focus near the center or the elliptical foci of the basin. This trapped wavefield is examined in different vertical cross-sections and time steps. The basin trends to trap the wave energy inside which results in much lengthened strong shacking motions. This trapped wave extends to the basin bottom with decreasing amplitude. (4) Site effects are strongly affected by the point-source-generated spherical wavefronts which interact with the basin boundary, causing very complex strong shacking patterns through wave focusing. This focusing and soft basin trapped wave energy essentially produce the large-amplitude and the long-duration strong ground shacking that usu ally infl icts severe damage to structures, such as in the case of the recent Northridge earthquake.