Seismic Wave Propagation in Basin Structure from Numerical Modeling

1 Institute of Seismology, National Chung Cheng University, Chia-Yi, Taiwan, ROC 2 Institute of Geophysics, National Central University, Chung-Li, Taiwan, ROC 3 Department of Natural Sciences, National Science Council, Taipei, Taiwan, ROC * Corresponding author address: Dr. Chau-Huei Chen, Institute of Seismology, National Chung Cheng University, Chia-Yi, Taiwan, ROC; E-mail: seichen@eq.ccu.edu.tw A sedimentary basin can significantly effect the amplification of earthquake ground motions; however, amplification at any given site is dependant upon the locality of the earthquake’s epicenter. A knowledge of basin response plays an important role in the study of probabilistic seismic hazard in basin areas. It is important to know average amplification and intrinsic variability of each site in a given basin area. To this end, we have applied the spectral element method to estimate three-dimensional seismic responses for various basin structures including: (1) a basin with heterogeneous sediments surrounded by bedrock; (2) a basin with complex topographic features at its surface; and (3) a basin with a thin near surface soft layer. Incorporating a double-couple point source, we simulate the ground motions of waves propagating through various numerical basin models. The results obtained in this study serve as a guide to expectations; in particular, they help in understanding what can occur when waves propagate through basin structures containing increasing amplification factors at sites located above the deeper parts of the basin. Another important task of this work is to test the accuracy and the stability of the spectral element method. By comparing waveforms generated by the spectral element method and the non-uniform grid spacing finite difference method developed by Pitarka (1999), we find the waveforms calculated by both methods are in remarkably good agreement; and more importantly, the main characteristics of the wavefield are well preserved. (


INTRODUCTION
Wave propagation simulation incorporated with 2D or 3D basin velocity structures, even using plane waves or line/point sources still can have significant basin effect (Frankel 1993;Yomogida and Etgen 1993;Scholler et al. 1994;Graves 1995;Olsen and Schuster 1995;Olsen et al. 1995a, b;Olsen and Archuleta 1996;Komatitsch et al. 2004).Basin effect is enhanced when seismic waves are unable to dispel quickly.The result of such energy resonance can cause a great deal of damage buildings and infrastructures (Vidale et al. 1991).Basin and topographic effects include the generation of surface waves at basin edges have strong correlation between the duration of ground shaking and site location within a basin (Olsen and Archuleta 1996;Komatitsch et al. 2004).
Besides the basin effect, a soft layer at shallow depth can also make energy difficult to dispel (Olsen et al. 1996).Liquefaction may occur when seismic energy passes through basin soil structures containing a good deal of loose sediments (Rollins et al. 2004).In addition, topographic scattering can affect peak particle amplitude and signal duration (Olsen et al. 1996).In order to better understand these effects, this research adopts the spectral element method (SEM) to simulate seismic wave propagation through basin structures.
The SEM was first put forward in the 1980s (Babuska et al. 1981;Babuska and Dorr 1981).Now the method is widely used for solving problems in hydrodynamics (Patera 1984;Maday and Patera 1989;Fisher 1990) and to model seismic-wave propagation in local and regional models (Cohen et al. 1993;Priolo et al. 1994;Komatitsch 1997; Komatitsch and Vilotte 1998;Komatitsch and Tromp 1999;Komatitsch et al. 2004).It was more recently introduced for simulation of global wave propagation by Chaljub (2000) and also extended for large-scale global wave propagation (Komatitsch and Tromp 2002a, b;Capdeville et al. 2003;Chaljub et al. 2003).This research utilizes the method to simulate seismic wave propagation through various basin structures, including a basin with complicated surface topography, and distributes through out the 3D medium and various geometric interfaces within the basin.The SEM is a kind of high-order differential equation in a physical field.This method is an efficient tool for the study of various types of fluctuation characteristics in a 3D structure.One important reason for adopting this method is that in order to exhibit complicated wave phases generated by seismic waves traveling through geological structures, a precise method is usually needed to simulate boundary conditions, and irregular and non-homogeneous interfaces, given the possibility of structural comprise strong velocity contrast and topographic effect.These effects cause the ground to shake in a complicated manner, and have a decisive influence on the structural design of buildings.In the past, the finite-difference method (FDM) was once widely applied in the calculation of seismic wave propagation through disparate media.However, when close to the big change of wave field, or coarse mesh operation will all cause "grid dispersion" (Alford et al. 1974;Kelly et al. 1976;Robertsson 1996).Under free surface conditions, one drawback of the FDM is that a lack of flexibility in geometry means that precision cannot be maintained.On the other hand, the application of the finite-element method (FEM), has its problems as well where low-order FEMs give serious dispersion (Marfurt 1984); and high-order FEMs can yield spurious waves (Sarma et al. 1998).A final disadvantage of using the FEM is that it is time consuming.In order to overcome the above problems, we adopt a SEM to perform simulation of ground motions, which has flexibility in geometry and precise fast time marching (Faccioli et al. 1997;Komatitsch 1997;Komatitsch and Vilotte 1998;Seriani 1998;Komatitsch and Tromp 1999;Paolucci et al. 1999).The SEM can also yield significant improvements in calculation efficiency, especially when ground motion simulation is carried out on parallel machines.Another significant benefit of using the SEM is that it converges (Fisher 1990) rapidly, which reduces computation time.Because of the above advantages, we applied the SEM to perform seismic ground motion simulations with four different basin structures.

NUMERICAL TESTING FOR ACCURACY AND STABILITY
When utilizing the SEM to simulate elastic wave propagation, Lagrange interpolation and Gauss-Lobatto Legendre quadrature formulas are selected (Komatitsch 1997;Komatitsch and Vilotte 1998;Komatitsch and Tromp 1999).Details about how the mathematical formula is constructed and the procedure of spatial and temporal discretization have been discussed fully in several papers (Komatitsch 1997;Komatitsch and Vilotte 1998;Komatitsch and Tromp 1999;Komatitsch and Tromp 2002a, b).One point that we need to emphasize is that due to the consistent integration scheme and the use of Gauss-Lobatto Legendre formulas, a significant characteristic of the method is that the mass matrix is always diagonal, thus making computation more efficient.
In order to demonstrate the accuracy and stability of the SEM, we performed several simulation tests to model wave propagation in laterally homogeneous and heterogeneous media.The results are compared with simulations obtained by 3D FDM with non-uniform spacing (Pitarka 1999).The reason for this comparison is that the testing models used in this research (see Fig. 1) are simple, which include: (1) the homogenous half-space model (denoted by TEST1); and (2) two-layered models (denoted by TEST2 and TEST3).The location of the source (star) and the stations (inverted triangles) are also shown in Fig. 1.For simplicity and reducing finite fault effect, we calculated velocity seismograms for a double-couple point source with a focal mechanism with strike = 30°, dip = 80°, and rake = 30°.A bell-shaped source-time function with duration time, T 0 = 1.0 sec, is applied.The source time function is written as: where T 0 is the duration time.The mesh is composed of 64 × 64 × 23 elements.With a polynomial approximation of order N = 5, the total number of collocation points is 144 -461 and the average distance between points is 0.125 km.The physical parameters of the model TEST1 are (see Fig. 1a): Vp = 4.0 km sec -1 , Vs = 2.3 km sec -1 , ρ =1.8 g cm -3 , and the depth of earthquake source is 3 km.For the discretization of the model, Pitarka (1999) used 0.125 km as the increment of horizontal grid spacing.However, for the vertical grid spacing, if the depth is less than 2 km, the increment is 0.125 km; otherwise the increment is 0.25 km.By using the above criteria, the model has approximately 6 and 12 grid points per shear wavelength.Based on this discretization, the resolution of the model is approximately 1.5 Hz.One of the important objectives of this study is to examine the accuracy and stability of the SEM.Thus, for model TEST2 and TEST3, we will carefully investigate the accuracy and stability of the SEM, but for model TEST1, we only examine the accuracy of the SEM.The synthetics are bandpass filtered between 0.1 and 1.5 Hz, using a Butterworth filter with 4 poles and the taper being 0.05. Figure 2 illustrates the results calculated from the SEM (lower solid line) and the FDM (upper solid line) for model TEST1.As we can see, the waveforms obtained by both methods fit remarkably well.For further testing, we double (dashed line) and quadruple (dotted line) the grid spacing.In other words, the number of elements is decreased.The waveforms calculated by the SEM and the FDM also exhibit high agreement, except for the micro oscillations in coda wave.
In the model TEST2 (Fig. 1b), Pitarka (1999) included a soft thin layer on the top of the half-space.Since we will compare results with those of Pitarka (1999), the synthetic seismograms are bandpass filtered between 0.1 and 1.3 Hz as in Pitarka (1999).In Fig. 3, we note that waveform agreement between these two techniques is still very good.According to the above comparison, the SEM can handle models with strong velocity contrast.And more importantly, the SEM also demonstrates excellent accuracy in waveform simulation.For the stability test, we also double and triple the grid spacing for spatial discretization of the model.The results indicate that for double grid spacing, the main phases generated by the SEM can be clearly observed (dashed line).By contrast, for quadruple grid spacing, fitting of the waveform in the later part of the seismogram is not good (Fig. 3: dotted line).Except for testing different grid spacing, we also change the thickness of the soft layer in model TEST2 to examine the stability of the SEM.Noticeably, the main characteristics of the waveform are well preserved, even when the thickness of soft layer is reduced by half (Fig. 4).Model TEST3 (Fig. 1c) tests the response of surface sediments with an earthquake occurring at shallow depth.The earthquake source is located at a depth of 0.5 km.The comparison between the synthetics generated by the    FDM and the SEM are shown in Fig. 5.Both synthetics are filtered between 0.1 and 1.3 Hz.In Fig. 5, the lower solid lines represent the waveform calculated by the SEM, while the upper solid lines indicate the synthetics of the FDM.The waveforms generated by both methods fit very well.One feature worth noting is that the surface waves remain trapped inside the softer surface layers which develop complicated waveforms with relatively large amplitudes.For the stability test, we double the grid spacing in the simulation for the SEM (Fig. 5: dashed line), the waveform fitting between the FDM and the SEM shows high consistency.However, when we quadruple grid spacing, the waveforms generated by the SEM become complicated and at the later parts of the waveform do not fit well between the two methods (Fig. 5: dotted line).
We also change the thickness of soft layer in model TEST3, when the thickness of soft layer becomes half, the amplitudes of waveforms are larger and the wave phases become more complicated (Fig. 6).

A Basin Model with Strong Velocity Contrast
In this section, we extend our analysis to a basin model with strong velocity contrast between the basin and the bedrock.The basin with semi-ellipsoidal shape, and the locations of the hypocenter (star) and stations (inverted triangles) are shown in Fig. 7.The velocity structure of the basin is: Vp = 1.9 km sec -1 , Vs = 0.8 km sec -1 , and density ρ = 2.0 g cm -3 , while the velocity structure of the bedrock is: Vp = 5.6 km sec -1 , Vs = 3.2 km sec -1 and density ρ = 2.2 g cm -3 .The synthetic velocity seismograms are calculated and filtered (frequency range: from 0.1 to 1.3 Hz) at receivers located on the free surface along a line across the basin.Figure 8 illustrates the synthetic velocity seismograms which are calculated by using a double-couple point source with focal mechanism of strike = 30°, dip = 90°, and rake = 30°.A bell-shaped source-time function with duration time of 0.8 sec is used.The seismic moment: M 0 = 10 25 dyne cm.The source is located at a depth of 10 km.The resulting waveforms clearly illustrate the develop- ment of complexity of ground motions as the waves propagate through the basin structure (Figs.8a, b).At stations 7 through 12, our results indicate that the most significant ground motion amplification is above the deepest part of the basin along the X-X' (Fig. 8a: phase is indicated by the arrow a).This is partially due to large-amplitude surface waves generated by the edges of the basin (Fig. 8a).The most prominent feature is that the direct S wave can be clearly identified on the Y component (phase is indicated by the arrow b).Due to the shape of the basin, the S wave is deflected.In Fig. 8a, the phase, indicated by arrow c (on the X component), demonstrates especially large amplification effect near the edge of the basin.This may have been caused by basin-induced diffracted waves and surface waves generated at the edge, that travel in phase causing constructive interference.However, in Fig. 8b, the phase, c, is not significant; this may have resulted from the effect of the basin edge not occurring at this range of distance along Y-Y'.Edge effect was a key factor damaging many buildings in Kobe during the Hyogo-ken Nanbu earthquake of 1995.

A Basin Model with Added 3D Hill
A 3D hill of height 500 m is added to the top of the basin model (Fig. 9).Waveforms are shown in Figs.9a and b.Compared with the results obtained by the previous basin model, phases, indicated by arrows a, b, and c, can also be observed along the X-X' and Y-Y' cross sections (see Fig. 9a).The amplification pattern on the hill exhibits strong variability.The maximum amplification occurs in a small area around the summit or at the mid-slope of the hill, while amplification occurs at the hill-top area.For one profile along the Y-Y', strong   amplification is observed around the top of hill (Fig. 9b).For the z component, arrow d indicates that the S wave converts to a P wave at the bottom of the basin (see Fig. 9a).After t = 9 sec, the wave field becomes more complicated (Figs.9a, b).This is caused by the coupling effect between the hill and the basin structure, or residence time for constructive interference of the diffracted waves within the hill structure.In conclusion, all significant energy seen in Figs. 8 and 9 is contributed by basin-generated phases.

A Basin Model with 3D Concave Valley
Contrary to the model used in section 3.2, we simulate ground motions with a basin model that has a 3D concave valley in the top section (Fig. 10).From our results, we can clearly observe the phases, indicated by arrows a, b, c, and d, described in section 3.1 (Fig. 10a).The results of wavefield simulations are associated with strong impedance contrast between the basin and the concave valley, which produce deamplification effects.The S-wave multiple arrivals created by reflection between the free surface and bottom of the basin are also observed.Since this model is has strong lateral variation, body waves are warped at these interfaces and have constructive interference near the surface.At the edge of the basin, the wave field is dominated by trapped Stonley type waves propagating along the basin boundary that are excited by incident waves to the basin as well as reflected waves from the basin (Figs.10a, b).Hill and Levander (1984) and Olsen et al. (1995a), among others, showed that reverberations occurring in a near-surface with low-velocity layer in a basin can significantly increase the duration time of a signal.Therefore, we add a near-surface thin soft layer with thickness of 500 m.The velocity structure for this thin layer is Vp = 1.65 km sec -1 , Vs = 0.41 km sec -1 , and ρ = 1.8 g cm -3 .We apply the same attenuation constant, Q = 0.05 × Vs, as in the study of Olsen et al. (2003).3-component synthetic seismograms along the X-X' and Y-Y' are shown in Fig. 11.We discover that the primary effect of adding this thin layer is to delay the arrival time of energy, which propagates through the low velocity sediments.The phases represented by arrows a, b, and d can be clearly seen.It is worth noting that at t = 6 sec, the complexity of the wave field begins to emerge.We also observe that the wave packet in the central part of the basin is trapped in the very shallow layers of sediment.As shown in the X-component of Fig. 11a, stations 11 ~ 14 exhibit very obvious wave packets, which are near the basin edge.Therefore, the effect of adding a near-surface low-velocity layer is to increase the peak of particle amplitude by as much as a factor of 3 (compared to a basin surface not including a thin layer) and also increase the duration time of the signal (see Figs. 11a, b).Another phenomenon is the surface wave generated by this model having a larger amplitude along the Y-Y' (Fig. 11b).Our result is consistent with the findings of Olsen et al. (1995a).According to their findings, mode conversion and resonance are probably the main factors that are responsible for generating larger amplitude ground motions when the basin model consists of a near-surface lowvelocity layer.Impedance effects and surface wave generation are other factors that may contribute to the amplification of ground motion.

DISCUSSIONS AND CONCLUSIONS
The main purpose of this study is to test the stability and accuracy of the SEM.In this study, we demonstrate that the SEM is an efficient and sufficiently accurate tool for modeling wave propagating in a 3D elastic media.We have applied the SEM to calculate ground motion with a flat-layered model and more complicated basin models with distinct topographic features.For accuracy and stability tests, by comparison with the FDM, the waveforms generated by the SEM and the FDM fit very well, and the main characteristics of the wavefield are well preserved.Temporal and spatial variations of amplification pattern are observed in SEM simulation with various basin models.Our results also clearly illustrate the complexity of ground motion response as waves propagate through the 3D basin structures.These limited responses already show several interesting features produced by the complicated basin structures used in this study.The important phenomena that we observed in this research include the presence of an S-wave multiple created by reflection at the bottom of the basin, small S-to-P mode conversion generated at the bottom of the basin, and obvious wave packets generated due to the near surface low-velocity layer located beneath the surface of the basin.In particularly, we find the basin-generated phases are strongly influenced by lateral variations of surface topography.Moreover, we also show that deep-basin reflection, reverberations in the near-surface lowvelocity layer, attenuation, and topographic scattering all have a strong effect on ground motion amplification.Another significant characteristic of basin response is that the trapping of waves inside the basin produces long wave duration.These results indicate that the SEM can catch the dominant characteristics of wave propagation.We have demonstrated that the SEM is a powerful tool in handling heterogeneous internal structures below the topography of a region and for predicting ground motion response, which leads to its appropriate applications in seismic risk assessment and microzonation studies.

Fig. 1 .
Fig. 1.Three flat-layered models for SEM simulations.The star indicates the dislocation point source and the inverted triangles are the locations of stations.

Fig. 2 .
Fig. 2. The comparison of three-component synthetic velocity seismograms generated by the SEM (indicated by solid, dashed and dotted lines, each line represents different grid spacing) and the FDM (indicated by solid line) for the TEST1 model.The number on the right hand side of each time series represents the value of peak velocity.

Fig. 4 .
Fig. 4. The synthetic seismograms for TEST2 model, we reduce the thickness of soft layer in half to examine the stability of the SEM.

Fig. 6 .
Fig. 6.The synthetic seismograms for model TEST3, we change the thickness of soft layer in half to examine stability with the SEM.

Fig. 7 .
Fig. 7.A half-ellipsoid basin structure embedded in a homogeneous half-space.Triangles indicate the locations of stations and the star indicates the location of epicenter at depth of 10 km.

Fig. 8 .
Fig. 8.The three-component synthetic velocity seismograms calculated by the SEM for the basin model with strong velocity contrast.The top panel displays the cross-section view of the basin.(a) The synthetics shown in the middle panel is for the profile X-X'; (b) while the synthetics for the profile Y-Y' are shown at the bottom panel.The inverse triangles indicate the locations of the site.

Fig. 9 .
Fig. 9. Same as Fig. 8, except using the basin model with a 3D hill.

Fig. 10 .
Fig. 10.Same as Fig. 8, except for using the basin model with a 3D concave valley.

Fig. 11 .
Fig. 11.Same as Fig. 8, except for using the basin model consisting of a thin near surface soft layer.