Staggered-grid Finite-difference Method Research Articles

Decoupled acoustoelasticity equations and stable FD modeling for wave propagation in prestressed media

Seismic exploration of deep oil/gas reservoirs involves wave propagation in high-pressure media. Acoustoelasticity with the third-order elastic constants can account for nonlinear strain responses due to stresses of finite magnitude. Since current seismic exploration depends primarily on P waves, further development in prospecting technology to handle deep reflection data will benefit from the ability to decouple the elastic waves from high-pressure deep formations. We decouple the first-order velocity-stress acoustoelasticity equations for wave propagation under isotropic (confining) and anisotropic (uniaxial and pure shear) prestress conditions. A rotated staggered-grid finite-difference (RSG-FD) method is used to solve the decoupled acoustoelasticity equations with an unsplit convolutional perfectly matched layer (CPML) absorbing boundary. Numerical examples demonstrate the significant impact of prestressed conditions on seismic responses in both phase and amplitude. The stress-induced seismic anisotropy with orthotropic characteristics is strongly related to the orientation of prestresses. We compare the coupled and the decoupled wavefields to validate the proposed decoupling method. Investigation of the influence of prestresses on the acoustoelastic decoupling illustrates that the isotropic prestress does not affect the decoupling, but the anisotropic prestresses significantly affect the decoupling accuracy.

Propagation Characteristics of Seismic Waves in Shallow Sea Sedimentary Layers

Zhao, S.; Gu, D.; Guo, C.; Wang, X., and Zhang, J., 2023. Propagation characteristics of seismic waves in shallow sea sedimentary layers. Journal of Coastal Research, 39(1), 32–44. Charlotte (North Carolina), ISSN 0749-0208. To explore new techniques for detecting and identifying targets in water, the propagation characteristics of seismic waves in shallow sea were deeply studied. Wave theory and staggered grid finite difference method were employed to analyze the propagation characteristics of seismic waves in a layer of underwater medium with a sedimentary layer. The results show that a Scholte wave is prone to disperse in an underwater medium with a sedimentary layer. The sediment thickness and acoustic parameters have great influence on the propagation characteristics of seismic waves. Moreover, the seismic wave signals were mainly concentrated below 20 Hz by collecting the vibration velocity of underwater interface. Furthermore, Scholte wave signals in underwater experiments were determined by utilizing four filtering methods and Scholte polarization characteristics, and the Butterworth filtering method can better preserve the integrity of vibration velocity signal. The theoretical, numerical, and experimental polarization characteristic results are basically consistent, which indicate that Scholte wave can be well extracted by the method of polarization characteristics.

Cold and Wet Mantle Transition Zone Beneath the Mediterranean Sea Revealed by the Electrical Image

A three-dimensional electrical conductivity model of the mantle beneath North Africa and Southern Europe is obtained by geomagnetic depth sounding. C-responses are estimated from geomagnetic data observed at observatories in and around the region and converted to the electrical structure of the mantle transition zone. The limited-memory quasi-Newton method is chosen to minimize the nonlinear objective function of inversion, while the forward modeling relies on a staggered-grid finite difference method in the spherical coordinate system. The data misfit term of the inversion function is measured using the L1-norm in order to suppress the response instability caused by the significant noise contained in the observed data. In order to adequately correct the ocean effect in observatories near the coast, a fixed shell comprised of ocean and land is incorporated in inversion. A banded zone with reduced conductivity is present in the three-dimensional model, primarily seen in the lower mantle transition zone and lower mantle beneath the Mediterranean Sea. Combining laboratory-measured conductivity models, we propose that subducted slabs causing reduced temperature and a water reservoir in the mantle transition zone should be responsible for the observed electrical model.

Dispersive staggered grid finite difference modelling of Rayleigh waves

