Finite-Difference Simulation of Elastic Wave with Separation in Pure P- and S-Modes

Abstract

Elastic wave equation simulation offers a way to study the wave propagation when creating seismic data. We implement an equivalent dual elastic wave separation equation to simulate the velocity, pressure, divergence, and curl fields in pure P- and S-modes, and apply it in full elastic wave numerical simulation. We give the complete derivations of explicit high-order staggered-grid finite-difference operators, stability condition, dispersion relation, and perfectly matched layer (PML) absorbing boundary condition, and present the resulting discretized formulas for the proposed elastic wave equation. The final numerical results of pure P- and S-modes are completely separated. Storage and computing time requirements are strongly reduced compared to the previous works. Numerical testing is used further to demonstrate the performance of the presented method.

1. Introduction

Elastic wave numerical simulation has proven to be very efficient for modeling seismic wave propagation and seismic acquisition, and it can also guide seismic processing and seismic interpretation. The traditional numerical modeling method using first-order elastic wave equation can only generate the synthetic seismograms of each component in isotropic medium, in which the P- and S-waves are coupled. To obtain wave field of the pure P- and S-modes, the general method is wave field separation processing of the coupled wave field of each component, but it is difficult to get completely separated seismograms. If we use P- and S-wave equations, respectively, to generate P- and S-wave, the converted P- and S-wave will not appear in wave field; thus, it is not equal to full elastic wave field simulation. Carrying out full separation of wave field modeling of pure P- and S-wave (in its reconstructed wave field, each component comprises P- and S-waves that are fully separated) makes no need of wave field separation in the following multiwave processing [1], and it is suitable not only for the isotropic medium but also for the anisotropic medium [2]. What is more, it is of great practical importance for us to study seismic wave propagation mechanism and structure of geology as well as oil reservoir characterization.

Ma and Zhu [3] presented an elastic wave numerical simulating method using second-order elastic wave equations to separate P- and S-waves with Fourier method. This leads to a new direction for elastic wave modeling study and inversion. But this method does not adapt to widely extended use because of its low efficiency and the difficulty to deal with the absorbing boundary condition. Jianlei et al. [4] proposed an equivalent first-order elastic wave equation, which separates velocity fields into pure P- and S-modes of each component with stress field unchanged. But it is not suitable for the fluid and can cause unstable computational problem. Li et al. [5] presented another equivalent first-order elastic wave equation, which separates both velocity fields and stress fields into pure P- and S-modes. However, when the number of variables of the proposed equation increases, the computational efficiency reduces. Then, we are facing on a problem how to improve efficiency and to be suitable for the fluid case.

At present, using first-order velocity-stress elastic wave equation is the most prevailing and appropriate modeling method [6, 7]. This can be solved using high-order staggered-grid finite-difference scheme. The first advantage is that it can avoid the derivatives to the elastic constant media (the known velocity and density). Thus, it is more efficient and accurate than that of traditional second-order elastic wave equation. The second advantage is that it is very convenient to apply the PML absorbing boundary condition [8]. The PML absorbing boundary condition has proven to be extremely efficient for absorbing all kinds of outgoing waves compared with traditional absorbing boundary condition. The idea is to add artificial absorbing layers with a certain thickness around the computational domain of the interest, in which the outgoing waves are trapped and attenuated [9, 10]. The remarkable advantage is of zero reflection for all angles of incidence and all frequencies with a certain thickness of the artificial boundary [11, 12].

The paper is arranged as follows. First, we derive the equivalent first-order dual elastic wave field separation equation and its corresponding PML absorbing boundary condition. Then, we present discretized formulations with high-order staggered-grid finite-difference scheme. Next, we discuss the stability condition and the dispersion relation. Finally, we present examples that demonstrate the desired effects of the proposed method in elastic wave field separation numerical simulation.

2. Basic Principle

Traditional 3D second-order elastic wave equation can be represented as follows: 𝜌𝜕^2⇀𝑆𝜕𝑡²=𝜆+𝜇∇∇⋅^⇀𝑆+𝜇∇^2⇀𝑆. (1)

Rewriting (1) as 𝜕^2⇀𝑆𝜕𝑡²=𝛼²∇∇⋅^⇀𝑆−𝛽²∇×∇×^⇀𝑆, (2) where ^⇀𝑆=(𝑢,𝑣,𝑤)^𝑇 is particle displacement vector in Cartesian coordinate (𝑥,𝑦,𝑧), 𝜌 is density, 𝜆 and 𝜇 are Lames’ coefficient, 𝛼 is P-wave velocity (𝛼=(𝜆+2𝜇)/𝜌), and 𝛽 is S-wave velocity (𝛽=𝜇/𝜌).

By introducing new variables ^⇀𝑆_𝑝=𝑢_𝑝,𝑣_𝑝,𝑤_𝑝^𝑇,^⇀𝑆_𝑠=𝑢_𝑠,𝑣_𝑠,𝑤_𝑠^𝑇, (3)

(2) becomes ^⇀𝑆=^⇀𝑆_𝑝+^⇀𝑆_𝑠,𝜕^2⇀𝑆_𝑝𝜕𝑡²=𝛼²∇∇⋅^⇀𝑆,𝜕^2⇀𝑆_𝑠𝜕𝑡²=−𝛽²∇×∇×^⇀𝑆. (4)

It is easy to prove the following relations [3]: 𝜕²∇×^⇀𝑆_𝑝𝜕𝑡²𝜕=0,²∇⋅^⇀𝑆_𝑠𝜕𝑡²=0. (5)

Equation (5) indicates that ^⇀𝑆_𝑝 is nonrotational pure P-wave field and ^⇀𝑆_𝑠 is nondispersed pure S-wave field. It means that this method can generate completely separated P- and S-wave by elastic wave numerical modeling.

For 2D case, (4) will be simplified to the following equations: 𝑢=𝑢_𝑝+𝑢_𝑠,𝑤=𝑤_𝑝+𝑤_𝑠,𝜕²𝑢_𝑝𝜕𝑡²=𝛼²𝜕²𝑢𝜕𝑥²+𝜕²𝑤𝜕𝜕𝑥𝜕𝑧₂𝑤_𝑝𝜕𝑡²=𝛼²𝜕²𝑢+𝜕𝜕𝑥𝜕𝑧²𝑤𝜕𝑧²,𝜕²𝑢_𝑠𝜕𝑡²=𝛽²𝜕²𝑢𝜕𝑧²−𝜕²𝑤𝜕𝜕𝑥𝜕𝑧₂𝑤_𝑠𝜕𝑡²=𝛽²−𝜕²𝑢+𝜕𝜕𝑥𝜕𝑧²𝑤𝜕𝑥². (6)

Because of the difficulty and complexity to solve the absorb boundary condition and to improve difference accuracy of the second-order elastic wave equation, especially the mixed partial derivative term 𝜕²/𝜕𝑥𝜕𝑧, therefore,we choose the first-order elastic wave equation with the advantages of avoiding the derivatives to the elastic constant media and applying the PML absorbing boundary condition conveniently.

