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

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 Pand Smodes, and apply it in full elastic wave numerical simulation. We give the complete derivations of explicit high-order staggeredgrid 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 Pand 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.


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 Swave 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 Pand 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 staggeredgrid 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 2 Journal of Computational Methods in Physics 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.
By introducing new variables (2) becomes It is easy to prove the following relations [3]: (5) Equation (5)  For 2D case, (4) will be simplified to the following equations: 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 /, 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: where V  and V  are particle velocity fields, V  and V  are pure P-and S-wave fields in -direction, and V  and V  are pure P-and S-wave fields in -direction. and  are auxiliary variables defined above.We can further derive the following relations:  where  ⇀  = (V  , V  )  .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 Pand S-modes.

Numerical Discretion
For typical grid cell (, ), take the first-order spatial and temporary derivatives in (7); we have [6,7] where  and  denote particle velocity fields and auxiliary variables, Δ and Δ are spatial steps, Δ is temporary step,  ()  denotes coefficients of -order staggered-grid finitedifference scheme,  denotes temporary grid, and  and  are coordinate systems illustrated in Figure 1.
Where P-velocity  is located in the integral grid, but Svelocity  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] Journal of Computational Methods in Physics where . ., .The finite-difference order  in the paper is 5, which can guarantee perfect numerical accuracy [15].

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 where F is the Fourier transformation of  with damping term exp[−  ⇀  ⋅  ⇀  ] and F 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: This implies the complex transformation of Cartesian variables: Or equivalently, upon difference: Therefore, the plane-wave solutions with PML in the direction can be rewritten as 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 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: where  and  denote particle velocity fields or auxiliary variables. and  denote variable such as P-velocity or Svelocity.
Following the PML splitting equations discussed in the conventional PML methods [18], the local coordinatesystem-based PML splitting equations are formulated in a local Cartesian coordinate system with the axis being normal to the local artificial boundary: 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 Equation ( 20) can be discretized with finite-difference scheme as follows: Finally, the resulting discretized formulations of the PML equation are derived:

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: We represent (23) with matrix form where  = (V , V , V , V , , )  and  is the matrix containing elastic constants and spatial difference in (23).Equation (24) indicates the following relation: Substitute the spatial derivatives in (25) with ( 9), becomes: where matrix  is finite-difference discretized formula of matrix .
We use Taylor series expansion technique and then (24) can be written as (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.
Solve (36) with (38); first two maximum eigenvalues can be written successively as By substituting (39) into (36), the dispersion relation can be derived as follows: According to (39) and (40), we can notice the upper part of ( 40) is P-wave's dispersion relation, the lower part is Swave's dispersion relation, and they are completely separated.
Substituting Nyquist wave number (  = /Δ and   = /Δ) into (38) yields Finally, we can derive the implicit stability condition as follows with (37), (39), and (41): 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 |; 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.

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 3 , 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.
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 Swave 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 Swave field are separated completely, which is also shown in    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.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 Pwave (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.

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 3 , 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.
Figure 6 shows Marmousi model and snapshots of elastic wave fields of horizontal component and vertical component  Figure 7 shows the numerical records of our method.Figures 7(a  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 staggeredgrid finite-difference method and PML absorbing boundary condition, no dispersion and spurious boundary reflection phenomenon appear, and entire modeling process is stable.

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 3 , respectively.Spatial step is 5 m in each direction, and temporary step is 0.5 ms.Dominant 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 where   and   are normal stress and   is shear stress.(b) Li et al. [5] equivalent elastic wave equation can be written as

Figure 1 :
Figure 1: Local mesh and spatial models for velocity and auxiliary variables updates.

Figure 2 :
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 3 .
Figure 3.The full elastic wave field equals the sum of pure Pwave 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.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
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.Figure4shows 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.

Figure 3 :
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 Pand S-wave fields, and (c) denotes the pure P-wave field and S-wave field based on divergence and curl operators, respectively.

Figure 4 :
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.

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

2 DistanceFigure 5 :
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 6 :
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.

Table 1 :
Comparison with different separating methods.

Table 2 :
Comparison with different difference orders of the proposed method.can be extended to heterogeneous media.PML absorbing boundary condition and the numerical discretion with highorder 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. and [4]nd  , are normal stress of pure P-wave,  , and  , are normal stress of pure shear wave, and  , is shear stress of pure shear wave.(c)Jianleietal.[4]equivalent elastic wave equation can be expressed as follows: ,