Non-destructive surface wave-based ground exploration techniques such as multichannel analysis of surface wave (MASW) require a rigorous numerical model to predict the soil properties accurately. The staggered grid finite difference (SGFD) method is widely used to simulate elastic wave propagation due to its stability and ease of implementation. However, in the presence of low-velocity layers, the SGFD method requires computationally expensive high spatial grid resolution to accurately model the surface waves. On the other hand, velocity dependent grid approaches become unstable for the highly irregular soil profile. In order to overcome these challenges, a vertically non-uniform dispersive staggered grid scheme is proposed that significantly reduces the computing time without sacrificing accuracy. The proposed non-uniform discretization strategy is based on the Rayleigh wave dispersion relation, where the penetration depth of each frequency is assumed to be roughly one and half of the wavelength. Based on the above assumption, an exponential relationship between grid size and grid number is developed. The proposed method calculates continuously increasing grid length while ensuring grid points at the exact layer interface. Five different soil profiles: a homogenous elastic half-space, an irregularly dispersive crustal level profile, a low velocity layer at a deeper depth model, a real-world model, and a lateral heterogeneous model are used to demonstrate the efficacy of the proposed approach. For all the cases, the results of the proposed dispersive grid match those of the uniform grid with a correlation coefficient value (R) > 0.99. Compared to the uniform and velocity dependent grid approaches, the proposed dispersive grid scheme requires ½ to 1/4th of elements and CPU time. Furthermore, unlike the velocity dependent grid scheme, the proposed method does not become unstable for a highly irregular profile. Results indicate that the proposed method is superior to the traditional uniform and velocity dependent grid approaches for modelling elastic waves in layered media.

Wavelet-Based Three-Dimensional Inversion for Geomagnetic Depth Sounding

The complexity of Earth’s structure poses a challenge to the multiscale detection capability of geophysics. In this paper, we present a new wavelet-based three-dimensional inversion method for geomagnetic depth sounding. This method is based on wavelet functions to transfer model parameters in the space domain into the wavelet domain. The model is represented by wavelet coefficients containing both large- and fine-scale information, enabling wavelet-based inversion to describe multiscale anomalies. L1-norm measurement is applied to measure the model roughness to accomplish the sparsity constraint in the wavelet domain. Meanwhile, a staggered-grid finite difference method in a spherical coordinate system is used to calculate the forward responses, and the limited-memory quasi-Newton method is applied to seek the solution of the inversion objective function. Inversion tests of synthetic data for multiscale models show that wavelet-based inversion is stable and has multiresolution. Although higher-order wavelets can lead to finer results, our tests present that a db6 wavelet is suitable for geomagnetic depth sounding inversion. The db6 inversion results of responses at 129 geomagnetic observatories around the world reveal a higher-resolution image of the mantle.

Analysis of Rayleigh wave response characteristics in complicated cave system and detection methods

The detection of caves in strata is of great importance in engineering geology, geophysics, and other fields. Accurate identification of caves by elastic waves requires a clear understanding of the effects of each important parameter of caves on a wave field. Given the fact that previous Rayleigh wave prediction studies confined themselves to individual caves in a homogeneous half-space, the present study employed high-precision finite difference numerical modeling to analyze the Rayleigh wave response characteristics in strata containing complex caves. The staggered-grid finite difference method was used to simulate the wave field response of complex cave strata; then, the τ-p transform was employed for effective wave field separation, while the S-transform time-frequency method was used to analyze the properties of caves through the scattered Rayleigh waves. We mainly analyzed the pattern of wave field variation in the horizontal and vertical directions in the presence of multiple caves. With respect to a multi-layer model, the effect of the overlying layers on the cave-Rayleigh wave interaction was discussed. Stochastic cave cluster models were designed and the scattering recognition characteristics of cave clusters on the wave field were summarized for the actual geological conditions. Ultimately, the basis for wave field judgment for cave filling materials was presented. The present study provided an effective guide to solving complex cave problems using the surface wave method.

Enhancement of In-Plane Seismic Full Waveform Inversion with CPU and GPU Parallelization

