Artificial Boundary Reflections Research Articles

Numerical simulation of 3-D seismic wave based on alternative flux finite-difference WENO scheme

SUMMARY High-frequency non-physical oscillations may occur due to shock waves in seismic wavefield and dynamic rupture simulation. In this study, we introduced the alternative flux finite-difference weighted essentially non-oscillatory scheme to address potential shock wave issues in computational seismology effectively. The wavefield of the body-fitted curvilinear domain was accurately computed through conservative grid mapping, ensuring accurate implementation of free surface boundary conditions on irregular surfaces using characteristic boundary conditions and minimizing artificial boundary reflections with exponential decay absorbing layers. Finally, we compared our scheme with the GRTM for flat surfaces and the CGFDM3D-EQR for irregular surfaces to demonstrate its correctness and accuracy, and validated its non-oscillatory characteristics. The aforementioned scheme is anticipated to assume a significant function in simulating more intricate seismic wavefields or dynamic ruptures.

Hybrid absorbing boundary condition based on transmitting boundary and its application in 3D fractional viscoacoustic modeling

An accurate numerical simulation for wave equations is essential for understanding of wave propagation in the earth's interior as well as full waveform inversion and reverse time migration. However, due to computational cost and hardware capability limitations, numerical simulations are often performed within a finite domain. Thus, an adequate absorbing boundary condition (ABC) is indispensable for obtaining accurate numerical simulation results. In this study, we develop a hybrid ABC based on a transmitting boundary, which is referred to as THABC, to eliminate artificial boundary reflections in 3D second-order fractional viscoacoustic numerical simulations. Furthermore, we propose an adaptive weighted coefficient to reconcile the transmitting and viscoacoustic wavefields in THABC. Through several numerical examples, we determine that the proposed THABC approach is characterized by the following benefits. First, with the same number of absorbing layers, THABC exhibits a better ability in eliminating boundary reflection than traditional ABC schemes. Second, THABC is more effective in computation, since it only requires the wavefields at the current and last time steps to solve the transmitting formula within the absorbing layers. Benefiting from a simple but effective combination between the transmitting equation and the second-order wave equation, our scheme performs well in the 3D fractional Laplacian viscoacoustic numerical simulation.

The Combined Compact Difference Scheme Applied to Shear-Wave Reverse-Time Migration

In this paper, the combined compact difference scheme (CCD) and the combined supercompact difference scheme (CSCD) are used in the numerical simulation of the shear-wave equation. According to the Taylor series expansion and shear-wave equation, the fourth-order discrete scheme of the displacement field is established; then, the CCD and CSCD schemes are used to calculate the spatial derivative of the displacement field. Additionally, the accuracy, dispersion, and stability of the CCD and CSCD are analyzed, and numerical simulation analyses are carried out using 1D uniform models. Lastly, based on the processing of artificial boundary reflection using PML boundary conditions, shear-wave reverse-time migrations are carried out using synthetic data. The results show that (1) CCD and CSCD have smaller truncation errors, higher simulation precision, and lower numerical dispersion than other normal difference schemes; (2) CCD and CSCD can use the coarse grid and larger time step to calculate, with less memory and high computational efficiency; (3) finally, the result of the shear-wave reverse-time migration of the 2D synthetic data model show that the reverse-time migration imaging is clear, and the proposed method for shear-wave reverse-time migration is practical and effective.

A new qP-wave approximation in tilted transversely isotropic media and its reverse time migration for areas with complex overburdens

