Simulation Of Seismic Waves Research Articles

Synthesizing ocean bottom pressure records including seismic wave and tsunami contributions: Toward realistic tests of monitoring systems

AbstractThe present study proposes a method for synthesizing the ocean bottom pressure records during a tsunamigenic earthquake. First, a linear seismic wave simulation is conducted with a kinematic earthquake fault model as a source. Then, a nonlinear tsunami simulation is conducted using the sea bottom movement calculated in the seismic wave simulation. By using these simulation results, this method can provide realistic ocean bottom pressure change data, including both seismic and tsunami contributions. A simple theoretical consideration indicates that the dynamic pressure change caused by the sea bottom acceleration can contribute significantly until the duration of ~90 s for a depth of 4000 m in the ocean. The performance of a tsunami monitoring system was investigated using the synthesized ocean bottom pressure records. It indicates that the system based on the hydrostatic approximation could not measure the actual tsunami height when the time does not elapse enough. The dynamic pressure change and the permanent sea bottom deformation inside the source region break the condition of a simple hydrostatic approximation. A tsunami source estimation method of tFISH is also examined. Even though the synthesized records contain a large dynamic pressure change, which is not considered in the algorithm, tFISH showed a satisfactory performance 5 min after the earthquake occurrence. The pressure records synthesized in this study, including both seismic wave and tsunami contributions, are more practical for evaluating the performance of our monitoring ability, whereas most tsunami monitoring tests neglect the seismic wave contribution.

Modeling of low-frequency seismic waves in a shallow sea using the staggered grid difference method

Elastic waves in the seabed generated by low-frequency noise radiating from ships are known as ship seismic waves and can be used to detect and identify ships. To obtain the propagation characteristics of ship seismic waves, an algorithm for calculating seismic waves at the seafloor is presented based on the staggered-grid finite difference method. The accuracy of the algorithm was tested by comparison with analytical solutions. Numerical simulation of seismic waves generated by a low-frequency point sound source in a typical shallow sea environment was carried out. Using various source frequencies and locations in the numerical simulation, we show that the seismic waves in the near field are composed mostly of transmitted S-waves and interface waves while transmitted P-waves are weak near the seafloor. However, in the far field, the wave components of the seismic wave are mainly normal modes and interface waves, with the latter being relatively strong in the waveforms. As the source frequency decreases, the normal modes become smaller and the interface waves dominate the time series of the seismic waves.

Discretization effects in the finite element simulation of seismic waves in elastic and elastic-plastic media

Presented here is a numerical investigation that (re-)appraises standard rules for space/time discretization in seismic wave propagation analyses. Although the issue is almost off the table of research, situations are often encountered where (established) discretization criteria are not observed and inaccurate results possibly obtained. In particular, a detailed analysis of discretization criteria is carried out for wave propagation through elastic and elastic-plastic media. The establishment of such criteria is especially important when accurate prediction of high-frequency motion is needed and/or in the presence of highly non-linear material models. Current discretization rules for wave problems in solids are critically assessed as a conditio sine qua non for improving verification/validation procedures in applied seismology and earthquake engineering. For this purpose, the propagation of shear waves through a 1D stack of 3D finite elements is considered, including the use of wide-band input motions in combination with both linear elastic and non-linear elastic-plastic material models. The blind use of usual rules of thumb is shown to be sometimes debatable, and an effort is made to provide improved discretization criteria. Possible pitfalls of wave simulations are pointed out by highlighting the dependence of discretization effects on time duration, spatial location, material model and specific output variable considered.

Multi-hole seismic modeling in 3-D space and cross-hole seismic tomography analysis for boulder detection

A boulder stone, a common geological feature in south China, is referred to the remnant of a granite body which has been unevenly weathered. Undetected boulders could adversely impact the schedule and safety of subway construction when using tunnel boring machine (TBM) method. Therefore, boulder detection has always been a key issue demanded to be solved before the construction. Nowadays, cross-hole seismic tomography is a high resolution technique capable of boulder detection, however, the method can only solve for velocity in a 2-D slice between two wells, and the size and central position of the boulder are generally difficult to be accurately obtained. In this paper, the authors conduct a multi-hole wave field simulation and characteristic analysis of a boulder model based on the 3-D elastic wave staggered-grid finite difference theory, and also a 2-D imaging analysis based on first arrival travel time. The results indicate that (1) full wave field records could be obtained from multi-hole seismic wave simulations. Simulation results describe that the seismic wave propagation pattern in cross-hole high-velocity spherical geological bodies is more detailed and can serve as a basis for the wave field analysis. (2) When a cross-hole seismic section cuts through the boulder, the proposed method provides satisfactory cross-hole tomography results; however, when the section is closely positioned to the boulder, such high-velocity object in the 3-D space would impact on the surrounding wave field. The received diffracted wave interferes with the primary wave and in consequence the picked first arrival travel time is not derived from the profile, which results in a false appearance of high-velocity geology features. Finally, the results of 2-D analysis in 3-D modeling space are comparatively analyzed with the physical model test vis-a-vis the effect of high velocity body on the seismic tomographic measurements.