We introduce new variables again and transform (6) into equivalent first-order elastic wave equations: 𝑣_𝑥=𝜕𝑢𝑣𝜕𝑡^𝑧=𝜕𝑤𝑣𝜕𝑡_𝑥=𝑣_𝑥𝑝+𝑣_𝑥𝑠𝑣_𝑧=𝑣_𝑧𝑝+𝑣_𝑧𝑠,𝜕𝑣_𝑥𝑝𝜕𝑡=𝛼²𝜕𝐴𝜕𝑥𝜕𝑣_𝑧𝑝𝜕𝑡=𝛼²𝜕𝐴𝜕𝑧𝜕𝑣_𝑥𝑠𝜕𝑡=𝛽²𝜕𝐵𝜕𝑧𝜕𝑣_𝑧𝑠𝜕𝑡=−𝛽²𝜕𝐵,𝜕𝑥𝐴=𝜕𝑢+𝜕𝑥𝜕𝑤𝜕𝑧𝐵=𝜕𝑢−𝜕𝑧𝜕𝑤𝜕𝑥𝜕𝐴=𝜕𝑡𝜕𝑣_𝑥+𝜕𝑥𝜕𝑣_𝑧𝜕𝑧𝜕𝐵=𝜕𝑡𝜕𝑣_𝑥−𝜕𝑧𝜕𝑣_𝑧,𝜕𝑥 (7) where 𝑣_𝑥 and 𝑣_𝑧 are particle velocity fields, 𝑣_𝑥𝑝 and 𝑣_𝑥𝑠 are pure P- and S-wave fields in 𝑥-direction, and 𝑣_𝑧𝑝 and 𝑣_𝑧𝑠 are pure P- and S-wave fields in 𝑧-direction. 𝐴 and 𝐵 are auxiliary variables defined above.

We can further derive the following relations: ∇×𝜕𝐴𝜕𝑡=∇×𝜕𝑣_𝑥+𝜕𝑥𝜕𝑣_𝑧𝜕𝑧=∇×∇⋅^⇀𝑢=∇⋅∇×^⇀𝑢=0,∇⋅𝜕𝐵𝜕𝑡=∇⋅𝜕𝑣_𝑥−𝜕𝑧𝜕𝑣_𝑧𝜕𝑥=∇⋅∇×^⇀𝑢=0, (8) where ^⇀𝑢=(𝑣_𝑥,𝑣_𝑧)^𝑇. Equation (8) indicates that 𝐴 and 𝐵 are nonrotational pure P-wave field and nondispersed pure S-wave field, respectively, and they are the divergence field and curl field embedded in (7). It means that our method can achieve dual elastic wave numerical simulation with pure P- and S-modes.

3. Numerical Discretion

For typical grid cell (𝑖,𝑗), take the first-order spatial and temporary derivatives in (7); we have [6, 7] 𝜕𝐺𝜕𝑡^𝑛+1/2_𝑖,𝑗=1𝐺|Δ𝑡^{𝑛+1𝑖,𝑗}−𝑔^𝑛_𝑖,𝑗+𝑜Δ𝑡²𝜕𝐹𝜕𝑥^𝑛_𝑖,𝑗=1Δ𝑥^𝑀_𝑚=1𝐶_𝑚^(𝑀)𝐹|^𝑛_{𝑖+𝑚−1/2,𝑗}−𝐹|^𝑛_{()𝑖−𝑚−1/2,𝑗}+𝑜Δ𝑥^2𝑀, (9) where 𝐹 and 𝐺 denote particle velocity fields and auxiliary variables, Δ𝑥 and Δ𝑧 are spatial steps, Δ𝑡 is temporary step, 𝐶_𝑚^(𝑀) denotes coefficients of 𝑀-order staggered-grid finite-difference scheme, 𝑛 denotes temporary grid, and 𝑖 and 𝑗 are coordinate systems illustrated in Figure 1.

Figure 1 

Local mesh and spatial models for velocity and auxiliary variables updates.

Where P-velocity 𝛼 is located in the integral grid, but S-velocity 𝛽 is in the half-grid, therefore, the average operation is taken to approximate no value grids. With this approximate treatment, the method can be used for heterogeneous media. Mitchell [13] and Boore [14] discussed various approximations and indicated that the following approximate treatment yields comparable results: 𝛼²_{𝑖+1/2,𝑗+1/2}=𝛼²_𝑖+1,𝑗+𝛼²_{𝑖+1,𝑗+1}+𝛼²_𝑖,𝑗+𝛼²_𝑖,𝑗+14𝛽²_𝑖,𝑗=𝛽²_{𝑖+1/2,𝑗+1/2}+𝛽²_{𝑖+1/2,𝑗−1/2}+𝛽²_{𝑖−1/2,𝑗+1/2}+𝛽²_{𝑖−1/2,𝑗−1/2}(4)⁻¹. (10)

Hence, (7) can be rewritten as 𝑣_𝑥^𝑛+1_𝑖,𝑗=𝑣_𝑥𝑝^𝑛+1_𝑖,𝑗+𝑣_𝑥𝑠^𝑛+1_𝑖,𝑗,𝑣_𝑧^𝑛+1_{𝑖+1/2,𝑗+1/2}=𝑣_𝑧𝑝^𝑛+1_{𝑖+1/2,𝑗+1/2}+𝑣_𝑧𝑠^𝑛+1_{𝑖+1/2,𝑗+1/2},𝑣_𝑥𝑝^𝑛+1_𝑖,𝑗=𝑣_𝑥𝑝^𝑛_𝑖,𝑗+𝛼²_𝑖,𝑗Δ𝑡Δ𝑥^𝑀_𝑚=1𝐶_𝑚^(𝑀)𝐴|^𝑛+1/2_{𝑖+𝑚−1/2,𝑗}−𝐴|^𝑛+1/2_{𝑖−(𝑚−1/2),𝑗},𝑣_𝑧𝑝^𝑛+1_{𝑖+1/2,𝑗+1/2}=𝑣_𝑧𝑝^𝑛_{𝑖+1/2,𝑗+1/2}+𝛼²_{𝑖+1/2,𝑗+1/2}Δ𝑡Δ𝑧^𝑀_𝑚=1𝐶_𝑚^(𝑀)𝐴|^𝑛+1/2_{𝑖+1/2,𝑗+𝑚}−𝐴|₍₎^𝑛+1/2_{𝑖+1/2,𝑗−𝑚−1},𝑣_𝑥𝑠^𝑛+1_𝑖,𝑗=𝑣_𝑥𝑠^𝑛_𝑖,𝑗+𝛽²_𝑖,𝑗Δ𝑡Δ𝑧^𝑀_𝑚=1𝐶_𝑚^(𝑀)𝐵|^𝑛+1/2_{𝑖,𝑗+𝑚−1/2}−𝐵|₍₎^𝑛+1/2_{𝑖,𝑗−𝑚−1/2},𝑣_𝑧𝑠^𝑛+1_{𝑖+1/2,𝑗+1/2}=𝑣_𝑧𝑠^𝑛_{𝑖+1/2,𝑗+1/2}−𝛽²_{𝑖+1/2,𝑗+1/2}Δ𝑡Δ𝑥^𝑀_𝑚=1𝐶_𝑚^(𝑀)𝐵|^𝑛+1/2_{𝑖+𝑚,𝑗+1/2}−𝐵|₍₎^𝑛+1/2_{𝑖−𝑚−1,𝑗+1/2},𝐴|^𝑛+1/2_{𝑖+1/2,𝑗}=𝐴|^𝑛−1/2_{𝑖+1/2,𝑗}+Δ𝑡Δ𝑥^𝑀_𝑚=1𝐶_𝑚^(𝑀)𝑣_𝑥^𝑛_{𝑖+𝑚,𝑗}−𝑣_𝑥^𝑛_{()𝑖−𝑚−1,𝑗}+Δ𝑡Δ𝑧^𝑀_𝑚=1𝐶_𝑚^(𝑀)𝑣_𝑧^𝑛_{()𝑖+1/2,𝑗+𝑚−1/2}−𝑣_𝑧^𝑛_{()𝑖+1/2,𝑗−𝑚−1/2},𝐵|^𝑛+1/2_{𝑖,𝑗+1/2}=𝐵|^𝑛−1/2_{𝑖,𝑗+1/2}+Δ𝑡Δ𝑧^𝑀_𝑚=1𝐶_𝑚^(𝑀)𝑣_𝑥^𝑛_{𝑖,𝑗+𝑚}−𝑣_𝑥^𝑛_{()𝑖,𝑗−𝑚−1}−Δ𝑡Δ𝑥^𝑀_𝑚=1𝐶_𝑚^(𝑀)𝑣_𝑧^𝑛_{()𝑖+𝑚−1/2,𝑗+1/2}−𝑣_𝑧^𝑛_{()𝑖−𝑚−1/2,𝑗+1/2}, (11) where 𝐶_𝑚^(𝑀)=(−1)^𝑚+1∏^𝑀_{𝑖=1,𝑖≠𝑚}(2𝑖−1)²/(2𝑚−1)∏^𝑀−1_𝑖=1[(2𝑚−1)²−(2𝑖−1)²], 𝑚=1,…,𝑀. The finite-difference order 𝑀 in the paper is 5, which can guarantee perfect numerical accuracy [15].

