Tilted Transversely Isotropic Research Articles

Optimal rotated staggered-grid finite-difference schemes for elastic wave modeling in TTI media

The rotated staggered-grid finite-difference (RSFD) is an effective approach for numerical modeling to study the wavefield characteristics in tilted transversely isotropic (TTI) media. But it surfaces from serious numerical dispersion, which directly affects the modeling accuracy. In this paper, we propose two different optimal RSFD schemes based on the sampling approximation (SA) method and the least-squares (LS) method respectively to overcome this problem. We first briefly introduce the RSFD theory, based on which we respectively derive the SA-based RSFD scheme and the LS-based RSFD scheme. Then different forms of analysis are used to compare the SA-based RSFD scheme and the LS-based RSFD scheme with the conventional RSFD scheme, which is based on the Taylor-series expansion (TE) method. The contrast in numerical accuracy analysis verifies the greater accuracy of the two proposed optimal schemes, and indicates that these schemes can effectively widen the wavenumber range with great accuracy compared with the TE-based RSFD scheme. Further comparisons between these two optimal schemes show that at small wavenumbers, the SA-based RSFD scheme performs better, while at large wavenumbers, the LS-based RSFD scheme leads to a smaller error. Finally, the modeling results demonstrate that for the same operator length, the SA-based RSFD scheme and the LS-based RSFD scheme can achieve greater accuracy than the TE-based RSFD scheme, while for the same accuracy, the optimal schemes can adopt shorter difference operators to save computing time.

Alternative stable qP wave equations in TTI media with their applications for reverse time migration

Numerical instabilities may arise if the spatial variation of symmetry axis is handled improperly when implementing P-wave modeling and reverse time migration in heterogeneous tilted transversely isotropic (TTI) media, especially in the cases where fast changes exist in TTI symmetry axis' directions. Based on the pseudo-acoustic approximation to anisotropic elastic wave equations in Cartesian coordinates, alternative second order qP (quasi-P) wave equations in TTI media are derived in this paper. Compared with conventional stable qP wave equations, the proposed equations written in stress components contain only spatial derivatives of wavefield variables (stress components) and are free from spatial derivatives involving media parameters. These lead to an easy and efficient implementation for stable P-wave modeling and imaging. Numerical experiments demonstrate the stability and computational efficiency of the presented equations in complex TTI media.

Reducing seismic positioning uncertainty with “simple” tilted transverse isotropy earth models — Gulf of Mexico Lobster field case study

The remaining potential in many mature fields in the Gulf of Mexico consists of smaller long-reach targets that require ever-more-accurate vertical and horizontal location for drilling precision. This is particularly true at the Lobster field in Ewing Bank, where, during the previous platform drilling campaigns, some midreach prospects were interpreted to have been missed due to mispositioning. Potential future drilling at the Lobster field would target primarily long-reach prospects in areas of significant geologic dip — increasing the risk of missing the prospects. This mispositioning risk could sometimes exceed the geologic risk, and it could have significant impact on the risked reserves and on the economics. Achieving the desired positioning accuracy requires the use of 3D seismic and anisotropic prestack depth migration with accurate velocity and anisotropy models. We have developed a case study of using a long offset walkaway checkshot and a geologic frame of reference to build a tilted transverse isotropy velocity model sufficiently accurate to support potential future drilling campaigns at the Lobster field.

Imaging complex geology through challenging surface terrain – a case study from West China

Huiwen Xie, Chao Wu, Caiming Luo, WenSheng Duan,Yingjie Zhao, Xin Ji, Sherman Yang, Xiaolong Yao, Dingxue Wang, Zhengguo Jia and Qinglin Liu describe a reprocessing project combining four 3D datasets that improved the imaging of a deep subsalt thrustfaulted target in an area of variable topography and complex near-surface conditions. Oil and gas exploration is moving into areas of increasingly complex geology, where depth imaging – along with the associated need to build accurate velocity models – is becoming an ever more important tool in efficiently finding and developing oilfields. This article presents a case study of structural imaging from an area where both the subsurface geology and surface conditions are extremely complex. The project combined several vintages of land 3D seismic data from the Dabei block, located in the Kuche foreland basin area of the northern Tarim Basin, West China. Several previous attempts had been made to reprocess these datasets, but the results were all considered inadequate to meet the imaging quality and accuracy required for reliable well placement. The goal of the new reprocessing project was to deliver accurate imaging of hydrocarbon-bearing strata beneath a large anticline in subsalt thrust-fault blocks at depths of around 7000 m. Tomographic static correction using carefully picked first breaks outperformed previous field static solutions based on upholes. A multi-domain noise attenuation sequence included a newly developed method that uses true shot and receiver coordinates, thereby handling irregularities in acquisition geometry. A tilted transverse isotropic (TTI) velocity model was built using a combination of diving wave tomography (DWT) in the shallow part; geologically constrained reflection tomography in the remainder of the overburden; and multi-azimuthal tomography in the offset vector tile (OVT) domain in the salt and target zone. In areas of thick high-velocity conglomerate deposits, a seismic joint inversion (SJI) model building technique was applied to include information from non-seismic data to better define strong lateral velocity changes.