Full waveform inversion is a widely used technique to estimate the subsurface parameters with the help of seismic measurements on the surface. Due to the amount of data, model size and non-linear iterative procedures, the numerical computation of Full Waveform Inversion are computationally intensive and time-consuming. This paper addresses the parallel computation of seismic full waveform inversion with Graphical Processing Units. Seismic full-waveform inversion of in-plane wave propagation in the finite difference method is presented here. The stress velocity formulation of the wave equation in the time domain is used. A four nodded staggered grid finite-difference method is applied to solve the equation, and the perfectly matched layers are considered to satisfy Sommerfeld’s radiation condition at infinity. The gradient descent method with conjugate gradient method is used for adjoined modelling in full-waveform inversion. The host code is written in C++, and parallel computation codes are written in CUDA C. The computational time and performance gained from CUDA C and OpenMP parallel computation in different hardware are compared to the serial code. The performance improvement is enhanced with increased model dimensions and remains almost constant after a certain threshold. A GPU performance gain of up to 90 times is obtained compared to the serial code.

Can we obtain the internal information of a surface crack from Rayleigh waves?

In this study, the Rayleigh wave responses of surface cracks with irregular planes are investigated utilizing novel numerical approaches. The effects of the equivalent width, filling situation and geometrical properties of the crack on the wave propagation are analyzed based on the scattered surface wave signals. The scattering process of the source due to an irregular surface crack is implemented by the staggered grid finite difference method (SG-FDM). The complex boundary conditions on the crack surface are classified by a modified stress image method (SIM) with a series of cells. A well comparison result is presented between SG-FDM and the finite element method. Three typical models are designed to analyze the effects of the crack width and filling material on Rayleigh waves. Furthermore, we focus on the influences of the effective width spatial variation characteristics of an irregular crack on Rayleigh wave propagation with different crack depth-wavelength ratios. The energy distribution of the scattered Rayleigh waves depends on the crack average width and the integrity of the filling material after damage repair. Our study reveals the specific wave phenomena (e.g., crack width-dependent and filling condition-dependent scattered waves) influenced by the crack interior structures, and this provides a method for detecting the subsurface structures (e.g., variation of crack width along the depth direction) of a breaking crack and evaluating the repair quality regarding local defects.

3D analysis of Rayleigh-type surface wave in fractured media: Insights into anisotropic characteristics

The propagation characteristics of Rayleigh waves on the surfaces of fracture-containing media have considerable theoretical and practical significance. Limitations are associated with using two-dimensional (2D) observations of Rayleigh waves to predict fracture features, because real fractures have a complex distribution inside media. Therefore, the three-dimensional (3D) propagation of Rayleigh waves in media with directional fractures is simulated in this study by combining the staggered-grid finite difference method with parallel computing. The wave velocity and amplitude of Rayleigh waves depend on the relative direction between fracture and wave path on the medium surface. That is, Rayleigh waves exhibit clearly anisotropic characteristics for a medium with directionally distributed fractures, where the degree of anisotropy is closely related to fracture parameters, such as the fracture density and fracture scale. Two new anisotropy parameters for Rayleigh waves are defined in this study: the anisotropy of the wave velocity and the anisotropy of attenuation (transmission coefficient). For a single fracture, the minor axis of the anisotropic quasi-semiellipse and the locations of discontinuities can indicate the normal direction of the fracture and the fracture distribution range, respectively. For multiple fractures, the major axis of the anisotropic quasi-semiellipse of the Rayleigh wave has almost the same direction as the main distribution of the fracture zone. This study provides theoretical support for the accurate prediction of complex fracture parameters inside media.

A hybrid explicit implicit staggered grid finite-difference scheme for the first-order acoustic wave equation modeling

