Elastic-wave-equation numerical simulation using a time–space domain dispersion-relation-based staggered-grid finite-difference method

Compared to the conventional high-order staggered-grid finite-difference method (C-SFD), the time–space domain dispersion-relation-based high-order staggered-grid finite-difference method (TS-SFD) can suppress the numerical dispersion more efficiently to achieve higher modeling accuracy when applied to numerically solving the velocity–stress acoustic equation. The enhanced modeling accuracy of TS-SFD is attributed to the difference-coefficient calculation approach based on the time–space domain dispersion relation, which makes the difference coefficients change adaptively with the propagation velocity of seismic waves in the medium. However, when numerical simulation of the velocity–stress elastic wave equation is conducted with TS-SFD, if the difference coefficients are calculated based on the time–space domain P-wave dispersion relation, the modeling accuracy of P-wave is high, while that of S-wave is low and vice versa. To address the limitation of TS-SFD in ensuring high modeling accuracy for both P- and S-wave, in this article, we propose a novel strategy to conduct elastic wave numerical simulation with TS-SFD, in which the decoupled P- and S-wave elastic wave equations are numerically solved with TS-SFD, and the difference coefficients are calculated based on the time–space domain P- and S-wave dispersion relations, respectively. Numerical dispersion analysis and numerical simulation examples show that the simulation strategy proposed in this article can ensure high simulation accuracy for both P- and S-wave. Stability analysis shows that this simulation strategy can effectively improve the stability of elastic wave numerical simulation with TS-SFD. In addition, the simulation strategy has the additional advantage of automatically separating the P- and S-wave, which provides a solid foundation for the subsequent analysis of P- and S-wave propagation characteristics in elastic media and elastic wave reverse time migration.

1 Introduction

Numerical simulation of the elastic wave equation is an important technique to study the characteristics of P- and S-waves in complex media (Ren and Liu, 2015; Chen et al., 2017), and it is an important research basis for elastic wave reverse time migration (Du et al., 2014; Su and Trad, 2018) and full-waveform inversion (Virieux and Operto, 2009; Singh et al., 2018). The finite-difference (FD) method, with its two important advantages of high computational efficiency and minimal memory requirements, has evolved into one of the most widely used numerical simulation methods for solving the elastic wave equation (Ren and Li, 2017; Cao and Chen, 2018; Wang et al., 2025). However, FD is affected by the inherent numerical dispersion, resulting in low simulation accuracy of both P- and S-waves. Therefore, suppressing the numerical dispersion to improve the simulation accuracy is an important issue in the finite-difference numerical simulation of elastic waves (Ren and Liu, 2015; Zhou et al., 2022). In addition, when analyzing the propagation characteristics of P- and S-waves in elastic media and conducting elastic wave reverse time migration, it is necessary to separate the P- and S-waves (Wang et al., 2015; Zhu, 2017), so it is of great significance to decouple P- and S-wave in the process of numerical simulation of elastic waves.

The finite-difference method was first proposed by Alterman and Karal (1968) for numerically solving the second-order elastic wave displacement equation in order to conduct elastic wave numerical simulations. Virieux (1984) and Virieux (1986) introduced a staggered-grid finite-difference method, first used in the numerical simulation of the electromagnetic wave equation (Kane, 1966), to the first-order velocity–stress elastic wave equation, realizing two-dimensional elastic wave numerical simulation using a staggered-grid finite-difference scheme. The staggered-grid finite-difference method is naturally associated with the first-order velocity–stress elastic wave equation, which facilitates the implementation of both stress sources and displacement sources. Consequently, it has emerged as the most prevalent numerical method for elastic wave simulation (Graves, 1996; Wang et al., 2022). The staggered-grid finite-difference method described by Virieux has only second-order temporal accuracy and second-order spatial accuracy, which cannot effectively suppress temporal and spatial numerical dispersion, resulting in low modeling accuracy. Levander (1988) developed a staggered-grid finite-difference scheme with second-order temporal accuracy and fourth-order spatial accuracy, significantly improving the modeling accuracy. Dong et al. (2000) further developed a staggered-grid finite-difference method with second-order temporal accuracy and 2Mth-order spatial accuracy and effectively improved the numerical simulation accuracy of elastic waves in VTI media. The temporal high-order accuracy has a significant impact on the amount of memory occupation and computation required. Consequently, scholars commonly adopt the staggered-grid finite-difference method with second-order temporal accuracy and 2$M$th-order spatial accuracy, which is referred to as the conventional high-order staggered-grid finite-difference method (C-SFD), for the numerical simulation of elastic waves. Liu and Sen (2011) indicated that the numerical simulation of the velocity–stress acoustic wave equation was conducted in both the spatial and temporal domains. Nevertheless, the difference coefficients of C-SFD were calculated based on only the dispersion relation in the spatial domain, which is unreasonable. Consequently, they proposed a time–space domain staggered-grid finite-difference method (TS-SFD), with the difference coefficients calculated based on the time–space domain dispersion relation. Drawing upon the difference-coefficient calculation method of TS-SFD, Tan and Huang (2014) proposed a novel TS-SFD based on a stencil that contains a few more grid points than the standard stencil, and Zhou et al. (2024) further developed a new multi-axial staggered-grid finite-difference method that further enhanced the simulation accuracy of the acoustic wave equation. Compared with C-SFD, TS-SFD can suppress the numerical dispersion more effectively, thereby enhancing the accuracy of the simulation. The primary feature of TS-SFD is that the difference coefficients change adaptively with the propagation velocity of the seismic wave in the medium. However, when TS-SFD is directly used in numerical simulation of the velocity–stress elastic wave equation—considering the existence of P- and S-wave velocity—if the difference coefficients change adaptively with the P-wave velocity, this would result in a high degree of P-wave simulation accuracy, accompanied by a lower degree of S-wave simulation. Conversely, if the difference coefficients change adaptively with the S-wave velocity, this would result in a high degree of S-wave simulation accuracy, accompanied by a lower degree of P-wave simulation accuracy.

The separation of the P- and S-waves constitutes a pivotal element in the realm of elastic wave simulation and migration. Among the array of methodologies used to address this separation, Helmholtz decomposition stands as a prominent approach (Sun et al., 2004; Zhu, 2017). However, the separated P-wave scalar wavefields and S-wave vector wavefields that are yielded from the divergence and curl operators in Helmholtz decomposition are subject to amplitude and phase distortion, necessitating further calibration (Sun et al., 2004; Du et al., 2012). Ma and Zhu (2003) derived the decoupled P- and S-wave elastic wave equation from the standard elastic wave equation, and they demonstrated that numerical simulation using this decoupled elastic wave equation can automatically achieve the separation of P- and S-waves. Wu et al. (2018) conducted a systematic comparison of elastic wavefield separation methods. It was highlighted that elastic wavefield separation based on the decoupled P- and S-wave elastic wave equations has a relatively high level of accuracy, with local energy anomalies only occurring at the interface, which can be eliminated by smoothing the velocity model. In this article, we attempt to introduce TS-SFD to the numerical simulation of the decoupled P- and S-wave elastic wave equation. The decoupled P- and S-wave equations are solved numerically with the difference coefficients calculated based on the time–space domain P- and S-wave dispersion relations, respectively. This ensures that both P- and S-waves are simulated with high accuracy. The structure of the present paper is as follows. First, the velocity–stress elastic wave equation and its difference discretization form are presented, and the difference-coefficient calculation approach is derived based on the time–space domain P- and S-wave dispersion relation. Second, an implementation scheme for solving the decoupled P- and S-wave elastic wave equations with TS-SFD is proposed. Third, dispersion analysis and stability analysis are conducted. Finally, two numerical simulation examples are performed on a layer model and the tectonic-reservoir model from the Sichuan Basin, respectively. The numerical simulation examples demonstrate that solving the decoupled P- and S-wave elastic wave equation with TS-SFD can enhance the simulation accuracy of both P- and S-waves, and P- and S-waves can be automatically separated.

2 Basic theory of TS-SFD

2.1 Finite-difference discrete elastic wave equation given by TS-SFD

The 2D velocity–stress elastic wave equation can be expressed as

$\begin{array}{c}\frac{\partial {\upsilon }_x}{\partial t}=\frac{1}{\rho }\left(\frac{\partial {\tau }_{xx}}{\partial x}+\frac{\partial {\tau }_{xz}}{\partial z}\right),\frac{\partial {\upsilon }_z}{\partial t}=\frac{1}{\rho }\left(\frac{\partial {\tau }_{xz}}{\partial x}+\frac{\partial {\tau }_{zz}}{\partial z}\right),\\ \frac{\partial {\tau }_{xx}}{\partial t}=\left(\lambda +2\mu \right)\frac{\partial {\upsilon }_x}{\partial x}+\lambda \frac{\partial {\upsilon }_z}{\partial z},\frac{\partial {\tau }_{zz}}{\partial t}=\lambda \frac{\partial {\upsilon }_x}{\partial x}+\left(\lambda +2\mu \right)\frac{\partial {\upsilon }_z}{\partial z},\\ \frac{\partial {\tau }_{xz}}{\partial t}=\mu \left(\frac{\partial {\upsilon }_x}{\partial z}+\frac{\partial {\upsilon }_z}{\partial x}\right),\end{array}(1)$

where ${\upsilon }_x={\upsilon }_x(x,z,t)$ and ${\upsilon }_z={\upsilon }_z(x,z,t)$ are the velocity fields of particle vibration; ${\tau }_{xx}={\tau }_{xx}(x,z,t)$, ${\tau }_{zz}={\tau }_{zz}(x,z,t)$, and ${\tau }_{xz}={\tau }_{xz}(x,z,t)$ are the stress fields; $\lambda =\lambda (x,z)$ and $\mu =\mu (x,z)$ are the Lamé constants; and $\rho =\rho (x,z)$ is the density.