Accurate simulations of pure quasi-P-waves in complex anisotropic media

Reverse time migration (RTM) in complex anisotropic media requires calculation of the propagation of a single-mode wave, the quasi-P-wave. This was conventionally realized by solving a [Formula: see text] system of second-order partial differential equations. The implementation of this [Formula: see text] system required at least twice the computational resources as compared with the acoustic wave equation. The S-waves, an introduced auxiliary function in this system, were treated as artifacts in the RTM. Furthermore, the [Formula: see text] system suffered numerical stability problems at the places in which abrupt changes of symmetric axis of anisotropy exist, which brings more challenges to real data implementation. On the other hand, the Alkhalifah’s equation, which governs the pure quasi-P-wave propagation, was hard to solve because it was a pseudodifferential equation. We proposed a pure quasi-P-wave equation that can be easily implemented within current imaging framework. Our new equation was obtained by decomposing the original pseudodifferential operator into two numerical solvable operators: a differential operator and a scalar operator. The combination of these two operators yielded an accurate phase of quasi-P-wave propagation. Our solution was S-wave free and numerically stable for very complicated models. Because only one equation was required to resolve numerically, the new proposed scheme was more efficient than those conventional schemes that solve the [Formula: see text] second-order differential equations system. For tilted transverse isotropy (TTI) RTM implementation, the required increase of numerical cost was minimal, and we could expect at least a factor of two of improvement of efficiency. We showed the effectiveness and robustness of our method with numerical examples with simple and very complicated TTI models, the SEG Advanced Modeling (SEAM) model. Further extension of our approach to orthorhombic anisotropy or tilted orthorhombic anisotropy was straightforward.

Formulation of the Riccati equation for the impedance operator in cylindrical coordinates for inhomogeneous anisotropic waveguides with the example of rectilinear anisotropy

The analysis of cylindrical waveguides with arbitrary anisotropy and inhomogeneities is more difficult than for isotropic or cylindrically anisotropic. The reason for this difficulty is because the independent consideration of different Fourier components of the wavefield is no longer possible. Taking into account the coupling between different Fourier components of the wavefield implies that the impedance tensor should be treated as the operator. Treating the impedance this way is the key difference between the isotropic and the cylindrical anisotropy case in which the independent consideration of the partial impedance matrices for each circumferential number n is possible. The natural approach to computing the impedance operator is to use the Riccati equation, which is presented by the authors in cylindrical coordinates. The matrix representations of the impedance operator and the Riccati equation over the basis of the time-harmonic cylindrical waves result in the infinite set of ordinary differential equations. The nondiagonal partial impedance matrices, which describe the coupling of different Fourier components, become nonzero. For the problems, which possess a symmetry, the formulation can be simplified by projecting the Fourier basis on that of the irreducible representations of the symmetry group. The numerical solution of the Riccati equation and its application to finding the dispersion curves of the eigenmodes of the waveguide have peculiarities due to the fact that the matrix representation is infinite. The feasibility of the proposed approach is demonstrated by computing the dispersion curves for the first monopole and dipole normal modes of a cylindrical waveguide in a homogeneous media possessing tilted transverse isotropy type of anisotropy.

Common-reflection-point migration velocity analysis of 2D P-wave data from TTI media