During seismic imaging in anisotropic media, undesired SV waves may be detected when conventional wave equations are used for wavefield simulation, which will contaminate the final seismic image. A pure qP-wave equation that is free of S-wave artifacts is desired, such that crosstalk introduced from the interference of different waves can be eliminated. We have developed a new fully decoupled P-wave equation in transversely isotropic media with titled symmetry axis that avoids unwanted energy during migration. The qP and qSV phase velocities for the new approximations are plotted and compared with the exact solution and some other approximations and indicate high accuracy for various anisotropic models (from low to high degrees of anisotropy). During wavefield simulation, the finite-difference and pseudospectral methods also are combined to accelerate the computation. Moreover, a hybrid perfectly matched layer (H-PML) that combines two different perfectly matched layers for second-order wave equations in the wavenumber domain also is proposed and applied to eliminate the artificial boundary reflections. Comparisons of different absorbing boundary layers also are illustrated to validate the wavenumber domain H-PML. Finally, the new proposed wave approximation together with the new absorbing layer is further incorporated in a new 3D anisotropic reverse time migration (RTM). We test the RTM on the 3D synthetic SEG Advanced Modeling Corporation data and offshore field data with a complicated overburden (salt and gas overburden) and obtain seismic images with high resolution. Benchmarks against commercial implementations also are demonstrated, which proves that the positioning of the reflectors is more reliable and accurate with the new algorithm.

An improved 25‐point finite‐difference scheme for frequency‐domain elastic wave modelling

ABSTRACTFrequency‐domain finite‐difference modelling is widely used in exploration geophysics. However, it involves solving a large system of linear algebraic equations, which may cause huge computational costs, especially for large‐scale models. To address this issue, we have proposed an improved 25‐point finite‐difference scheme for wavefield modelling of two‐dimensional frequency‐domain elastic wave equations. The biggest difference between the improved 25‐point scheme with other 25‐point schemes is the finite‐difference formulas for spatial derivatives. The proposed 25‐point scheme applies to equal and unequal grid intervals, and its optimization coefficients depend on the grid‐spacing ratio and Poisson's ratio. The dispersion analysis indicates that within the phase velocity errors of 1% and 2%, the improved 25‐point scheme only requires approximately 2.3 and 2.2 grid points per shear wavelength, which has greater accuracy than the existing schemes. To eliminate artificial boundary reflections, we apply the perfectly matched layer boundary conditions to the model edges. Several numerical examples are presented to prove the validity of the improved 25‐point scheme.

Trapezoid-Grid Finite-Difference Time-Domain Method for 3D Seismic Wavefield Modeling Using CPML Absorbing Boundary Condition

The large computational memory requirement is an important issue in 3D large-scale wave modeling, especially for GPU calculation. Based on the observation that wave propagation velocity tends to gradually increase with depth, we propose a 3D trapezoid-grid finite-difference time-domain (FDTD) method to achieve the reduction of memory usage without a significant increase of computational time or a decrease of modeling accuracy. It adopts the size-increasing trapezoid-grid mesh to fit the increasing trend of seismic wave velocity in depth, which can significantly reduce the oversampling in the high-velocity region. The trapezoid coordinate transformation is used to alleviate the difficulty of processing ununiform grids. We derive the 3D acoustic equation in the new trapezoid coordinate system and adopt the corresponding trapezoid-grid convolutional perfectly matched layer (CPML) absorbing boundary condition to eliminate the artificial boundary reflection. Stability analysis is given to generate stable modeling results. Numerical tests on the 3D homogenous model verify the effectiveness of our method and the trapezoid-grid CPML absorbing boundary condition, while numerical tests on the SEG/EAGE overthrust model indicate that for comparable computational time and accuracy, our method can achieve about 50% reduction on memory usage compared with those on the uniform-grid FDTD method.

PML and CFS-PML boundary conditions for a mesh-free finite difference solution of the elastic wave equation

Mesh-free finite difference (FD) methods can improve the geometric flexibility of modeling without the need for lattice mapping or complex meshing process. Radial-basis-function-generated FD is among the most commonly used mesh-free FD methods and can accurately simulate seismic wave propagation in the non-rectangular computational domain. In this paper, we propose a perfectly matched layer (PML) boundary condition for a mesh-free FD solution of the elastic wave equation, which can be applied to the boundaries of the non-rectangular velocity model. The performance of the PML is, however, severely reduced for near-grazing incident waves and low-frequency waves. We thus also propose the complex-frequency-shifted PML (CFS-PML) boundary condition for a mesh-free FD solution of the elastic wave equation. For two PML boundary conditions, we derive unsplit time-domain expressions by constructing auxiliary differential equations, both of which require less memory and are easy for programming. Numerical experiments demonstrate that these two PML boundary conditions effectively eliminate artificial boundary reflections in mesh-free FD simulations. When compared with the PML boundary condition, the CFS-PML boundary condition results in better absorption for near-grazing incident waves and evanescent waves.

