Simulation Of Seismic Waves Research Articles

Numerical Simulation of Seismic Transverse Waves on the Interface Between an Elastic Medium and Saturated Frozen Soil

Numerical Simulation of Seismic Transverse Waves on the Interface Between an Elastic Medium and Saturated Frozen Soil

Simulation of seismic waves in fluid-solid coupled thermoelastic media

This research advances seismic wave propagation analysis by incorporating temperature effects within fluid-solid coupled media. Leveraging the Lord-Shulman theory, we develop propagation equations that apply thermoacoustic and thermoelastic concepts to the intricacies of fluid-solid-thermal coupled media. This method facilitates the energy exchange between different regions by using distinct wave propagation equations in various media and implementing appropriate boundary conditions at the interfaces. Subsequently, we use high-order staggered grids and formulate corresponding perfectly matched layer absorbing boundaries to enhance simulation accuracy and reduce computational demands. This approach is validated through two numerical examples and rock-physics experiments, highlighting the crucial influence of thermal variations on subsurface wave dynamics and the method’s effectiveness in synthesizing seismic records and simulating converted wave information. The results indicate that our method provides a new approach for analyzing seismic wave propagation characteristics in fluid-solid coupled environments with temperature variations.

Acoustic wave simulation in strongly heterogeneous models using a discontinuous Galerkin method

In recent years, the discontinuous Galerkin method (DGM) has been rapidly developed for the numerical simulation of seismic waves. For wavefield propagation between two adjacent elements, it is common practice to apply a numerical flux to the boundary of each element to propagate waves between adjacent elements. Several fluxes, including the center, penalty, local Lax-Friedrich (LLF), upwind, and Rankine-Hugoniot jump condition-based (RH condition) fluxes, are widely used in numerical seismic wave simulation. However, some fluxes do not account for media differences between adjacent elements. Although different fluxes have been successfully used in DGM for many velocity models, it is unclear whether they can produce sufficiently accurate or stable results for strongly heterogeneous models, such as checkerboard models commonly used in tomographic studies. We test different fluxes using the acoustic wave equation. We analyze the accuracy of the penalty, LLF, upwind, and RH-condition fluxes based on the results of the numerical simulations of the homogeneous and two-layer models. We conduct simulations using checkerboard models, and the results indicate that the LLF, penalty, and upwind fluxes may have instability problems in heterogeneous models with long-time simulations. We observe instability issues in the LLF, penalty, and upwind fluxes when the wave-impedance contrast is high at the media interface. However, the results of the RH-condition flux remain consistently stable. The series of numerical examples presented in this work provide insights into the characteristics and application of fluxes for seismic wave modeling.

How Accurate Numerical Simulation of Seismic Waves in a Heterogeneous Medium Can Be?

ABSTRACT Analysis of equations of motion by Moczo et al. (2022) led to the conclusion that the discrete (grid) representation of the heterogeneous medium must be wavenumber bandlimited up to the Nyquist frequency. This is a consequence of the spatial discretization. Mittet (2021a) reported that if the discrete grid model of medium coincides with the true medium up to some wavenumber, the simulated wavefield is accurate only up to a half of this wavenumber. Here, we present results of the systematic and comprehensive analysis focused on the principal limits of accuracy of numerically simulated wavefields. First, we analyze wavenumber spectra of (1) exact wavefields in a heterogeneous elastic medium, (2) wavenumber bandlimited wavefields, and (3) spatially discretized wavefields. Then, we derive spatial dependence of the frequency spectrum of waves generated by a finite source, and perturbing wavefields due to a small perturbation of the medium and due to a small wavenumber bandlimited perturbation of the medium. We analyze an interaction of an incoming wave with the medium perturbation through a change of phase difference and through wavenumber spectra. We draw conclusions on the wavenumber limitation of wavefields in the wavenumber bandlimited heterogeneous medium. We numerically verify the fundamental finding using exact solutions. The main consequence for the finite-difference (FD) modeling based on spatial discretization of the computational domain is: Due to spatial sampling, the medium must be wavenumber limited up to the Nyquist frequency. Then, the wavefield should not be sampled by less than four spatial grid spacings per shortest wavelength to obtain sufficiently accurate results. This applies to any heterogeneous FD scheme.