We have developed a common-reflection-point (CRP)-based kinematic migration velocity analysis for 2D P-wave reflection data to estimate the four transversely isotropic (TI) parameters [Formula: see text], [Formula: see text], and [Formula: see text], and the tilt angle [Formula: see text] of the symmetry axis in a TI medium. In each iteration, the tomographic parameter was updated alternately with prestack anisotropic ray-based migration. Iterations initially used layer stripping to reduce the number of degrees of freedom; after convergence was reached, a couple of more iterations over all parameters and all CRPs ensured global interlayer coupling and parameter interaction. The TI symmetry axis orientation was constrained to be locally perpendicular to the reflectors. The [Formula: see text] dominated the inversion, and so it was weighted less than [Formula: see text] and [Formula: see text] in the parameter updates. Estimates of [Formula: see text] and [Formula: see text] were influenced if the error in [Formula: see text] was [Formula: see text]; estimates of [Formula: see text] were also influenced if the error in [Formula: see text] was [Formula: see text]. Examples included data for a simple model with a homogeneous TI layer whose dips allowed recovery of all anisotropy parameters from noise-free data, and a more realistic model (the BP tilted transversely isotropic (TTI) model) for which only [Formula: see text], [Formula: see text], and [Formula: see text] were recoverable. The adequacy of the traveltimes predicted by the inverted anisotropic models was tested by comparing migrated images and common image gathers, with those produced using the known velocity models.

Performance of CPU/GPU compiler directives on ISO/TTI kernels

GPUs are slowly becoming ubiquitous devices in High Performance Computing, as their capabilities to enhance the performance per watt of compute intensive algorithms as compared to multicore CPUs have been identified. The primary shortcoming of a GPU is usability, since vendor specific APIs are quite different from existing programming languages, and it requires a substantial knowledge of the device and programming interface to optimize applications. Hence, lately a growing number of higher level programming models are targeting GPUs to alleviate this problem. The ultimate goal for a high-level model is to expose an easy-to-use interface for the user to offload compute intensive portions of code (kernels) to the GPU, and tune the code according to the target accelerator to maximize overall performance with a reduced development effort. In this paper, we share our experiences of three of the notable high-level directive based GPU programming models--PGI, CAPS and OpenACC (from CAPS and PGI) on an Nvidia M2090 GPU. We analyze their performance and programmability against Isotropic (ISO)/Tilted Transversely Isotropic (TTI) finite difference kernels, which are primary components in the Reverse Time Migration (RTM) application used by oil and gas exploration for seismic imaging of the sub-surface. When ported to a single GPU using the mentioned directives, we observe an average 1.5---1.8x improvement in performance for both ISO and TTI kernels, when compared with optimized multi-threaded CPU implementations using OpenMP.

The impact of 3D tilted resistivity anisotropy on marine CSEM measurements

Marine CSEM modeling and data-interpretation algorithms currently used by the E&P industry are based on application of vertical transverse isotropic (VTI) resistivity models, which in many practical cases cannot adequately describe complex tilted anisotropic geologic structures typical for accumulation of hydrocarbons (e.g., syncline or anticline types). Use of VTI models is a serious limitation of CSEM technology and might be one of the core reasons for some false-positive and false-negative outcomes of marine CSEM surveys. To investigate this phenomenon in more detail and assess the impact of arbitrary tilted resistivity anisotropy on CSEM data, we performed simulations for a set of 3D tilted transverse isotropic (TTI) and VTI benchmark models. The 3D modeling results show that ignoring the TTI effect, if it is present, might lead to incorrect assessment of the feasibility of a potential CSEM survey as well as erroneous assessments and/or mispositioning of hydrocarbon targets during interpretation of field data and, in turn, might compromise the reliability of marine CSEM technology. Therefore, there is a strong practical need to start using TTI models, which provide more accurate descriptions of complex anisotropic resistivity structures than VTI models do. Application of 3D TTI modeling will improve the reliability of CSEM technology and its derisking potential.

Ray tracing of multiple transmitted/reflected/converted waves in 2-D/3-D layered anisotropic TTI media and application to crosswell traveltime tomography