When the staggered-grid finite-difference method is used to solve the elastic wave equation, the wavefield variables $\left({\tau }_{xx},{\tau }_{zz},{\tau }_{xz},{\upsilon }_x,{\upsilon }_z\right)$ and the elastic parameters $\left(\lambda ,\mu ,\rho \right)$ are defined on the staggered-grid locations, as illustrated in Figure 1. ${\upsilon }_z$ is defined at spatial half-grid and temporal integer grid, and its expression is ${\upsilon }_{z(m-1/2,n-1/2)}^j={\upsilon }_z\left[x+(m-1/2)h,z+(n-1/2)h,t+j\mathrm{\Delta }t\right]$, where $m$, $n$, and $j$ are integers; $h$ is the spatial sampling interval; and $\mathrm{\Delta }t$ is the temporal sampling interval. ${\upsilon }_{x(m,n)}^j$, ${\tau }_{xx(m-1/2,n)}^{j-1/2}$, ${\tau }_{zz(m-1/2,n)}^{j-1/2}$, and ${\tau }_{xz(m,n-1/2)}^{j-1/2}$ have similar expressions.

Figure 1

Figure 1. Schematic diagram of the relative positions of the wavefield variables and elastic parameters in a 2D elastic wave staggered-grid finite-difference method.

Both C-SFD and TS-SFD use a second-order difference operator to approximate the partial derivatives of the wavefield variables with respect to $t$, and $\partial {\upsilon }_x/\partial t$, $\partial {\upsilon }_z/\partial t$, $\partial {\tau }_{xx}/\partial t$, $\partial {\tau }_{zz}/\partial t$, and $\partial {\tau }_{xz}/\partial t$ can be approximated as

$\begin{array}{c}{\left.\frac{\partial {\upsilon }_x}{\partial t}\right|}_{\left(\mathrm{0,0}\right)}^{\frac{1}{2}}\approx \frac{{\upsilon }_{x\left(\mathrm{0,0}\right)}^1-{\upsilon }_{x\left(\mathrm{0,0}\right)}^0}{\mathrm{\Delta }t},{\left.\frac{\partial {\upsilon }_z}{\partial t}\right|}_{\left(\frac{1}{2},\frac{1}{2}\right)}^{\frac{1}{2}}\approx \frac{{\upsilon }_{z\left(\frac{1}{2},\frac{1}{2}\right)}^1-{\upsilon }_{z\left(\frac{1}{2},\frac{1}{2}\right)}^0}{\mathrm{\Delta }t},\\ {\left.\frac{\partial {\tau }_{xx}}{\partial t}\right|}_{\left(\frac{1}{2},0\right)}^0\approx \frac{{\tau }_{xx\left(\frac{1}{2},0\right)}^{\frac{1}{2}}-{\tau }_{xx\left(\frac{1}{2},0\right)}^{-\frac{1}{2}}}{\mathrm{\Delta }t},{\left.\frac{\partial {\tau }_{zz}}{\partial t}\right|}_{\left(\frac{1}{2},0\right)}^0\approx \frac{{\tau }_{zz\left(\frac{1}{2},0\right)}^{\frac{1}{2}}-{\tau }_{zz\left(\frac{1}{2},0\right)}^{-\frac{1}{2}}}{\mathrm{\Delta }t},\\ {\left.\frac{\partial {\tau }_{xz}}{\partial t}\right|}_{\left(0,\frac{1}{2}\right)}^0\approx \frac{{\tau }_{xz\left(0,\frac{1}{2}\right)}^{\frac{1}{2}}-{\tau }_{xz\left(0,\frac{1}{2}\right)}^{-\frac{1}{2}}}{\mathrm{\Delta }t},\end{array}(2)$

where ${\upsilon }_{z(1/\mathrm{2,1}/2)}^0$ stands for ${\upsilon }_{z(m-1/2,n-1/2)}^j$ with $m=1$, $n=1$, and $j=1$.

Both C-SFD and TS-SFD use a 2Mth-order difference operator to approximate the partial derivatives of the wavefield variables with respect to $x$ and $z$, and $\partial {\tau }_{xx}/\partial x$, $\partial {\tau }_{xz}/\partial z$, $\partial {\tau }_{xz}/\partial x$, $\partial {\tau }_{zz}/\partial z$, $\partial {\upsilon }_x/\partial x$, $\partial {\upsilon }_z/\partial z$, $\partial {\upsilon }_x/\partial z$, and $\partial {\upsilon }_z/\partial x$ can be approximated as

$\begin{array}{c}{\left.\frac{\partial {\tau }_{xx}}{\partial x}\right|}_{\left(\mathrm{0,0}\right)}^{\frac{1}{2}}\approx \frac{1}{h}{\displaystyle \sum _{m=1}^M}{a}_m\left[{\tau }_{xx\left(m-\frac{1}{2},0\right)}^{\frac{1}{2}}-{\tau }_{xx\left(-m+\frac{1}{2},0\right)}^{\frac{1}{2}}\right],{\left.\frac{\partial {\tau }_{xz}}{\partial z}\right|}_{\left(\mathrm{0,0}\right)}^{\frac{1}{2}}\approx \frac{1}{h}{\displaystyle \sum _{m=1}^M}{a}_m\left[{\tau }_{xz\left(0,m-\frac{1}{2}\right)}^{\frac{1}{2}}-{\tau }_{xz\left(0,-m+\frac{1}{2}\right)}^{\frac{1}{2}}\right],\\ {\left.\frac{\partial {\tau }_{xz}}{\partial x}\right|}_{\left(\frac{1}{2},\frac{1}{2}\right)}^{\frac{1}{2}}\approx \frac{1}{h}{\displaystyle \sum _{m=1}^M}{a}_m\left[{\tau }_{xz\left(m,\frac{1}{2}\right)}^{\frac{1}{2}}-{\tau }_{xz\left(-m+1,\frac{1}{2}\right)}^{\frac{1}{2}}\right],{\left.\frac{\partial {\tau }_{zz}}{\partial z}\right|}_{\left(\frac{1}{2},\frac{1}{2}\right)}^{\frac{1}{2}}\approx \frac{1}{h}{\displaystyle \sum _{m=1}^M}{a}_m\left[{\tau }_{zz\left(\frac{1}{2},m\right)}^{\frac{1}{2}}-{\tau }_{zz\left(\frac{1}{2},-m+1\right)}^{\frac{1}{2}}\right],\\ {\left.\frac{\partial {\upsilon }_x}{\partial x}\right|}_{\left(\frac{1}{2},0\right)}^0\approx \frac{1}{h}{\displaystyle \sum _{m=1}^M}{a}_m\left[{\upsilon }_{x\left(m,0\right)}^0-{\upsilon }_{x\left(-m+\mathrm{1,0}\right)}^0\right],{\left.\frac{\partial {\upsilon }_z}{\partial z}\right|}_{\left(\frac{1}{2},0\right)}^0\approx \frac{1}{h}{\displaystyle \sum _{m=1}^M}{a}_m\left[{\upsilon }_{z\left(\frac{1}{2},m-\frac{1}{2}\right)}^0-{\upsilon }_{z\left(\frac{1}{2},-m+\frac{1}{2}\right)}^0\right],\\ {\left.\frac{\partial {\upsilon }_x}{\partial z}\right|}_{\left(0,\frac{1}{2}\right)}^0\approx \frac{1}{h}{\displaystyle \sum _{m=1}^M}{a}_m\left[{\upsilon }_{x\left(0,m\right)}^0-{\upsilon }_{x\left(0,-m+1\right)}^0\right],{\left.\frac{\partial {\upsilon }_z}{\partial x}\right|}_{\left(0,\frac{1}{2}\right)}^0\approx \frac{1}{h}{\displaystyle \sum _{m=1}^M}{a}_m\left[{\upsilon }_{z\left(m-\frac{1}{2},\frac{1}{2}\right)}^0-{\upsilon }_{z\left(-m+\frac{1}{2},\frac{1}{2}\right)}^0\right],\end{array}(3)$

where $a_1,a_2,\dots ,a_M$ are the difference coefficients.

Substituting Equations 2, 3 into Equation 1 yields

$\begin{array}{c}\frac{{\upsilon }_{x\left(\mathrm{0,0}\right)}^1-{\upsilon }_{x\left(\mathrm{0,0}\right)}^0}{\mathrm{\Delta }t}\approx \frac{1}{\rho h}{\displaystyle \sum _{m=1}^M}{a}_m\left[{\tau }_{xx\left(m-\frac{1}{2},0\right)}^{\frac{1}{2}}-{\tau }_{xx\left(-m+\frac{1}{2},0\right)}^{\frac{1}{2}}+{\tau }_{xz\left(0,m-\frac{1}{2}\right)}^{\frac{1}{2}}-{\tau }_{xz\left(0,-m+\frac{1}{2}\right)}^{\frac{1}{2}}\right],\\ \frac{{\upsilon }_{z\left(\frac{1}{2},\frac{1}{2}\right)}^1-{\upsilon }_{z\left(\frac{1}{2},\frac{1}{2}\right)}^0}{\mathrm{\Delta }t}\approx \frac{1}{\rho h}{\displaystyle \sum _{m=1}^M}{a}_m\left[{\tau }_{xz\left(m,\frac{1}{2}\right)}^{\frac{1}{2}}-{\tau }_{xz\left(-m+1,\frac{1}{2}\right)}^{\frac{1}{2}}+{\tau }_{zz\left(\frac{1}{2},m\right)}^{\frac{1}{2}}-{\tau }_{zz\left(\frac{1}{2},-m+1\right)}^{\frac{1}{2}}\right],\\ \frac{{\tau }_{xx\left(\frac{1}{2},0\right)}^{\frac{1}{2}}-{\tau }_{xx\left(\frac{1}{2},0\right)}^{-\frac{1}{2}}}{\mathrm{\Delta }t}\approx \frac{\lambda +2\mu }{h}{\displaystyle \sum _{m=1}^M}{a}_m\left[{\upsilon }_{x\left(m,0\right)}^0-{\upsilon }_{x\left(-m+\mathrm{1,0}\right)}^0\right]+\frac{\lambda }{h}{\displaystyle \sum _{m=1}^M}{a}_m\left[{\upsilon }_{z\left(\frac{1}{2},m-\frac{1}{2}\right)}^0-{\upsilon }_{z\left(\frac{1}{2},-m+\frac{1}{2}\right)}^0\right],\\ \frac{{\tau }_{zz\left(\frac{1}{2},0\right)}^{\frac{1}{2}}-{\tau }_{zz\left(\frac{1}{2},0\right)}^{-\frac{1}{2}}}{\mathrm{\Delta }t}\approx \frac{\lambda }{h}{\displaystyle \sum _{m=1}^M}{a}_m\left[{\upsilon }_{x\left(m,0\right)}^0-{\upsilon }_{x\left(-m+\mathrm{1,0}\right)}^0\right]+\frac{\lambda +2\mu }{h}{\displaystyle \sum _{m=1}^M}{a}_m\left[{\upsilon }_{z\left(\frac{1}{2},m-\frac{1}{2}\right)}^0-{\upsilon }_{z\left(\frac{1}{2},-m+\frac{1}{2}\right)}^0\right],\\ \frac{{\tau }_{xz\left(0,\frac{1}{2}\right)}^{\frac{1}{2}}-{\tau }_{xz\left(0,\frac{1}{2}\right)}^{-\frac{1}{2}}}{\mathrm{\Delta }t}\approx \frac{\mu }{h}{\displaystyle \sum _{m=1}^M}{a}_m\left[{\upsilon }_{x\left(0,m\right)}^0-{\upsilon }_{x\left(0,-m+1\right)}^0+{\upsilon }_{z\left(m-\frac{1}{2},\frac{1}{2}\right)}^0-{\upsilon }_{z\left(-m+\frac{1}{2},\frac{1}{2}\right)}^0\right].\end{array}(4)$

