Simulation Of Seismic Wave Propagation Research Articles

Local time–space mesh refinement for simulation of elastic wave propagation in multi-scale media

This paper presents an original approach to local time–space grid refinement for the numerical simulation of wave propagation in models with localized clusters of micro-heterogeneities. The main features of the algorithm are–the application of temporal and spatial refinement on two different surfaces;–the use of the embedded-stencil technique for the refinement of grid step with respect to time;–the use of the Fast Fourier Transform (FFT)-based interpolation to couple variables for spatial mesh refinement. The latter makes it possible to perform filtration of high spatial frequencies, which provides stability in the proposed finite-difference schemes.In the present work, the technique is implemented for the finite-difference simulation of seismic wave propagation and the interaction of such waves with fluid-filled fractures and cavities of carbonate reservoirs. However, this approach is easy to adapt and/or combine with other numerical techniques, such as finite elements, discontinuous Galerkin method, or finite volumes used for approximation of various types of linear and nonlinear hyperbolic equations.

Numerical earthquake model of the 31 October 2013 Ruisui, Taiwan, earthquake: Source rupture process and seismic wave propagation

We build a numerical earthquake model, including numerical source and wave propagation models, to understand the rupture process and the ground motion time history of the 2013 ML 6.4 Ruisui earthquake in Taiwan. This moderately large event was located in the Longitudinal Valley, a suture zone of the Philippine Sea Plate and the Eurasia Plate. A joint source inversion analysis by using teleseismic body wave, GPS coseismic displacement and near field ground motion data was performed first. The inversion results derived from a western dipping fault plane indicate that the slip occurred in depths between 10 and 20km. The rupture propagated from south to north and two asperities were resolved. The largest one was located approximately 15km north of the epicenter with a maximum slip about 1m. A 3D seismic wave propagation simulation based on the spectral-element method was then carried out by using the inverted source model. A strong rupture directivity effect in the northern area of the Longitudinal Valley was found, which was due to the northward rupture process. Forward synthetic waveforms could explain most of the near-field ground motion data for frequencies between 0.05 and 0.2Hz. This numerical earthquake model not only helps us confirm the detailed rupture processes on the Central Range Fault but also gives contribution to regional seismic hazard mitigation for future large earthquakes.

Numerical Simulation of Seismic Wave Propagation in Different Media

Theory of elastic waves layered homogeneous medium or even medium for the study can not meet the actual demand for seismic exploration, especially for fine rock to construct reservoirs for the study, had to consider small-scale heterogeneity of seismic wave propagation effects. In this thesis, multi-scale model of a complex medium, in ensuring the premise to further improve simulation accuracy simulation efficiency issue, the introduction of a variable grid numerical simulation techniques, and were analyzed for different types of grid difference, establish a different media model simulation results verify the validity of simulation and analyzes its efficiency and accuracy problems.

Numerical study of the interface errors of finite-difference simulations of seismic waves

Numerical simulations of wave propagation produce different errors and the most well known is numerical dispersion, which is only valid for homogeneous media. However, there is a lack of error studies for heterogeneous media or even for the canonical case of media that have two constant velocity layers. The error associated with media that have two layers is called an interface error, and it typically converges to zero with a lower order of convergence compared to the theoretical convergence rate of the finite-difference schemes (FDS) for homogeneous media. We evaluated a detailed numerical study of the interface error for three staggered-grid FDS that are commonly used in the simulation of seismic-wave propagation. We determined that a standard staggered-grid scheme (SSGS) (also known as the Virieux scheme), a rotated staggered-grid scheme (RSGS), and a Lebedev scheme (LS) preserve the second order of convergence at horizontal/vertical solid-solid interfaces when the medium parameters have been properly modified, such as by harmonic averaging of finely layered media for the stiffness tensor and arithmetic mean for the density. However, for a fluid-solid interface aligned with the grid line, a second-order convergence can only be achieved by an SSGS. In addition, the presence of a fluid-solid interface reduces the order of convergence for the LS and the RSGS to a first order of convergence. The presence of inclined interfaces makes high-order (second and more) convergence impossible.

Scattering of teleseismic P-waves by the Japan Trench: A significant effect of reverberation in the seawater column