4. PML Construction

The PML can be viewed as an analytical continuation of the real coordinates in the complex space [16, 17]. The idea is to construct a new wave equation that permits plane-wave solutions with the form 𝐹^′=Φexp𝑖^⇀𝑘⋅^⇀𝑥−𝜔𝑡−^⇀𝛾⋅^⇀𝑥=Φexp−^⇀𝛾⋅^⇀𝑥exp𝑖^⇀𝑘⋅^⇀𝑥−𝜔𝑡=exp−^⇀𝛾⋅^⇀𝑥𝐹, (12) where 𝐹^′ is the Fourier transformation of 𝐹 with damping term exp[−^⇀𝛾⋅^⇀𝑥] and 𝐹 is Fourier transformation of 𝐹. Here, 𝐹 denotes variables, such as particle velocity fields and auxiliary variables. Φ represents plane wave’s amplitude. ^⇀𝑘=(𝑘_𝑥,𝑘_𝑦,𝑘_𝑧)^𝑇 and ^⇀𝛾=(𝛾_𝑥,𝛾_𝑦,𝛾_𝑧)^𝑇 are wave vector and damping coefficient vector with Cartesian components ^⇀𝑥=(𝑥,𝑦,𝑧)^𝑇, respectively. Consider 𝛾_𝑥,𝛾_𝑦,𝛾_𝑧≥0; this means that the plane waves decay exponentially in each direction of increasing ^⇀𝑥; it can ensure reflection coefficient between the effective region and the PML region is exactly zero for all angles of incidence and all frequencies. Take 𝑥-direction as example; PML formulation consists of the following simple substitution: 𝜕⟹𝜕𝜕𝑥𝜕𝑥^′=𝑖𝜔(𝑥)𝜕𝑖𝜔+𝛾.𝜕𝑥 (13)

This implies the complex transformation of Cartesian variables: 𝑥^′𝑖=𝑥−𝜔^𝑥₀𝛾(𝑠)𝑑𝑠. (14)

Or equivalently, upon difference: 𝜕𝑥𝜕𝑥^′=𝑖𝜔(𝑥).𝑖𝜔+𝛾 (15)

Therefore, the plane-wave solutions with PML in the 𝑥-direction can be rewritten as 𝐹^′−𝑘=Φexp^𝑥𝜔⋅^𝑥₀𝛾(𝑠)𝑑𝑠⋅exp𝑖^⇀𝑘⋅^⇀𝑥−𝜔𝑡. (16)

The PML absorbing boundary condition satisfies the following. If plane waves propagate within the effective region (𝑥≤0), there is no reflection at the interface; it means the PML region (𝑥>0) is matched perfectly with the effective region. If the outgoing waves propagate within the PML region, they are damped. Therefore, the damping coefficient in the PML region can be written as 𝜁_𝑥=𝐹^′𝐹−𝑘=exp^𝑥𝜔⋅^𝑥₀𝛾(𝑠)𝑑𝑠. (17)

Equation (17) implies that the outgoing waves decrease exponentially in each direction of increasing 𝑥, and the damping coefficient 𝜁_𝑥 depends on the direction of wave’s propagation and the magnitude of PML attenuation coefficient 𝛾. More precisely, 𝜁_𝑥 decreases very quickly for waves propagating normally to the interface and decreases more and more slowly as the direction of propagation approaches the parallel to the interface [16].

Equation (7) can be written with a universal format: 𝜕𝐹𝜕𝑡=𝑎𝜕𝐺𝜕𝑥+𝑏𝜕𝐺,𝜕𝑧 (18) where 𝐹 and 𝐺 denote particle velocity fields or auxiliary variables. 𝑎 and 𝑏 denote variable such as P-velocity or S-velocity.

Following the PML splitting equations discussed in the conventional PML methods [18], the local coordinate-system-based PML splitting equations are formulated in a local Cartesian coordinate system with the axis being normal to the local artificial boundary: 𝜕𝐹₌𝜕𝑡=𝑎𝜕𝐺,𝜕𝑥𝜕𝐹_⟂𝜕𝑡=𝑏𝜕𝐺,𝜕𝑧𝐹=𝐹₌+𝐹_⟂, (19) where subscript = means that only the derivative perpendicular to the boundary is considered, and subscript ⟂ means that only the derivatives parallel to the boundary are considered.

We transform (19) into frequency domain with Fourier transformation, then carry out complex substitution using (13), and finally transform it back into time domain, yielding 𝜕𝐹₌𝜕𝑡+𝛾_𝑥𝐹₌=𝑎𝜕𝐺,𝜕𝑥𝜕𝐹_⟂𝜕𝑡+𝛾_𝑧𝐹_⟂=𝑏𝜕𝐺.𝜕𝑧 (20)