Multi-point shaking table test for long tunnels subjected to non-uniform seismic loadings - part II: Application to the HZM immersed tunnel

A series of multi-point shaking table tests are conducted on a long immersed tunnel designed for the Hongkong-Zhuhai-Macau linkage (HMZ linkage) under non-uniform seismic excitations. Details of experimental setup are first presented with particular focuses on: shaking table array of the four independent shaking tables; dynamic similitude design between scale model of soil and tunnel structure and prototype model; design and fabrication of the scaled model tunnel and model soil; and procedure for simulation of non-uniform seismic waves. A series of testing cases of the shaking tables are carried out on the model tunnel using input motions with different seismic amplitudes. Dynamic responses measured from the experiment include acceleration of soil stratum and model tunnel, as well as extension and closure of immersion joints. Acceleration and deformation responses in the tunnel segments and their connecting joints are investigated and analyzed in the test cases under both uniform and non-uniform seismic loadings, and the capacity of critical structural components, such as extension of joints is evaluated. Results show that the non-uniform seismic excitation significantly aggravates the extension of immersion joints compared to the uniform excitation, and thus the non-uniform seismic effect should be considered in the design of immersed tunnels.

A dynamic lattice method for elastic seismic modeling in anisotropic media

Basing on a particle network model, we have developed a dynamic lattice approach to simulate seismic wave propagation in transversely isotropic (TI) media with tilted symmetry axis (TTI). Different from other wave-equation-based numerical methods, the dynamic lattice approach calculates the micromechanical interactions between particles in the lattice instead of solving the wave equation. We implement Lagrange’s equations to transform these interactions into elastic forces acting upon each particle. By solving the equations of motion, we obtain the disturbances of particles. Therefore, seismic wave motions in continua are approximated by displacements of these particles. Elastic features of the continuum are represented by properties of the particle lattice, including the physical properties of particles and the micromechanical interactions between particles. In the case of TI media, it is a challenge to find the correct particle lattice model that can reflect the anisotropic nature of TI media. We have determined the theoretical connection between the TI medium and the particle lattice model, allowing us to model elastic seismic waves in heterogeneous TI media. Furthermore, we have linearized the propagator to improve the CPU efficiency of our method for seismic wave simulation. We have applied the method to reverse time migration on a TI model to test its usefulness in complex media imaging.

Using time filtering to control the long-time instability in seismic wave simulation

Long-time instabilities have been observed in various scenarios of numerical simulation for seismic wave propagation. They appear as slowly magnifying spurious oscillations in the seismograms at the late stage of the simulation. Their magnifying speed is typically much slower than the magnifying speed observed when the Courant–Friedrichs–Lewy condition is violated. The simulations can therefore continue to proceed for a relatively long period without floating-point overflow. The impact of the long-time instabilities on the simulation accuracy at the early stage can be negligible in some cases. In existing literatures, spatial-filtering techniques that, in principle, average the solution within certain spatial range at the same time level are typically utilized to control the long-time instability. In this paper, we present an alternative time-filtering approach that, in principle, averages the solution at different time levels of the same spatial location to control the long-time instability. Comparing with the spatial filtering, the advantages of this time-filtering approach lie in its flexibility, particularly when boundaries or interfaces are involved, its simplicity and low additional arithmetic operations, at the expense of extra memory cost. When application of the time filtering is localized to regions where long-time instabilities are emitted from, for example, a boundary or an interface layer, the additional cost is negligible when compared with the cost of wave simulation. For linear wave equations, this time-filtering approach can be understood as the introduction of artificial diffusion. Its application has impact on the accuracy of the solution and the restriction of the time step size. We present an indicator-based approach to adjust the filtering parameters both spatially and temporally, in order to provide the best trade-off between accuracy and stability. The indicator is calculated heuristically by monitoring the spurious oscillation as the simulation evolves in time.

Asymmetric chirplet transform for sparse representation of seismic data

