Journal of Computational Methods in Physics

Volume 2014 (2014), Article ID 108713, 14 pages

http://dx.doi.org/10.1155/2014/108713

Research Article

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

Ke-Yang Chen

Exploration and Development Research Institute of Daqing Oilfield Company Limited, Daqing, Heilongjiang 163712, China

Received 19 January 2014; Revised 15 April 2014; Accepted 15 April 2014; Published 18 May 2014

Academic Editor: Ivan D. Rukhlenko

Copyright © 2014 Ke-Yang Chen. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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:

Rewriting (1) as

where is particle displacement vector in Cartesian coordinate , is density, and are Lames’ coefficient, is P-wave velocity (), and is S-wave velocity ().

By introducing new variables

(2) becomes

It is easy to prove the following relations [3]:

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:

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:

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:

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]

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:

Hence, (7) can be rewritten as

where , . 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

where is the Fourier transformation of with damping term 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 ; 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 (), there is no reflection at the interface; it means the PML region () 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 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:

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:

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:

We represent (23) with matrix form

where 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

Substituting (26) into (27) yields

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

is unit matrix. Complex matrixes and are the Fourier transform of the matrix and the vector , respectively, where

The Frequency domain expression of matrix is , where

Because is hyperbolic equation system, there must be a nonsingular matrix satisfying the following relation:

where matrix is diagonal matrix containing six eigenvalues .

Substitute (32) into (29), and use similar matrix transform property:

Equation (29) becomes

Equation (34) can be rewritten as

We let . 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

In order to maintain stability of the numerical iteration computation, eigenvalues must satisfy the following relation:

According to the property of Fourier transform, spatial wave number and derivative can be expressed as

where and . represents Fourier transform.

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 S-wave’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: ; let variable ; we can compute the second- and fourth-order spatial difference’s variable and the second- and tenth-order spatial difference’s variable .

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.

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
8

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.

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.

Appendix

Previous Works

(a) The traditional elastic wave equation can be represented as

where and are normal stress and is shear stress.

(b) Li et al. [5] equivalent elastic wave equation can be written as

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:

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.

References

K. Y. Chen, “First-order velocity-stress elastic wave field separation scheme for biot two-phase isotropic medium,” Chinese Journal of Computational Physics, vol. 28, no. 3, pp. 404–412, 2011.

J. Dellinger and J. Etgen, “Wave-field separation in two-dimensional anisotropic media,” Geophysics, vol. 55, no. 7, pp. 914–919, 1990.

D. Ma and G. Zhu, “Numerical modeling of P-wave and S-wave separation in elastic wavefield,” Oil Geophysical Prospecting, vol. 38, no. 5, pp. 482–486, 2003.

J.-L. Zhang, Z.-P. Tian, and C.-X. Wang, “P- and S-wave separated elastic wave equation numerical modeling using 2D staggered-grid,” in Proceedings of the SEG/San Antonio Annual Meeting, pp. 2104–2109, 2007.

Z.-C. Li, H. Zhang, Q.-M. Liu, and W.-G. Han, “Numeric simulation of elastic wavefield separation by staggering grid high-order finite-difference algorithm,” Oil Geophysical Prospecting, vol. 42, no. 5, pp. 510–515, 2007.

L.-G. Dong, Z.-T. Ma, J.-Z. Cao et al., “A staggered-grid high-order difference method of one-order elastic wave equation,” Chinese Journal of Geophysics, vol. 43, no. 3, pp. 411–419, 2000.

L.-G. Dong, Z.-T. Ma, and J.-Z. Cao, “A study on stability of the staggered-grid high-order difference method of first-order elastic wave equation,” Chinese Journal of Geophysics, vol. 43, no. 6, pp. 856–864, 2000.

J.-P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” Journal of Computational Physics, vol. 114, no. 2, pp. 185–200, 1994.

K. Y. Chen, “Study on perfectly matched layer absorbing boundary condition,” Geophysical Prospecting for Petroleum, vol. 49, no. 5, pp. 472–477, 2010.

K. Y. Chen, “Anisotropic analysis of difference approximation in seismic wave numerical modeling,” Geophysical Prospecting for Petroleum, vol. 49, no. 1, pp. 19–22, 2010.

W. C. Chew and Q. H. Liu, “Perfectly matched layers for elastodynamics: a new absorbing boundary condition,” Journal of Computational Acoustics, vol. 4, no. 4, pp. 341–359, 1996.

J.-P. Berenger, “Three-dimensional perfectly matched layer for the absorption of electromagnetic waves,” Journal of Computational Physics, vol. 127, no. 2, pp. 363–379, 1996.

A. R. Mitchell, Computational Methods in Partial Differential Equations, John Wiley & Sons, New York, NY, USA, 1969.

D. M. Boore, “Finite-difference methods for seismic wave propagation in heterogeneous materials,” in Methods in Computational Physics, pp. 21–22, Academic Press, New York, NY, USA, 1972.

P. Moczo, J. O. A. Robertsson, and L. Eisner, “The finite-difference time-domain method for modeling of seismic wave propagation,” Advances in Geophysics, vol. 48, pp. 421–516, 2007.

F. Collino and C. Tsogka, “Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media,” Geophysics, vol. 66, no. 1, pp. 294–307, 2001.

D. Komatitsch and R. Martin, “An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation,” Geophysics, vol. 72, no. 5, pp. SM155–SM167, 2007.

C. Marcinkovich and K. Olsen, “On the implementation of perfectly matched layers in a three-dimensional fourth-order velocity-stress finite difference scheme,” Journal of Geophysical Research B: Solid Earth, vol. 108, no. 5, pp. 2276–2292, 2003.