An Optimized Choice of UCPML to Truncate Lattices With Rotated Staggered Grid Scheme for Ground Penetrating Radar Simulation

Efficient and accurate simulation of ground penetrating radar (GPR) in the open region helps immensely in both grasping the features of echoes and facilitating the interpretation of real GPR data. Due to the limitation of the computer model, however, the strong artificial boundary reflections, especially the low-frequency propagating waves encountered at the late stage of simulation greatly affect the simulation accuracy of GPR. This paper presents an innovative optimized unsplit-field convolutional perfectly matched layer (UCPML) based on rotated staggered grid (RSG) scheme to truncate the finite-difference time-domain (FDTD) lattices. Rather than obey the sharp variation based on an ${m}$ th-order polynomial, the optimized approach employs a novel optimized term and an adjustment factor to seek a gentle variation on optimal constitutive coefficients. This guarantees that the determination of optimal constitutive coefficients can be less influenced by the order of polynomial and especially, to improve the absorptive performance on low-frequency propagating waves. The calculating efficiency and accuracy of the RSG-FDTD scheme, as well as the absorbing performance of the optimized UCPML, are verified by two numerical examples. In particular, the analysis of the amplitude-frequency features of low-frequency clutters at steady state of the electromagnetic (EM) field and the corresponding global reflection error in the time–frequency domain is also presented.

Trapezoid coordinate finite difference modeling of acoustic wave propagation using the CPML boundary condition

Due to compaction of clastic sedimentary rocks caused by gravitational forces, the wave propagation velocity tends to gradually increase with depth. Based on this observation, a trapezoid coordinate system is built to integrate the general velocity variation, and an acoustic wave equation is derived in the new coordinate system. Numerically, a conventional uniform rectangular grid FD scheme can be evoked to solve the new equation. Physically, a variable grid mesh is used: that is, fine grids for shallow areas with low velocity and coarse grids for deep areas with high velocity. This method is free of artificial reflections caused by the transition of different grid sizes. Furthermore, the Convolutional Perfectly Matched Layer (CPML) boundary condition is adapted to eliminate artificial boundary reflections. A homogeneous velocity model is used to verify the effectiveness of the CPML boundary. Tests using the Marmousi benchmark model show that with the demonstrated comparable accuracy, the proposed method is 2.9 times as fast as the conventional physical uniform grid FD algorithm while saving as much as 75% of the computer memory.

Perfectly matched layer boundary conditions for frequency‐domain acoustic wave simulation in the mesh‐free discretization

ABSTRACTMesh‐free discretization, flexibly distributing nodes without computationally expensive meshing process, is able to deal with staircase problem, oversampling and undersampling problems and saves plenty of nodes through distributing nodes suitably with respect to irregular boundaries and model parameters. However, the time‐domain mesh‐free discretization usually exhibits poorer stability than that in regular grid discretization. In order to reach unconditional stability and easy implementation in parallel computing, we develop the frequency‐domain finite‐difference method in a mesh‐free discretization, incorporated with two perfectly matched layer boundary conditions. Furthermore, to maintain the flexibility of mesh‐free discretization, the nodes are still irregularly distributed in the absorbing zone, which complicates the situation of artificial boundary reflections. In this paper, we implement frequency‐domain acoustic wave modelling in a mesh‐free system. First, we present the perfectly matched layer boundary condition to suppress spurious reflections. Moreover, we develop the complex frequency shifted–perfectly matched layer boundary condition to improve the attenuation of grazing waves. In addition, we employ the radial‐basis‐function‐generated finite difference method in the mesh‐free discretization to calculate spatial derivatives. The numerical experiment on a rectangle homogeneous model shows the effectiveness of the perfectly matched layer boundary condition and the complex frequency shifted–perfectly matched layer boundary condition, and the latter one is better than the former one when absorbing large angle incident waves. The experiment on the Marmousi model suggests that the complex frequency shifted–perfectly matched layer boundary condition works well for complicated models.