Self-adaptive numerical dispersion suppressed method for shallow seismic simulation

Near-surface in seismic exploration generally refers to the low deceleration zone and a section of stratum below it, which can be described in geophysical language as “low velocity”, “low Q”, and “Free surface”, etc. Due to the presence of low-velocity layers in near-surface, shallow seismic simulation is plagued by numerical dispersion. Numerical dispersion affects the accuracy of seismic wave simulation, especially severe in shallow areas containing low-velocity layers. To improve the simulated quality, denser grids or higher-order difference schemes are used, but both of which would seriously reduce computational efficiency, especially for frequency-domain numerical modeling. Differentiation is replaced by difference in wave field simulation, which would inevitably result in numerical errors. These numerical errors are manifested as the difference between the phase velocity of the discretized wave field and the medium velocity. The waves with different wavenumber components have different phase velocities, and the larger the wavenumber, the more that its phase velocity lags the group velocity. Based on this theoretical analysis, a regularization factor is added to the frequency-domain acoustic wave equation to correct the phase velocity of high wavenumber components. According to the Von Neumann stability requirements, a regularization factor that adaptively changes with simulation parameters is derived. Compared to the fixed regularization factor, adaptive regularization factor can better match the velocity field and protect the effective wave field to the maximum extent while suppressing numerical dispersion. The improved acoustic wave equation can effectively suppress numerical dispersion without increasing the number of grids. Different numerical models demonstrate the effectiveness and efficiency of the improved acoustic equation with the adaptive regularization factor for suppressing numerical dispersion.

Cube2sph : A toolkit enabling flexible and accurate continental-scale seismic wave simulations using the SPECFEM3D_Cartesian package

To enable flexible and accurate seismic wave simulations at continental scales (10°−60°) based on the spectral-element method using the open-source SPECFEM3D_Cartesian package, we develop a toolkit, Cube2sph, that allows the generation of customized spherical meshes that account for the Earth’s curvature. This toolkit enables the usage of the perfectly matched layer (PML) absorbing boundary condition even when the artificial boundaries do not align with the coordinate axes. A series of numerical experiments are presented to validate the effectiveness of this toolkit. From these numerical experiments, we conclude that (1) continental-scale seismic wave simulations, especially surface wave simulations, can be more efficiently performed without the loss of accuracy by truncating the mesh at an appropriate depth, (2) curvilinear-grid PML can be used to effectively suppress artificial reflections for seismic wave simulations at continental scales, and (3) the Earth’s spherical geometry needs to be accurately meshed in order to obtain accurate simulation results for study regions larger than 8°.

Error Propagation and Control in 2D and 3D Hybrid Seismic Wave Simulations for Box Tomography

ABSTRACT To enhance the local resolution of global waveform tomography models, particularly in areas of interest within the Earth’s deep structures, a higher resolution localized tomography approach (referred to as “box tomography”) is crucial for a more detailed understanding of the Earth’s internal structure and geodynamics. Because the small-scale features targeted by box tomography are finer than those in global reference models, distinct spatial meshes are necessary for global and local (hybrid) forward simulations. Within the spectral element method (SEM) framework, we employ the intrinsic Lagrangian spatial interpolation to compute and store hybrid inputs (displacement/potential) in the global numerical simulation. These hybrid inputs are subsequently imposed into the localized domain during the iterative box tomography. However, inaccurate spatial Lagrange interpolation can lead to imprecise hybrid inputs, and this error can propagate from the global simulation to the hybrid simulation. It is essential to quantitatively analyze this error propagation and control it to ensure the credibility of box tomography. We introduce a unique spatial window function into the conventional “direct discrete differentiation” hybrid method. When the local mesh and structure align with those in the global simulation, the synthetic hybrid waveforms match the global ones, serving as a reference for quantitatively assessing error propagation stemming from changes in the local spatial mesh during hybrid simulation. Significantly, the relative waveform error arising due to spatial Lagrange interpolation is around 5% when employing the traditional SEM with five Gauss–Lobatto–Legendre points per minimum wavelength in the 3D global simulation through SPECFEM3D_GLOBE. Ultimately, we achieve hybrid waveforms with an accuracy of about 1.5% by increasing the spectral elements by about 1.5 times in the standard global simulation.