Equation (20) can be discretized with finite-difference scheme as follows: 𝐹₌^𝑛+1/2_{𝑖+1/2,𝑗}=1−𝛾_𝑥Δ𝑡/21+𝛾_𝑥𝐹Δ𝑡/2^{=𝑛−1/2𝑖+1/2,𝑗}+Δ𝑡𝑎_{𝑖+1/2,𝑗}Δ𝑥1+𝛾_𝑥Δ𝑡/2^𝑁_𝑚=1𝐶_𝑚^(𝑁)𝐺|^𝑛_{𝑖+𝑚,𝑗}−𝐺^𝑛_{𝑖−𝑚+1,𝑗},𝐹_⟂^𝑛+1/2_{𝑖,𝑗+1/2}=1−𝛾_𝑧Δ𝑡/21+𝛾_𝑧𝐹Δ𝑡/2^{⟂𝑛−1/2𝑖,𝑗+1/2}+𝑎_{𝑖,𝑗+1/2}Δ𝑡Δ𝑥1+𝛾_𝑧Δ𝑡/2^𝑁_𝑚=1𝐶_𝑚^(𝑁)𝐺|^𝑛_{𝑖,𝑗+𝑚}−𝐺|^𝑛_{𝑖,𝑗−𝑚+1}. (21)

Finally, the resulting discretized formulations of the PML equation are derived: 𝑣_𝑥𝑝=^𝑛+1_𝑖,𝑗=1−𝛾_𝑥Δ𝑡/21+𝛾_𝑥𝑣Δ𝑡/2^{𝑥𝑝=𝑛𝑖,𝑗}+Δ𝑡𝛼²_𝑖,𝑗Δ𝑥1+𝛾_𝑥Δ𝑡/2^𝑀_𝑚=1𝐶_𝑚^(𝑀)𝐴|^𝑛+1/2_{𝑖+𝑚−1/2,𝑗}−𝐴|₍₎^𝑛+1/2_{𝑖−𝑚−1/2,𝑗},𝑣_𝑧𝑝⟂^𝑛+1_{𝑖+1/2,𝑗+1/2}=1−𝛾_𝑧Δ𝑡/21+𝛾_𝑧𝑣Δ𝑡/2^{𝑧𝑝⟂𝑛𝑖+1/2,𝑗+1/2}+Δ𝑡𝛼²_{𝑖+1/2,𝑗+1/2}Δ𝑧1+𝛾_𝑧Δ𝑡/2^𝑀_𝑚=1𝐶_𝑚^(𝑀)𝐴|^𝑛+1/2_{𝑖+1/2,𝑗+𝑚}−𝐴|₍₎^𝑛+1/2_{𝑖+1/2,𝑗−𝑚−1},𝑣_𝑥𝑠⟂^𝑛+1_𝑖,𝑗=1−𝛾_𝑧Δ𝑡/21+𝛾_𝑧𝑣Δ𝑡/2^{𝑥𝑠⟂𝑛𝑖,𝑗}+Δ𝑡𝛽²_𝑖,𝑗Δ𝑧1+𝛾_𝑧Δ𝑡/2^𝑀_𝑚=1𝐶_𝑚^(𝑀)𝐵|^𝑛+1/2_{𝑖,𝑗+𝑚−1/2}−𝐵|₍₎^𝑛+1/2_{𝑖,𝑗−𝑚−1/2},𝑣_𝑧𝑠=^𝑛+1_{𝑖+1/2,𝑗+1/2}=1−𝛾_𝑥Δ𝑡/21+𝛾_𝑥𝑣Δ𝑡/2^{𝑧𝑠=𝑛𝑖+1/2,𝑗+1/2}−Δ𝑡𝛽²_{𝑖+1/2,𝑗+1/2}Δ𝑥1+𝛾_𝑥Δ𝑡/2^𝑀_𝑚=1𝐶_𝑚^(𝑀)𝐵|^𝑛+1/2_{𝑖+𝑚,𝑗+1/2}−𝐵|₍₎^𝑛+1/2_{𝑖−𝑚−1,𝑗+1/2},𝐴₌^𝑛+1/2_{𝑖+1/2,𝑗}=1−𝛾_𝑥Δ𝑡/21+𝛾_𝑥𝐴Δ𝑡/2^{=𝑛−1/2𝑖+1/2,𝑗}+Δ𝑡Δ𝑥1+𝛾_𝑥Δ𝑡/2^𝑀_𝑚=1𝐶_𝑚^(𝑀)𝑣_𝑥^𝑛_{𝑖+𝑚,𝑗}−𝑣_𝑥^𝑛_{()𝑖−𝑚−1,𝑗},𝐴_⟂^𝑛+1/2_{𝑖+1/2,𝑗}=1−𝛾_𝑧Δ𝑡/21+𝛾_𝑧𝐴Δ𝑡/2^{⟂𝑛−1/2𝑖+1/2,𝑗}+Δ𝑡Δ𝑧1+𝛾_𝑧Δ𝑡/2^𝑀_𝑚=1𝐶_𝑚^(𝑀)𝑣_𝑧^𝑛_{()𝑖+1/2,𝑗+𝑚−1/2}−𝑣_𝑧^𝑛_{()𝑖+1/2,𝑗−𝑚−1/2},𝐵₌^𝑛+1/2_{𝑖,𝑗+1/2}=1−𝛾_𝑥Δ𝑡/21+𝛾_𝑥𝐵Δ𝑡/2^{=𝑛−1/2𝑖,𝑗+1/2}−Δ𝑡Δ𝑥1+𝛾_𝑥Δ𝑡/2^𝑀_𝑚=1𝐶_𝑚^(𝑀)𝑣_𝑧^𝑛_{()𝑖+𝑚−1/2,𝑗+1/2}−𝑣_𝑧^𝑛_{()𝑖−𝑚−1/2,𝑗+1/2},𝐵_⟂^𝑛+1/2_{𝑖,𝑗+1/2}=1−𝛾_𝑧Δ𝑡/21+𝛾_𝑧𝐵Δ𝑡/2^{⟂𝑛−1/2𝑖,𝑗+1/2}+Δ𝑡Δ𝑧1+𝛾_𝑧Δ𝑡/2^𝑀_𝑚=1𝐶_𝑚^(𝑀)𝑣_𝑥^𝑛_{𝑖,𝑗+𝑚}−𝑣_𝑥^𝑛_{()𝑖,𝑗−𝑚−1},𝐴|^𝑛+1/2_{𝑖+1/2,𝑗}=𝐴₌^𝑛+1/2_{𝑖+1/2,𝑗}+𝐴_⟂^𝑛+1/2_{𝑖+1/2,𝑗}𝐵|^𝑛+1/2_{𝑖,𝑗+1/2}=𝐵₌^𝑛+1/2_{𝑖,𝑗+1/2}+𝐵_⟂^𝑛+1/2_{𝑖,𝑗+1/2}. (22)

5. Stable Condition

The process of elastic wave numerical simulation contains multiple times of wave field iteration. If each iteration process is divergence, then the results are incorrect and cannot be accepted [9, 10]. Therefore, it is necessary to derive the stability condition to avoid the problem. Its idea is to construct a relationship among the discretized parameters.