We detected a scattered wave train in data from the high-sensitivity seismograph network in Japan (Hi-net) following the arrival of the near-vertically incident P-wave generated by the 2009 earthquake (Mw 7.8) off the South Island of New Zealand. The scattered wave train represented predominantly vertical ground motion at a period of 20 to 50 s and with an apparent velocity of 3.5 km/s; it propagated cylindrically westward through the Kanto area of central Japan. Array analysis showed that the scattered wave train developed beneath the Pacific Ocean near the Boso triple junction, southeast of the Kanto area. A 3D finite-difference simulation of seismic wave propagation using a high-resolution model incorporating subsurface structure, topography, and bathymetry revealed that the strong scattered waves that were generated along the Japan Trench and propagated normal to the trench axis represented multiple reverberations of seismic waves between the seafloor and the Pacific plate boundary. In addition, strong reverberation of acoustic waves in the seawater column above the Boso triple junction causes elongated scattered waves, which reasonably explains our observations.

Numerical modeling of 3D zero-offset laboratory data by a discretized Kirchhoff integral method

Accurate simulation of seismic wave propagation in complex geologic structures is of particular interest nowadays. However, difficulties arise for complex geologic structures with great and rapid structural changes, due, for instance, to the presence of shadow zones, head waves, diffractions and/or edge effects. Different methods have thus been developed and are typically tested on synthetic configurations against analytical solutions for simple canonical problems, reference methods, or via direct comparison with real data acquired in situ. Such approaches have limitations, especially if the propagation occurs in a complex environment with strong-contrast reflectors and surface irregularities because it can be difficult to determine the method that gives the best approximation of the “real” solution or to interpret the results obtained without an a priori knowledge of the geologic environment. An alternative approach for seismics consists in comparing the synthetic data with data obtained in laboratory experiments. In contrast to in situ experiments, high-quality data are collected under controlled conditions for a known configuration. In contrast with numerical experiments, laboratory data possess many of the characteristics of field data because real waves propagate through models with no numerical approximations. Our main purpose was to test the approach of using laboratory data as reference data for benchmarking 3D numerical methods and techniques using the setup that we have designed for this study. We performed laboratory-scaled measurements of zero-offset reflection of broadband pulses from a strong topographic environment immersed in a water tank. We compared these measurements with numerical data simulated by means of a discretized Kirchhoff integral method. The comparisons of synthetic and laboratory data indicated a good quantitative fit in terms of time arrivals and acceptable fit in amplitudes. Thus, the first step of the approach was successfully applied.

Strong seismic wave scattering in the low-velocity anomaly associated with subduction of oceanic plate

Analyses of dense seismic records in Kanto, Japan, revealed distinct pulse broadening and peak delay of high-frequency S waves at central Chiba. These phenomena are observed at frequency range of 1–8 Hz and exist only for ray paths passing through the low-velocity (LV) zone at depth of 20–40 km beneath northwestern Chiba. To obtain a more detailed understanding of these phenomena, we conducted 2-D and 3-D finite difference method simulations of seismic wave propagation using a realistic heterogeneous structure model. Through numerous simulations we demonstrated that strong seismic scattering, due to localized strong small-scale heterogeneities in the LV zone and in the oceanic crust, is a major cause of strong pulse broadening and peak delay of high-frequency S waves. After comparing simulation results with observations, the most preferable small-scale velocity heterogeneity in the LV zone is characterized by a Gaussian power spectral density function (PSDF) with correlation distance a of 1–2 km and rms value ε = 0.07–0.09, superposed on a background exponential PSDF (a = 3 km, ε = 0.07). Assuming strong velocity heterogeneities, observed amplitude decay at Chiba is also well explained by strong scattering attenuation in the LV zone. Because the LV zone, which has been reported by seismic tomography studies, is interpreted as being constructed by the dehydration of the subducting oceanic crust of the Philippine Sea Plate, strong small-scale velocity heterogeneity in the LV zone may be related to the random distribution of fluid in this volume.

Strong Ground Motion Prediction Using Virtual Earthquakes