2D FDM Simulation of Seismic Waves and Tsunamis Based on Improved Coupling Equations Under Gravity

To understand the characteristics of seismic waves and tsunamis recorded simultaneously by the ocean-bottom observation networks, the coupling between the solid Earth and the ocean has to be modeled in the presence of gravity. However, previous coupled simulations adopted approximate equations that did not fully incorporate the effects of gravity. In this study, we derived correctly linearized governing equations under gravity and compared them with those of previous studies. Numerical experiments were performed for a two-dimensional P-SV wavefield, using the finite difference method (FDM). To validate the accuracy of the calculated tsunamis, we computed the theoretical tsunami dispersion relation using a propagator matrix and compared it with our results and those of previous studies. We found that our proposed method provided more accurate results than those of previous studies, particularly in the short-period band. We also investigated the applicability of the proposed method to distant tsunamis by examining the difference between calculated and theoretical tsunami phase velocities in the long-period band. The proposed formulation provides accurate results that properly incorporate gravity into the simultaneous simulation of seismic waves and tsunamis.

Comparative Study of 2D Lattice Boltzmann Models for Simulating Seismic Waves

The simulation of seismic wavefields holds paramount significance in understanding subsurface structures and seismic events. The lattice Boltzmann method (LBM) provides a computational framework adept at capturing detailed wave interactions, offering a new approach to improve seismic wavefield simulations. Our study involves a novel comparative analysis of wavefields using different lattice Boltzmann models, focusing on how relaxation times, discrete velocity models, and collision operators affect simulation accuracy and efficiency. We explore the impacts of distinct relaxation times and evaluate their effects on wave propagation speed and fidelity. By incorporating four discrete velocity models of LBM, we innovatively investigate the trade-off between spatial resolution and computational complexity. Additionally, we delve into the implications of employing three collision operators—single relaxation time (SRT), two relaxation times (TRT), and multiple relaxation times (MRT). By comparing their accuracy and stability, we provide insights into selecting the most suitable collision operator for capturing complex wave interactions. Our research provides a comprehensive framework to optimize the LBM parameters, enhancing both accuracy and efficiency in seismic wave simulations, and offers valuable insights to benefit wave simulation across diverse disciplines.

Double-pole unsplit complex-frequency-shifted multiaxial perfectly matched layer combined with strong-stability-preserved Runge-Kutta time discretization for seismic wave equation based on the discontinuous Galerkin method

The classical complex-frequency-shifted perfectly matched layer (CFS-PML) technique has attracted widespread attention for seismic wave simulations. However, few studies have addressed the double-pole variant of the CFS-PML scheme. The double-pole CFS-PML has a stronger capacity to absorb near-grazing incident waves and evanescent waves than the classical CFS-PML. Using the discontinuous Galerkin (DG) method, we derive a double-pole unsplit auxiliary ordinary differential equation CFS-multiaxial PML (AODE CFS-MPML) formulation, which combines a fourth-order strong-stability-preserved Runge-Kutta time discretization for wavefield simulation on an unstructured grid. The double-pole unsplit CFS-MPML formulations are obtained by introducing auxiliary memory variables and AODEs. The original stress-velocity equations and the double-pole unsplit AODE CFS-MPML equations are all first-order hyperbolic systems and suitable for the DG method. The attenuative variables are added directly to the original seismic wave equations without changing their formats. In contrast to the split perfectly matched layer (PML), we avoid reformulating PML equations in the nonattenuative modeling region. The original seismic wave equation is solved in the nonattenuative modeling domain, whereas the double-pole unsplit AODE CFS-MPML equation is implemented in the PML absorbing region. Three numerical examples validate the performance of the double-pole unsplit AODE CFS-MPML technique. The isotropic and anisotropic experiments demonstrate that our developed double-pole unsplit AODE CFS-MPML is more stable and obtains more accurate solutions than the classical CFS-PML. The second example indicates the flexibility of the combination of the DG method with the double-pole CFS-MPML on undulating topography. The final example displays the applicability and effectiveness of our method in a 3D situation.

Reducing the low-wavenumber dispersion error by building the Lagrange dual problem with a powerful local restriction