Rewrite (7) as follows: 𝑣_𝑥𝑝𝑣_𝑥𝑠𝑣_𝑧𝑝𝑣_𝑧𝑠𝐴𝐵_𝑡=0000𝑣_𝑝𝜕0𝜕𝑥00000𝑣_𝑠𝜕𝜕𝑧0000𝑣_𝑝𝜕0𝜕𝑧00000−𝑣_𝑠𝜕𝑣𝜕𝑥_𝑝𝜕𝑣𝜕𝑥^𝑝𝜕𝑣𝜕𝑥^𝑝𝜕𝑣𝜕𝑧^𝑝𝜕𝑣𝜕𝑧00_𝑠𝜕𝑣𝜕𝑧^𝑠𝜕𝜕𝑧−𝑣_𝑠𝜕𝜕𝑥−𝑣_𝑠𝜕×𝑣𝜕𝑥00_𝑥𝑝𝑣_𝑥𝑠𝑣_𝑧𝑝𝑣_𝑧𝑠𝐴𝐵. (23)

We represent (23) with matrix form 𝜕𝑈(𝑥,𝑧,𝑡)(),𝜕𝑡=𝑄𝑈𝑥,𝑧,𝑡 (24) where 𝑈=(𝑣_𝑥𝑝,𝑣_𝑥𝑠,𝑣_𝑧𝑝,𝑣_𝑧𝑠,𝐴,𝐵)^𝑇 and 𝑄 is the matrix containing elastic constants and spatial difference in (23).

Equation (24) indicates the following relation: 𝜕^2𝑘−1𝑈(𝑡)𝜕𝑡^2𝑘−1=𝑄^2𝑘−1𝑈(𝑡). (25)

Substitute the spatial derivatives in (25) with (9), becomes: 𝜕^2𝑘−1𝑈(𝑥,𝑧,𝑡)𝜕𝑡^2𝑘−1=𝑅^2𝑘−1𝑈(𝑥,𝑧,𝑡), (26) where matrix 𝑅 is finite-difference discretized formula of matrix 𝑄.

We use Taylor series expansion technique and then (24) can be written as 𝑈()()𝑥,𝑧,𝑡+Δ𝑡/2=𝑈𝑥,𝑧,𝑡−Δ𝑡/2+2^𝐾_𝑘=1(Δ𝑡/2)^2𝑘−1()!𝜕2𝑘−1^2𝑘−1𝑈(𝑥,𝑧,𝑡)𝜕𝑡^2𝑘−1+𝑂Δ𝑡^2𝐾. (27)

Substituting (26) into (27) yields 𝑈(𝑥,𝑧,𝑡+Δ𝑡/2)=𝑈(𝑥,𝑧,𝑡−Δ𝑡/2)+2^𝐾_𝑘=1(Δ𝑡/2)^2𝑘−1()!𝑅2𝑘−1^2𝑘−1𝑈(𝑥,𝑧,𝑡)+𝑂Δ𝑡^2𝐾+Δ𝑥^2𝑀+Δ𝑧^2𝑀. (28)

Equation (28) is the arbitrary order finite-difference numerical formulation of the proposed method with staggered-grid technique in both temporary and spatial domain.

Transforming (28) into frequency domain yields𝐸sin²𝜔Δ𝑡2+^𝐾_𝑘=1Δ𝑡^2𝑘2^2𝑘()!⋅𝐺2𝑘−1^𝑘⋅𝑈𝑘_𝑥,𝑘_𝑧,𝜔=0. (29)

𝐸 is unit matrix. Complex matrixes 𝐺 and 𝑈 are the Fourier transform of the matrix 𝐺(𝑥,𝑧)=𝑅²(𝑥,𝑧) and the vector 𝑈(𝑥,𝑧,𝑡), respectively, where 𝐺()=𝐷𝑥,𝑧¹¹00𝐷₂₂,𝐷₁₁=𝑣²_𝑝𝜕²𝜕𝑥²𝑣²_𝑝𝜕²𝜕𝑥²𝑣²_𝑝𝜕²𝑣𝜕𝑥𝜕𝑧^2𝑝𝜕²𝑣𝜕𝑥𝜕𝑧_2𝑠𝜕²𝜕𝑧²𝑣²_𝑠𝜕²𝜕𝑧²−𝑣²_𝑠𝜕²𝜕𝑥𝜕𝑧−𝑣²_𝑠𝜕²𝑣𝜕𝑥𝜕𝑧_2𝑝𝜕²𝑣𝜕𝑥𝜕𝑧^2𝑝𝜕²𝑣𝜕𝑥𝜕𝑧^2𝑝𝜕²𝜕𝑧²𝑣²_𝑝𝜕²𝜕𝑧²−𝑣²_𝑠𝜕²𝜕𝑥𝜕𝑧−𝑣²_𝑠𝜕²𝑣𝜕𝑥𝜕𝑧^2𝑠𝜕²𝜕𝑥²𝑣²_𝑠𝜕²𝜕𝑥²,𝐷₂₂=𝑣²_𝑝𝜕²𝜕𝑥²+𝜕²𝜕𝑧²00𝑣²_𝑠𝜕²𝜕𝑥²+𝜕²𝜕𝑧². (30)

The Frequency domain expression of matrix 𝐺(𝑥,𝑧) is 𝐺(𝑥,𝑧)=^𝐷110_0𝐷22, where 𝐷₁₁=−𝑣²_𝑝𝑘²_𝑥−𝑣²_𝑝𝑘²_𝑥−𝑣²_𝑝𝑘_𝑥𝑘_𝑧−𝑣²_𝑝𝑘_𝑥𝑘_𝑧−𝑣²_𝑠𝑘²_𝑧−𝑣²_𝑠𝑘²_𝑧𝑣²_𝑠𝑘_𝑥𝑘_𝑧𝑣²_𝑠𝑘_𝑥𝑘_𝑧−𝑣²_𝑝𝑘_𝑥𝑘_𝑧−𝑣²_𝑝𝑘_𝑥𝑘_𝑧−𝑣²_𝑝𝑘²_𝑧−𝑣²_𝑝𝑘²_𝑧𝑣²_𝑠𝑘_𝑥𝑘_𝑧𝑣²_𝑠𝑘_𝑥𝑘_𝑧−𝑣²_𝑠𝑘²_𝑥−𝑣²_𝑠𝑘²_𝑥,𝐷₂₂=−𝑣²_𝑝𝑘²_𝑥+𝑘²_𝑧00−𝑣²_𝑠𝑘²_𝑥+𝑘²_𝑧. (31)

Because 𝜕²𝑈(𝑥,𝑧,𝑡)/𝜕𝑡²=𝑄²𝑈(𝑥,𝑧,𝑡) is hyperbolic equation system, there must be a nonsingular matrix 𝑆 satisfying the following relation: 𝐺=𝑆⁻¹Λ𝑆, (32) where matrix Λ is diagonal matrix containing six eigenvalues 𝜆_𝑖(𝑖=1,2,…,6).

Substitute (32) into (29), and use similar matrix transform property: 𝑆⁻¹Λ𝑆^𝑘=𝑆⁻¹Λ^𝑘𝑆. (33)