For a fast and accurate extraction of important information in seismic signals, a sparse representation based on the physical parameters of the given data is crucial. We have used the asymmetric Gaussian chirplet model (AGCM) and established a dictionary-free variant of the orthogonal matching pursuit, a greedy algorithm for sparse approximation. The atoms of AGCM, so-called chirplets, display asymmetric oscillation attenuation properties, which make the AGCM very suitable for sparse representation of absorption decay seismic signals. Unlike the Fourier transform, which assumes that the seismic signals consist of plane waves, the AGCM assumes that the seismic signal consists of nonstationary compressed plane waves, i.e., symmetric or asymmetric chirplets. Thus, AGCM is a general model for seismic wave simulation, and its model parameters include an envelope part and a phase part. We have mainly concentrated on the parameters of the envelope part, such as the envelope amplitude and arrival time. We provide numerical examples using the algorithm for seismic signal approximation and arrival-time detection. The results indicated promising performance but also may be improved considering the spatial correlations of seismic data.

An immersed free-surface boundary treatment for seismic wave simulation

Finite-difference methods are popular for wave simulation within the seismic exploration community, thanks to their efficiency. However, difficulties arise when encountering complex topography due to the regular grid pattern of the finite-difference schemes. Despite alternatives that can handle the free surface with little effort, such as the spectral element or discontinuous Galerkin’s methods, incorporating a free-surface boundary condition within the finite-difference framework is still appealing, even at the cost of extra algorithm complexity and stronger requirement of computational resources. We present a free-surface boundary treatment within the finite-difference framework, belonging to the family of the immersed-boundary methods. Inherently, the presented boundary treatment is separated from the rest of the wave simulation, which makes it easy to be integrated in existing finite-difference codes. Specifically, we construct an extrapolation operator for each grid point above the free surface, if requested by the finite-difference stencil, to estimate its fictitious wavefield value at each time step. These operators are constructed only once and remain unchanged for all the time steps and source locations. The memory requirement of these operators is significant. Fortunately, grouping together multiple simulations concerning different source locations makes it possible to dilute the memory burden to a negligible level. Additionally, applying these operators incurs numerical noise, which may lead to long time instabilities. In such a scenario, additional numerical procedures, for instance, introducing artificial diffusion, are necessary to control the instability and obtain sensible simulation results. Successful applications of the presented boundary treatment to elastic-wave equations on domains with nontrivial topographies, in 2D and 3D, are presented. Robust and efficient numerical techniques to control high-frequency numerical noise remain to be investigated.

A staggered discontinuous Galerkin method for the simulation of seismic waves with surface topography

Accurate simulation of seismic waves is of critical importance in a variety of geophysical applications. Based on recent works on staggered discontinuous Galerkin methods, we have developed a new method for the simulations of seismic waves, which has energy conservation and extremely low grid dispersion, so that it naturally provided accurate numerical simulations of wave propagation useful for geophysical applications and was a generalization of classical staggered-grid finite-difference methods. Moreover, it could handle with ease irregular surface topography and discontinuities in the subsurface models. Our new method discretized the velocity and the stress tensor on this staggered grid, with continuity imposed on different parts of the mesh. The symmetry of the stress tensor was enforced by the Lagrange multiplier technique. The resulting method was an explicit scheme, requiring the solutions of a block diagonal system and a local saddle point system in each time step, and it was, therefore, very efficient. To tailor our scheme to Rayleigh waves, we developed a mortar formulation of our method. Specifically, a fine mesh was used near the free surface and a coarse mesh was used in the rest of the domain. The two meshes were in general not matching, and the continuity of the velocity at the interface was enforced by a Lagrange multiplier. The resulting method was also efficient in time marching. We also developed a stability analysis of the scheme and an explicit bound for the time step size. In addition, we evaluated some numerical results and found that our method was able to preserve the wave energy and accurately computed the Rayleigh waves. Moreover, the mortar formulation gave a significant speed up compared with the use of a uniform fine mesh, and provided an efficient tool for the simulation of Rayleigh waves.

Second-order seismic wave simulation in the presence of a free-surface by pseudo-spectral method

In this paper second-order seismic wave equations are adopted to simulate the wave propagation in the presence of a free-surface. In the simulations, horizontal and vertical spatial derivatives are solved by Fourier and Chebyshev pseudo-spectral algorithms respectively. In addition, Fourier–Chebyshev pseudo-spectral method is also used to solve the boundary conditions of the free-surface, which are Neumann boundaries. That makes results at the free-surface keep same high precision with inner domain. The new scheme is tested against Lamb's problem in a uniform elastic half space. Agreement between numerical and analytical results is every good. Compared with the former simulation scheme based first-order velocity–stress equations, the new scheme based on second-order equations is more efficient because it has fewer controlling equations.