Sedimentary basins increase the damaging effects of earthquakes by trapping and amplifying seismic waves. Simulations of seismic wave propagation in sedimentary basins capture this effect; however, there exists no method to validate these results for earthquakes that have not yet occurred. We present a new approach for ground motion prediction that uses the ambient seismic field. We apply our method to a suite of magnitude 7 scenario earthquakes on the southern San Andreas fault and compare our ground motion predictions with simulations. Both methods find strong amplification and coupling of source and structure effects, but they predict substantially different shaking patterns across the Los Angeles Basin. The virtual earthquake approach provides a new approach for predicting long-period strong ground motion.

Numerical modeling of zero-offset laboratory data in a strong topographic environment: results for a spectral-element method and a discretized Kirchhoff integral method

Accurate simulation of seismic wave propagation in complex geological structures is of particular interest nowadays. However conventional methods may fail to simulate realistic wavefields in environments with great and rapid structural changes, due for instance to the presence of shadow zones, diffractions and/or edge effects. Different methods, developed to improve seismic modeling, are typically tested on synthetic configurations against analytical solutions for simple canonical problems or reference methods, or via direct comparison with real data acquired in situ. Such approaches have limitations, especially if the propagation occurs in a complex environment with strong-contrast reflectors and surface irregularities, as it can be difficult to determine the method which gives the best approximation of the “real” solution, or to interpret the results obtained without an a priori knowledge of the geologic environment. An alternative approach for seismics consists in comparing the synthetic data with high-quality data collected in laboratory experiments under controlled conditions for a known configuration. In contrast with numerical experiments, laboratory data possess many of the characteristics of field data, as real waves propagate through models with no numerical approximations. We thus present a comparison of laboratory-scaled measurements of 3D zero-offset wave reflection of broadband pulses from a strong topographic environment immersed in a water tank with numerical data simulated by means of a spectral-element method and a discretized Kirchhoff integral method. The results indicate a good quantitative fit in terms of time arrivals and acceptable fit in amplitudes for all datasets.

A Hamiltonian Particle Method with a Staggered Particle Technique for Simulating Seismic Wave Propagation

We present a Hamiltonian particle method (HPM) with a staggered particle technique for simulating seismic wave propagation. In the conventional HPM, physical variables, such as particle displacement and stress, are defined at the center, i.e., at the same position, of each particle. As most seismic simulations using finite difference methods (FDM) are practiced with staggered grid techniques, we know the staggered alignment of space variables could improve the numerical accuracy. In the present study, we hypothesized that staggered technique could improve the numerical accuracy also in the HPM and tested the hypothesis. First, we conducted a plane wave analysis for the HPM with the staggered particles in order to verify the validity of our strategy. The comparison of grid dispersion in our strategy with that in the conventional one suggests that the accuracy would be improved dramatically by use of the staggered technique. It is also observed that the dispersion of waves is dependent on the propagation direction due to the difference in the average spacing of the neighboring two particles for the same parameters, as is usually observed in FDM with a rotated staggered grid. Next, we compared the results from the conventional Lamb’s problem using our HPM with those from an analytical approach in order to demonstrate the effectiveness of the staggered particle technique. Our results showed better agreement with the analytical solutions than those from HPM without the staggered particles. We conclude that the staggered particle technique would be a method to improve the calculation accuracy in the simulation of seismic wave propagation.

Lg wave propagation in the area around Japan: observations and simulations

Regional wavefields are strongly influenced by crustal structure heterogeneity variations along their propagation paths. Observations of such effects between the Asian continent and the Japanese subduction zone across the Sea of Japan (East Sea) have been strongly assisted recently by the development of high-density seismic networks in Japan and Korea, as well as the supercomputer-based three-dimensional finite difference method seismic wave propagation simulations using detailed heterogeneous crustal models. In this study, an Lg propagation map derived from 289,000 ray paths connecting sources to observation stations reveals efficient Lg wave propagation from continental Asia to Kyushu through the Korean Peninsula, and to Hokkaido, thus indicating a laterally consistent crustal structure extending from continental Asia to Japan. However, the Lg wave propagation to the Japanese main island of Honshu is totally blocked as it crosses the continental-oceanic boundary surrounding the Sea of Japan. Three-dimensional (3-D) finite-difference method seismic wave propagation simulations performed using a detailed crustal structural model allow us to clearly visualize the way in which an Lg wave develops from a shallow source in the crust and its propagation in the crustal waveguide by means of multiple post-critical S wave reflections in the continental structure. The sudden thinning of the continental crust at the edge of the Asian continent adjacent to the oceanic crust in the Sea of Japan, which involves a thickness change from 30 to 10 km within a 100-km distance (with the thinner crust extending over 600 km), decreases the Lg wave energy by 10%. It has been confirmed that 50% of this Lg wave energy loss occurs during wave passage along the thinner crust and that the other 50% results from conversion into P energy in the overlying seawater.