Equation (29) becomes 𝐸sin²𝜔Δ𝑡2+𝑆^−1𝐾_𝑘=1Δ𝑡^2𝑘2^2𝑘()!⋅Λ2𝑘−1^𝑘𝑆⋅𝑈𝜔,𝑘_𝑥,𝑘_𝑧=0. (34)

Equation (34) can be rewritten as 𝑆⁻¹𝐸sin²𝜔Δ𝑡2+^𝐾_𝑘=1Δ𝑡^2𝑘2^2𝑘()!⋅Λ2𝑘−1^𝑘𝑆⋅𝑈𝜔,𝑘_𝑥,𝑘_𝑧=0. (35)

We let Ω=𝑆⁻¹𝑇𝑆=𝑆⁻¹{𝐸sin²(𝜔Δ𝑡/2)+∑^𝐾_𝑘=1(Δ𝑡^2𝑘/2^2𝑘(2𝑘−1)!)⋅Λ^𝑘}𝑆. It means the matrix Ω is similar to diagonal matrix 𝑇. Given the vector 𝑈(𝜔,𝑘_𝑥,𝑘_𝑧) has nonzero solution, (35)’s coefficient matrix must be equal to zero, and it becomes sin²𝜔Δ𝑡2+^𝐾_𝑘=1Δ𝑡^2𝑘2^2𝑘()!⋅𝜆2𝑘−1^𝑘𝑖=0,𝑖=1,2…,6. (36)

In order to maintain stability of the numerical iteration computation, eigenvalues 𝜆_𝑖 must satisfy the following relation: 0≤−sin²𝜔Δ𝑡2=^𝐾_𝑘=1Δ𝑡^2𝑘2^2𝑘()!⋅𝜆2𝑘−1^𝑘𝑖≤1. (37)

According to the property of Fourier transform, spatial wave number and derivative can be expressed as 𝑘_𝑥=2Δ𝑥𝑝𝑘_𝑧=2𝜕Δ𝑧𝑞,FFT2𝜕𝑥=𝑗𝜕Δ𝑥𝑝FFT2𝜕𝑧=𝑗𝜕Δ𝑧𝑞,FFT²𝜕𝑥²2=−𝑝Δ𝑥²𝜕FFT²𝜕𝑧²2=−𝑞Δ𝑧², (38) where 𝑝=∑^𝑀_𝑚=1𝐶_𝑚^(𝑀)sin[Δ𝑥(𝑚−0.5)𝑘_𝑥] and 𝑞=∑^𝑀_𝑚=1𝐶_𝑚^(𝑀)sin[Δ𝑧(𝑚−0.5)𝑘_𝑧]. FFT(⋅) represents Fourier transform.

Solve (36) with (38); first two maximum eigenvalues can be written successively as 𝜆₁=−4𝑣²_𝑝𝑝²Δ𝑥²+𝑞²Δ𝑧²,𝜆₂=−4𝑣²_𝑠𝑝²Δ𝑥²+𝑞²Δ𝑧². (39)

By substituting (39) into (36), the dispersion relation can be derived as follows: sin²𝜔Δ𝑡2+^𝐾_𝑘=1Δ𝑡^2𝑘2^2𝑘()!⋅𝜆2𝑘−1^𝑘1=0,sin²𝜔Δ𝑡2+^𝐾_𝑘=1Δ𝑡^2𝑘2^2𝑘()!⋅𝜆2𝑘−1^𝑘2=0. (40)

According to (39) and (40), we can notice the upper part of (40) is P-wave’s dispersion relation, the lower part is S-wave’s dispersion relation, and they are completely separated.

Substituting Nyquist wave number (𝑘_𝑥=𝜋/Δ𝑥 and 𝑘_𝑧=𝜋/Δ𝑧) into (38) yields 𝑝=𝑞=^𝑀_𝑚=1𝐶_𝑚^(𝑀)(−1)^𝑚−1. (41)

Finally, we can derive the implicit stability condition as follows with (37), (39), and (41): ^𝐾_𝑘=1Δ𝑡^2𝑘𝑣²_𝑝∑^𝑀_𝑚=1𝐶_𝑚^(𝑀)21/Δ𝑥²+1/Δ𝑧^2𝑘(2𝑘−1)!≤1. (42)

Given the discretized parameters of spatial step, medium constant, and difference order, the maximum allowed temporal step Δ𝑡 can be computed with (42). For the second-order temporal difference, the corresponding stability condition is as follows: Δ𝑡≤(Δ𝑥Δ𝑧/𝑣_𝑝∑^𝑀_𝑚=1|𝐶_𝑚^(𝑀)|)1/(Δ𝑥²+Δ𝑧²); let variable 𝜁=∑^𝑀_𝑚=1|𝐶_𝑚^(𝑀)|; we can compute the second- and fourth-order spatial difference’s variable 𝜁=1.1667 and the second- and tenth-order spatial difference’s variable 𝜁=1.317.

6. Numerical Experiments

6.1. Homogeneous Model

To show the application and validity of the proposed method, first we use our code to simulate the waveform response in isotropic homogeneous media. We use Ricker wavelet with dominant frequency of 20 Hz, which is added in horizontal component of velocity field and located in the center of the model (Figure 2(a)). This source can produce both P- and S-wave. P-wave velocity is 2000 m/s; Poisson ratio and density are 0.25 and 1.0 g/cm³, respectively. Temporal step is 0.5 ms and spatial grid size is 5 m in both 𝑥- and 𝑧-direction. We locate a receiver at 200 m, 200 m. The computational results are compared with the analytical solution.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

Figure 2 

Collection design of isotropic model. Consider (a) and snapshots at 0.2 s, (b) is full elastic wave field of horizontal component, (c) is full elastic wave field of vertical component, (d) is the separated pure P-wave of horizontal component, (e) is the separated pure P-wave of vertical component, (f) is the pure P-wave based on divergence operator, (g) is the separated pure S-wave of horizontal component, (h) is the separated pure S-wave of vertical component, and (i) is the pure S-wave based on curl operator.

Figure 2 shows velocity model and snapshots of elastic wave fields of horizontal component and vertical component at 0.2 s. Figures 2(b) and 2(c) are full elastic wave fields. Figures 2(d) and 2(e) are pure P-wave fields. Figures 2(g) and 2(h) are pure S-wave fields. Figures 2(f) and 2(i) are pure P-wave field based on divergence operator and pure S-wave field based on curl operator, respectively. In snapshots of full elastic wave fields, there are P- and S-waves, and they are coupled. By introducing our method, P-wave field and S-wave field are separated completely, which is also shown in Figure 3. The full elastic wave field equals the sum of pure P-wave field and S-wave field of the corresponding elastic wave component, and the difference between analytical solution and numerical solution is almost negligible. What is more, pure P- and S-wave fields based on divergence and curl operators are also computed with perfect numerical accuracy. In summary, the proposed method performs well.

(a)

(b)

(c)

Figure 3 

