Peter Moczo

3D Heterogeneous Staggered-grid Finite-difference Modeling of Seismic Motion with Volume Harmonic and Arithmetic Averaging of Elastic Moduli and Densities

We analyze the problem of a heterogeneous formulation of the equation of motion and propose a new 3D fourth-order staggered-grid finite-difference (FD) scheme for modeling seismic motion and seismic-wave propagation. We first consider a 1D problem for a welded planar interface of two half-spaces. A simple physical model of the contact of two media and mathematical considerations are shown to give an averaged medium representing the contact of two media. An exact heterogeneous formulation of the equation of motion is a basis for constructing the corresponding heterogeneous FD scheme. In a much more complicated 3D problem we analyze a planar-interface contact of two isotropic media (both with interface parallel to a coordinate plane and interface in general position in the Cartesian coordinate system) and a nonplanar-interface contact of two isotropic media. Because in the latter case 21 elastic coefficients at each point are necessary to describe the averaged medium, we consider simplified boundary conditions for which the averaged medium can be described by only two elastic coefficients. Based on the simplified approach we construct the explicit heterogeneous 3D fourth-order displacement-stress FD scheme on a staggered grid with the volume harmonic averaging of the shear modulus in grid positions of the stress-tensor components, volume harmonic averaging of the bulk modulus in grid positions of the normal stress-tensor components, and volume arithmetic averaging of density in grid positions of the displacement components. Our displacement-stress FD scheme can be easily modified into the velocity-stress or displacement-velocity-stress FD schemes. The scheme allows for an arbitrary position of the material discontinuity in the spatial grid. Numerical tests for 12 configurations in four types of models show that our scheme is more accurate than the staggered-grid schemes used so far. Numerical examples also show that differences in thickness of a soft surface or interior layer smaller than one grid spacing can cause considerable changes in seismic motion. The results thus underline the importance of having a FD scheme with sufficient sensitivity to heterogeneity of the medium. Manuscript received 21 May 2001.

The Finite-Difference Time-Domain Method for Modeling of Seismic Wave Propagation

We present a review of the recent development in finite-difference time-domain modeling of seismic wave propagation and earthquake motion. The finite-difference method is a robust numerical method applicable to structurally complex media. Due to its relative accuracy and computational efficiency it is the dominant method in modeling earthquake motion and it also is becoming increasingly more important in the seismic industry and for structural modeling. We first introduce basic formulations and properties of the finite-difference schemes including promising recent advances. Then we address important topics as material discontinuities, realistic attenuation, anisotropy, the planar free surface boundary condition, free-surface topography, wavefield excitation (including earthquake source dynamics), non-reflecting boundaries, and memory optimization and parallelization.

3D Fourth-Order Staggered-Grid Finite-Difference Schemes: Stability and Grid Dispersion

We investigated stability and grid dispersion in the 3D fourth-order in space, second-order in time, displacement-stress staggered-grid finite-difference scheme. Though only displacement-stress scheme is explicitly treated, all results also apply to the velocity-stress and displacement-velocity-stress finite-difference schemes. We derived independent stability conditions for the P and S waves by exact separation of equations for the two types of waves. Since the S -wave group velocity can differ from the actual velocity as much as 5% for the sampling ratio 1/5 (that is usually used in modeling), we recommend to sample a minimum S wavelength by six grid spacings. Grid dispersion is strongest for a wave propagating in the direction of a coordinate axis and weakest for a wave propagating along a body diagonal. Grid dispersion in the fourth-order scheme for the sampling ratios s = 1/5 and s = 1/6 is smaller than grid dispersion in the second-order scheme for s = 1/10 and s = 1/12, respectively.

Seismic-Wave Propagation in Viscoelastic Media with Material Discontinuities: A 3D Fourth-Order Staggered-Grid Finite-Difference Modeling