An improved hybrid absorbing boundary condition for wave equation modeling

A recently developed hybrid absorbing boundary condition (ABC) can significantly suppress artificial boundary reflections by inserting a transition area between the boundary and the inner area. The wavefield in the transition area is computed by linearly weighting the wavefield from one-way wave equation (OWWE) and that from two-way wave equation. This leads to a smooth variation from the computational area to the boundary via the transition area and thus greatly reduces boundary reflections. In this paper, we propose two techniques to further enhance the absorbing effect. First, we widen the boundary from one point to several points. Second, we adopt nonlinear weighting coefficients instead of linear ones. Numerical simulations with 2D acoustic finite-difference modeling, combined with Clayton–Engquist OWWE, demonstrate that the proposed techniques can significantly improve the absorption effect without increasing computational cost. Tests using Higdon OWWE-based ABC suggest that the second technique can greatly enhance the absorption effect. Experiments with frequency-domain finite-difference modeling demonstrate similar improvement in absorption effect.

A new high-order scheme based on numerical dispersion analysis of the wave phase velocity for semidiscrete wave equations

In the numerical computation of wave equations, numerical dispersion is a persistent problem arising from inadequate discretization of the continuous wave equation. To thoroughly understand the mechanism of numerical dispersion, we separately analyze the numerical dispersion relations of time-stepping and spatial discretization schemes by Fourier analysis. The relevant results show that the numerical dispersion errors of time-marching schemes depend on the time step length or the Courant number, whereas the numerical dispersion errors of spatial discretization schemes are determined by the error between the eigenvalues of the numerical spatial differential operator and the continuous spatial differential operator. We also find that the much better numerical dispersion accuracy of the stereo-modeling discrete (SMD)-type operator can be attributed to the inclusion of diversified basis functions for its eigenvalue. Based on these findings, we combine the optimal four-stage symplectic partitioned Runge-Kutta and eighth-order SMD as a new fully discrete scheme. The subsequent analysis of its normalized phase velocity is consistent with the numerical dispersion analysis in semidiscrete forms. This is followed by an acoustic wave simulation in a homogeneous model and corresponding computational efficiency comparison. The results show that the new scheme is much more accurate and antidispersive on a coarse grid. In the final two numerical experiments, we use the new scheme to model the acoustic wave propagation in a three-layer model and the Marmousi model. The convolutional perfectly matched layer is applied to eliminate artificial boundary reflections. Our semidiscrete numerical dispersion analysis provides an efficient tool to quantitatively evaluate the time-stepping and spatial discretization schemes. It can facilitate the development of more accurate fully discrete numerical schemes for solving seismic wave equations.

A multi-block finite difference method for seismic wave equation in auxiliary coordinate system with irregular fluid–solid interface

PurposeThe purpose of this paper is to develop a new finite difference method for solving the seismic wave propagation in fluid-solid media, which can be described by the acoustic and viscoelastic wave equations for the fluid and solid parts, respectively.Design/methodology/approachIn this paper, the authors introduced a coordinate transformation method for seismic wave simulation method. In the new method, the irregular fluid–solid interface is transformed into a horizontal interface. Then, a multi-block coordinate transformation method is proposed to mesh every layer to curved grids and transforms every interface to horizontal interface. Meanwhile, a variable grid size is used in different regions according to the shape and the velocity within each region. Finally, a Lebedev-standard staggered coupled grid scheme for curved grids is applied in the multi-block coordinate transformation method to reduce the computational cost.FindingsThe instability in the auxiliary coordinate system caused by the standard staggered grid scheme is resolved using a curved grid viscoelastic wave field separation strategy. Several numerical examples are solved using this new method. It has been shown that the new method is stable, efficient and highly accurate in solving the seismic wave equation defined on domain with irregular fluid–solid interface.Originality/valueFirst, the irregular fluid–solid interface is transformed into a horizontal interface by using the coordinate transformation method. The conversion between pressures and stresses is easy to implement and adaptive to different irregular fluid–solid interface models, because the normal stress and shear stress vanish when the normal angle is 90° in the interface. Moreover, in the new method, the strong false artificial boundary reflection and instability caused by ladder-shaped grid discretion are resolved as well.