Receiver waveform comparisons of horizontal component at 200 m, 200 m with temporal length of 0.5 s. (a) denotes the comparison between full elastic wave field and analytical solution, (b) denotes the comparison between full elastic wave field and the separated pure P- and S-wave fields, and (c) denotes the pure P-wave field and S-wave field based on divergence and curl operators, respectively.

Figures 4(a) and 4(b) are the snapshots of the horizontal component with dominant frequency of 30 Hz at 0.2 s with fourth-order and tenth-order spatial difference, respectively. Figures 4(c) and 4(d) are the snapshots of the horizontal component with dominant frequency of 30 Hz with tenth-order spatial difference at 0.5 s and 2 s, respectively. Figure 4 shows the numerical stability of the PML absorbing boundary and the negligible boundary reflection error, and the larger spatial difference order is used, the weaker numerical dispersion occurs.

(a)

(b)

(c)

(d)

Figure 4 

Snapshots of horizontal component with Figure 2(a) survey. (a) is the snapshot at 0.2 s with fourth-order spatial difference, (b) is the snapshot at 0.2 s with tenth-order spatial difference, and (c) and (d) are the snapshot at 0.5 s and 2 s with tenth-order spatial difference.

6.2. Planar Horizontal Model

A model with two layers is tested (Figure 5(a)). P-wave velocity of the upper layer and the lower layer is 2000 m/s and 3000 m/s, respectively. Poisson ratio is 0.25 and density is 1.0 g/cm³. We use Ricker wavelet as seismic source with the dominant frequency of 25 Hz, which is added in the normal stress component and located at 1000 m, 900 m. This source can only produce P-wave. The temporal step is 0.5 ms and spatial grid size is 5 m in both 𝑥- and 𝑧-direction.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

Figure 5 

Collection design of the velocity model with two layers (a) and snapshots at 0.3 s. (a) is the survey of theoretical model, (b) is full elastic wave field of horizontal component, (c) is full elastic wave field of vertical component, (d) is the separated pure P-wave of horizontal component, (e) is the separated pure P-wave of vertical component, (f) is the pure P-wave based on divergence operator, (g) is the separated pure S-wave of horizontal component, (h) is the separated pure S-wave of vertical component, and (i) is the pure S-wave based on curl operator.

Figure 5 shows the designed model and snapshots of elastic wave fields of horizontal component and vertical component at 0.3 s. Figures 5(b) and 5(c) are full elastic wave fields. Figures 5(d) and 5(e) are pure P-wave fields. Figures 5(g) and 5(h) are pure S-wave fields. Figures 5(f) and 5(i) are pure P-wave field based on divergence operator and pure S-wave field based on curl operator, respectively. In snapshots of full elastic wave fields, there are P-wave (HP), reflection P-wave (RPP), reflective S-wave (RPS), sliding P-wave (HPP), slide S-wave (HPS), transmitted P-wave (TPP), and transmitted S-wave (TPS). Head waves (such as HPP and HPS) are formed with incident angle greater than critical angle, and their characteristic can be observed obviously. With our method, all kinds of P-waves and S-waves in the full elastic wave field are completely separated, and the numerical accuracy is satisfactory.

6.3. Marmousi Model

In order to further investigate our method, Marmousi model (Figure 6(a)) is tested. The truncated domain size is 3.4 km in 𝑥-direction and 1.4 km in 𝑧-direction. P-velocity ranges from 1028 m/s to 4670 m/s; Poisson ratio and density are 0.25 and 1.0 g/cm³, respectively. The explosive P-wave source is used and located at 1.7 km, 0.05 km, and it is a Ricker wavelet with the dominant frequency of 25 Hz. In this modeling test, the temporal step is 0.2 ms and spatial grid size is 5 m in both 𝑥- and 𝑧-direction. Trace interval is 5 m and the maximum offset is 1.7 km. Horizontal receiver line is located at the depth of 50 m below the surface.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

Figure 6 

Marmousi model (a) and snapshots at 0.4 s. (a) is the survey of theoretical model, (b) is full elastic wave field of horizontal component, (c) is full elastic wave field of vertical component, (d) is the separated pure P-wave of horizontal component, (e) is the separated pure P-wave of vertical component, (f) is the pure P-wave based on divergence operator, (g) is the separated pure S-wave of horizontal component, (h) is the separated pure S-wave of vertical component, and (i) is the pure S-wave based on curl operator.

Figure 6 shows Marmousi model and snapshots of elastic wave fields of horizontal component and vertical component at 0.4 s. Figures 6(b) and 6(c) are full elastic wave fields. Figures 6(d) and 6(e) are pure P-wave fields. Figures 6(g) and 6(h) are pure S-wave fields. Figures 6(f) and 6(i) are pure P-wave field based on divergence operator and pure S-wave field based on curl operator, respectively.

Figure 7 shows the numerical records of our method. Figures 7(a) and 7(b) are full elastic wave fields. Figures 7(c) and 7(d) are pure P-wave fields. Figures 7(f) and 7(g) are pure S-wave fields. Figures 7(e) and 7(h) are pure P-wave field based on divergence operator and pure S-wave field based on curl operator, respectively. The full elastic wave fields are very complex in this heterogeneous medium, and it is difficult to identify pure P-waves and S-waves. With the proposed method, pure P- and S-waves in full elastic wave fields are separated automatically and completely. It proves that our method can simulate separated pure P- and S-wave field in complex media and achieve the paper’s destination. What is more, with application of high-order staggered-grid finite-difference method and PML absorbing boundary condition, no dispersion and spurious boundary reflection phenomenon appear, and entire modeling process is stable.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 7 

Numerical records of dual separated elastic wave method. (a) is full elastic wave field of horizontal component, (b) is full elastic wave field of vertical component, (c) is the separated pure P-wave of horizontal component, (d) is the separated pure P-wave of vertical component, (e) is the pure P-wave based on divergence operator, (f) is the separated pure S-wave of horizontal component, (g) is the separated pure S-wave of vertical component, and (h) is the pure S-wave based on curl operator.

7. Computational Cost Comparison

Numerical simulating cost depends mostly on the number of variables allocated for memory storage and the number of equations for iteration. For a case study, we use 2D homogenous model with same computer condition (CPU) to do the actual cost comparison. The model parameters are as follows. Grid size is 800 grids × 800 grids and the PML thickness is 50 grids in the damping boundary; spatial grid size is 5 m in both 𝑥- and 𝑧-direction. P-velocity is 2000 m/s; Poisson ratio and density are 0.25 and 1.0 g/cm³, respectively. Spatial step is 5 m in each direction, and temporary step is 0.5 ms. Dominant frequency is 40 Hz, and recording time is 5 s. Our method has only 8 equations and 10 variables and it is with brief form and lowest computational complexity compared with Zhang and Li’s methods. Computational results in Table 1 show the best efficiency of the proposed method. What is more, additional tests show the similar accuracy of above separating methods with same finite-difference order.

Table 1 

Comparison with different separating methods.