Multi-scale Coupling Simulation of Seismic Waves and Building Vibrations Using ppOpen-HPC

Abstract In order to simulate an earthquake shock originating from the earthquake source and the damage it causes to buildings, not only the seismic wave that propagates over a wide region of several 100 km2, but also the building vibrations that occur over a small region of several 10 m2 must be resolved concurrently. Such a multi-scale simulation is difficult because such kind of modeling and implementation by only a specific application are limited. To overcome these problems, a multi-scale weak-coupling simulation of seismic wave and building vibrations using “ppOpen-HPC” libraries is conducted. The ppOpen-HPC, wherein “pp” stands for “post-peta scale”, is an open source infrastructure for development and execution of optimized and reliable simulation codes on large-scale parallel computers. On the basis of our evaluation, we confirm that an acceptable result can be achieved that ensures that the overhead cost of the coupler is negligible and it can work on large-scale computational resources.

Solution of Large-scale Seismic Modeling Problems

Abstract The grid-characteristic method is commonly used to accurately model elastic waves propagation within the application of mechanical stress tests, collision of solid bodies and others. In this note, we discuss the application of the grid-characteristic method to numerical simulation of seismic waves which is an important problem in geophysics. A computational code for this application has been developed in our team. We studied its low-level optimization with streaming SIMD instructions (SSE and AVX). Beyond, we parallelized our code with MPI and reached parallel efficiency of up to 70% when using 16’000 computational cores.

Seismic attenuation in fractured media

The prime objective of this paper is to quantitatively estimate seismic attenuation caused by fractures with different physical parameters. In seismic wave simulation, the fractured media are treated as the anisotropic media and fractures are represented by frequency-dependent elastic constants. Based on numerical experiments with three different parameters, namely viscosity, porosity and the Lame parameters, this paper has the following observations. First, seismic attenuation is not affected by the viscosity within fractures, although it increases with the increase of porosity and decreases with the increase of the Lame parameters within fractures. Among the latter two parameters, seismic attenuation is more sensitive to the Lame parameters than to the porosity. Second, for the attenuation anisotropy, low frequencies have more anisotropic effect than high frequencies. For example, a 50 Hz wavefield has the strongest anisotropy effect if compared to 100 and 150 Hz wavefields. The attenuation anisotropy for low frequency (say 50 Hz) is more sensitive to the viscosity than the porosity and the Lame parameters have the weakest effect among these three parameters. These observations suggest that low-frequency seismic attenuation, and especially the attenuation anisotropy in low frequency, would have great potential for fluid discrimination within fractured media.

Decoupling of Pacific subduction zone guided waves beneath central Japan: Evidence for thin slab

AbstractThe fine‐scale seismic structure of the northeast Japan subduction zone is studied based on waveform analyses of moderate‐sized (M4.5–6), deep‐focus earthquakes (h >350 km) and the finite difference method (FDM) simulation of high‐frequency (up to 8 Hz) wave propagation. Strong regional S wave attenuation anomalies for specific source‐receiver paths connecting the cluster of events occurring in central part of the Sea of Japan recorded at fore arc stations in northern and central Japanese Islands (Honshu) are used to model the deeper structure of the subducting Pacific Plate, where recent teleseismic tomography has shown evidence for a possible slab tear westward beneath the Sea of Japan. The character of the observed anomalous S wave attenuation and the following high‐frequency coda can be captured with the two‐dimensional (2‐D) FDM simulation of seismic waves in heterogeneous plate model, incorporating the thinning of the plate at depth, which is also compared with other possible causes of dramatic attenuation of high‐frequency S wave due to low‐Q or much weaker heterogeneities in the slab. The results of simulation clearly demonstrate that the dramatic loss of high‐frequency S wavefield from the plate into the surrounding mantle occurred due to the variation in the plate geometry (i.e., thinning of the plate) at depth near the source rather than due to variation in physical properties, such as due to the lowered‐Q and weaker heterogeneities in the plate. The presence of such a thin zone defocuses the high‐frequency slab‐guided S wave energy from the subducting plate into the surrounding mantle and acts as a geometric antiwaveguide. Based on the sequence of simulation results obtained, we propose thinning of Pacific Plate at depth subducting beneath northeastern Japan, localized to central part of Honshu, in agreement with the observations.

Simulation of seismic waves at the earth's crust (brittle–ductile transition) based on the Burgers model