Implicit staggered-grid finite-difference (SGFD) methods are widely used for the first-order acoustic wave-equation modeling. The identical implicit SGFD operator is commonly used for all of the first-order spatial derivatives in the first-order acoustic wave-equation. In this paper, we propose a hybrid explicit implicit SGFD (HEI-SGFD) scheme which could simultaneously preserve the wave-equation simulation accuracy and increase the wave-equation simulation speed. We use a second-order explicit SGFD operator for half of the first-order spatial derivatives in the first-order acoustic wave-equation. At the same time, we use the implicit SGFD operator with added points in the diagonal direction for the other first-order spatial derivatives in the first-order acoustic wave-equation. The proposed HEI-SGFD scheme nearly doubles the wave-equation simulation speed compared to the implicit SGFD schemes. In essence, the proposed HEI-SGFD scheme is equivalent to the second-order FD scheme with ordinary grid format. We then determine the HEI-SGFD coefficients in the time–space domain by minimizing the phase velocity error using the least-squares method. Finally, the effectiveness of the proposed method is demonstrated by dispersion analysis and numerical simulations.

Poro-acoustoelasticity finite-difference simulation of elastic wave propagation in prestressed porous media

Insights into wave propagation in prestressed porous media are important in geophysical applications, such as monitoring changes in geo-pressure. This can be addressed by poro-acoustoelasticity theory, which extends the classical acoustoelasticity of solids to porous media. The relevant poro-acoustoelasticity equations can be derived from anisotropic poroelasticity equations by replacing the poroelastic stiffness matrix with an acoustoelastic stiffness matrix consisting of second-order and third-order elastic constants. The theory considers the poroelasticity equations to be nonlinear due to the cubic strain-energy function with linear strains under finite-magnitude prestresses. A rotated staggered-grid finite-difference method with an unsplit convolutional perfectly matched layer absorbing boundary is used to solve a first-order velocity-stress formulation of poro-acoustoelasticity equations for elastic wave propagation in prestressed porous media. Numerical solutions are partially verified by computing the velocities of fast P wave, slow P wave, and S wave as a function of hydrostatic prestress and are compared with the exact values. Numerical simulations of wave propagation are carried out for the model of poro-acoustoelastic homogeneous space under three states—prestress confining (hydrostatic), uniaxial, and pure shear—and for the model of two poro-acoustoelastic homogeneous half-spaces in the planar contact under confining (hydrostatic) prestress. The resulting wavefield snapshots show fast P-wave, slow P-wave, and S-wave propagations in poro-acoustoelastic media under loading prestresses, which illustrate that the stress-induced velocity anisotropy is of orthotropy strongly related to the orientation of prestresses. These examples demonstrate the significant impact of prestressing conditions on seismic responses in velocity and anisotropy.

A Study on Integrity Testing for Platform-Pile Systems considering Wave Propagation

Platform-pile systems are typically used in bridge foundations and seaport structures. These systems are difficult to analyze using a conventional low-strain integrity test (LST) due to periodic oscillations. A staggered-grid finite-difference (SFD) method is proposed for a platform-pile system in soil and verified with ABAQUS and measured data. Different impact locations on the top surface of the platform and the lateral surface of the pile are analyzed, and parallel velocities are obtained at different locations. Choosing an appropriate impact location reduces the influence of the periodic oscillation signal; the impact locations at the center of symmetry are the best choices. At the same impact energy, the impact location should not be outside the cross section of the pile when the impact occurs on the top surface of the platform to achieve a meaningful damage signal. When the impact locations are on the lateral surface of the pile, the interferences are relatively small, and the damage vibration signal can be detected. The material and size of the hammer determine the impact duration and should be chosen to reduce the high-frequency influences and obtain reasonable data for analysis.

Study on propagation characteristics of Scholte wave in shallow sea at non-horizontal underwater interface

Theoretical and simulation studies are carried out to study the propagation characteristics of Scholte waves in non-horizontal underwater interface of shallow sea. The dispersion and polarization characteristics of Scholte waves in horizontal layered underwater interface are analyzed by using wave theory and staggered grid finite difference method. The effectiveness of the method used in this paper is proved from the perspective of theory and simulation. The propagation characteristics and polarization characteristics of Scholte waves in non-horizontal underwater interface are also analyzed. At the wedge-shaped underwater interface, the particle trajectory of Scholte wave is related to the angle between the propagation direction and the underwater interface, and with the increase of the angle, the deflection angle of the particle trajectory increases. At the bulged underwater interface, the particle trajectory of seismic wave has nothing to do with the propagation direction and the bulged underwater interface, but the bulged underwater interface has a great influence on the amplitude of seismic wave. The research results are of great significance to underwater target detection and geological exploration.