We further carry out the numerical experiments on the fourth-order and tenth-order spatial difference cases with the above acquisition parameters. Take fourth-order spatial difference case as example, the spatial grid size is reduced by half, the total number of grids in 2D model is increased by a factor of four, and the number of temporal steps with the same recording length is also doubled, while the computation amount with the spatial difference order is reduced approximately by half; therefore, the total computational time is increased by a factor of four compared with tenth-order spatial difference case. The results are shown in Table 2, which verify the better computational efficiency with high-order spatial difference under the same error bound.

Table 2 

Comparison with different difference orders of the proposed method.

8. Conclusions

Equivalent dual elastic wave separating method is presented. It is reformulated from the traditional elastic wave equation and can be extended to heterogeneous media. PML absorbing boundary condition and the numerical discretion with high-order stagger-grid finite-difference scheme are detailed, as well as stability condition and dispersion relation. With three numerical examples, we can conclude that the pure P- and S-waves are naturally separated in the numerical modeling instead of the separation of wave field afterwards. The separated pure P- and S-waves in each elastic wave component are comparable with the pure P- or S-waves based on divergence and curl operators. Compared with the previous works, the computational cost of the proposed method is the smallest, and it is suitable for fluid case. The complex model demonstrates that our modeling method works well.

Appendix

Previous Works

(a) The traditional elastic wave equation can be represented as 𝜕𝑣_𝑥=𝜕𝑡𝜕𝜏_𝑥𝑥+𝜕𝑥𝜕𝜏_𝑥𝑧𝜕𝑧𝜕𝑣_𝑧=𝜕𝑡𝜕𝜏_𝑥𝑧+𝜕𝑥𝜕𝜏_𝑧𝑧,𝜕𝑧𝜕𝜏_𝑥𝑥𝜕𝑡=𝑣²_𝑝𝜕𝑣_𝑥+𝑣𝜕𝑥^2𝑝−2𝑣²_𝑠𝜕𝑣_𝑧𝜕𝑧𝜕𝜏_𝑧𝑧=𝑣𝜕𝑡^2𝑝−2𝑣²_𝑠𝜕𝑣_𝑥𝜕𝑥+𝑣²_𝑝𝜕𝑣_𝑧𝜕𝑧𝜕𝜏_𝑥𝑧𝜕𝑡=𝑣²_𝑠𝜕𝑣_𝑧+𝜕𝑥𝜕𝑣_𝑥,𝜕𝑧 (A.1) where 𝜏_𝑥𝑥 and 𝜏_𝑧𝑧 are normal stress and 𝜏_𝑥𝑧 is shear stress.

(b) Li et al. [5] equivalent elastic wave equation can be written as 𝑣_𝑥=𝑣_𝑥𝑝+𝑣_𝑥𝑠𝑣_𝑧=𝑣_𝑧𝑝+𝑣_𝑧𝑠,𝜕𝑣_𝑥𝑝=𝜕𝑡𝜕𝜏_{𝑝,𝑥𝑥}𝜕𝑥𝜕𝑣_𝑧𝑝=𝜕𝑡𝜕𝜏_{𝑝,𝑧𝑧}𝜕𝑧𝜕𝜏_{𝑝,𝑥𝑥}𝜕𝑡=𝑣²_𝑝𝜕𝑣_𝑥+𝜕𝑥𝜕𝑣_𝑧𝜕𝑧𝜕𝜏_{𝑝,𝑧𝑧}𝜕𝑡=𝑣²_𝑝𝜕𝑣_𝑥+𝜕𝑥𝜕𝑣_𝑧,𝜕𝑧𝜕𝑣_𝑥𝑠=𝜕𝑡𝜕𝜏_{𝑠,𝑥𝑥}+𝜕𝑥𝜕𝜏_{𝑠,𝑥𝑧}𝜕𝑧𝜕𝑣_𝑧𝑠=𝜕𝑡𝜕𝜏_{𝑠,𝑥𝑧}+𝜕𝑥𝜕𝜏_{𝑠,𝑧𝑧}𝜕𝑧𝜕𝜏_{𝑠,𝑥𝑥}𝜕𝑡=−2𝑣²_𝑠𝜕𝑣_𝑧𝜕𝑧𝜕𝜏_{𝑠,𝑧𝑧}𝜕𝑡=−2𝑣²_𝑠𝜕𝑣_𝑥𝜕𝑥𝜕𝜏_{𝑠,𝑥𝑧}𝜕𝑡=𝑣²_𝑠𝜕𝑣_𝑥+𝜕𝑧𝜕𝑣_𝑧,𝜕𝑥 (A.2) where 𝜏_{𝑝,𝑥𝑥} and 𝜏_{𝑝,𝑧𝑧} are normal stress of pure P-wave, 𝜏_{𝑠,𝑥𝑥} and 𝜏_{𝑠,𝑥𝑥} are normal stress of pure shear wave, and 𝜏_{𝑠,𝑥𝑧} is shear stress of pure shear wave.

(c) Jianlei et al. [4] equivalent elastic wave equation can be expressed as follows: 𝑣_𝑥=𝑣_𝑥𝑝+𝑣_𝑥𝑠𝑣_𝑧=𝑣_𝑧𝑝+𝑣_𝑧𝑠,𝜕𝑣_𝑥𝑝=𝑣𝜕𝑡2𝜆+^2𝑝2𝜇𝜕𝜏_𝑥𝑥+𝜕𝑥𝜕𝜏_𝑧𝑧𝜕𝑥𝜕𝑣_𝑧𝑝=𝑣𝜕𝑡2𝜆+^2𝑝2𝜇𝜕𝜏_𝑥𝑥+𝜕𝑧𝜕𝜏_𝑧𝑧𝜕𝑧𝜕𝑣_𝑥𝑠𝜕𝑡=𝑣²_𝑠1𝜇𝜕𝜏_𝑥𝑧−1𝜕𝑧2𝜆𝜇+2𝜇²𝜆+2𝜇𝜕𝜏_𝑧𝑧𝜕𝑥−𝜆𝜕𝜏_𝑥𝑥𝜕𝑥𝜕𝑣_𝑧𝑠𝜕𝑡=𝑣²_𝑠1𝜇𝜕𝜏_𝑥𝑧−1𝜕𝑥2𝜆𝜇+2𝜇²𝜆+2𝜇𝜕𝜏_𝑥𝑥𝜕𝑧−𝜆𝜕𝜏_𝑧𝑧,𝜕𝑧𝜕𝜏_𝑥𝑥=𝜕𝑡𝜆+2𝜇𝜕𝑣_𝑥𝜕𝑥+𝜆𝜕𝑣_𝑧𝜕𝑧𝜕𝜏_𝑧𝑧=𝜕𝑡𝜆+2𝜇𝜕𝑣_𝑧𝜕𝑧+𝜆𝜕𝑣_𝑥𝜕𝑥𝜕𝜏_𝑥𝑧𝜕𝑡=𝜇𝜕𝑣_𝑥+𝜕𝑧𝜕𝑣_𝑧,𝜕𝑥 (A.3) where 𝜏_𝑥𝑥 and 𝜏_𝑧𝑧 are normal stress and 𝜏_𝑥𝑧 is shear stress.

Conflict of Interests

The author declares that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

Thanks are due to the anonymous reviewer for constructive comments that improved the paper.