Abstract The broadband finite-difference (FD) coefficients computed by a cost function have been widely applied to suppress of numerical dispersion. Under the same conditions, the FD coefficients with small low-wavenumber dispersion error will produce a more accurate numerical solution in the long-time seismic wave simulation. Thus, how to reduce the low-wavenumber dispersion error becomes a crucial problem. According to the research into the zero-point position at the dispersion curve for three types of cost function, we found that the more zero points concentrate on the low-wavenumber region, the less the dispersion error. Therefore, the concentration of zero points is a good way to reduce dispersion error, which can be implemented by modified wavenumbers of zero points. Then, we design a Lagrange dual problem with a restriction based on the modified wavenumbers. The requirements for constructing the Lagrange dual problem are the optimization function and restricted condition, where the former is based on the dispersion relation, and the latter comprises the modified wavenumbers. The solution of this optimization problem, calculated by the dual ascent method, affords a less low-wavenumber dispersion error than the solution yielded by the alternating direction method of multipliers (ADMM). The theoretical analysis and numerical modeling suggest that the proposed method is superior to the existing optimal FD coefficients in reducing numerical error accumulation in low-frequency simulation.

Helmholtz-equation solution in nonsmooth media by a physics-informed neural network incorporating quadratic terms and a perfectly matching layer condition

Frequency-domain simulation of seismic waves plays an important role in seismic inversion, but it remains challenging in large models. The recently proposed physics-informed neural network (PINN), as an effective deep-learning method, has achieved successful applications in solving a wide range of partial differential equations (PDEs), and there is still room for improvement on this front. For example, PINN can lead to inaccurate solutions when PDE coefficients are nonsmooth and describe structurally complex media. Thus, we solve the acoustic and visco-acoustic scattered-field (Lippmann-Schwinger) wave equation in the frequency domain with PINN instead of the wave equation to remove the source singularity. We first illustrate that nonsmooth velocity models lead to inaccurate wavefields when no boundary conditions are implemented in the loss function. Then, we add the perfectly matched layer (PML) conditions in the loss function to better couple the real and imaginary parts of the wavefield. Moreover, we design new neurons by replacing the classical affine function with a quadratic function in the argument of the activation function to better capture the nonsmooth features of the wavefields. We find that the PML condition and the quadratic functions improve the results including handling attenuation and discuss the reason for this improvement. We also illustrate that a network trained to predict a wavefield for a specific medium can be used as an initial model of the neural network for predicting other wavefields corresponding to PDE-coefficient alterations and improving the convergence speed accordingly. This pretraining strategy should find applications in iterative full-waveform inversion and time-lag target-oriented imaging when the model perturbation between two consecutive iterations or two consecutive experiments is small.

Semi-analytical simulation of seismic waves generated by a point source in 1-D layered half-space of double-porosity media

SUMMARYIn this study, we present a semi-analytical method to simulate the propagation of seismic waves in horizontally stratified double-porosity media. We solve the governing equations in the frequency–wavenumber domain and then compute the time–space domain solutions by Hankel and fast Fourier transforms. We conduct numerical simulations to study the properties of the seismic waves in the double-porosity media. The results show the existence of three P waves, one fast P wave and two slow P waves. The two slow P waves are highly attenuated and can be observed by assuming a low fluid viscosity. By comparing the fast P wave in a single-porosity medium and that in a double-porosity medium, we find that the double-porosity model predicts much higher attenuation of the fast P wave due to the local fluid flow between the background medium and the inclusions. We also study the parameters affecting the attenuation, such as the radius of the spherical inclusions, the viscosity of the pore fluid, the permeability of the background medium and the porosity of the inclusions. Finally, we present a cross-well survey model to investigate the seismic responses in the double-porosity media and compare them with the responses in the single-porosity media. We find that both the amplitude and velocity of the fast P wave decrease with the increase of the volume ratio of the inclusions.

Unified framework based parallel FEM code for simulating marine seismoacoustic scattering