We address the basic theoretical and algorithmic aspects of memory-efficient implementation of realistic attenuation in the staggered-grid finite-difference modeling of seismic-wave propagation in media with material discontinuities. We show that if averaging is applied to viscoelastic moduli in the frequency domain, it is possible to determine anelastic coefficients of the averaged medium representing a material discontinuity. We define (1) the anelastic functions in a new way, as being independent of anelastic coefficients, and (2) a new coarse spatial distribution of the anelastic functions in order to properly account for material discontinuities and, at the same time, to have it memory efficient. Numerical tests demonstrate that our approach enables more accurate viscoelastic modeling than other approaches.

Misfit Criteria for Quantitative Comparison of Seismograms

We have developed and numerically tested quantitative misfit criteria for comparison of seismograms. The misfit criteria are based on the time-frequency representation of the seismograms obtained as the continuous wavelet transform with the analyzing Morlet wavelet. The misfit criteria include time-frequency envelope and phase misfits, time-dependent envelope and phase misfits, frequency-dependent envelope and phase misfits, and single-valued envelope and phase misfits. We tested properties of the misfit criteria using canonical signals. The canonical signals, taken as the reference signals, were specifically amplitude, phase shift, time shift, and frequency modified to demonstrate the ability of the misfit criteria to prop- erly quantify the misfits and recognize the character and cause of the misfits between the reference and modified signals. In all cases the misfit criteria properly quantified and characterized the misfits. The misfit criteria were also calculated for four different numerical solutions for a single layer over half-space (the SCEC LOH.3 problem) and the reference FK so- lution. The misfit criteria provided useful insight into the misfits between individual numerical solutions and the reference solution. The standard RMS misfit matches the single-valued envelope misfit only in the case of a pure amplitude modification of the signal. In all other cases RMS consid- erably overestimates the misfits and does not characterize them.

Quantitative Comparison of Four Numerical Predictions of 3D Ground Motion in the Grenoble Valley, France

This article documents a comparative exercise for numerical simulation of ground motion, addressing the seismic response of the Grenoble site, a typical Alpine valley with complex 3D geometry and large velocity contrasts. Predictions up to 2 Hz were asked for four different structure wave-field configurations (point source and extended source, with and without surface topography). This effort is part of a larger exercise organized for the third international symposium on the effects of surface geology (ESG 2006), the complete results of which are reported elsewhere (Tsuno et al., 2009). While initial, blind computations significantly differed from one another, a remarkable fit was obtained after correcting for some nonmethodological errors for four 3D methods: the arbitrary high-order derivative discontinuous Galerkin method (ADER-DGM), the velocity-stress finite-difference scheme on an arbitrary discontinuous staggered grid (FDM), and two implementations of the spectral-element method (SEM1 and SEM2). Their basic formulation is briefly recalled, and their implementation for the Grenoble Valley and the corresponding requirements in terms of computer resources are detailed. Besides a visual inspection of PGV maps, more refined, quantitative comparisons based on time-frequency analysis greatly help in understanding the origin of differences, with a special emphasis on phase misfit. The match is found excellent below 1 Hz, and gradually deteriorates for increasing frequency, reflecting differences in meshing strategy, numerical dispersion, and implementation of damping properties. While the numerical prediction of ground motion cannot yet be considered a mature, push-button approach, the good agreement reached by four participants indicates that, when used properly, numerical simulation is actually able to handle correctly wave radiation from extended sources in complex 3D media. The main recommendation to obtain reliable numerical predictions of earthquake ground motion is to use at least two different but comparably accurate methods, for instance the present formulations and implementations of the FDM, SEM, and ADER-DGM.

Hybrid seismic modeling based on discrete-wave number and finite-difference methods

Any calculation of seismic wave propagation comprising the seismic source, the travel path, and the receiver site in a single finite-difference (FD) model requires a considerable amount of computer time and memory. Moreover, the methods currently available for including point sources in the 2D FD calculations are far-field approximations only. Therefore we have developed a new hybrid method for treating the seismic wave fields at localized 2D near-surface structures embedded in a 1D background medium, and excited by a point source. The source radiation and propagation in the background model is solved by the discrete-wave number (DW) method, while the propagation in the local 2D structure is calculated by the FD method. The coupling between the two sets of calculations is performed on a rectangular excitation box surrounding the local structure. We show the usefulness of the method in ground-motion studies where both near-field source effects and local site effects are important. Technical problems connected with the inconsistency between the 3D source radiation and the 2D FD calculation are minor for the relatively distant in-plane point explosive sources, but are more serious for the in-plane dislocation sources.