Equation 4 represents the finite-difference discrete elastic wave equation corresponding to Equation 1. Both C-SFD and TS-SFD realize numerical simulations of elastic waves by iterating Equation 4, and the fundamental difference between them lies in the method of calculating the difference coefficients $a_1,a_2,\dots ,a_M$. As C-SFD and TS-SFD use the same difference discrete elastic wave equation, their computational demands and memory requirements are essentially equivalent.

2.2 Difference-coefficient calculation

C-SFD calculates the difference coefficients only based on the space domain dispersion relation, with a detailed derivation of this algorithm provided by Peng et al. (2022). For the sake of brevity, the detailed derivation will not be reproduced here. Instead, the generic solution of the difference coefficients of C-SFD is presented as follows:

$a_m=\frac{{\left(-1\right)}^{M-1}}{2m-1}{\displaystyle \prod _{1\le k\le M,k\ne m}}\frac{{\left(2k-1\right)}^2}{{\left(2m-1\right)}^2-{\left(2k-1\right)}^2}\left(m=\mathrm{1,2},\dots ,M\right).(5)$

From Equation 5, it can be observed that the difference coefficients of C-SFD are only related to the difference order $M$ and not to the propagation velocity of the seismic wave in the medium.

The algorithm for calculating the difference coefficients of TS-SFD based on the time–space domain dispersion relation can be summarized as follows: first, the time–space domain P- and S-wave dispersion relations are derived based on plane-wave theory. Second, the triangular function in the dispersion relation is expanded by the Taylor series, and a system of equations on the difference coefficients is established. Finally, the difference coefficients can be obtained by solving this system of equations. The following part of this section provides a comprehensive description of the algorithm for calculating the difference coefficients of TS-SFD, as outlined in the preceding steps.

In a homogeneous medium, the velocity–stress elastic wave equation has the following discrete plane wave solution:

$\begin{array}{cc}\hfill {\upsilon }_{x\left(m,n\right)}^j& =A_{{\upsilon }_x}{\mathrm{e}}^{\text{i}\left[k_x\left(x + mh\right) + k_z\left(z + nh\right) - \omega \left(t + j\mathrm{\Delta }t\right)\right]},\hfill \\ \hfill {\upsilon }_{z\left(m-\frac{1}{2},n-\frac{1}{2}\right)}^j& =A_{{\upsilon }_z}{\mathrm{e}}^{\text{i}\left[k_x\left(x + \left(m - \frac{1}{2}\right)h\right) + k_z\left(z + \left(n - \frac{1}{2}\right)h\right) - \omega \left(t + j\mathrm{\Delta }t\right)\right]},\hfill \\ \hfill {\tau }_{xx\left(m-\frac{1}{2},n\right)}^{j-\frac{1}{2}}& =A_{{\tau }_{xx}}{\mathrm{e}}^{\text{i}\left[k_x\left(x + \left(m - \frac{1}{2}\right)h\right) + k_z\left(z + nh\right) - \omega \left(t + \left(j - \frac{1}{2}\right)\mathrm{\Delta }t\right)\right]},\hfill \\ \hfill {\tau }_{zz\left(m-\frac{1}{2},n\right)}^{j-\frac{1}{2}}& =A_{{\tau }_{zz}}{\mathrm{e}}^{\text{i}\left[k_x\left(x + \left(m - \frac{1}{2}\right)h\right) + k_z\left(z + nh\right) - \omega \left(t + \left(j - \frac{1}{2}\right)\mathrm{\Delta }t\right)\right]},\hfill \\ \hfill {\tau }_{xz\left(m,n-\frac{1}{2}\right)}^{j-\frac{1}{2}}& =A_{{\tau }_{xz}}{\mathrm{e}}^{\text{i}\left[k_x\left(x + mh\right) + k_z\left(z + \left(n - \frac{1}{2}\right)h\right) - \omega \left(t + \left(j - \frac{1}{2}\right)\mathrm{\Delta }t\right)\right]},\hfill \end{array}(6)$

where $A_{{\upsilon }_x}$, $A_{{\upsilon }_z}$, $A_{{\tau }_{xx}}$, $A_{{\tau }_{zz}}$, and $A_{{\tau }_{xz}}$ are the amplitudes of the plane wave; $k_x=k\cos \theta$, $k_z=k\sin \theta$, and $k$ are the wave numbers; $\theta$ is the angle between the propagation direction of the plane wave and the positive direction of the $x$ axis; $\omega$ is the angular frequency; and $\mathrm{i}$ the imaginary unit, i.e., ${\mathrm{i}}^2=-1$.

Substituting Equation 6 into Equation 4 yields

$\begin{array}{c}g{A}_{{\upsilon }_x}\approx \frac{A_{{\tau }_{xx}}}{\rho }{f}_x+\frac{A_{{\tau }_{xz}}}{\rho }{f}_z,g{A}_{{\upsilon }_z}\approx \frac{A_{{\tau }_{xz}}}{\rho }{f}_x+\frac{A_{{\tau }_{zz}}}{\rho }{f}_z,\\ g{A}_{{\tau }_{xx}}\approx \left(\lambda +2\mu \right)A_{{\upsilon }_x}{f}_x+\lambda A_{{\upsilon }_z}{f}_z,\\ g{A}_{{\tau }_{zz}}\approx \lambda A_{{\upsilon }_x}{f}_x+\left(\lambda +2\mu \right)A_{{\upsilon }_z}{f}_z,g{A}_{{\tau }_{xz}}\approx \mu A_{{\upsilon }_x}{f}_x+\mu A_{{\upsilon }_z}{f}_z,\end{array}(7)$

where

$\begin{array}{c}g=\frac{-\sin \left(\frac{\omega \mathrm{\Delta }t}{2}\right)}{\mathrm{\Delta }t},f_x=\frac{1}{h}{\displaystyle \sum _{m=1}^M}{a}_m\sin \left[\left(m-\frac{1}{2}\right)k_{x}h\right],\\ f_z=\frac{1}{h}{\displaystyle \sum _{m=1}^M}{a}_m\sin \left[\left(m-\frac{1}{2}\right)k_{z}h\right].\end{array}(8)$

Eliminating $A_{{\upsilon }_x}$, $A_{{\upsilon }_z}$, $A_{{\tau }_{xx}}$, $A_{{\tau }_{zz}}$, and $A_{{\tau }_{xz}}$ from Equation 7 yields Equation 9:

$g^2\left[g^2-\frac{\lambda +2\mu }{\rho }\left(f_x^2+f_z^2\right)\right]\left[g^2-\frac{\mu }{\rho }\left(f_x^2+f_z^2\right)\right]\approx 0,(9)$

where $g^2$ is not constantly zero, so we obtain

$\left[g^2-\frac{\lambda +2\mu }{\rho }\left(f_x^2+f_z^2\right)\right]\approx 0,\left[g^2-\frac{\mu }{\rho }\left(f_x^2+f_z^2\right)\right]\approx 0.(10)$

Substituting Equation 8 into Equation 10 and considering $v_{\text{p}}^2=(\lambda +2\mu )/\rho$ and $v_{\text{s}}^2=\mu /\rho$, we obtain

$\begin{array}{ll}\hfill \frac{1}{{\left(v_{\text{p}}\mathrm{\Delta }t\right)}^2}{\sin }^2\left(\frac{r_{\text{p}}kh}{2}\right)& \approx \frac{1}{h^2}{\left\{{\displaystyle \sum _{m=1}^M}{a}_m\sin \left[\left(m-\frac{1}{2}\right)k_{x}h\right]\right\}}^2\hfill \\ \hfill & +\frac{1}{h^2}{\left\{{\displaystyle \sum _{m=1}^M}{a}_m\sin \left[\left(m-\frac{1}{2}\right)k_{z}h\right]\right\}}^2,\hfill \end{array}(11)$

$\begin{array}{ll}\hfill \frac{1}{{\left(v_{\text{s}}\mathrm{\Delta }t\right)}^2}{\sin }^2\left(\frac{r_{\text{s}}kh}{2}\right)& \approx \frac{1}{h^2}{\left\{{\displaystyle \sum _{m=1}^M}{a}_m\sin \left[\left(m-\frac{1}{2}\right)k_{x}h\right]\right\}}^2\hfill \\ \hfill & +\frac{1}{h^2}{\left\{{\displaystyle \sum _{m=1}^M}{a}_m\sin \left[\left(m-\frac{1}{2}\right)k_{z}h\right]\right\}}^2,\hfill \end{array}(12)$