Wavefield simulation of 3D borehole dipole radiation

Reflections in acoustic borehole-logging data can be used to image near-borehole geologic structures, using either monopole or dipole measurements, the latter of which resolve azimuthal ambiguities. Dipole methods in particular can characterize small and subtle fractures in reservoir environments. Accurate tools for wavefield simulation are required to account for wave modes in and out of the borehole, and for the receiver responses to these modes. We have evaluated a 3D elastic staggered-grid finite-difference method for isotropic and anisotropic media in the context of the borehole acoustic imaging problem. Based on the convolutional perfectly matched layer and the multiaxial perfectly matched layer schemes, a hybrid perfectly matched layer method was developed and implemented to reduce artificial boundary reflections. The superiority of this hybrid perfectly matched layer beyond the convolutional perfectly matched layer and the multiaxial perfectly matched layer was proven based on numerical benchmarks in isotropic and anisotropic media. Numerical simulation of radiation, reflection, and multipole reception of elastic waves, excited by a dipole source, was carried out. To optimize azimuthal detection, the relationships between S-wave polarization as well as source-receiver offset and source-reflector angle was analyzed. Results indicated that the S-S reflection is most sensitive to the angle between the incident ray and the normal to the reflector. Its maximum amplitude occurs as the incident angle reaches its critical value, a fact that can be used to calculate the total propagation distance of the S-S wave. The critical angle as well as the SH-wave velocities of geologic structures outside the borehole can thus be determined.

Acoustic and elastic modeling by optimal time-space-domain staggered-grid finite-difference schemes

Staggered-grid finite-difference (SFD) methods are widely used in modeling seismic-wave propagation, and the coefficients of finite-difference (FD) operators can be estimated by minimizing dispersion errors using Taylor-series expansion (TE) or optimization. We developed novel optimal time-space-domain SFD schemes for acoustic- and elastic-wave-equation modeling. In our schemes, a fourth-order multiextreme value objective function with respect to FD coefficients was involved. To yield the globally optimal solution with low computational cost, we first used variable substitution to turn our optimization problem into a quadratic convex one and then used least-squares (LS) to derive the optimal SFD coefficients by minimizing the relative error of time-space-domain dispersion relations over a given frequency range. To ensure the robustness of our schemes, a constraint condition was imposed that the dispersion error at each frequency point did not exceed a given threshold. Moreover, the hybrid absorbing boundary condition was applied to remove artificial boundary reflections. We compared our optimal SFD with the conventional, TE-based time-space-domain, and LS-based SFD schemes. Dispersion analysis and numerical simulation results suggested that the new SFD schemes had a smaller numerical dispersion than the other three schemes when the same operator lengths were adopted. In addition, our LS-based time-space-domain SFD can obtain the same modeling accuracy with shorter spatial operator lengths. We also derived the stability condition of our schemes. The experiment results revealed that our new LS-based SFD schemes needed a slightly stricter stability condition.

A mixed-grid finite element method with PML absorbing boundary conditions for seismic wave modelling

We have developed a mixed-grid finite element method (MGFEM) to simulate seismic wave propagation in 2D structurally complex media. This method divides the physical domain into two subdomains. One subdomain covering the major part of the physical domain is divided by regular quadrilateral elements, while the other subdomain uses triangular elements to correctly fit a rugged free surface topography. The local stiffness matrix of any quadrilateral element is identical and matrix-vector production is calculated using an element-by-element technique, which avoids assembling a huge global stiffness matrix. As only a few triangular elements exist in the subdomain containing the rugged free surface topography, the memory requirements for storing the assembled subdomain global stiffness matrix are significantly reduced. To eliminate artificial boundary reflections, the MGFEM is also implemented to solve the system equations of PML absorbing boundary conditions (PML ABC). The accuracy and efficiency of the MGFEM is tested in numerical experiments by comparing it with conventional methods, and numerical comparisons also indicate its tremendous ability to describe rugged surfaces.