To overcome the deficiency of some current grid-/cell-based ray tracing algorithms, which are only able to handle first arrivals or primary reflections (or conversions) in anisotropic media, we have extended the functionality of the multistage irregular shortest-path method to 2-D/3-D tilted transversely isotropic (TTI) media. The new approach is able to track multiple transmitted/reflected/converted arrivals composed of any kind of combinations of transmissions, reflections and mode conversions. The basic principle is that the seven parameters (five elastic parameters plus two polar angles defining the tilt of the symmetry axis) of the TTI media are sampled at primary nodes, and the group velocity values at secondary nodes are obtained by tri-linear interpolation of the primary nodes across each cell, from which the group velocities of the three wave modes (qP, qSV and qSH) are calculated. Finally, we conduct grid-/cell-based wave front expansion to trace multiple transmitted/reflected/converted arrivals from one region to the next. The results of calculations in uniform anisotropic media indicate that the numerical results agree with the analytical solutions except in directions of SV-wave triplications, at which only the lowest velocity value is selected at the singularity points by the multistage irregular shortest-path anisotropic ray tracing method. This verifies the accuracy of the methodology. Several simulation results show that the new method is able to efficiently and accurately approximate situations involving continuous velocity variations and undulating discontinuities, and that it is suitable for any combination of multiple transmitted/reflected/converted arrival tracking in TTI media of arbitrary strength and tilt. Crosshole synthetic traveltime tomographic tests have been performed, which highlight the importance of using such code when the medium is distinctly anisotropic.

Seismic wavefield simulation in 2D elastic and viscoelastic tilted transversely isotropic media: comparisons between four different kinds of finite-difference grid schemes

In this paper, we use the staggered grid, the auxiliary grid, the rotated staggered grid and the non-staggered grid finite-difference methods to simulate the wavefield propagation in 2D elastic tilted transversely isotropic (TTI) and viscoelastic TTI media, respectively. Under the stability conditions, we choose different spatial and temporal intervals to get wavefront snapshots and synthetic seismograms to compare the four algorithms in terms of computational accuracy, CPU time, phase shift, frequency dispersion and amplitude preservation. The numerical results show that: (1) the rotated staggered grid scheme has the least memory cost and the fastest running speed; (2) the non-staggered grid scheme has the highest computational accuracy and least phase shift; (3) the staggered grid has less frequency dispersion even when the spatial interval becomes larger.

SEAM Update

SEAM Phase I-RPSEA is nearing the completion of its work. The project was initiated with the goal of building and conducting simulations on an acoustic model for a deep-water resource that is based on what might be found in the U.S. Gulf of Mexico. With support from the Research Partnership to Secure Energy for America (RPSEA) (Subcontract no. 07121-2001 under US DOE prime contract no. DE-AC 26-07NT42677) and the SEAM Phase I industrial participants, other geophysical parameters have been subsequently defined for the model and several additional surveys have been simulated on the model. We are now reaching the conclusion of our work. Surveys conducted to date include: variable density acoustic, variable density acoustic with an absorbing upper surface, gravity, controlled-source electromagnetic, tilted transverse isotropy acoustic, elastic, and magnetotelluric.

An efficient hybrid pseudospectral/finite-difference scheme for solving the TTI pure P-wave equation

The pure P-wave equation for modelling and migration in tilted transversely isotropic (TTI) media has attracted more and more attention in imaging seismic data with anisotropy. The desirable feature is that it is absolutely free of shear-wave artefacts and the consequent alleviation of numerical instabilities generally suffered by some systems of coupled equations. However, due to several forward–backward Fourier transforms in wavefield updating at each time step, the computational cost is significant, and thereby hampers its prevalence. We propose to use a hybrid pseudospectral (PS) and finite-difference (FD) scheme to solve the pure P-wave equation. In the hybrid solution, most of the cost-consuming wavenumber terms in the equation are replaced by inexpensive FD operators, which in turn accelerates the computation and reduces the computational cost. To demonstrate the benefit in cost saving of the new scheme, 2D and 3D reverse-time migration (RTM) examples using the hybrid solution to the pure P-wave equation are carried out, and respective runtimes are listed and compared. Numerical results show that the hybrid strategy demands less computation time and is faster than using the PS method alone. Furthermore, this new TTI RTM algorithm with the hybrid method is computationally less expensive than that with the FD solution to conventional TTI coupled equations.

Generalized migration in frequency-wavenumber domain (MGF-K) in anisotropic media

In this paper, the background of MGF-K migration in dual domain (wavenumber-frequency K-F and space-time) in anisotropic media is presented. Algorithms for poststack (zero-offset) and prestack migration are based on downward extrapolation of acoustic wavefield by shift-phase with correction filter for lateral variability of medium’s parameters. In anisotropic media, the vertical wavenumber was determined from full elastic wavefield equations for two dimensional (2D) tilted transverse isotropy (TTI) model. The method was tested on a synthetic wavefield for TTI anticlinal model (zero-offset section) and on strongly inhomogeneous vertical transverse isotropy (VTI) Marmousi model. In both cases, the proper imaging of assumed media was obtained.

