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

Authors: Ke-Yang Chen
Journal: Journal of Computational Methods in Physics
 (2014)

---

## 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.


---

## Body

## 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:(1)ρ∂2S⃑∂t2=(λ+μ)∇(∇·S⃑)+μ∇2S⃑.Rewriting (1) as
(2)∂2S⃑∂t2=α2∇(∇·S⃑)-β2∇×(∇×S⃑),
where S⃑=(u,v,w)T is particle displacement vector in Cartesian coordinate (x,y,z), ρ is density, λ and μ are Lames’ coefficient, α is P-wave velocity (α=(λ+2μ)/ρ), and β is S-wave velocity (β=μ/ρ).By introducing new variables(3)S⃑p=(up,vp,wp)T,S⃑s=(us,vs,ws)T,(2) becomes
(4)S⃑=S⃑p+S⃑s,∂2S⃑p∂t2=α2∇(∇·S⃑),∂2S⃑s∂t2=-β2∇×(∇×S⃑).It is easy to prove the following relations [3]:
(5)∂2∇×S⃑p∂t2=0,∂2∇·S⃑s∂t2=0.Equation (5) indicates that S⃑p is nonrotational pure P-wave field and S⃑s 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)u=up+us,w=wp+ws,∂2up∂t2=α2(∂2u∂x2+∂2w∂x∂z)∂2wp∂t2=α2(∂2u∂x∂z+∂2w∂z2),∂2us∂t2=β2(∂2u∂z2-∂2w∂x∂z)∂2ws∂t2=β2(-∂2u∂x∂z+∂2w∂x2).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∂2/∂x∂z, 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)vx=∂u∂tvz=∂w∂tvx=vxp+vxsvz=vzp+vzs,∂vxp∂t=α2∂A∂x∂vzp∂t=α2∂A∂z∂vxs∂t=β2∂B∂z∂vzs∂t=-β2∂B∂x,A=∂u∂x+∂w∂zB=∂u∂z-∂w∂x∂A∂t=∂vx∂x+∂vz∂z∂B∂t=∂vx∂z-∂vz∂x,
where vx and vz are particle velocity fields, vxp and vxs are pure P- and S-wave fields in x-direction, and vzp and vzs are pure P- and S-wave fields in z-direction. A and B are auxiliary variables defined above.We can further derive the following relations:(8)∇×(∂A∂t)=∇×(∂vx∂x+∂vz∂z)∇×(∂A∂t)=∇×∇·u⃑=∇·∇×u⃑=0,ii∇·(∂B∂t)=∇·(∂vx∂z-∂vz∂x)=∇·∇×u⃑=0,
where u⃑=(vx,vz)T. Equation (8) indicates that A and B 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(i,j), take the first-order spatial and temporary derivatives in (7); we have [6, 7]
(9)∂G∂t|i,jn+1/2=1Δt{G|i,jn+1-g|i,jn}+o(Δt2)IIII∂F∂x|i,jn=1Δx∑m=1MCm(M){F|i+m-1/2,jn-F|i-(m-1/2),jn}∂F∂x|i,jn=IIII+o(Δx2M),
where F and G denote particle velocity fields and auxiliary variables, Δx and Δz are spatial steps, Δt is temporary step, Cm(M) denotes coefficients of M-order staggered-grid finite-difference scheme, n denotes temporary grid, and i and j 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:
(10)α2|i+1/2,j+1/2=α2|i+1,j+α2|i+1,j+1+α2|i,j+α2|i,j+14β2|i,j=(β2|i+1/2,j+1/2+β2|i+1/2,j-1/2β2|i,jii=+β2|i-1/2,j+1/2+β2|i-1/2,j-1/2)(4)-1.Hence, (7) can be rewritten as
(11)vx|i,jn+1=vxp|i,jn+1+vxs|i,jn+1,vz|i+1/2,j+1/2n+1=vzp|i+1/2,j+1/2n+1+vzs|i+1/2,j+1/2n+1,vxp|i,jn+1=vxp|i,jnvxp|i,jn+1+α2|i,jΔtΔx∑m=1MCm(M){A|i+m-1/2,jn+1/2-A|i-(m-1/2),jn+1/2},vzp|i+1/2,j+1/2n+1=vzp|i+1/2,j+1/2n+α2|i+1/2,j+1/2ΔtΔz∑m=1MCm(M){A|i+1/2,j+mn+1/2iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii-A|i+1/2,j-(m-1)n+1/2},vxs|i,jn+1=vxs|i,jn+β2|i,jΔtΔz∑m=1MCm(M){B|i,j+m-1/2n+1/2-B|i,j-(m-1/2)n+1/2},vzs|i+1/2,j+1/2n+1=vzs|i+1/2,j+1/2n-β2|i+1/2,j+1/2ΔtΔx∑m=1MCm(M){B|i+m,j+1/2n+1/2iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii-B|i-(m-1),j+1/2n+1/2},A|i+1/2,jn+1/2=A|i+1/2,jn-1/2+ΔtΔx∑m=1MCm(M){vx|i+m,jn-vx|i-(m-1),jn}+ΔtΔz∑m=1MCm(M){vz|i+1/2,j+(m-1/2)n-vz|i+1/2,j-(m-1/2)n},B|i,j+1/2n+1/2=B|i,j+1/2n-1/2+ΔtΔz∑m=1MCm(M){vx|i,j+mn-vx|i,j-(m-1)n}-ΔtΔx∑m=1MCm(M){vz|i+(m-1/2),j+1/2n-vz|i-(m-1/2),j+1/2n},
where Cm(M)=(-1)m+1∏i=1,i≠mM(2i-1)2/(2m-1)∏i=1M-1[(2m-1)2-(2i-1)2], m=1,…,M. The finite-difference order M 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
(12)F~′=Φexp⁡[i(k⃑·x⃑-ωt)-γ⃑·x⃑]=Φexp⁡[-γ⃑·x⃑]exp⁡[i(k⃑·x⃑-ωt)]=exp⁡[-γ⃑·x⃑]F~,
where F~′ is the Fourier transformation of F with damping term exp⁡[-γ⃑·x⃑] and F~ is Fourier transformation of F. Here, F denotes variables, such as particle velocity fields and auxiliary variables. Φ represents plane wave’s amplitude. k⃑=(kx,ky,kz)T and γ⃑=(γx,γy,γz)T are wave vector and damping coefficient vector with Cartesian components x⃑=(x,y,z)T, respectively. Consider γx,γy,γz≥0; this means that the plane waves decay exponentially in each direction of increasing x⃑; 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 x-direction as example; PML formulation consists of the following simple substitution:
(13)∂∂x⟹∂∂x′=iωiω+γ(x)∂∂x.This implies the complex transformation of Cartesian variables:(14)x′=x-iω∫0xγ(s)ds.Or equivalently, upon difference:(15)∂x∂x′=iωiω+γ(x).Therefore, the plane-wave solutions with PML in thex-direction can be rewritten as
(16)F~′=Φexp⁡[-kxω·∫0xγ(s)ds]·exp⁡[i(k⃑·x⃑-ωt)].The PML absorbing boundary condition satisfies the following. If plane waves propagate within the effective region (x≤0), there is no reflection at the interface; it means the PML region (x>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
(17)ζx=∥F~′∥∥F~∥=exp⁡[-kxω·∫0xγ(s)ds].Equation (17) implies that the outgoing waves decrease exponentially in each direction of increasing x, and the damping coefficient ζx depends on the direction of wave’s propagation and the magnitude of PML attenuation coefficient γ. More precisely, ζx 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)∂F∂t=a∂G∂x+b∂G∂z,
where F and G denote particle velocity fields or auxiliary variables. a and b 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)∂F=∂t=a∂G∂x,∂F⊥∂t=b∂G∂z,F=F=+F⊥,
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)∂F=∂t+γxF==a∂G∂x,∂F⊥∂t+γzF⊥=b∂G∂z.Equation (20) can be discretized with finite-difference scheme as follows:
(21)F=|i+1/2,jn+1/2=(1-γxΔt/2)(1+γxΔt/2)F=|i+1/2,jn-1/2F=|i+1/2,jn+1/2+Δtai+1/2,jΔx(1+γxΔt/2)∑m=1NCm(N)(G|i+m,jn-Gi-m+1,jn),F⊥|i,j+1/2n+1/2=(1-γzΔt/2)(1+γzΔt/2)F⊥|i,j+1/2n-1/2F=|i+1/2,jn+1/2+ai,j+1/2ΔtΔx(1+γzΔt/2)∑m=1NCm(N)(G|i,j+mn-G|i,j-m+1n).Finally, the resulting discretized formulations of the PML equation are derived:(22)vxp=|i,jn+1=(1-γxΔt/2)(1+γxΔt/2)vxp=|i,jn+Δtαi,j2Δx(1+γxΔt/2)∑m=1MCm(M){A|i+m-1/2,jn+1/2-A|i-(m-1/2),jn+1/2},vzp⊥|i+1/2,j+1/2n+1=(1-γzΔt/2)(1+γzΔt/2)vzp⊥|i+1/2,j+1/2n+Δtαi+1/2,j+1/22Δz(1+γzΔt/2)∑m=1MCm(M){A|i+1/2,j+mn+1/2-A|i+1/2,j-(m-1)n+1/2},vxs⊥|i,jn+1=(1-γzΔt/2)(1+γzΔt/2)vxs⊥|i,jn+Δtβi,j2Δz(1+γzΔt/2)∑m=1MCm(M){B|i,j+m-1/2n+1/2-B|i,j-(m-1/2)n+1/2},vzs=|i+1/2,j+1/2n+1=(1-γxΔt/2)(1+γxΔt/2)vzs=|i+1/2,j+1/2n-Δtβi+1/2,j+1/22Δx(1+γxΔt/2)∑m=1MCm(M){B|i+m,j+1/2n+1/2-B|i-(m-1),j+1/2n+1/2},A=|i+1/2,jn+1/2=(1-γxΔt/2)(1+γxΔt/2)A=|i+1/2,jn-1/2+ΔtΔx(1+γxΔt/2)∑m=1MCm(M){vx|i+m,jn-vx|i-(m-1),jn},A⊥|i+1/2,jn+1/2=(1-γzΔt/2)(1+γzΔt/2)A⊥|i+1/2,jn-1/2+ΔtΔz(1+γzΔt/2)∑m=1MCm(M){vz|i+1/2,j+(m-1/2)niiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii-vz|i+1/2,j-(m-1/2)n},B=|i,j+1/2n+1/2=(1-γxΔt/2)(1+γxΔt/2)B=|i,j+1/2n-1/2-ΔtΔx(1+γxΔt/2)∑m=1MCm(M){vz|i+(m-1/2),j+1/2niiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii-vz|i-(m-1/2),j+1/2n},B⊥|i,j+1/2n+1/2=(1-γzΔt/2)(1+γzΔt/2)B⊥|i,j+1/2n-1/2+ΔtΔz(1+γzΔt/2)∑m=1MCm(M){vx|i,j+mn-vx|i,j-(m-1)n},A|i+1/2,jn+1/2=A=|i+1/2,jn+1/2+A⊥|i+1/2,jn+1/2B|i,j+1/2n+1/2=B=|i,j+1/2n+1/2+B⊥|i,j+1/2n+1/2.