Three-dimensional magnetotelluric inversion considering electrical anisotropy with synthetic and real data

Electrical anisotropy increases the non-uniqueness of magnetotelluric (MT) inversions. However, it is still useful to develop MT anisotropic inversion to provide new knowledge for this field. In this study, we developed a three-dimensional (3-D) MT inversion code considering electrical anisotropy. The forward procedure is based on a staggered-grid finite-difference method and solves the governing equation for the electric fields with a direct method. The values of a traditional objective function and its gradient could then be calculated. The inversion procedure was completed with a well-developed subroutine of non-linear conjugate gradient method. Synthetic tests for different isotropic models indicate the validity of this inversion code. Furthermore, we performed 3-D anisotropic inversions for synthetic and real data. The results indicate that vertical resistivity is still indistinguishable in the anisotropic inversion results. In our synthetic tests, it becomes more difficult to recover low-resistivity components of anisotropic anomalies as the depth increases, and the vertical sizes of anisotropic anomalies are reduced, as characterized by a much shallower bottom boundary. The anisotropic strike angle is not as stable as the principal resistivities in the inversion process. Under the same changes in model parameters, the strike angle may cause larger deviations of the responses than the principal horizontal resistivities. With this code, we revisited the field data collected at the northern margin of the Tibetan Plateau. The inversion results indicate a possible anisotropic anomaly exists below the boundary between the plateau and its northern neighbor. This anisotropic anomaly could represent an area of concentration of shear deformation.

An optimized staggered-grid finite-difference operator for seismic wave simulation in poroelastic media

Most of the finite-difference methods for simulating seismic wave propagation in poroelastic media compute the spatial derivatives using fixed-length, which is the 2Mth-order operator. An optimized staggered-grid finite-difference (SGFD) method, which uses mixed operator length to approximate spatial derivatives, is derived from Biot’s equations to improve accuracy and stability. The numerical stability and dispersion of our optimized scheme are derived using the dispersion relation for the poroelastic equations. We have compared the dispersion curve plots for the developed and the conventional method. We introduce evidence that the developed scheme is comparatively more stable and accurate than the existing conventional method. The numerical solutions obtained using our optimized scheme closely approximate the analytical solution, and all wave types predicted by Biot’s theory are generated. Our scheme is used to simulate seismic wave propagation in a 2D homogeneous model, a two-layered model, a moderately complex model derived from the Utsira Formation, Sleipner North Sea, where a large amount of CO2 is being sequestrated, and a 3D homogeneous model. The developed scheme takes less computational time than the conventional SGFD method without compromising the accuracy or the stability.

A novel equivalent staggered-grid finite-difference scheme and its optimization strategy for variable-density acoustic wave modelling

Compared with the standard staggered-grid finite-difference (FD) methods, equivalent staggered-grid (ESG) ones can significantly reduce the computational memory for acoustic wave modelling in the variable-density media. To further enhance the simulation efficiency and accuracy, one way is to optimize the FD coefficients, another way is to design new FD stencils. In this paper, we propose a modified ESG (M-ESG) scheme which can significantly accelerate the wavefield simulation process while preserving or even improving the modelling accuracy. We calculate the FD coefficients by approximating the temporal and spatial derivatives simultaneously based on time–space domain (TS-D) dispersion relation of the discrete wave equation. Our M-ESG scheme in the TS-D can maintain basically the same accuracy as the conventional ESG (C-ESG) one when the FD coefficients are derived by the Taylor-series expansion (TE) approach. Note that the TS-D dispersion relation is nonlinear with respect to the FD coefficients of the C-ESG scheme, so it is difficult to obtain the optimized FD coefficients for the discrete wave equation. However, we can minimize the L2-norm error of the dispersion relation based on our M-ESG scheme to implement a linear FD coefficients optimization strategy, which is easy and efficient. Comparisons with TE- and optimization-based C-ESG schemes demonstrate the accuracy, stability, and efficiency superiorities of our TE- and optimization-based M-ESG ones.