where $v_{\text{p}}$ and $v_{\text{s}}$ represent the propagation velocity of P- and S-wave in the medium, respectively, and $r_{\text{p}}=v_{\text{p}}\mathrm{\Delta }t/h$ and $r_{\text{s}}=v_{\text{s}}\mathrm{\Delta }t/h$ represent the Courant number of P- and S-wave, respectively. Equation 11 and Equation 12 represent the time–space domain P- and S-wave dispersion relation, respectively.

The principle of calculating the difference coefficients based on the time–space domain dispersion relation of the P- or S-wave is the same, so we will only explain the algorithm of calculating the difference coefficients based on the time–space domain dispersion relation of the P-wave in detail below. Performing a Taylor series expansion of the sine function in Equation 11, we obtain

$\begin{array}{ll}\hfill {\left[{\displaystyle \sum _{j=0}^{\infty }}{r}_{\text{p}}^{2j}{\beta }_j{\left(\frac{k}{2}\right)}^{2j+1}{h}^{2j}\right]}^2& \approx {\left[{\displaystyle \sum _{j=0}^{\infty }}{d}_j{\beta }_j{\left(\frac{k\cos \theta }{2}\right)}^{2j+1}{h}^{2j}\right]}^2\hfill \\ \hfill & +{\left[{\displaystyle \sum _{j=0}^{\infty }}{d}_j{\beta }_j{\left(\frac{k\sin \theta }{2}\right)}^{2j+1}{h}^{2j}\right]}^2\hfill \end{array}(13)$

where

${\beta }_j=\frac{{\left(-1\right)}^j}{\left(2j+1\right)!},d_j={\displaystyle \sum _{m=1}^M}{\left(2m-1\right)}^{2j+1}{a}_m.(14)$

Equalizing the coefficients of $k^{2j+2}{h}^{2j}(j=\mathrm{0,1,2},\dots ,M-1)$ on the left and right sides of Equation 13 yields

$d_0^2\left({\cos }^2\theta +{\sin }^2\theta \right)=1\left(j=0\right),(15)$

$\begin{array}{ll}\hfill & {\displaystyle \sum _{q=0}^j}{d}_q{d}_{j-q}{\beta }_q{\beta }_{j-q}\left({\cos }^{2j+2}\theta +{\sin }^{2j+2}\theta \right)\hfill \\ \hfill & ={\displaystyle \sum _{q=0}^j}{\beta }_q{\beta }_{j-q}{r}_{\text{p}}^{2j}\left(j=\mathrm{1,2},\dots ,M-1\right).\hfill \end{array}(16)$

According to Equation 15, we obtain $d_0=\pm 1$. When $d_0$ varies from 1 to $-1$, the difference coefficients $a_1,a_2,\dots ,a_M$ change to their opposite, which has no effect on the final result. Here, we take $d_0=1$, and according to Equation 16, we obtain

$d_j=\frac{{\displaystyle \sum _{q=0}^j}{\beta }_q{\beta }_{j-q}{r}_{\text{p}}^{2j}-{\displaystyle \sum _{q=1}^{j-1}}{d}_q{d}_{j-q}{\beta }_q{\beta }_{j-q}\left({\cos }^{2j+2}\theta +{\sin }^{2j+2}\theta \right)}{2d_0{\beta }_0{\beta }_j\left({\cos }^{2j+2}\theta +{\sin }^{2j+2}\theta \right)}\left(j=\mathrm{1,2},\dots ,M-1\right).(17)$

In Equation 17, ${\beta }_j$ is given by Equation 14, and by choosing a certain value for $\theta$, $d_1,d_2,\dots ,d_{M-1}$ can be calculated in turn. In this way, $d_0,d_1,d_2,\dots ,d_{M-1}$ are all worded out, and then according to Equation 14, we obtain the following matrix equation:

$\begin{array}{ll}\hfill & \left[\begin{array}{ccccc}\hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill \cdots \hfill & \hfill 1\hfill \\ \hfill 1^2\hfill & \hfill 3^2\hfill & \hfill 5^2\hfill & \hfill \cdots \hfill & \hfill {\left(2M-1\right)}^2\hfill \\ \hfill 1^4\hfill & \hfill 3^4\hfill & \hfill 5^4\hfill & \hfill \cdots \hfill & \hfill {\left(2M-1\right)}^4\hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill ⋮\hfill \\ \hfill 1^{2M-2}\hfill & \hfill 3^{2M-2}\hfill & \hfill 5^{2M-2}\hfill & \hfill \cdots \hfill & \hfill {\left(2M-1\right)}^{2M-2}\hfill \end{array}\right]\left[\begin{array}{c}\hfill 1a_1\hfill \\ \hfill 3a_2\hfill \\ \hfill 5a_3\hfill \\ \hfill ⋮\hfill \\ \hfill \left(2M-1\right)a_M\hfill \end{array}\right]\hfill \\ \hfill & =\left[\begin{array}{c}\hfill d_0\hfill \\ \hfill d_1\hfill \\ \hfill d_2\hfill \\ \hfill ⋮\hfill \\ \hfill d_{M-1}\hfill \end{array}\right].\hfill \end{array}(18)$

Equation 18 is a Vandermonde matrix equation. From the progress of calculating $d_1,d_2,\dots ,d_{M-1}$, it is clear that $d_j$ is related to the value of $\theta$, and then, the difference coefficients $a_1,a_2,\dots ,a_M$ are also related to the value of $\theta$. If $\theta$ takes an arbitrary value, a simple expression for $d_j$ may not exist. As a result, the generic solution of the difference coefficients $a_1,a_2,\dots ,a_M$ cannot be derived, and then, the difference coefficients $a_1,a_2,\dots ,a_M$ must be calculated by solving Equation 18.

Taking $\theta =0$, we can derive $d_j=r_{\text{p}}^{2j}(j=\mathrm{0,1,2},\dots ,M-1)$, and the matrix Equation 18 can be rewritten as

$\begin{array}{ll}\hfill & \left[\begin{array}{ccccc}\hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill \cdots \hfill & \hfill 1\hfill \\ \hfill 1^2\hfill & \hfill 3^2\hfill & \hfill 5^2\hfill & \hfill \cdots \hfill & \hfill {\left(2M-1\right)}^2\hfill \\ \hfill 1^4\hfill & \hfill 3^4\hfill & \hfill 5^4\hfill & \hfill \cdots \hfill & \hfill {\left(2M-1\right)}^4\hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill ⋮\hfill \\ \hfill 1^{2M-2}\hfill & \hfill 3^{2M-2}\hfill & \hfill 5^{2M-2}\hfill & \hfill \cdots \hfill & \hfill {\left(2M-1\right)}^{2M-2}\hfill \end{array}\right]\left[\begin{array}{c}\hfill 1a_1\hfill \\ \hfill 3a_2\hfill \\ \hfill 5a_3\hfill \\ \hfill ⋮\hfill \\ \hfill \left(2M-1\right)a_M\hfill \end{array}\right]\hfill \\ \hfill & =\left[\begin{array}{c}\hfill 1\hfill \\ \hfill r_{\text{p}}^2\hfill \\ \hfill r_{\text{p}}^4\hfill \\ \hfill ⋮\hfill \\ \hfill r_{\text{p}}^{2M-2}\hfill \end{array}\right].\hfill \end{array}(19)$

Solving Equation 19 yields

$a_m=\frac{{\left(-1\right)}^{M-1}}{2m-1}{\displaystyle \prod _{\begin{array}{c}1\le k\le M,\\ k\ne m\end{array}}}\frac{r_{\text{p}}^2-{\left(2k-1\right)}^2}{{\left(2m-1\right)}^2-{\left(2k-1\right)}^2}\left(m=\mathrm{1,2},\dots ,M\right).(20)$

Taking $\theta =\mathrm{\pi }/4$, we can derive $d_j={(\sqrt{2}{r}_{\text{p}})}^{2j}(j=\mathrm{0,1,2},\dots ,M-1)$. Substituting this expression into matrix Equation 18 and solving it yields

$a_m=\frac{{\left(-1\right)}^{M-1}}{2m-1}{\displaystyle \prod _{\begin{array}{c}1\le k\le M,\\ k\ne m\end{array}}}\frac{{\left(\sqrt{2}{r}_{\text{p}}\right)}^2-{\left(2k-1\right)}^2}{{\left(2m-1\right)}^2-{\left(2k-1\right)}^2}\left(m=\mathrm{1,2},\dots ,M\right).(21)$

Equation 20 and Equation 21 provide the generic solution of the difference coefficients of TS-SFD calculated based on the time–space domain P-wave dispersion relation by taking $\theta =0$ and $\theta =\mathrm{\pi }/4$, respectively. If $\theta =\mathrm{\pi }/8$, the corresponding generic solution of the difference coefficients is difficult to derive; at this point, it is necessary to calculate the value of $d_0,d_1,d_2,\dots ,d_{M-1}$ using Equation 15 and Equation 17 and then to calculate the difference coefficients by solving Equation 18.

The algorithm for calculating the difference coefficients of TS-SFD based on the time–space domain P-wave dispersion relation is given above. To obtain the algorithm for calculating the difference coefficients based on the time–space domain S-wave dispersion relation, it is only necessary to replace $r_{\text{p}}$ with $r_{\text{s}}$ in the equations arising in the derivation process. Consequently, by replacing $r_{\text{p}}$ with $r_{\text{s}}$, Equation 20 and Equation 21 become the generic solution of the difference coefficients of TS-SFD calculated based on the time–space domain S-wave dispersion relation by taking $\theta =0$ and $\theta =\mathrm{\pi }/4$, respectively.

Comparing the generic solution of the difference coefficients of C-SFD (Equation 5) and TS-SFD (Equation 20; Equation 21), we can observe that the generic solution of the difference coefficients of C-SFD is a special case of the generic solution of the difference coefficients of TS-SFD with $r_{\text{p}}=0$. The difference coefficients of C-SFD are related only to the value of M but are unrelated to the velocity of the seismic wave propagating in the medium. The difference coefficients of TS-SFD are not only related to the value of M and $\theta$ but also vary with the value of the Courant number $r_{\text{p}}=v_{\text{p}}\mathrm{\Delta }t/h$ or $r_{\text{s}}=v_{\text{s}}\mathrm{\Delta }t/h$. The temporal sampling interval $\mathrm{\Delta }t$ and the spatial sampling interval $h$ are usually fixed in the numerical simulation, so the difference coefficients vary adaptively with P-wave velocity $v_{\text{p}}$ or S-wave velocity $v_{\text{s}}$, which is the main reason why TS-SFD can achieve higher modeling accuracy than C-SFD.