The simulation of seismic wave propagation in marine areas, which belongs to seismoacoustic scattering problem, is complicated due to the fluid-solid interaction between seawater and seabed, especially when the seabed is saturated with fluid. Meanwhile, huge computation resources are required for large-scale marine seismic wave simulation. In the paper, an efficient parallel simulation code is developed to solve the near-field seismoacoustic scattering problem. The method and technologies used in this code includes: 1) Unified framework for acoustic-solid-poroelastic interaction analysis, in which seawater and dry bedrock are considered as generalized saturated porous media with porosity equals to one and zero, respectively, and the coupling between seawater, saturated seabed and dry bedrock can be analyzed in the unified framework of generalized saturated porous media and avoid interaction between solvers of different differential equation; 2) Element-by-element strategy and voxel finite-element method (VFEM), with these strategies, it only needs to calculate several classes of element matrix and avoid assembling and storing the global system matrices, which significantly reduces the amount of memory required; 3) Domain-partitioning procedure and parallel computation technology, it performs 3D and 2D model partitioning for the 3D and 2D codes respectively, sets up the velocity structure model for the partitioned domain on each CPU or CPU core, and calculates the seismic wave propagation in the domain using Message Passing Interface data communication at each time step; 4) Local transmitting boundary condition, we adopt multi-transmitting formula, which is independent of specific wave equations, to minimize reflections from the boundaries of the computational model. A horizontal layered model with the plane P-wave incident vertically from bottom is used to demonstrate the computational efficiency and accuracy of our code. Then, the code is used to simulate the wave propagation in Tokyo Bay. All codes were written following to the standards of Fortran 95.

Numerical Solver-Independent Seismic Wave Simulation Using Task-Decomposed Physics-Informed Neural Networks

Numerical Solver-Independent Seismic Wave Simulation Using Task-Decomposed Physics-Informed Neural Networks

Analysis on the seismic wave caused by low frequency sound source in shallow sea based on multi-transmitting formula artificial boundary

The elastic waves propagating in seabed caused by sailing ships are called ship seismic waves, which can be used to identify ship targets. The wave components and the influences of source frequency, source depth and depth of seawater to the seismic waves are important to the application of seismic wave to detect ship targets. Thus, a forward numerical simulation of seismic waves in shallow sea excited by low frequency sources in time domain was carried out with finite element method based on the Multi-Transmitting Formula artificial boundary. The numerical results show that the body waves and conical waves decay faster in space, but the interface wave decays relatively slow. Multi-transmitting Formula (MTF) artificial boundary has achieved good transmitting effect for longitudinal wave, upstream and downstream acoustic wave and transverse wave. When the source frequency is very low and the seawater layer is shallow, due to the effect of low-frequency cutoff, there is no normal mode wave which can propagate without attenuation in the seawater waveguide, and only the interface wave is propagating near the sea bottom surface. When the frequency of the source or the depth of the seawater layer increases, the effect of low frequency cutoff in the shallow water waveguide is weakened, and the normal mode waves propagating in seawater layer are gradually excited. The intensity of the interface wave caused by low-frequency point sound source mainly depends on the distance between the source and the seafloor.

Basin effect quantitation of four deep basins in Japan to explain the gap between empirical and numerical basin‐depth effect models

AbstractTo incorporate basin effects into seismic hazard assessment of urban areas and infrastructure, basin‐depth effect models that are the dependence of basin amplification factor on basin depth have been derived via numerical and empirical approaches. However, commonly observed quantitative and qualitative differences between numerical and empirical models remain unresolved. In particular, empirical models tend to predict smaller basin amplification factors. In this study, a modified empirical approach from Choi et al. was applied to the derivation of basin‐depth effect models for four deep basins in Japan using a database consisting of 71MJ > 6 (Mw > 5.7) earthquakes and 13,562 records from strong‐motion seismograph networks (K‐NET and KiK‐net). The obtained basin‐depth models vary among basins and crustal‐subduction earthquake pairs, which confirms the recent trend of accounting for regional‐seismotectonic basin amplification in seismic hazard analysis. The basin amplification observed for crustal earthquakes is generally larger than that for subduction earthquakes. Moreover, the differences between the numerical and empirical models can be explained by changing the options of basin or outside‐basin stations for computing the inter‐event residuals. These results bridge a knowledge gap in seismic hazard analysis between ground‐motion prediction equations and those using large‐scale seismic wave simulations, which are now the main tools for regions lacking seismic records.

High-frequency wavefield extrapolation using the Fourier neural operator