Accelerating the discontinuous Galerkin method for seismic wave propagation simulations using multiple GPUs with CUDA and MPI

We have successfully ported an arbitrary high-order discontinuous Galerkin method for solving the three-dimensional isotropic elastic wave equation on unstructured tetrahedral meshes to multiple Graphic Processing Units (GPUs) using the Compute Unified Device Architecture (CUDA) of NVIDIA and Message Passing Interface (MPI) and obtained a speedup factor of about 28.3 for the single-precision version of our codes and a speedup factor of about 14.9 for the double-precision version. The GPU used in the comparisons is NVIDIA Tesla C2070 Fermi, and the CPU used is Intel Xeon W5660. To effectively overlap inter-process communication with computation, we separate the elements on each subdomain into inner and outer elements and complete the computation on outer elements and fill the MPI buffer first. While the MPI messages travel across the network, the GPU performs computation on inner elements, and all other calculations that do not use information of outer elements from neighboring subdomains. A significant portion of the speedup also comes from a customized matrix–matrix multiplication kernel, which is used extensively throughout our program. Preliminary performance analysis on our parallel GPU codes shows favorable strong and weak scalabilities.

An Advanced Geophysical Software Tool for Seismic Modeling Based on the Finite Element Method

We demonstrate an application of an advanced geophysical software tool to simulate seismic wave propagation and generate shot records. Simulation of seismic wave propagation was performed using eScript, a partial differential equation solver based on the finite element method. We solved the 2D acoustic wave equation for P-wave propagation by applying the partial differential equation to eScript. A traditional absorbing boundary condition was applied to reducing artificial wave reflections at the boundary. We tested several different models by employing meshes for layered media and other complex geological structures. Seismometers were placed equidistantly on the surface to record the wave field and generate shot records. Model tests prove that eScript is a user-friendly software platform which allows rapid development of numerical seismic modeling using python and XML.

Finite difference simulations of seismic wave propagation for understanding earthquake physics and predicting ground motions: Advances and challenges

Seismic waves radiated from an earthquake propagate in the Earth and the ground shaking is felt and recorded at (or near) the ground surface. Understanding the wave propagation with respect to the Earth's structure and the earthquake mechanisms is one of the main objectives of seismology, and predicting the strong ground shaking for moderate and large earthquakes is essential for quantitative seismic hazard assessment. The finite difference scheme for solving the wave propagation problem in elastic (sometimes anelastic) media has been more widely used since the 1970s than any other numerical methods, because of its simple formulation and implementation, and its easy scalability to large computations. This paper briefly overviews the advances in finite difference simulations, focusing particularly on earthquake mechanics and the resultant wave radiation in the near field. As the finite difference formulation is simple (interpolation is smooth), an easy coupling with other approaches is one of its advantages. A coupling with a boundary integral equation method (BIEM) allows us to simulate complex earthquake source processes.

Numerical analysis of the relationship between time-variant coda-Q and the variation in crustal stress