Abstract. The earth's crust presents two dissimilar rheological behaviors depending on the in situ stress-temperature conditions. The upper, cooler part is brittle, while deeper zones are ductile. Seismic waves may reveal the presence of the transition but a proper characterization is required. We first obtain a stress–strain relation, including the effects of shear seismic attenuation and ductility due to shear deformations and plastic flow. The anelastic behavior is based on the Burgers mechanical model to describe the effects of seismic attenuation and steady-state creep flow. The shear Lamé constant of the brittle and ductile media depends on the in situ stress and temperature through the shear viscosity, which is obtained by the Arrhenius equation and the octahedral stress criterion. The P and S wave velocities decrease as depth and temperature increase due to the geothermal gradient, an effect which is more pronounced for shear waves. We then obtain the P−S and SH equations of motion recast in the velocity-stress formulation, including memory variables to avoid the computation of time convolutions. The equations correspond to isotropic anelastic and inhomogeneous media and are solved by a direct grid method based on the Runge–Kutta time stepping technique and the Fourier pseudospectral method. The algorithm is tested with success against known analytical solutions for different shear viscosities. A realistic example illustrates the computation of surface and reverse-VSP synthetic seismograms in the presence of an abrupt brittle–ductile transition.

Seismic wave simulation by velocity–stress wave equations in two-phase anisotropic media

In two-phase anisotropic media, fast P-waves, slow P-waves and SV-waves are coupled with each other. Fast P-wavefronts show an elliptical anisotropic feature. Although the wavelet phase of the fast P-wavefront in a solid and that in a fluid are the same, the wavelet phases of the slow P-waves are opposite. In cases where there is an obvious slow P-wave in the wavefield, the slow P-wave in fluid is stronger than that in the solid. The attenuation mechanism introduced in Biot's theory is of secondary importance for the fast P-wave and the SV-wave, but dissipation coefficients have a significant effect on the slow P-wave. If there is a large dissipation in the media, the slow P-wave will be attenuated very quickly. When the dissipation coefficients have strong anisotropy, the slow P-wave is attenuated not just along the direction that has a big dissipation coefficient, but for the entire slow P-wavefield. Hence the slow P-wave is hardly observable in practice on surface seismic records. It is clearly a limitation of conventional seismic analysis which we should be aware of in reservoir geophysics.

Equivalent Boundary Conditions for Heterogeneous Acoustic Media

In this work, we have worked on possibilities to model artificial boundaries needed in the simulation of wave propagation in acoustic heterogeneous media. Our motivation is to restrict the computational domain in the simulation of seismic waves that are propagated from the earth and transmitted to the stratified heterogeneous media composed by ocean and atmosphere. Two possibilities were studied and compared in computational tests: the use of absorbing boundary conditions on an artificial boundary in the atmosphere layer and the elimination of the atmosphere layer using an equivalent boundary condition that mimics the propagation of waves through the atmosphere. <br />

A new kind of optimal second-order symplectic scheme for seismic wave simulations

Here we introduce generalized momentum and coordinate to transform seismic wave displacement equations into Hamiltonian system. We define the Lie operators associated with kinetic and potential energy, and construct a new kind of second order symplectic scheme, which is extremely suitable for high efficient and long-term seismic wave simulations. Three sets of optimal coefficients are obtained based on the principle of minimum truncation error. We investigate the stability conditions for elastic wave simulation in homogeneous media. These newly developed symplectic schemes are compared with common symplectic schemes to verify the high precision and efficiency in theory and numerical experiments. One of the schemes presented here is compared with the classical Newmark algorithm and third order symplectic scheme to test the long-term computational ability. The scheme gets the same synthetic surface seismic records and single channel record as third order symplectic scheme in the seismic modeling in the heterogeneous model.

A study of damping factors in perfectly matched layers for the numerical simulation of seismic waves

When simulating seismic wave propagation in free space, it is essential to introduce absorbing boundary conditions to eliminate reflections from artificially truncated boundaries. In this paper, a damping factor referred to as the Gaussian damping factor is proposed. The Gaussian damping factor is based on the idea of perfectly matched layers (PMLs). This work presents a detailed analysis of the theoretical foundations and advantages of the Gaussian damping factor. Additionally, numerical experiments for the simulation of seismic waves are presented based on two numerical models: a homogeneous model and a multi-layer model. The results show that the proposed factor works better. The Gaussian damping factor achieves a higher Signal-to-Noise Ratio (SNR) than previously used factors when using same number of PMLs, and requires less PMLs than other methods to achieve an identical SNR.