3 High-accuracy modeling strategy with the decoupled elastic wave equation

Directly using TS-SFD to solve the elastic wave Equation 1, the P- and S-waves are coupled together, and when the difference coefficients are calculated based on the time–space domain P-wave dispersion relation, the P-wave modeling accuracy is high, while the S-wave modeling accuracy is low; on the contrary, when the difference coefficients are calculated based on the time–space domain S-wave dispersion relation, the S-wave modeling accuracy is high, while the P-wave modeling accuracy is low. In this article, we propose to solve the decoupled P- and S-wave elastic wave equation, ensuring that both P- and S-waves are simulated with high accuracy.

The first-order velocity–stress elastic wave equation with decoupled P- and S-waves (Li et al., 2007) can be expressed as follows:

$\begin{array}{c}\hfill {\upsilon }_x={\upsilon }_x^{\text{p}}+{\upsilon }_x^{\text{s}},{\upsilon }_z={\upsilon }_z^{\text{p}}+{\upsilon }_z^{\text{s}},\end{array}(\mathrm{22a})$

$\begin{array}{ll}\hfill & \begin{array}{cc}\hfill & \frac{\partial {\upsilon }_x^{\text{p}}}{\partial t}=\frac{1}{\rho }\frac{\partial {\tau }_{xx}^{\text{p}}}{\partial x},\frac{\partial {\tau }_{xx}^{\text{p}}}{\partial t}=\left(\lambda +2\mu \right)\left(\frac{\partial {\upsilon }_x}{\partial x}+\frac{\partial {\upsilon }_z}{\partial z}\right),\hfill \\ \hfill & \frac{\partial {\upsilon }_z^{\text{p}}}{\partial t}=\frac{1}{\rho }\frac{\partial {\tau }_{zz}^{\text{p}}}{\partial z},\frac{\partial {\tau }_{zz}^{\text{p}}}{\partial t}=\left(\lambda +2\mu \right)\left(\frac{\partial {\upsilon }_x}{\partial x}+\frac{\partial {\upsilon }_z}{\partial z}\right),\hfill \end{array}\hfill \end{array}(\mathrm{22b})$

$\begin{array}{ll}\hfill & \begin{array}{cc}\hfill & \frac{\partial {\upsilon }_x^{\text{s}}}{\partial t}=\frac{1}{\rho }\left(\frac{\partial {\tau }_{xx}^{\text{s}}}{\partial x}+\frac{\partial {\tau }_{xz}^{\text{s}}}{\partial z}\right),\frac{\partial {\upsilon }_z^{\text{s}}}{\partial t}=\frac{1}{\rho }\left(\frac{\partial {\tau }_{xz}^{\text{s}}}{\partial x}+\frac{\partial {\tau }_{zz}^{\text{s}}}{\partial z}\right),\hfill \\ \hfill & \frac{\partial {\tau }_{xx}^{\text{p}}}{\partial t}=-2\mu \frac{\partial {\upsilon }_z}{\partial z},\frac{\partial {\tau }_{zz}^{\text{p}}}{\partial t}=-2\mu \frac{\partial {\upsilon }_x}{\partial x},\frac{\partial {\tau }_{xz}^{\text{s}}}{\partial t}=\mu \left(\frac{\partial {\upsilon }_x}{\partial z}+\frac{\partial {\upsilon }_z}{\partial x}\right),\hfill \end{array}\hfill \end{array}(\mathrm{22c})$

where ${\upsilon }_x^{\text{p}}$ and ${\upsilon }_z^{\text{p}}$, along with ${\upsilon }_x^{\text{s}}$ and ${\upsilon }_z^{\text{s}}$, represent the P- and S-wave components of the velocity filed of the particle vibration, respectively, and ${\tau }_{xx}^{\text{p}}$ and ${\tau }_{zz}^{\text{p}}$, along with ${\tau }_{xx}^{\text{s}}$, ${\tau }_{zz}^{\text{s}}$ and ${\tau }_{xz}^{\text{s}}$, represent the P- and S-wave components of the stress field, respectively. Equation 22 constitutes the elastic wave equation with decoupled P- and S-waves, where Equation 22b represents the decoupled P-wave equation and Equation 22c represents the decoupled S-wave equation. The main steps of numerical simulation of this equation using TS-SFD are as follows:

1. Derive the corresponding difference-discrete P- and S-wave equations by performing difference discretization on the decoupled P-wave (Equation 22b) and S-wave (Equation 22c).

2. Iteratively solve the difference-discrete P- and S-wave equations, respectively, with the difference coefficients calculated based on the time–space domain P- and S-wave dispersion relation.

3. Calculate ${\upsilon }_x$ and ${\upsilon }_z$ using Equation 22a at the current moment.

4. Repeat steps 2 and 3 until the set maximum recording time is reached.

As observed from the above steps, in the process of solving the decoupled P- and S-wave elastic wave equations with TS-SFD, the decoupled P-wave and S-wave equations are solved with the difference coefficients calculated with the time–space domain P-wave and S-wave dispersion relations, respectively, ensuring that both the P- and S-waves achieve high modeling accuracy. Furthermore, by solving the decoupled P- and S-wave elastic wave equations, the separation of the P- and S-waves can be realized automatically.

4 Numerical dispersion analysis and stability analysis

4.1 Numerical dispersion analysis

The magnitude of the numerical dispersion directly reflects the accuracy of the staggered-grid finite-difference method for numerical simulation of the elastic wave equation. In this article, the normalized P-wave phase velocity error function ${ϵ}_{\text{ph}\left(\text{p}\right)}(kh,\theta )=v_{\text{ph}\left(\text{p}\right)}/v_{\text{p}}-1$ is used to quantitatively describe the numerical dispersion of the P-wave phase velocity. According to the definition of the phase velocity $v_{\text{ph}\left(\text{p}\right)}=\omega /k$ and the time–space domain P-wave dispersion relation, i.e., Equation 11, the expression of the normalized P-wave phase velocity error function ${ϵ}_{\text{ph}\left(\text{p}\right)}(kh,\theta )=v_{\text{ph}\left(\text{p}\right)}/v_{\text{p}}-1$ can be rewritten as

$\begin{array}{ll}\hfill {ϵ}_{\text{ph}\left(\text{p}\right)}\left(kh,\theta \right)& =\frac{2}{r_{\text{p}}kh}si{n}^{-1}\left(r_{\text{p}}\sqrt{q}\right)-1,\hfill \\ \hfill q& ={\left\{{\displaystyle \sum _{m=1}^M}{a}_m\sin \left[\left(m-\frac{1}{2}\right)kh\cos \theta \right]\right\}}^2\hfill \\ \hfill & +{\left\{{\displaystyle \sum _{m=1}^M}{a}_m\sin \left[\left(m-\frac{1}{2}\right)kh\sin \theta \right]\right\}}^2.\hfill \end{array}(23)$

The expression of the normalized P-wave phase velocity error function ${ϵ}_{\text{ph}\left(\text{p}\right)}(kh,\theta )$ given in Equation 23 is suitable for both C-SFD and TS-SFD; only the value of the difference coefficients $a_1,a_2,\dots ,a_M$ is different for them. In the same way, the expression of the normalized S-wave phase velocity error function ${ϵ}_{\text{ph}\left(\text{p}\right)}(kh,\theta )$ can be derived using the time–space domain S-wave dispersion relation, which has the same format as Equation 23, except that $r_{\text{p}}$ should be replaced by $r_{\text{s}}$.

When ${ϵ}_{\text{ph}\left(\text{p}\right)}(kh,\theta )$ or ${ϵ}_{\text{ph}\left(\text{s}\right)}(kh,\theta )$ is equal to 0, i.e., the phase velocity is equal to the real velocity, there is no numerical dispersion; when ${ϵ}_{\text{ph}\left(\text{p}\right)}(kh,\theta )$ or ${ϵ}_{\text{ph}\left(\text{s}\right)}(kh,\theta )$ is greater than 0, i.e., the phase velocity is greater than the real velocity, there is temporal dispersion, which is manifested as phase lead; when ${ϵ}_{\text{ph}\left(\text{p}\right)}(kh,\theta )$ or ${ϵ}_{\text{ph}\left(\text{s}\right)}(kh,\theta )$ is less than 0, i.e., the phase velocity is smaller than the real velocity, there is spatial dispersion, which is manifested as phase lag. According to the expressions of ${ϵ}_{\text{ph}\left(\text{p}\right)}(kh,\theta )$ and ${ϵ}_{\text{ph}\left(\text{s}\right)}(kh,\theta )$, the P- and S-wave phase velocity dispersion curves of C-SFD and TS-SFD can be plotted, which facilitates the analysis of the P- and S-wave phase velocity dispersion characteristics.

Figure 2 presents P-wave phase velocity dispersion curves of TS-SFD $(M=8)$, with the difference coefficients calculated based on the time–space domain P-wave dispersion relation. When calculating the value of the P-wave phase velocity error function for plotting the dispersion curves in this figure, $r_{\text{p}}=0.36$ is fixed, and the value of $\theta$ for calculating the difference coefficients is alternately set as 0, $\mathrm{\pi }/16$, $\mathrm{\pi }/8$, $\mathrm{3}\mathrm{\pi }/16$, and $\mathrm{\pi }/4$.

Figure 2

Figure 2. P-wave phase velocity dispersion curves of TS-SFD $(M=8)$, with the difference coefficients calculated based on the time–space domain P-wave dispersion relation $(r_{\text{p}}=0.36)$. The value of $\theta$ taken for calculating the difference coefficients is (A) $\theta =0$, (B) $\theta =\mathrm{\pi }/16$, (C) $\theta =\mathrm{\pi }/8$, (D) $\theta =\mathrm{3}\mathrm{\pi }/16$, and (E) $\theta =\mathrm{\pi }/4$.