SUMMARY Coda-Q is a stochastic parameter reflecting the heterogeneities of medium that seismic waves travel through. We confirmed that coda-Q would vary with the stress loaded to an elastic medium using numerical simulations of seismic wave propagation. When the stress is loaded, cracks in the crust could either close or newly open. The closure and opening of the cracks are not random but depending on the magnitude and the direction of the stress and the crack aspect ratio. The cracks in the medium after loading stress could be aligned in a specific orientation, and elastic wave velocity field would become anisotropic due to the alignment of specific crack orientations. Elastic wave velocity is in general faster along the direction corresponding with the crack orientation while slower along the perpendicular direction. In the numerical simulation, the effect of anisotropy in elastic wave velocity field due to the selective closure and opening of the cracks is calculated using a 2-D finite difference method assuming elastic wave velocity to be a function of the magnitude of loaded stress. The coda-Q calculated from seismic waves simulated for a model varies when the averaged normal stress changes. Our simulation indicated that the sensitivity of coda-Q −1 , that is the reciprocal of the coda-Q, would be 1.0 × 10 −2 (1.0MPa ‐1 ) against the magnitude of the confining pressure and 1.0 × 10 −3 (1.0 deg ‐1 ) against the direction of principal stress. We would like to conclude that coda-Q, a stochastic parameter reflecting heterogeneities of subsurface medium, could become a quantitative state indicator of the stress field of the medium where seismic waves propagate through. Spatiotemporal variation of coda-Q reflects change in the stress field in the crust.

Full waveform tomography of the upper mantle in the South Atlantic region: Imaging a westward fluxing shallow asthenosphere?

A prominent feature of the South Atlantic region is its strongly asymmetric residual bathymetry across the ocean basin. It has been suggested that the residual bathymetry is dynamic in nature, arising from the large slow velocity seismic anomaly located in the lower mantle beneath the African plate. Unfortunately, the pattern of mantle heterogeneity particularly in the upper mantle is not well known owing to the sparsity of seismic stations and the existence of large aseismic regions on the African and South American plates. Here we present a new seismic tomographic study of the South Atlantic upper mantle. Our model is based on a full seismic waveform inversion of ≈ 4000 high-quality seismograms for isotropic 3-D seismic structure using a powerful adjoint methodology capable of extracting maximum information from each seismogram. The theory requires simulation of seismic wave propagation in 3-D heterogeneous earth models computed with a spectral-element method where the differences between observed and synthetic seismograms are quantified using phase misfits obtained through a time–frequency transform. The model images a continuous channel of pronounced slow seismic velocity in the shallow sublithospheric mantle between ≈ 150 and ≈ 300 km depth that branches in between the cratonic roots under the African and South American continents. At greater depth, below 300–350 km, the slow anomalies are less pronounced and a change in the convective planform is indicated by isolated, round shaped patches in an overall faster mantle. It is possible that the depthwise change of the convective planform from vertical to horizontal advection of hot buoyant material in a low viscosity asthenosphere can reconcile the anomalous residual bathymetry in the region, and we present a simple fluid dynamic model of pressure driven flow to assess the feasibility of this scenario.

Small‐scale heterogeneities in the oceanic lithosphere inferred from guided waves

We analyze seismic waveforms from deep‐focus earthquakes occurring in the subducting slab beneath Japan, recorded by broadband ocean bottom seismometers (BBOBSs) installed on the northwestern Pacific Ocean seafloor. The data reveal waveforms with a low‐frequency direct P onset, followed by large‐amplitude, high‐frequency, long‐duration Po and So waves. From the analysis of the BBOBS records and a numerical finite‐difference method simulation of seismic wave propagation, we elucidate the generation and propagation processes of such guided waves. We demonstrate that the low‐frequency direct P and S waves propagate in the asthenosphere and that the following high‐frequency, long‐duration Po and So waves are developed by multiple forward scattering of P and S waves. The scattering occurs due to laterally elongated heterogeneities in both the subducting and horizontal parts of the oceanic lithosphere, with the apparent velocities (Vp = 8.1 km/s, Vs = 4.6 km/s) being close to the velocities of oceanic lithosphere.

Quantitative comparison between simulations of seismic wave propagation in heterogeneous poro-elastic media and equivalent visco-elastic solids for marine-type environments