## 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:
(23)(vxpvxsvzpvzsAB)|t=(0000vp∂∂x000000vs∂∂z0000vp∂∂z000000-vs∂∂xvp∂∂xvp∂∂xvp∂∂zvp∂∂z00vs∂∂zvs∂∂z-vs∂∂x-vs∂∂x00)×(vxpvxsvzpvzsAB).We represent (23) with matrix form
(24)∂U(x,z,t)∂t=QU(x,z,t),
where U=(vxp,vxs,vzp,vzs,A,B)T and Q is the matrix containing elastic constants and spatial difference in (23).Equation (24) indicates the following relation:
(25)∂2k-1U(t)∂t2k-1=Q2k-1U(t).Substitute the spatial derivatives in (25) with (9), becomes: (26)∂2k-1U(x,z,t)∂t2k-1=R2k-1U(x,z,t),
where matrix R is finite-difference discretized formula of matrix Q.We use Taylor series expansion technique and then (24) can be written as
(27)U(x,z,t+Δt/2)=U(x,z,t-Δt/2)+2∑k=1K(Δt/2)2k-1(2k-1)!∂2k-1U(x,z,t)∂t2k-1+O(Δt2K).Substituting (26) into (27) yields
(28)U(x,z,t+Δt/2)=U(x,z,t-Δt/2)+2∑k=1K(Δt/2)2k-1(2k-1)!R2k-1U(x,z,t)+O(Δt2K+Δx2M+Δz2M).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(29)[Esin2(ωΔt2)+∑k=1KΔt2k22k(2k-1)!·G~k]·U~(kx,kz,ω)=0.E is unit matrix. Complex matrixes G~ and U~ are the Fourier transform of the matrix G(x,z)=R2(x,z) and the vector U(x,z,t), respectively, where
(30)G(x,z)=(D1100D22),D11=(vp2∂2∂x2vp2∂2∂x2vp2∂2∂x∂zvp2∂2∂x∂zvs2∂2∂z2vs2∂2∂z2-vs2∂2∂x∂z-vs2∂2∂x∂zvp2∂2∂x∂zvp2∂2∂x∂zvp2∂2∂z2vp2∂2∂z2-vs2∂2∂x∂z-vs2∂2∂x∂zvs2∂2∂x2vs2∂2∂x2),D22=(vp2(∂2∂x2+∂2∂z2)00vs2(∂2∂x2+∂2∂z2)).The Frequency domain expression of matrixG(x,z) is G~(x,z)=(D~1100D~22), where
(31)D~11=(-vp2kx2-vp2kx2-vp2kxkz-vp2kxkz-vs2kz2-vs2kz2vs2kxkzvs2kxkz-vp2kxkz-vp2kxkz-vp2kz2-vp2kz2vs2kxkzvs2kxkz-vs2kx2-vs2kx2),D~22=(-vp2(kx2+kz2)00-vs2(kx2+kz2)).Because∂2U(x,z,t)/∂t2=Q2U(x,z,t) is hyperbolic equation system, there must be a nonsingular matrix S~ satisfying the following relation:
(32)G~=S~-1Λ~S~,
where matrix Λ~ is diagonal matrix containing six eigenvalues λi(i=1,2,…,6).Substitute (32) into (29), and use similar matrix transform property:
(33)(S~-1Λ~S~)k=S~-1Λ~kS~.Equation (29) becomes
(34){Esin2(ωΔt2)+S~-1[∑k=1KΔt2k22k(2k-1)!·Λ~k]S~}·U~(ω,kx,kz)=0.Equation (34) can be rewritten as
(35)S~-1{Esin2(ωΔt2)+∑k=1KΔt2k22k(2k-1)!·Λ~k}S~·U~(ω,kx,kz)=0.We letΩ=S~-1TS~=S~-1{Esin2(ωΔt/2)+∑k=1K(Δt2k/22k(2k-1)!)·Λ~k}S~. It means the matrix Ω is similar to diagonal matrix T. Given the vector U~(ω,kx,kz) has nonzero solution, (35)’s coefficient matrix must be equal to zero, and it becomes
(36)sin2(ωΔt2)+∑k=1KΔt2k22k(2k-1)!·λ~ik=0,i=1,2…,6.In order to maintain stability of the numerical iteration computation, eigenvaluesλ~i must satisfy the following relation:
(37)0≤|-sin2(ωΔt2)|=|∑k=1KΔt2k22k(2k-1)!·λ~ik|≤1.According to the property of Fourier transform, spatial wave number and derivative can be expressed as(38)kx=2Δxpkz=2Δzq,FFT(∂∂x)=j2ΔxpFFT(∂∂z)=j2Δzq,FFT(∂2∂x2)=-(2Δxp)2FFT(∂2∂z2)=-(2Δzq)2,
where p=∑m=1MCm(M)sin[Δx(m-0.5)kx] and q=∑m=1MCm(M)sin[Δz(m-0.5)kz]. FFT(·) represents Fourier transform.Solve (36) with (38); first two maximum eigenvalues can be written successively as
(39)λ~1=-4vp2(p2Δx2+q2Δz2),λ~2=-4vs2(p2Δx2+q2Δz2).By substituting (39) into (36), the dispersion relation can be derived as follows:
(40)sin2(ωΔt2)+∑k=1KΔt2k22k(2k-1)!·λ~1k=0,sin2(ωΔt2)+∑k=1KΔt2k22k(2k-1)!·λ~2k=0.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 (kx=π/Δx and kz=π/Δz) into (38) yields
(41)p=q=∑m=1MCm(M)(-1)m-1.Finally, we can derive the implicit stability condition as follows with (37), (39), and (41):
(42)∑k=1KΔt2k{vp2(∑m=1MCm(M))2(1/Δx2+1/Δz2)}k(2k-1)!≤1.Given the discretized parameters of spatial step, medium constant, and difference order, the maximum allowed temporal stepΔt can be computed with (42). For the second-order temporal difference, the corresponding stability condition is as follows: Δt≤(ΔxΔz/vp∑m=1M|Cm(M)|)1/(Δx2+Δz2); let variable ζ=∑m=1M|Cm(M)|; 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 (Figure2(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/cm3, respectively. Temporal step is 0.5 ms and spatial grid size is 5 m in both x- and z-direction. We locate a receiver at 200 m, 200 m. The computational results are compared with the analytical solution.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.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Figure2 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.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.
(a)
(b)
(c)Figures4(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.Snapshots of horizontal component with Figure2(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.
(a)
(b)
(c)
(d)


### 6.2. Planar Horizontal Model

A model with two layers is tested (Figure5(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/cm3. 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 x- and z-direction.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.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Figure5 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 (Figure6(a)) is tested. The truncated domain size is 3.4 km in x-direction and 1.4 km in z-direction. P-velocity ranges from 1028 m/s to 4670 m/s; Poisson ratio and density are 0.25 and 1.0 g/cm3, 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 x- and z-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.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.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Figure6 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.Figure7 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.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.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)


## 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 (Figure2(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/cm3, respectively. Temporal step is 0.5 ms and spatial grid size is 5 m in both x- and z-direction. We locate a receiver at 200 m, 200 m. The computational results are compared with the analytical solution.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.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Figure2 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.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.
(a)
(b)
(c)Figures4(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.Snapshots of horizontal component with Figure2(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.
(a)
(b)
(c)
(d)


## 6.2. Planar Horizontal Model

A model with two layers is tested (Figure5(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/cm3. 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 x- and z-direction.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.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Figure5 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 (Figure6(a)) is tested. The truncated domain size is 3.4 km in x-direction and 1.4 km in z-direction. P-velocity ranges from 1028 m/s to 4670 m/s; Poisson ratio and density are 0.25 and 1.0 g/cm3, 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 x- and z-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.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.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Figure6 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.Figure7 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.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.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)


## 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 x- and z-direction. P-velocity is 2000 m/s; Poisson ratio and density are 0.25 and 1.0 g/cm3, 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.
Methods
Zhang’s method
Li’s method
Proposed method
CPU computational time
4260 s
6614 s
3923 s
Numbers of variables
11
13
10
Numbers of equations
9
11
8We 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 Table2, 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.
Spatial difference order
Fourth-order
Tenth-order
CPU computational time
16534 s
3923 s
Spatial grid size
2.5 m
5 m
Memory cost
103.8 MB
27.5 MB


## 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.



---

*Source: 108713-2014-05-18.xml*