A comparison of the P-wave phase velocity dispersion curves in Figure 2 reveals the following observations: when taking $\theta =0$, it shows mainly temporal dispersion; increasing $\theta$ from 0 to $\mathrm{\pi }/4$, the temporal dispersion decreases, while the spatial dispersion increases. When taking $\theta =\mathrm{\pi }/4$, it shows mainly spatial dispersion. Further comparison of the amplitude of the phase velocity dispersion shows that when taking $\theta =\mathrm{\pi }/4$, the amplitude of the phase velocity dispersion is the largest, while when taking $\theta =\mathrm{\pi }/8$, the amplitude of the phase velocity dispersion is the smallest. Therefore, when calculating the difference coefficients of TS-SFD based on the time–space domain P-wave dispersion relation, taking $\theta =\mathrm{\pi }/8$ ensures that the amplitude of the P-wave phase velocity dispersion reaches the minimum.

Similarly, if we plot the S-wave phase velocity dispersion curves of TS-SFD $(M=8)$, with the difference coefficients calculated based on the time–space domain S-wave dispersion relation, and then analyze the characteristic of the numerical dispersion, a similar conclusion can be obtained: when calculating the difference coefficients of TS-SFD based on the time–space domain S-wave dispersion relation, taking $\theta =\mathrm{\pi }/8$ ensures that the amplitude of the S-wave phase velocity dispersion reaches the minimum. Therefore, $\theta =\mathrm{\pi }/8$ is the optimal choice for the calculation of the difference coefficients of TS-SFD based on the time–space domain P- or S-wave dispersion relation. In the following sections, $\theta =\mathrm{\pi }/8$ is used as the default choice for the calculation of the difference coefficients of TS-SFD based on the time–space domain P- or S-wave dispersion relation, unless specifically stated.

Figure 3 presents P- and S-wave phase velocity dispersion curves of C-SFD $(M=\mathrm{2,5,8})$ and TS-SFD $(M=\mathrm{2,5,8})$. When calculating the value of the P- and S-wave phase velocity error function for plotting the dispersion curves in this figure, $r_{\text{p}}=0.45$ and $r_{\text{p}}=0.25$ are fixed. An analysis and comparison of the P- and S-wave phase velocity dispersion curves in Figure 3 reveal the following findings:

1. For C-SFD$(M=2)$ and TS-SFD$(M=2)$, both P- and S-waves exhibit spatial dispersion, and thus, the modeling accuracy of both P- and S-waves is low.

2. For C-SFD $(M=\mathrm{5,8})$, both P- and S-waves exhibit temporal dispersion, and thus, the modeling accuracy of both P- and S-waves remains low. For TS-SFD $(M=\mathrm{5,8})$, when the difference coefficients are calculated based on the time–space domain P-wave dispersion relation, the P-wave has relatively small numerical dispersion, but the S-wave exhibits obvious spatial dispersion; when the difference coefficients are calculated based on the time–space domain S-wave dispersion relation, the S-wave has relatively small numerical dispersion, but the P-wave exhibits obvious temporal dispersion; when the values of the normalized P- and S-wave phase velocity error functions are computed with the difference coefficients calculated based on the time–space domain P- and S-wave dispersion relations, respectively, both P- and S-waves have relatively small numerical dispersion. In other words, performing elastic numerical simulation by solving the decoupled P- and S-wave elastic wave equation can ensure that both the P and S-waves exhibit relatively small numerical dispersion.

3. The amplitude of both the normalized P- and S-wave phase velocity error functions of TS-SFD $(M=\mathrm{5,8})$, with the difference coefficients calculated based on the time–space domain S-wave dispersion relation, is smaller than that of C-SFD $(M=\mathrm{5,8})$.

Figure 3

Figure 3. P- and S-wave phase velocity dispersion curves of C-SFD $(M=\mathrm{2,5,8})$ and TS-SFD $(M=\mathrm{2,5,8})$ ($r_{\mathrm{p}}=0.45$ and $r_{\mathrm{s}}=0.25$). (A–C) C-SFD $(M=\mathrm{2,5,8})$; (D–F) TS-SFD $(M=\mathrm{2,5,8})$, with the difference coefficients calculated based on the time–space domain P-wave dispersion relation; (G–I) TS-SFD $(M=\mathrm{2,5,8})$, with the difference coefficients calculated based on the time–space domain S-wave dispersion relation; (J–L) TS-SFD $(M=\mathrm{2,5,8})$; the P-wave and S-wave phase velocity curves are plotted with the FD coefficients calculated based on the time–space domain P- and S-wave dispersion relations, respectively.

The above numerical dispersion analysis shows that when M takes a small value, such as $M=2$, TS-SFD shows no significant advantage over C-SFD in suppressing phase velocity numerical dispersion. When M takes a relatively large value, such as $M>=5$, performing elastic numerical simulation by directly solving the velocity–stress elastic wave equation, TS-SFD with the difference coefficients calculated based on the time–space domain S-wave dispersion relation can achieve higher modeling accuracy of both P- and S-waves than C-SFD. When TS-SFD with a relatively large $M$ value is used to conduct elastic numerical simulation by solving the decoupled P- and S-wave elastic wave equations and the difference coefficients for solving the decoupled P- and S-wave equations are calculated based on the time–space domain P- and S-wave dispersion relations, respectively, it can ensure that both P- and S-waves can achieve relatively high modeling accuracy.

4.2 Stability analysis

According to the time–space domain P-wave dispersion relation, i.e., Equation 11, we obtain Equation 24:

$\begin{array}{ll}\hfill {\sin }^2\left(\frac{r_{\text{p}}kh}{2}\right)& \approx r_{\text{p}}^2{\left\{{\displaystyle \sum _{m=1}^M}{a}_m\sin \left[\left(m-\frac{1}{2}\right)k_{x}h\right]\right\}}^2\hfill \\ \hfill & +r_{\text{p}}^2{\left\{{\displaystyle \sum _{m=1}^M}{a}_m\sin \left[\left(m-\frac{1}{2}\right)k_{z}h\right]\right\}}^2\hfill \end{array}(24)$

Taking the maximum space wave number $k_x=k_z=\mathrm{\pi }$ and considering $0\le {\sin }^2(r_{\text{p}}kh/2)\le 1$ yields

$r_{\text{p}}\le S=\frac{1}{\sqrt{2}\left|{\displaystyle \sum _{m=1}^M}{\left(-1\right)}^{m-1}{a}_m\right|},(25)$

where $S$ is the stability factor. Equation 25 is the P-wave stability condition of C-SFD and TS-SFD. Similarly, according to the time–space domain S-wave dispersion relation, i.e., Equation 12, the S-wave stability condition can be derived, which is identical in form to Equation 25, except that $r_{\text{p}}$ is replaced by $r_{\text{s}}$. Considering the general case $r_{\text{p}}>r_{\text{s}}$, the stability of C-SFD and TS-SFD is contingent on the P-wave stability condition.

According to the expression of the P-wave stability condition, the curve of the maximum $r_{\text{p}}$ under the constraint of the stability condition with $M$ can be drawn, which is referred to as the P-wave stability curve. A larger value of the maximum $r_{\text{p}}$ under the constraint of the P-wave stability condition indicates stronger stability.

Figure 4 presents the P-wave stability curves of C-SFD and TS-SFD for directly solving the velocity–stress elastic wave equation. From this figure, we can observe the following: the stability of C-SFD is the weakest; the stability of TS-SFD with the difference coefficients calculated based on the time–space domain S-wave dispersion relation is stronger than that of C-SFD, and the stability of TS-SFD with the difference coefficients calculated based on the time–space domain P-wave dispersion relation is the strongest. In addition, it must be emphasized that when performing elastic numerical simulation by solving the decoupled P- and S-wave elastic wave equations with TS-SFD, with the difference coefficients for solving the decoupled P- and S-wave equations calculated based on the time–space domain P- and S-wave dispersion relation, respectively, the stability is the same as that of TS-SFD for directly solving the velocity–stress elastic wave equation, with the difference coefficients calculated based on the time–space domain P-wave dispersion relation.

Figure 4

Figure 4. P-wave stability curves of C-SFD and TS-SFD for directly solving the velocity–stress elastic wave equation. TS-SFD (S-wave) represents the P-wave stability curve of TS-SFD, with the difference coefficients calculated based on the time–space domain S-wave dispersion relation, and TS-SFD (P-wave) represents the P-wave stability curve of TS-SFD, with the difference coefficients calculated based on the time–space domain P-wave dispersion relation.

5 Numerical simulation example

5.1 Homogeneous model

Figure 5 presents the wavefield snapshots of the ${\upsilon }_z$ component at 2.1 s simulated by C-SFD and TS-SFD for a homogeneous model. The model size is $8\mathrm{k}\mathrm{m}\times 8\mathrm{k}\mathrm{m}$, with the spatial sampling interval $h=10\mathrm{m}$, and thus, the grid number is $801\times 801$. The P- and S-wave velocities of the model are 3,200 and 2080 m/s, respectively. A horizontally oriented force function characterized by a Ricker wavelet with a dominant frequency of 20 Hz located at (0.5 km, 0.5 km) is used as the source wavelet.

Figure 5

Figure 5. Wavefield snapshots of the ${\upsilon }_z$ component at 2.1 s simulated by C-SFD and TS-SFD for a homogeneous model. (A) Snapshot without numerical dispersion simulated by C-SFD $(M=20)$ with a very small temporal sampling interval $\mathrm{\Delta }t=0.1\mathrm{m}\mathrm{s}$ serving as the reference solution. (B–D) Snapshots simulated by TS-SFD $(M=8)$ solving the decoupled P- and S-wave elastic wave equations by taking $\theta =0,\mathrm{\pi }/8,\mathrm{\pi }/4$ for calculating the difference coefficients based on the time–space domain P- and S-wave dispersion relations, respectively.