Earthquake Ground Motion in the Mygdonian Basin, Greece: The E2VP Verification and Validation of 3D Numerical Simulation up to 4 Hz

n a low‐seismicity context, the use of numerical simulations becomes essential due to the lack of representative earthquakes for empirical approaches. The goals of the EUROSEISTEST Verification and Validation Project (E2VP) are to provide (1) a quantitative analysis of accuracy of the current, most advanced numerical methods applied to realistic 3D models of sedimentary basins (verification) and (2) a quantitative comparison of the recorded ground motions with their numerical predictions (validation). The target is the EUROSEISTEST site located within the Mygdonian basin, Greece. The site is instrumented with surface and borehole accelerometers, and a 3D model of the medium is available. The simulations are performed up to 4 Hz, beyond the 0.7 Hz fundamental frequency, thus covering a frequency range at which ground motion undergoes significant amplification. The discrete representation of material heterogeneities, the attenuation model, the approximation of the free surface, and nonreflecting boundaries are identified as the main sources of differences among the numerical predictions. The predictions well reproduce some, but not all, features of the actual site effect. The differences between real and predicted ground motions have multiple origins: the accuracy of source parameters (location, hypocentral depth, and focal mechanism), the uncertainties in the description of the geological medium (damping, internal sediment layering structure, and shape of the sediment‐basement interface). Overall, the agreement reached among synthetics up to 4 Hz despite the complexity of the basin model, with code‐to‐code differences much smaller than predictions‐to‐observations differences, makes it possible to include the numerical simulations in site‐specific analysis in the 3D linear case and low‐to‐intermediate frequency range.

Simulation of the Planar Free Surface with Near-Surface Lateral Discontinuities in the Finite-Difference Modeling of Seismic Motion

Kristek et al. (2002) developed a technique for simulating the planar free surface in the 3D fourth-order staggered-grid finite-difference (FD) modeling of seismic motion. The technique is based on (1) explicit application of zero values of the stress-tensor components at the free surface and (2) adjusted FD approximations (AFDAs) to vertical derivatives at and near the free surface. The technique was shown to be more accurate and efficient than the standard stress-imaging technique in 1D models. In this study, we tested accuracy of the AFDA technique in media with lateral material discontinuities reaching the free surface. We compared the FD synthetics with synthetics calculated by the standard finite-element (FE) method because the FE method naturally and sufficiently accurately satisfies the boundary conditions at the free surface and the traction interface continuity conditions at internal material discontinuities. The comparison showed a very good level of accuracy of the AFDA technique. We also demonstrated the very good sensitivity of our FD modeling to different positions of the same physical model in the spatial FD grid.

Stability and Grid Dispersion of the P-SV 4th-Order Staggered-Grid Finite-Difference Schemes

Stability and grid dispersion in the P-SV 4th-order in space, 2nd-order in time, displacement-stress staggered-grid finite-difference scheme is investigated in the case of a homogeneous unbounded medium. All results, however, also apply to the velocity-stress and displacement- velocity-stress finite-difference schemes.Independent stability conditions for the P and S waves are obtained by exact separation of equations for the two types of waves.Since the S-wave group velocity can differ from the actual velocity as much as 5% for the sampling ratio 1/5, commonly used in numerical modelling, the sampling of the minimum S wavelength by 6 grid spacings (with the velocity difference not larger than 2.5%) is recommended.Grid dispersion is strongest for a wave propagating in a direction of a coordinate axis and weakest for a wave propagating along a plane diagonal.Grid dispersion in the 4th-order scheme for the sampling ratios s = 1/5 and s = 1/6 is smaller than grid dispersion in the 2nd-order scheme for s = 1/10 and s = 1/12, respectively.