Numerical Modeling for Elastic Wave Propagation With a Hybrid of the Pseudo-Spectrum and Finite-Element Methods

A numerical hybrid method was developed to model elastic wave propa­ gation. This algorithm was implemented with both the pseudo-spectrum and the finite-element methods. The pseudo-spectrum is currently a popu­ lar numerical method in earthquake seismology studies due to its high effi­ ciency and accuracy. On the other hand, its most significant drawback is the difficulty of implementing a free surface or absorption boundary owing to the nature of its periodic boundary. In addition, since the grid space must be defined globally within a model to prevent grid dispersion depend­ ing on the region of strong velocity contrasts, computations may become very expensive. However, these drawbacks can be overcome with a hybrid of the pseudo-spectrum and the finite-element techniques. With the imple­ mentation based on the finite-element formulation, grid spacing can be de­ termined according to local velocity within a velocity model. In so doing, the coding of the boundary conditions becomes much easier as well. The advantages of this proposed hybrid method consist of both reducing the amount of computational time and memory needed and obtaining both ac­ curate and stable results during calculation. Some examples are shown to demonstrate the advantages of the hybrid method. This method can also be easily expanded to 3-D situations with minor modifications. (

the difficulty of implementing a free surface or absorption boundary owing to the nature of its periodic boundary. In addition, since the grid space must be defined globally within a model to prevent grid dispersion depend ing on the region of strong velocity contrasts, computations may become very expensive. However, these drawbacks can be overcome with a hybrid of the pseudo-spectrum and the finite-element techniques. With the imple mentation based on the finite-element formulation, grid spacing can be de (Key words: Hybrid method, Finite-element method, Pseudo-spectrum method, Wave propagation)

INTRODUCTION
Modeling wave propagation in the earth plays a very important role in the study of the solid earth. Numerical modeling especially for seismic data has long been recognized as a powerful technique to enable seismologists to better understand complicated wave propaga tion and interactions and to use as a tool for data interpretation. In modeling seismic wave propagation globally and computing synthetic seismograms analytically, the earth structure is 1 Institute of Earth Sciences, Academia Sinica, Nankang, Taipei, Taiwan, R.O.C.
usually treated as a layered medium, for example with radial or vertical inhomogeneous struc tures. Rowever, compiex structures wn1Ch 'mC1uae varymg 1aterll i 1y· umomogeneous s ' lrucunes and elastic parameters are necessary for regional and local scale waveform modeling. Usually, numerical methods need to be employed. To solve the wave equation numerically, the models of a more realistic world should therefore be simulated.
Many different numerical algorithms have been proposed to solve wave propagation in a complex medium. Limited by the approximation of the wave equation, different algorithms are only suitable for a limited class of models. In the past two decades, different numerical methods have been developed to model seismic wave propagation with the commonly used techniques consisting of the finite-element method (FEM) (Smith, 1975), the finite-difference method (FDM) (Boore, 1972), the boundary integral equation method (BIEM) (Dravinski, 1983) and the pseudo-spectrum method (PSM) (Gazdag, 1981;Huang, 1992). Among these, the PSM requires the smallest required memory for computing wave propagation and has been well tested for scalar wave equations (Johnson, 1984).
Restrictions on all numerical methods are of a practical nature. With limitations imposed by computer memory size and processor speed, the earth models usually have to be truncated, and the frequency contents of a propagation wave are band limited. In this case, numerical modeling generally faces problems of artificial boundary and signal resolution restrictions. In the case of elastic waves, the instability of the free-surface boundary condition was found while the PSM solution was being used (Kosloff et al., 1984;Reshef and Kosloff, 1985). This is because, unlike the case of acoustic waves, the boundary conditions involve derivatives of the solution itself (R eshef et al., 1988 a, b ).
A hybrid technique uses one numerical algorithm in combination with another numerical or analytic method to compensate for the shortcomings of the individual numerical method. Zienkiewicz et al. ( 1977), Van den Berg ( 1988) and Emmerich ( 1989) have shown that hybrid methods are able to improve upon the flexible description with individual limitations. In this study, due to their robustness, the wave equation is solved with a hybrid approach using the FEM in the vertical direction and the PSM in the horizontal direction. This proposed hybrid method yields an algorithm which is non-periodic as well as which has variable grid spaces in the vertical direction. Furthermore, it reduces the unstable free boundary condition of the pure PSM.

Elastic Wave Equation Formulation
The basic equations.governing elastic wave propagation in the 2-D medium are described by the following: d2 Ux dTxx dTxz f" where (x,z) are the Cartesian coordinates, and d is the differential operator. (Ux, Uz) is the displacement vector, ( Txx, Tzz, Txz) is the stress tensor, t denotes time, p denotes the density, A, µ are Lame coefficients, and ft, f are body forces.
In this study, time marching is implemented based on the central-difference approximation. From the central-difference definition (Wylie, 1975), the physical quantity a1u I at2 of Equation (1) can be expressed as: where !it denotes the increment in time in numerical computations. Hence, the wave field at time t +!it is estimated from the values at times t and t -!it as follows:

Spatial Discretization
In order to numerically solve Equation (1), in this study, the spatial derivatives of the wavefield are calculated using the Fourier transform method in the horizontal direction and the finite-element approximation in the vertical direction. Following the Fourier method, the spatial derivative of the displacement field U(x) can be solved by a simple multiplication operation in the wave number domain as follows: According to finite-element discretization, the generalized displacement of the 1-D ele ment can be expressed as: where the terms N1=( 1-z/L) and N2=( z!L) are the shape functions of the displacement field. U1, U2 are the node displacements of element e. z is the Cartesian coordinate, and L is the length of the element. The derivative of ue with respect to z can be shown as: Finally, the global derivative matrix of displacement fields is obtained by assembling the element vectors in Equation (7). For an equal grid space of mesh, this formulation is equiva lent to the FDM.

Absorption Boundary Condition
In order to suppress numerical reflections from the artificial boundaries in the vertical direction, the absorbing boundary method proposed by Smith ( 197 4) for the finite-element method was used. This involved the superposition of the independently calculated wave equa tion solutions with the Neumann or Dirichlet condition. In theory, the summed wavefield completely eliminates any artificial reflections. Using the absorbing boundary method (Smith, 1974), the wave propagation allowed for the computation on a small vertical size by which the observations were distorted by numerical reflections without this boundary absorption. To avoid the wraparound at the horizontal artificial boundaries, an absorbing boundary condition based on the gradual reduction of wave amplitudes in the vicinity of the boundary (Cerjan et al., 1985;Sochacki et al., 1987) was used. The purpose of this reduction was to gradually taper the wavefield through the use of a weighting function W(x,z). For the case of absorption boundaries on the four edges of a model, the weighting function is in the form: where a is the absorption coefficient, d is the width of the tapered zone (expressed in the number of grid points), and L,: Lz are the model sizes of the x and z components.
In this section, four examples of elastic wave propagation are presented to show the flex ibility, efficiency and accuracy of the proposed hybrid method. With these cases, the calcula tion procedures are demonstrated, and the results are compared with other numerical solu tions. The advantages of the proposed hybrid method are shown as well.

Example 1: Comparison With Other Methods
In this example, the elastic wave propagation through a simple model of a homogeneous medium was computed. The velocities for P-and S-waves were chosen to be 2700 mis and 1500 mis, respectively. Density was taken to be 2.0 g/cm3• A Ricker-like time history with a maximum frequency of 50 Hz and a peak frequency of 25 HZ was used as the source wave form. The source, acting as an explosive type mechanism, was located at the center of the numerical grid. Figure 1 shows three snapshots computed by the FEM (Figure la), the PSM (Figure lb) and the hybrid method of this study (Figure le), respectively. It is seen that the three different numerical algorithms produced the same accurate results. In other words, the proposed hybrid method was able to generate the wavefields just as accurately as the other two well established numerical methods.

Example 2: Absorption Boundary Conditions
To simulate an infinite earth structure in a finite model, it is necessary to apply boundary conditions that make the edge of the computational grid appear transparent. In this way, artifi cial reflections introduced by the edges of the computational region can be minimized. In this example, a model with the same elastic parameters as those in the previous example was employed; however, instead of using a full-space model, absorption boundaries were imple mented on the bottom and top sides. To better absorb boundary reflections, the absorption boundary condition of Smith (1974) was applied to both the FEM and the hybrid methods. An ab sorption coefficient of a = 0.02 and the absorbing layer width d = 25 grid points of Equation (8) were selected to reduce the wavefield wraparound for the PSM. The results of this example are shown in Figure 2. The hybrid method (Figure 2a) shows just as stable results as the FEM (Figure 2b). However, in Figure 2c, the boundary reflections were not absorbed completely when the PSM was employed. To reduce the artificial boundary reflections of Figure 2c, the grid points in the taper bands had to be increased. To maintain the most stable results, it was found that 40-grid points per wavelength was necessary to attenuate boundary reflections similar to those assessable levels for waveform modeling used by Huang et al. (1995).

Example 3: Free-surface Boundary Conditions
In this example, the implementation for the free-boundary condition of the proposed method was demonstrated. It is known that the most significant drawback of the pseudo-spectrum algorithm is its failure to implement a free surface or absorption boundary due to the nature of its periodic boundary. In this example, the robustness of the hybrid method in dealing with the Snapshot computed by the pseudo-spectrum method. No special bound ary condition is considered in Fig. 3b.
free-surface boundary condition was shown. Herein, the same model as that of Example 2 was selected, but a free-surface boundary condition on the top side and absorption boundary condi tion on the bottom side of model were implemented. The results for the test of the free-surface boundary are shown in Figure 3. After careful examination, it is found that the results from the hybrid method (Figure 3a) have the same accuracy as those computed by the FEM. Because of the nature of its periodic boundary, no boundary conditions were considered in the computa tion of the snapshot of Figure 3b using the PSM. From Figure 3b, the same wavefront, which propagated away from the top surface, propagated into the model from the bottom surface again. To avoid the boundary transmission wave, usually, a thick zone of zero velocity was implemented by the PSM to simulate the free surface boundary (Reshef et al., 1988 a, b;Huang et al., 1995). Employing this method, however, requires extra computing costs and the simulation for free surface reflection is not exact.

Example 4: Non-uniform Grid Spacing
Grid dispersion limits the usefulness of the discretization schemes for the wave equation. In forward modeling problems, spatial grid spacing must be determined according to the me dium velocity so that the accuracy criterion can be satisfied. In general, the grid spacing for the PSM is defined globally to prevent grid dispersion. If there exist strong velocity contrasts in a model, then the grid spacing should be greatly reduced. In such a case, computation becomes very expensive. Nevertheless, the above drawbacks can be overcome with a hybrid of the ps eudo-spectrum and the finite-element techniques. Because the grid spacing of the model is allowed to change. When the finite element method is employed, the local grid spacing is estimated on the basis of the media velocities within the velocity model. The variable-size griding would help to reduce the computer time and memory needed. In this example, a two layer model was used to illustrate the usefulness of non-uniform grid spacing in the calcula tion. The velocity model is defined as a two-layer model in which the velocities of the P-and S-waves for the upper layer are 1500 mis and 1100 mis, respectively. The velocities of the P and S-waves for the bottom layer are 3500 mis and 2000 mis, respectively. Here, the densities on both layers are assumed to be 2.0 g/cm3• The source is located in the lower layer. Figure 4a shows the results of the hybrid method by using non-uniform grid spacing. In this model, the grid spacing is determined according to the velocities of the two respective layers. The wave propagation in both half spaces is stable and accurate. In comparison to the one in Figure 4a, the snapshot of Figure 4b was computed by the pseudo-spectrum method using uniform grid spacing which was determined according to the velocity of the lower layer. It is found that unstable wavefields are present near the velocity boundary (solid line). To keep the computa tions as accurate as those in Figure 4a, the uniform fine grids are necessary when the pseudo spectrum method is employed. However, because the grid spacing can be adjusted by the hy brid method, computer time and memory can be reduced.

DISCUSSION
As for the limitations of the Fast Fourier Transform (FFT), the input grid sizes of the PSM (both in the horizontal and vertical directions) are limited to a length that can be FFT'ed, i.e., 2". Although this can be partially improved by using the prime factor FFT algorithm (Temperton, 1985), in which the size of the FFT is the product of different prime numbers, the selection of the model size is strongly restricted by this requirement. Snapshot computed by the pseudo-spectrum method using uni form grids of the same size as those in the lower part of Fig. 4 (a). and FEM, eliminates artificial reflections and reducing computational costs. Compared to the conventional stand with the PSM and the FEM alone, the hybrid scheme greatly reduces the memory requirements and somewhat saves computer time. This is particularly true for earth models with strong velocity contrasts between layers, where the grid spacing in the vertical direction may be largely reduced. However, the hybrid method described in this paper is an initial study which can be used in lots of applications. For instance, the hybrid method may be used to model seismic wave propagation in the study of near surface structures, where a more complex model is needed to describe the real earth and where strong velocity contrasts exist between the weathering layer and structures beneath it. Additionally, low velocity guide wave modeling (Huang et al., 1995) may also be computed by this method. Additionally, the hybrid method proposed herein can be easily expanded to the 3-D case by using a hybrid of either the 1-D FEM and 2-D PSM or the 1-D PSM and 2-D FEM. The proposed hybrid method could be refined with further work. The computational algo rithm can be implemented with a parallel computing technique  to improve its computational efficiency. Also, the incorporation of a different hybrid computational method into the scheme would enhance its power, a valuable feature for processing wide-angle seismic reflection data sets which generally have sparse receiver points. For instance, the Kirchhoff integral can compute the passage of seismic energy through the low velocity surface layers, while the PSM or FEM method can be used to compute the propa gation in the structure beneath the low velocity medium.

CONCLUSIONS
A robust numerical modeling method for wave propagation has been proposed. This method hybrids the FEM and PSM to compute the spatial derivatives in the vertical and horizontal directions, respectively. The advantages of the proposed hybrid method consist of reducing the computational cost, reducing the memory needed and obtaining both accurate and stable results during calculation.