To analyze the influence of the value of $\theta$, which is taken for calculating the difference coefficients of TS-SFD based on the time–space domain P- and S-wave dispersion relation, on the modeling accuracy, Figure 5A presents the snapshot without numerical dispersion simulated by C-SFD $(M=20)$ with a very small temporal sampling $\mathrm{\Delta }t=0.1\mathrm{m}\mathrm{s}$, which serves as the reference solution. Figures 5B–D presents the snapshots simulated by TS-SFD $(M=8)$, with $\mathrm{\Delta }t=1.5\mathrm{m}\mathrm{s}$ solving the decoupled P- and S-wave elastic wave equation by taking $\theta =0,\mathrm{\pi }/8,\mathrm{\pi }/4$, respectively, for calculating the difference coefficients based on the time–space domain P- and S-wave dispersion relation. Comparing the snapshots simulated by TS-SFD $(M=8)$ with the reference solution, it can be observed that when taking $\theta =0$, both P- and S-wave show slight temporal numerical dispersion, as illustrated by the red arrows in Figure 5B; when taking $\theta =\mathrm{\pi }/8$, both P- and S-waves show no visible numerical dispersion. When taking $\theta =\mathrm{\pi }/4$, both P- and S-wave show slight spatial numerical dispersion, as illustrated by the green arrows in Figure 5D.

The above analysis demonstrates that taking $\theta =\mathrm{\pi }/8$ for calculating the difference coefficients of TS-SFD based on the time–space domain P- and S-wave dispersion relations can ensure that P- and S-waves achieve relatively high modeling accuracy, which is consistent with the conclusion drawn from the numerical dispersion based on Figure 2.

5.2 Layer model

Figure 6 provides a layer model. The model size is $6\mathrm{k}\mathrm{m}\times 6\mathrm{k}\mathrm{m}$, with the spatial sampling interval $h=10\mathrm{m}$, and thus, the grid number is $601\times 601$. Figure 7 provides the wavefield snapshots of the ${\upsilon }_z$ component at 2.4 s simulated by C-SFD $(M=8)$ and TS-SFD $(M=8)$, with $\mathrm{\Delta }t=1.5\mathrm{m}\mathrm{s}$ for the layer model presented in Figure 6. In the simulations, a horizontally oriented force function characterized by a Ricker wavelet with a dominant frequency of 20 Hz located at (0.5 km, 0.5 km) is used as the source wavelet. Figure 8 presents the single-trace waveform at $x=3.2\mathrm{k}\mathrm{m}$ and $x=5.0\mathrm{k}\mathrm{m}$ extracted from the unseparated elastic wavefield snapshots of the ${\upsilon }_z$ component presented in Figure 7.

Figure 6

Figure 6. Layer model. (A) P-wave velocity model; (B) S-wave velocity model.

Figure 7

Figure 7. Wavefield snapshots of the ${\upsilon }_z$ component at 2.4 s simulated by C-SFD and TS-SFD for the layer model shown in Figure 6. Left: unseparated elastic wave wavefield snapshots; middle: separated P-wave wavefield snapshots; right: separated S-wave wavefield snapshots. (A) C-SFD $(M=8)$. (B) TS-SFD $(M=8)$, with the difference coefficients calculated based on the time–space domain P-wave dispersion relation. (C) TS-SFD $(M=8)$, with the difference coefficients calculated based on the time–space domain S-wave dispersion relation. (D) TS-SFD $(M=8)$ solving the decoupled P- and S-wave elastic wave equation, in which the difference coefficients for solving the decoupled P- and S-wave equations are calculated based on the time–space domain P- and S-wave dispersion relations, respectively.

Figure 8

Figure 8. Single-trace waveform at (A) $x=3.2\mathrm{k}\mathrm{m}$ and (B) $x=5.0\mathrm{k}\mathrm{m}$ extracted from the unseparated elastic wavefield snapshots of the ${\upsilon }_z$ component shown in Figure 7, in which the blue curve serving as the reference solution is simulated by C-SFD $(M=20)$, with a very small temporal sampling interval $\mathrm{\Delta }t=0.1\mathrm{m}\mathrm{s}$, and the red curves represent the simulated waveforms.① C-SFD $(M=8)$; ② TS-SFD $(M=8)$, with the difference coefficients calculated based on the time–space domain P-wave dispersion relation; ③ TS-SFD $(M=8)$, with the difference coefficients calculated based on the time–space domain S-wave dispersion relation; ④ TS-SFD $(M=8)$ solving the decoupled P- and S-wave elastic wave equations, in which the difference coefficients for solving the decoupled P- and S-wave equation are calculated based on the time–space domain P- and S-wave dispersion relation, respectively.

From Figures 7, 8, it can be observed that in the wavefield snapshots simulated by C-SFD $(M=8)$, both P- and S-waves show obvious temporal numerical dispersion; in the wavefield snapshots simulated by TS-SFD $(M=8)$, with the difference coefficients calculated based on the time–space domain P-wave dispersion relation, P-wave shows no visible numerical dispersion, but S-wave shows obvious spatial dispersion; in the wavefield snapshots simulated by TS-SFD $(M=8)$ with the difference coefficients calculated based on the time–space domain S-wave dispersion relation, S-wave shows no visible numerical dispersion, but P-wave shows some temporal numerical dispersion; and in the wavefield snapshots simulated by TS-SFD $(M=8)$ solving the decoupled P- and S-wave elastic wave equations, in which the difference coefficients for the decoupled P- and S-wave equation are calculated based on the time–space domain P- and S-wave dispersion relation, respectively, both P-wave and S-wave show no visible numerical dispersion.

The above layer model numerical simulation example demonstrates that directly using TS-SFD to solve the velocity–stress elastic wave equation, in which the P and S-waves are coupled together, irrespective of the difference coefficients calculated based on the time–space domain P- or S-wave dispersion relation, cannot ensure that both P- and S-waves reach relatively high modeling accuracy; performing elastic wave simulation by TS-SFD solving the decoupled P- and S-wave elastic wave equations, with the difference coefficients for solving the decoupled P- and S-wave equations calculated based on the time–space domain P- and S-wave dispersion relations, respectively, can ensure that both P- and S-waves achieve relatively high modeling accuracy. This conclusion is consistent with that drawn from the numerical dispersion based on Figure 3.

5.3 Tectonic-reservoir model from the Sichuan basin

Figure 9 presents a tectonic-reservoir model from the Sichuan basin. The model size is $27.3\mathrm{k}\mathrm{m}\times 7.0\mathrm{k}\mathrm{m}$, with the spatial sampling interval $h=10\mathrm{m}$, and thus, the grid number is $2731\times 701$. To perform elastic wave numerical simulation for this model with C-SFD $(M=8)$ and TS-SFD $(M=8)$ with an temporal sampling interval $\mathrm{\Delta }t=0.75\mathrm{m}\mathrm{s}$, an explosion source characterized by a Ricker wavelet, with a dominant frequency of 20 Hz located at (0.5 km, 0.5 km), is used. In the numerical simulation, C-SFD $(M=8)$ directly solves the velocity–stress elastic wave equation, and TS-SFD $(M=8)$ solves the decoupled P- and S-wave elastic wave equations, with the difference coefficients for solving the decoupled P- and S-wave equation calculated based on the time–space domain P- and S-wave dispersion relations, respectively.

Figure 9

Figure 9. Tectonic-reservoir model from the Sichuan basin. (A) P-wave velocity model; (B) S-wave velocity model.

Figure 10 presents the seismic record of the ${\upsilon }_z$ component simulated by TS-SFD $(M=8)$. It shows that TS-SFD $(M=8)$ performs elastic wave simulation by solving the decoupled P- and S-wave elastic wave equation, thereby automatically and effectively separating the P- and S-waves. In addition, as an explosion source is used, the P-wave energy is dominant in the seismic record, and the S-wave of low energy belongs mainly to the converted wave.

Figure 10

Figure 10. Seismic record of the ${\upsilon }_z$ component simulated by TS-SFD $(M=8)$ solving the decoupled P- and S-wave elastic wave equations for the tectonic-reservoir model from the Sichuan basin presented in Figure 9. (A) Unseparated seismic record. (B) Separated P-wave seismic record. (C) Separated S-wave seismic record.

Figure 11 presents the 1000th trace waveform of the unseparated ${\upsilon }_z$ component simulated by C-SFD $(M=8)$ and TS-SFD $(M=8)$. It shows that the waveform simulated by C-SFD $(M=8)$ exhibits a marked discrepancy from the reference waveform, indicating that both the P- and S-waves simulated by C-SFD $(M=8)$ have significant numerical dispersion, while the waveform simulated by TS-SFD (M = 8) exhibits a high degree of coincidence with the reference waveform, suggesting that the P- and S-waves simulated by TS-SFD $(M=8)$ have no visible numerical dispersion.

Figure 11

Figure 11. The 1000th trace waveform of the unseparated ${\upsilon }_z$ component simulated by C-SFD $(M=8)$ and TS-SFD $(M=8)$ for the tectonic-reservoir model from the Sichuan basin presented in Figure 9, in which the blue curve serving as the reference solution is simulated by C-SFD (M = 20), with a very small temporal sampling $\mathrm{\Delta }t=0.1\mathrm{m}\mathrm{s}$, and the red curves represent the simulated waveforms.① C-SFD $(M=8)$, ② TS-SFD $(M=8)$ solving the decoupled P- and S-wave elastic wave equation, in which the difference coefficients for solving the decoupled P- and S-wave equation are calculated based on the time–space domain P- and S-wave dispersion relations, respectively. (A) Time range 1.5–3.0 s, in which the reflected P-wave dominates; (B) time range 3.0–4.5 s, in which the converted S-wave dominates.

Numerical simulation examples on this tectonic-reservoir model from the Sichuan basin further demonstrate that performing elastic wave numerical simulation with TS-SFD $(M=8)$ by solving the decoupled P- and S-wave elastic wave equation, with the difference coefficients for solving the decoupled P- and S-wave equation calculated based on the time–space domain P- and S-wave dispersion relations, respectively, can ensure that both P- and S-waves achieve relatively high modeling accuracy and can automatically and effectively separate the P- and S-waves. Moreover, this example also demonstrates that this simulation strategy has good applicability to complex models.