There is increasing evidence to suggest that the presence of mesoscopic heterogeneities constitutes an important seismic attenuation mechanism in porous rocks. As a consequence, centimetre-scale perturbations of the rock physical properties should be taken into account for seismic modelling whenever detailed and accurate responses of specific target structures are desired, which is, however, computationally prohibitive. A convenient way to circumvent this problem is to use an upscaling procedure to replace each of the heterogeneous porous media composing the geological model by corresponding equivalent visco-elastic solids and to solve the visco-elastic equations of motion for the inferred equivalent model. While the overall qualitative validity of this procedure is well established, there are as of yet no quantitative analyses regarding the equivalence of the seismograms resulting from the original poro-elastic and the corresponding upscaled visco-elastic models. To address this issue, we compare poro-elastic and visco-elastic solutions for a range of marine-type models of increasing complexity. We found that despite the identical dispersion and attenuation behaviour of the heterogeneous poro-elastic and the equivalent visco-elastic media, the seismograms may differ substantially due to diverging boundary conditions, where there exist additional options for the poro-elastic case. In particular, we observe that at the fluid/porous-solid interface, the poro- and visco-elastic seismograms agree for closed-pore boundary conditions, but differ significantly for open-pore boundary conditions. This is an important result which has potentially far-reaching implications for wave-equation-based algorithms in exploration geophysics involving fluid/porous-solid interfaces, such as, for example, wavefield decomposition.

Scattering of high-frequency P wavefield derived by dense Hi-net array observations in Japan and computer simulations of seismic wave propagations

SUMMARY We studied the scattering properties of high-frequency seismic waves due to the distribution of small-scale velocity fluctuations in the crust and upper mantle beneath Japan based on an analysis of three-component short-period seismograms and comparison with finite difference method (FDM) simulation of seismic wave propagation using various stochastic random velocity fluctuation models. Using a large number of dense High-Sensitivity Seismograph network waveform data of 310 shallow crustal earthquakes, we examined the P-wave energy partition of transverse component (PEPT), which is caused by scattering of the seismic wave in heterogeneous structure, as a function of frequency and hypocentral distances. At distance of less than D = 150km, the PEPT increases with increasing frequency and is approximately constant in the range of from D = 50 to 150 km. The PEPT was found to increase suddenly at a distance of over D = 150 km and was larger in the high-frequency band (f > 4 Hz). Therefore, strong scattering of P wave may occur around the propagation path (upper crust, lower crust andaroundMohodiscontinuity)oftheP-wavefirstarrivalphaseatdistancesoflargerthanD = 150km.WealsofoundaregionaldifferenceinthePEPT value,wherebythePEPT valueislarge at the backarc side of northeastern Japan compared with southwestern Japan and the forearc side of northeastern Japan. These PEPT results, which were derived from shallow earthquakes, indicate that the shallow structure of heterogeneity at the backarc side of northeastern Japan is stronger and more complex compared with other areas. These hypotheses, that is, the depth andregionalchangeofsmall-scalevelocityfluctuations,areexaminedby3-DFDMsimulation using various heterogeneous structure models. By comparing the observed feature of the PEPT with simulation results, we found that strong seismic wave scattering occurs in the lower crust due to relatively higher velocity and stronger heterogeneities compared with that in the upper crust. To explain the observed regional difference, the velocity fluctuation model with 3‐4 per cent stronger fluctuation and smaller κ is required at the backarc side of northeastern Japan.

Probabilistic full waveform inversion based on tectonic regionalization—development and application to the Australian upper mantle

We present a first study to investigate the feasibility of a probabilistic 3-D full waveform inversion based on spectral-element simulations of seismic wave propagation and Monte Carlo exploration of the model space. Through a tectonic regionalization we reduce the dimension of the model space to 12, and we incorporate complete seismograms in order to exploit all available information in the period range from 60 to 200 s. S-wave velocity variations in the Australian Archean and Proterozoic lithospheres are generally well-constrained and strongly positive, in agreement with previous inferences from deterministic tomography. The maximum likelihood model reveals significantly elevated P velocities. While consistent with body wave studies, they are, however, not well constrained by our data, thereby providing an interesting example of a comparatively insignificant maximumlikelihood model. Our data are notably affected by 3-D density variations. The effect, however, appears to be misleading. Both the maximum-likelihood model and the posterior probability densities strongly prefer unrealistically large positive density variations that are inconsistent with independent information from geodynamics and mineral physics. This suggests that highly probable and less extreme density models may be hidden in small and hardly detectable subvolumes of the 12-D model space. It follows that deterministic full waveform inversions for density may require particularly accurate initial models. From a methodological perspective we must conclude that a transition to significantly higher dimensions would currently be difficult. Available computing power clearly imposes restrictions. However, even when Moore’s law continues to hold, the largest obstacle appears to be our inability to efficiently map small high-dimensional subvolumes with high probability.