Anisotropic parsimonious prestack depth migration

An efficient Kirchhoff-style prestack depth migration, called “parsimonious” migration, was developed a decade ago for isotropic 2D and 3D media by using measured slownesses to reduce the amount of ray tracing by orders of magnitude. It is conceptually similar to “map” migration, but its implementation has some differences. We have extended this approach to 2D tilted transversely isotropic (TTI) media and illustrated it with synthetic P-wave data. Although the framework of isotropic parsimonious may be retained, the extension to TTI media requires redevelopment of each of the numerical components, calculation of the phase and group velocity for TTI media, development of a new two-point anisotropic ray tracer, and substitution of an initial-angle isotropic shooting ray-trace algorithm for an anisotropic one. The model parameterization consists of Thomsen’s parameters ([Formula: see text], [Formula: see text], [Formula: see text]) and the tilt angle of the symmetry axis of the TI medium. The parsimonious anisotropic migration algorithm is successfully applied to synthetic data from a TTI version of the Marmousi2 model. The quality of the image improves by weighting the impulse response by the calculation of the anisotropic Fresnel radius. The accuracy and speed of this migration makes it useful for anisotropic velocity model building. The elapsed computing time for 101 shots for the Marmousi2 TTI model is 35 s per shot (each with 501 traces) in 32 Opteron cores.

Wavefield propagation characteristics in fracture-induced TTI double-porosity medium

The direction of symmetry axis of parallel fracture set in fractured hydrocarbon reservoir affects the transmission of seismic waves markedly, so a medium named fracture-induced TTI (tilted transverse isotropy) double-porosity medium is studied here to discuss the effect of different dip and azimuth angles of a fracture system. Based on the theories of fracture-induced HTI (horizental transverse isotropy) double-porosity medium, the softness and dispersion matrixes of fracture-induced TTI double-porosity medium are derived with the application of Bond transform, and finally, single-order velocity-stress equations are obtained. Furthermore, numerical simulations in xoz plane of 2.5 dimensional vector wavefield are carried out by the method of high-order staggered-grid finite-difference under perfect matched layer (PML) boundary conditions. The results show that the dipand azimuth angles of fractures have great impacts on seismic wave propagation, since the angles can cause the phenomena of shear wave splitting and, in the two-layer model of fracture-induced TTI double-porosity, converted shear wave splitting and shear wave sub-splitting. All of these increase the complexity of seismic wavefield and will lay a foundation of further studies on seismic wave propagation in actual earth layers.

Acceleration of stable TTI P-wave reverse-time migration with GPUs

When a pseudo-acoustic TTI (tilted transversely isotropic) coupled wave equation is used to implement reverse-time migration (RTM), shear wave energy is significantly included in the migration image. Because anisotropy has intrinsic elastic characteristics, coupling P-wave and S-wave modes in the pseudo-acoustic wave equation is inevitable. In RTM with only primary energy or the P-wave mode in seismic data, the S-wave energy is regarded as noise for the migration image. To solve this problem, we derive a pure P-wave equation for TTI media that excludes the S-wave energy. Additionally, we apply the rapid expansion method (REM) based on a Chebyshev expansion and a pseudo-spectral method (PSM) to calculate spatial derivatives in the wave equation. When REM is incorporated with the PSM for the spatial derivatives, wavefields with high numerical accuracy can be obtained without grid dispersion when performing numerical wave modeling. Another problem in the implementation of TTI RTM is that wavefields in an area with high gradients of dip or azimuth angles can be blown up in the progression of the forward and backward algorithms of the RTM. We stabilize the wavefields by applying a spatial-frequency domain high-cut filter when calculating the spatial derivatives using the PSM. In addition, to increase performance speed, the graphic processing unit (GPU) architecture is used instead of traditional CPU architecture. To confirm the degree of acceleration compared to the CPU version on our RTM, we then analyze the performance measurements according to the number of GPUs employed.

On the instability in second-order systems for acoustic VTI and TTI media