Abstract In seismic wave simulation, solving the wave equation in the frequency domain requires calculating the inverse of the impedance matrix. The total cost strictly depends on the number of frequency components that are considered, if using a finite-difference method. For the applications such as seismic imaging and inversion, high-frequency information is always required and thus the wave simulation is always a challenging task as it demands tremendous computational cost for obtaining dispersion-free high-frequency wavefields for large subsurface models. This paper demonstrates that a data-driven machine learning method, called the Fourier neural operator (FNO), is capable of predicting high-frequency wavefields, based on a limited number of low-frequency components. As the FNO method is for the first time applied to seismic wavefield extrapolation, the experiment reveals three attractive features with FNO: high efficiency, high accuracy and, importantly, the predicted high-frequency wavefields are dispersion free.

Local wavefield refinement using Fourier interpolation and boundary extrapolation for finite-element method based on domain reduction method

The domain reduction method based on the finite-element method (FEM) is promising for multiscale seismic wave simulations for local complex structures and rough topography. However, a transition region is required to refine the wavefields from coarse to fine grids, which may lead to strong artifacts. Traditional wavefield interpolation methods can only handle a small upsampling ratio due to their low accuracy. Here, we propose a high accuracy upsampling method using Fourier interpolation and boundary extrapolation. It requires that the wavefield in the transition region is evenly sampled, which may slightly reduce the flexibility of FEMs. Fortunately, this constraint only exists in two layers of nodes, in which the coarse and fine grids are collocated. To further apply our interpolation method to the wavefield upsampling on the nonuniform grid, we adopt cubic interpolation within the local region of interest, which is adequate in accuracy for regularizing the wavefield on the fine nonuniform grid. For a topographic surface, we perform wavefield extrapolation on the ground surface to avoid abnormal amplitude variations below and beyond the surface, which can greatly reduce potential artifacts caused by the Fourier interpolation. Numerical experiments find that our method is superior to traditional interpolation methods in accuracy, especially for large upsampling ratios. Its relative error is less than 2% even for a 40-time upsampling; in contrast, the relative error of linear and cubic interpolation is up to 90% and 41% in the same situation, respectively. We further verify the feasibility of our method for 3D heterogeneous models with rough topography, multilayer heterogeneity, and random perturbations of background models. The proposed method yields a significant speeding of the multiscale seismic simulation using the FEM for 3D models with complex local structures and topographies while being highly accurate.

Theoretical and quantitative evaluation of hybrid PML-ABCs for seismic wave simulation

A good artificial boundary treatment in a seismic wave grid-based numerical simulation can reduce the size of the computational region and increase the computational efficiency, which is becoming increasingly important for seismic migration and waveform inversion tasks requiring hundreds or thousands of simulations. Two artificial boundary techniques are commonly used: perfectly matched layers (PMLs), which exhibit the excellent absorption performance but impose a greater computational burden by using finite layers to gradually reduce wave amplitudes; and absorbing boundary conditions (ABCs), which have the high computational efficiency but are less effective in absorption because they employ the one-way wave equation at the exterior boundary. Naturally, PMLs have been combined with ABCs to reduce the number of PMLs, thus improving the computational efficiency; many studies have proposed such hybrid PMLs. Depending on the equations from which the ABCs are derived, there are two hybrid PML variants: the PML+unstretched ABC (UABC), in which the ABC is derived from a physical equation; or the PML+stretched ABC (SABC), in which the ABC is derived from the PML equation. Even though all the previous studies concluded that hybrid PMLs can improve the absorption performance, none of them quantified how many PMLs can be removed by combining the PML with the ABC compared with the pure PML. In this paper, we systematically study the absorption performance of the two hybrid PML variants. We develop a method to distinguish the artificial reflections from the PML-interior interface and those caused by the PML exterior boundary to accurately approximate the additional absorption achieved by using the UABC and the SABC. The reflection coefficients based on a theoretical derivation and numerical tests both show that the UABC amplifies most reflections and is not recommended in any situation; conversely, the SABC can always diminish reflections, but the additional absorption achieved by the SABC is relatively poor and cannot effectively reduce the number of PMLs. In contrast, we find that simply increasing the damping parameter improves absorption better than the PML+SABC. Our results show that the improvement in absorption achieved by combining the PML with either the SABC or the UABC is not better than that obtained by simply adjusting the damping profile of the PML; thus, combining the PML with the ABC is not recommended in practice.