6 Conclusion

In this article, we propose a novel strategy to perform elastic wave numerical simulation. This strategy is characterized by TS-SFD solving the decoupled P- and S-wave elastic wave equations, with the difference coefficients for solving the decoupled P- and S-wave equations calculated based on the time–space domain P- and S-wave dispersion relations, respectively. Through the numerical dispersion analysis, stability analysis, and numerical simulation experiments, the following conclusions can be drawn:

1. When calculating the difference coefficients for TS-SFD based on the time–space domain P- and S-wave dispersion relations, the differential coefficients are related to the value of $\theta$, and taking $\theta =\mathrm{\theta }/8$ can ensure that P- and S-waves achieve relatively high modeling accuracy.

2. Numerical dispersion analysis and stability analysis demonstrate that the novel strategy proposed in this article can not only ensure that both P- and S-waves achieve relatively high modeling accuracy but also improve the numerical stability.

3. Numerical simulation experiments demonstrate that the novel strategy proposed in this paper can ensure that both P- and S-waves achieve relatively high modeling accuracy and can automatically and effectively separate the P- and S-waves. Moreover, this strategy has good applicability to complex models.

Despite the distinct advantages offered by the elastic wave numerical simulation method proposed in this article, there are several limitations that require further investigation and resolution:

1. The algorithm for calculating the difference coefficients of TS-SFD is derived based on square grids. A variant of TS-SFD that can be used with rectangular grids should be developed to extend its range of applications.

2. The method described in this article is applicable solely to isotropic elastic wave equations; further development may be pursued to establish simulation techniques suitable for anisotropic elastic wave equations.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

XS: Methodology, Visualization, Writing – original draft. CW: Conceptualization, Formal analysis, Resources, Writing – review and editing. WL: Investigation, Methodology, Validation, Writing – review and editing. DZ: Data curation, Software, Supervision, Writing – review and editing. ZH: Conceptualization, Data curation, Formal analysis, Investigation, Writing – original draft. YF: Investigation, Project administration, Validation, Visualization, Writing – review and editing. YG: Investigation, Methodology, Resources, Supervision, Writing – review and editing. ZZ: Conceptualization, Data curation, Resources, Supervision, Writing – original draft.

Funding

The author(s) declared that financial support was not received for this work and/or its publication.

Conflict of interest

Authors XS, CW, DZ, YF, YG, and ZZ were employed by Shale Gas Research Institute, PetroChina Southwest Oil & Gas Field Company.

Authors CW and YG were employed by BGP Inc., China National Petroleum Corporation.

The remaining author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Alterman, Z., and Karal, F. C. (1968). Propagation of elastic waves in layered media by finite difference methods. Bull. Seismol. Soc. Am. 58, 367–398.

Google Scholar

Cao, J., and Chen, J.-B. (2018). A parameter-modified method for implementing surface topography in elastic-wave finite-difference modeling. Geophysics 83, T313–T332. doi:10.1190/geo2018-0098.1

CrossRef Full Text | Google Scholar

Chen, H., Zhou, H., Zhang, Q., and Chen, Y. (2017). Modeling elastic wave propagation using k-space operator-based temporal high-order staggered-grid finite-difference method. IEEE Trans. Geoscience Remote Sens. 55, 801–815. doi:10.1109/TGRS.2016.2615330

CrossRef Full Text | Google Scholar

Dong, L., Ma, Z., Cao, J., Wang, H., Geng, J., Bing, L., et al. (2000). A staggered-grid high-order difference method of one-order elastic wave equation (in Chinese). Chin. J. Geophys. 43, 411–419.

Google Scholar

Du, Q., Zhu, Y., and Ba, J. (2012). Polarity reversal correction for elastic reverse time migration. Geophysics 77, S31–S41. doi:10.1190/geo2011-0348.1

CrossRef Full Text | Google Scholar

Du, Q., Gong, X., Zhang, M., Zhu, Y., and Fang, G. (2014). 3d ps-wave imaging with elastic reverse-time migration. Geophysics 79, S173–S184. doi:10.1190/geo2013-0253.1

CrossRef Full Text | Google Scholar

Graves, R. W. (1996). Simulating seismic wave propagation in 3d elastic media using staggered-grid finite differences. Bull. Seismol. Soc. Am. 86, 1091–1106. doi:10.1785/bssa0860041091

CrossRef Full Text | Google Scholar

Kane, Y. (1966). Numerical solution of initial boundary value problems involving maxwell’s equations in isotropic media. IEEE Trans. Antennas Propag. 14, 302–307. doi:10.1109/TAP.1966.1138693

CrossRef Full Text | Google Scholar

Levander, A. R. (1988). Fourth-order finite-difference p-sv seismograms. Geophysics 53, 1425–1436. doi:10.1190/1.1442422

CrossRef Full Text | Google Scholar

Li, Z., Zhang, H., Liu, Q., and Han, W. (2007). Numerical simulation of elastic wavefield separation by staggering grid high-order finite-difference algorithm (in Chinese). Oil Geophys. Prospect. 42, 510–515.

Google Scholar

Liu, Y., and Sen, M. K. (2011). Scalar wave equation modeling with time–space domain dispersion-relation-based staggered-grid finite-difference schemes. Bull. Seismol. Soc. Am. 101, 141–159. doi:10.1785/0120100041

CrossRef Full Text | Google Scholar

Ma, D., and Zhu, G. (2003). Numerical modeling of p-wave and s-wave separation in elastic wavefield (in Chinese). Oil Geophys. Prospect. 38, 482–486.

Google Scholar

Peng, G., Liu, W., Guo, N., Hu, Z., Xu, K., and Pei, G. (2022). A time-space domain dispersion-relationship-based staggered-grid finite-difference scheme for modeling the stress-velocity acoustic wave equation (in Chinese). Geophys. Prospect. Petroleum 61, 156–165. doi:10.3969/j.issn.1000-1441.2022.01.016

CrossRef Full Text | Google Scholar

Ren, Z., and Li, Z. C. (2017). Temporal high-order staggered-grid finite-difference schemes for elastic wave propagation. Geophysics 82, T207–T224. doi:10.1190/geo2017-0005.1

CrossRef Full Text | Google Scholar

Ren, Z., and Liu, Y. (2015). Acoustic and elastic modeling by optimal time-space-domain staggered-grid finite-difference schemes. Geophysics 80, T17–T40. doi:10.1190/geo2014-0269.1

CrossRef Full Text | Google Scholar

Singh, S., Tsvankin, I., and Naeini, E. Z. (2018). Bayesian framework for elastic full-waveform inversion with facies information. Lead. Edge 37, 924–931. doi:10.1190/tle37120924.1

CrossRef Full Text | Google Scholar

Su, Z., and Trad, D. (2018). Elastic modeling and reverse time migration. CREWES Res. Rep. 30.

Google Scholar

Sun, R., McMechan, G. A., Hsiao, H.-H., and Chow, J. (2004). Separating p- and s-waves in prestack 3d elastic seismograms using divergence and curl. Geophysics 69, 286–297. doi:10.1190/1.1649396

CrossRef Full Text | Google Scholar

Tan, S., and Huang, L. (2014). An efficient finite-difference method with high-order accuracy in both time and space domains for modelling scalar-wave propagation. Geophys. J. Int. 197, 1250–1267. doi:10.1093/gji/ggu077

CrossRef Full Text | Google Scholar

Virieux, J. (1984). Sh-wave propagation in heterogeneous media: velocity-stress finite-difference method. Geophysics 49, 1933–1942. doi:10.1190/1.1441605

CrossRef Full Text | Google Scholar

Virieux, J. (1986). P-sv wave propagation in heterogeneous media: velocity-stress finite-difference method. Geophysics 51, 889–901. doi:10.1190/1.1442147

CrossRef Full Text | Google Scholar

Virieux, J., and Operto, S. (2009). An overview of full-waveform inversion in exploration geophysics. Geophysics 74, WCC1–WCC26. doi:10.1190/1.3238367

CrossRef Full Text | Google Scholar

Wang, W., McMechan, G. A., and Zhang, Q. (2015). Comparison of two algorithms for isotropic elastic p and s vector decomposition. Geophysics 80, T147–T160. doi:10.1190/geo2014-0563.1

CrossRef Full Text | Google Scholar

Wang, H., He, B., and Shao, X. (2022). Numerical simulation of first-order velocity-dilatation-rotation elastic wave equation with staggered grid (in Chinese). Oil Geophys. Prospect. 57, 1352–1361. doi:10.13810/j.cnki.issn.1000-7210.2022.06.010

CrossRef Full Text | Google Scholar

Wang, M., Li, Y., Chen, Y., Fu, Z., Wang, D., and Tang, F. (2025). An average-derivative based second-order staggered-grid finite-difference scheme for 3d frequency-domain elastic wave modeling. Geophysics 90, T79–T97. doi:10.1190/geo2024-0574.1

CrossRef Full Text | Google Scholar

Wu, X., Liu, Y., and Cai, X. (2018). Elastic wavefield separation methods and their applications in reverse time migration (in Chinese). Oil Geophys. Prospect. 53, 710–721. doi:10.13810/j.cnki.issn.1000-7210.2018.04.007

CrossRef Full Text | Google Scholar

Zhou, H., Liu, Y., and Wang, J. (2022). Elastic wave modeling with high-order temporal and spatial accuracies by a selectively modified and linearly optimized staggered-grid finite-difference scheme. IEEE Trans. Geoscience Remote Sens. 60, 1–22. doi:10.1109/TGRS.2021.3078626

CrossRef Full Text | Google Scholar

Zhou, H., Zhang, L., and Zhang, Y. (2024). Simulation of scalar wave propagation with high-order temporal and spatial accuracy by a new multi-axial staggered-grid finite-difference scheme. Bull. Seismol. Soc. Am. 115, 1–21. doi:10.1785/0120240148

CrossRef Full Text | Google Scholar

Zhu, H. (2017). Elastic wavefield separation based on the helmholtz decomposition. Geophysics 82, S173–S183. doi:10.1190/geo2016-0419.1

CrossRef Full Text | Google Scholar