We studied second-order wave propagation systems for vertical transversely isotropic (VTI) and tilted transversely isotropic (TTI) acoustic media with variable axes of symmetry that have their shear-wave speeds set to zero. Acoustic TTI systems are commonly used in reverse-time migration, but these second-order systems are susceptible to instablities appearing as nonphysical stationary noise growing linearly in time, particularly in variable-tilt TTI media. We found an explanation of the cause of this phenomenon. The instabilities are not caused only by the numerical schemes; they are inherent to the differential equations. These instabilities are present even in homogeneous VTI media. These instabilities are caused by zero wave speeds at a wide variety of wavenumbers — a direct consequence of setting the shear-wave speeds to zero — coupled with the second time derivative in these systems. Although the second-order isotropic wave equation allows smooth time-growing solutions, a larger class of time-growing solutions exists for the second-order acoustic TI systems, including nonsmooth solutions. Boundary conditions appear to be less effective in controlling these time-growing solutions than they are for the isotropic wave equation. These systems conserve an incomplete energy that does not prevent the instabilities. The corresponding steady-state systems are no longer elliptic differential equations and can have nonsmooth solutions that are related to the instabilities. We started initially with homogeneous VTI media, and then extended these results to heterogeneous variable-tilt TTI media. We also developed a second-order acoustic system for heterogeneous variable-tilt TTI media derived directly from the full-elastic system for heterogeneous variable-tilt TTI media. All second-order systems with a dispersion relation obtained by setting the shear-wave speeds to zero in the elastic dispersion relation allowed these nonphysical time-growing solutions; however, knowing the cause of these instabilities, it may be possible to prevent or control the activation of these solutions.

First-order systems for elastic and acoustic variable-tilt TI media

We derive and compare first-order wave propagation systems for variable-tilt elastic and acoustic tilted transversely isotropic (TTI) media. Acoustic TTI systems are commonly used in reverse-time migration. Starting initially with homogeneous vertical transversely isotropic (VTI) media, and then extending to heterogeneous variable-tilt TI media, we derive a pseudoacoustic [Formula: see text] first-order system of differential equations by setting the shear-wave speeds to zero and simplifying the full-elastic system accordingly. This [Formula: see text] system conserves a complete energy, but only when the anelliptic anisotropy parameter [Formula: see text]. For [Formula: see text] (including isotropic media), the system allows linearly time-growing and spatially nonpropagating nonphysical solutions frequently taken for numerical noise. We modified this [Formula: see text] acoustic first-order system by changing the stress variables to obtain a system that stays stable for [Formula: see text]. This system for homogeneous VTI media is generalized to heterogeneous variable-tilt TI media by rotating the stress and strain variables in the full elastic system before setting the shear-wave speeds to zero; the system obtained can be greatly simplified by combining the rotational terms, resulting in only one rotation and extra lower-order terms compared to the [Formula: see text] first-order acoustic system for VTI media. This new system can be simplified further by neglecting the lower-order terms. Both systems (with and without lower-order terms) conserve the same complete energy. Finally, the corresponding [Formula: see text] full elastic system for variable-tilt acoustic TI media can be used for the purposes of benchmarking.

Elastic modelling in tilted transversely isotropic media with convolutional perfectly matched layer boundary conditions

To simulate wave propagation in a tilted transversely isotropic (TTI) medium with a tilting symmetry-axis of anisotropy, we develop a 2D elastic forward modelling algorithm. In this algorithm, we use the staggered-grid finite-difference method which has fourth-order accuracy in space and second-order accuracy in time. Since velocity-stress formulations are defined for staggered grids, we include auxiliary grid points in the z-direction to meet the free surface boundary conditions for shear stress. Through comparisons of displacements obtained from our algorithm, not only with analytical solutions but also with finite element solutions, we are able to validate that the free surface conditions operate appropriately and elastic waves propagate correctly. In order to handle the artificial boundary reflections efficiently, we also implement convolutional perfectly matched layer (CPML) absorbing boundaries in our algorithm. The CPML sufficiently attenuates energy at the grazing incidence by modifying the damping profile of the PML boundary. Numerical experiments indicate that the algorithm accurately expresses elastic wave propagation in the TTI medium. At the free surface, the numerical results show good agreement with analytical solutions not only for body waves but also for the Rayleigh wave which has strong amplitude along the surface. In addition, we demonstrate the efficiency of CPML for a homogeneous TI medium and a dipping layered model. Only using 10 grid points to the CPML regions, the artificial reflections are successfully suppressed and the energy of the boundary reflection back into the effective modelling area is significantly decayed.