Elastic modelling in tilted transversely isotropic media with convolutional perfectly matched layer boundary conditions

To simulate wave propagation in a tilted transversely isotropic (TTI) medium with a tilting symmetry-axis of anisotropy, we develop a 2D elastic forward modelling algorithm. In this algorithm, we use the staggered-grid finite-difference method which has fourth-order accuracy in space and second-order accuracy in time. Since velocity-stress formulations are defined for staggered grids, we include auxiliary grid points in the z-direction to meet the free surface boundary conditions for shear stress. Through comparisons of displacements obtained from our algorithm, not only with analytical solutions but also with finite element solutions, we are able to validate that the free surface conditions operate appropriately and elastic waves propagate correctly. In order to handle the artificial boundary reflections efficiently, we also implement convolutional perfectly matched layer (CPML) absorbing boundaries in our algorithm. The CPML sufficiently attenuates energy at the grazing incidence by modifying the damping profile of the PML boundary. Numerical experiments indicate that the algorithm accurately expresses elastic wave propagation in the TTI medium. At the free surface, the numerical results show good agreement with analytical solutions not only for body waves but also for the Rayleigh wave which has strong amplitude along the surface. In addition, we demonstrate the efficiency of CPML for a homogeneous TI medium and a dipping layered model. Only using 10 grid points to the CPML regions, the artificial reflections are successfully suppressed and the energy of the boundary reflection back into the effective modelling area is significantly decayed.

2-D elastic wave modeling with frequency-space 25-point finite-difference operators

Numerical simulation in the frequency-space domain has inherent advantages, such as: it is possible to simulate wave propagation from multiple sources simultaneously; there are no cumulative errors; only the interesting frequencies can be selected; and it is more suitable for wave propagation in viscoelastic media. The only obstacle to using the method is the requirement of huge computer storage. We extend the compressed format for storing the coefficient matrix. It can reduce the required computer storage dramatically. We get the optimal coefficients by least-squares method to suppress the numerical dispersion and adopt the perfectly matched layer (PML) boundary conditions to eliminate the artificial boundary reflections. Using larger grid intervals decreases computer storage requirements and provides high computational efficiency. Numerical experiments demonstrate that these means are economic and effective, providing a good basis for elastic wave imaging and inversion.

2.5D modelling of elastic waves in transversely isotropic media using the spectral element method

The spectral-element method provides an accurate alternative to the finite-element method for modelling elastic waves in anisotropic media. With the aim of reducing the high computational overheads of 3D modelling, we have implemented 2.5D spectral-element modelling of elastic waves, initially for vertically transversely isotropic media, and then extended to a tilted transversely isotropic medium with a dipping symmetry-axis. We have investigated different categories of absorbing boundaries to minimise artificial boundary reflections, including viscous boundary conditions and two distinct 2.5D formulations of perfectly matched layers (PMLs). Both PML methods use complex coordinate stretching, but the first method is only applicable in the frequency-domain. The other method uses the decomposed gradient operator and can be applied in both the time-domain and frequency-domain. Traditional 2D modelling exhibits cylindrical wave propagation characteristics from a line source. Our 2.5D modelling results, presented as images of seismic wavefields, illustrate the more realistic spherical propagation from a point source.

2.5D modelling of elastic waves in anisotropic media using the spectral-element method – a preliminary investigation

The spectral-element method provides an accurate alternative to the finite element method for modelling elastic waves in anisotropic media. However its high computational overheads make it impractical for inversion. With the aim of reducing the overheads, we have implemented 2.5D spectral-element modelling of elastic waves, initially for vertically transversely isotropic (VTI) media. We have incorporated perfectly matched layers to minimise artificial boundary reflections, using the split gradient operator for the time domain and complex coordinate stretching for the frequency domain. Traditional 2D modelling does not reflect the true coupling between wave modes. Our preliminary results, presented as images of seismic waveforms, illustrate the expected coupling. We also present an intuitive explanation of our unique implementation of the spectral-element method.