A coupled geomechanics and reservoir simulator with a staggered grid finite difference method

For a stress-sensitive reservoir, the constant rock compressibility term used in a conventional reservoir simulator (CRS), does not account for change of porosity and permeability. This paper develops a coupled geomechanics and reservoir simulator (CGRS) which accounts for changes in porosity and permeability related to deformation.

Time–space domain scalar wave modeling by a novel hybrid staggered-grid finite-difference method with high temporal and spatial accuracies

Staggered-grid finite-difference (SFD) scheme is favored in wave equation simulation due to its superior accuracy and stability to center-grid finite-difference (CFD) scheme. However, for the scalar wave equation (SWE) modeling, the conventional SFD (CSFD) scheme only reaches second-order accuracy in space and time for the discrete SWE; the recently developed modified SFD (MSFD) scheme improves the accuracy to (2N)th-order (N<4) but is costly because the MSFD scheme is coined by adding many extra grid points to the CSFD scheme. To tackle these issues, we develop a cost-effective hybrid SFD (HSFD) scheme, which combines the features of the CSFD and MSFD schemes; we prove that the new HSFD scheme can simultaneously reach (2N)th-order accuracy in space and time for the discrete wave equation. In addition, to deal with the optimization difficulties due to the nonlinear dispersion relation of the SFD schemes, we propose a two-step linear optimization method to improve the accuracy of the new HSFD scheme. The analyses on dispersion, stability properties and numerical simulation examples demonstrate that the new HSFD scheme owns better accuracy and stability than the CSFD scheme. Computational cost analysis shows that the HSFD scheme can be more efficient than the MSFD scheme because it only requires half the additional grid points.

Optimal staggered-grid finite-difference method for wave modeling based on artificial neural networks

High-order staggered-grid finite-difference methods have been widely used in wave equation modeling. Imprecise difference coefficients and inappropriate difference order may result in numerical dispersion. We extend least square method (LSM) objective function to a flexible wave band to obtain its difference coefficients. Our method (model-driven LSM) is based on flexible model parameters to obtain corresponding optimal difference coefficients, rather than a uniform set of coefficients applied to all models. To reduce the computation cost in parameter estimation, we introduce the architecture, training and application of artificial neural networks (ANNs) to evaluate the staggered-grid finite-difference order and coefficients efficiently. In ANNs, double loss functions are adopted to ensure the stable and efficient convergence of training. Several numerical experiments for both two-dimensional and three-dimensional wave modeling are provided to verify the efficiency of the proposed scheme.

The Characteristics of Seismic Rotations in VTI Medium

Under the assumptions of linear elasticity and small deformation in traditional elastodynamics, the anisotropy of the medium has a significant effect on rotations observed during earthquakes. Based on the basic theory of the first-order velocity-stress elastic wave equation, this paper simulates the seismic wave propagation of the translational and rotational motions in two-dimensional isotropic and VTI (transverse isotropic media with a vertical axis of symmetry) media under different source mechanisms with the staggered-grid finite-difference method with respect to nine different seismological models. Through comparing the similarities and differences between the translational and rotational components of the wave fields, this paper focuses on the influence of anisotropic parameters on the amplitude and phase characteristics of the rotations. We verify that the energy of S waves in the rotational components is significantly stronger than that of P waves, and the response of rotations to the anisotropic parameters is more sensitive. There is more abundant information in the high-frequency band of the rotational components. With the increase of Thomsen anisotropic parameters ε and δ, the energy of the rotations increases gradually, which means that the rotational component observation may be helpful to the study of anisotropic parameters.