Downloaded 09/19/16 to 106.68.109.41. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

GEOPHYSICS, VOL. 81, NO. 5 (SEPTEMBER-OCTOBER 2016); P. C265–C278, 7 FIGS. 10.1190/GEO2015-0311.1

Acoustic wave propagation in tilted transversely isotropic media: Incorporating topography

Jeffrey Shragge1

a domain boundary conformal to free-surface topography. A generalized curvilinear transformation is used to specify a system of equations governing 3D acoustic TTI wave propagation in the “topographic” coordinate system. The developed finitedifference time-domain numerical solution adapts existing Cartesian TTI operators to this more generalized geometry with little additional computational overhead. Numerical evaluations illustrate that 2D and 3D impulse responses are wellmatched to those simulated on Cartesian meshes and analytic traveltimes for homogeneous elliptical TTI media. Accordingly, these generalized acoustic TTI propagators and their numerical adjoints are useful for undertaking most RTM or FWI applications using computational domains conforming to free-surface topography.

ABSTRACT Simulating two-way acoustic wavefield propagation directly from a free-surface boundary in the presence of topography remains a computational challenge for applications of reverse time migration (RTM) or full-waveform inversion (FWI). For land-seismic settings involving heavily reworked geology (e.g., fold and thrust belts), two-way wavefield propagation operators should also handle commonly observed complex anisotropy including tilted transversely isotropic (TTI) media. To address these issues, I have extended a system of coupled partial differential equations used to model 3D acoustic TTI wave propagation in Cartesian coordinates to more generalized 3D geometries, including a deformed computational mesh with

equation migration (Ursenbach and Bale, 2009). Although partial-wavefield approaches (and similarly prestack time-migration methods) are less computationally intensive and have lower sensitivity to irregularities in shooting geometries than full-wavefield approaches (Vestrum et al., 2011), their main drawback is an inherent inability to model all TTI wavefield arrivals (e.g., wavefield multipathing, turning waves, and free-surface scattering), which can lead to degraded imaging and inversion performance. This observation helps to motivate the study of full-wavefield 3D TTI propagation operators able to handle TTI anisotropy in the presence of complex free-surface topography and strong lateral velocity contrast. Ideally, these modeling engines would form the main computational kernels for 3D (acoustic) reverse time migration (RTM) (Baysal et al., 1983; McMechan, 1983; Whitmore, 1983) or full-waveform inversion (FWI; Tarantola, 1984; Pratt et al., 1998) algorithms for the types of land seismic scenarios described above. Cartesian-based RTM using time-evolution solutions of the (twoway) acoustic wave equation has become an industry standard for

INTRODUCTION Obtaining high-quality 3D seismic images in areas of complex geology remains a challenge for the land seismic exploration community. Many geologic regions (e.g., the Canadian Foothills, Turkish Anatolia, Pakistani Foothills, Colombia, and Western China) exhibit irregular surface topography and heavily reworked geology that result in strong lateral velocity contrasts and moderate-to-strong tilted transversely isotropic (TTI) media. Generating high-quality prestack depth migration images in these complex scenarios could benefit from full-wavefield anisotropic seismic imaging (and corresponding model building) techniques capable of addressing these challenges. The past decade has seen the development of more accurate partial-wavefield TTI migration algorithms and their application in complex land imaging scenarios. Reported methods include TTI Kirchhoff migration (Vestrum, 2002), TTI Gaussian beam migration (Vestrum, 2002; Gittins and Vestrum, 2004), and TTI wave-

Manuscript received by the Editor 12 June 2015; revised manuscript received 1 April 2016; published online 02 September 2016. 1 The University of Western Australia, Centre for Energy Geoscience, School of Earth and Environment and School of Physics, Crawley, Australia. E-mail: [email protected]. © 2016 Society of Exploration Geophysicists. All rights reserved. C265

Downloaded 09/19/16 to 106.68.109.41. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

C266

Shragge

high-end 3D marine seismic imaging applications. Extensions of the two-way operators used in isotropic acoustic RTM to TTI media largely originate from the vertical transverse isotropy dispersion relationship developed in Alkhalifah (2000), which introduces a pseudoacoustic approximation by setting the S-wave velocity to zero along the symmetry axis. Fletcher et al. (2009) and Zhang and Li (2013) present alternate formulations of 3D acoustic TTI wave equations. In particular, Fletcher et al. (2009) use a P-SV dispersion relation based on plane-wave solutions of the Christoffel equation (Tsvankin, 2001) to derive an “asymptotically exact” pseudoacoustic TTI wave equation. For implementation and efficiency reasons, Fletcher et al. (2009) reduce the resulting fourth-order partial differential equation (PDE) into a system of coupled second-order PDEs, and provide a recipe for stabilizing the simulation of quasi-P waves by introducing a nonzero S-wave velocity component (that also can generate unwanted S-wave artifacts). Zhang et al. (2009) and Guan et al. (2011) present different strategies for attenuating unwanted S-wave artifacts based on eigenvalue analysis and physical operator-based filtering, respectively. Xu and Zhou (2014) provide an alternate approach for P-wave TTI propagation applicable in general TTI scenarios that is free of S-wave artifacts and demonstrably more efficient to compute. Bin Waheed and Alkhalifah (2015) present a method for calculating an effective tilted elliptically anisotropic model that, although has more complex coefficients, leads to a single TTI propagation equation whose numerical solutions are free of S-wave artifacts and more efficient to calculate. Although most of TTI RTM applications reported in the literature predominantly focus on marine applications (e.g., subsalt imaging), questions remain regarding how best to extend 3D TTI RTM to complex land geologic scenarios, such as those described above. Shragge (2014a) introduces a methodology for performing isotropic acoustic wave propagation directly from topographic surfaces (Lan et al. (2014) present a related approach). This approach uses a 2D/3D coordinate transformation between Cartesian and a “topographic” coordinate system (Tessmer et al., 1992; Hestholm, 1999) formed by smooth linear interpolation between an irregular topographic surface and a user-specified horizon at depth. Wavefield solutions, single shot-record images, or complete image volumes computed in the topographic computational domain are then interpolated back to Cartesian for interpretation or further analysis. This topographic mapping facilitates the implementation of free-surface boundary conditions because the computational grid is conformal to the physical domain boundary. However, for such an RTM or FWI strategy to be more generally applicable, it must be extended to handle scenarios involving anisotropic TTI media commonly found in complex areas, where these full-wavefield methods are more likely to make a significant impact. The theory presented here aims to address this issue by formulating an approach for performing 2D/3D acoustic TTI wave propagation directly from a topographic freesource boundary using a deformed coordinate system. The paper begins by reviewing the system of coupled second-order PDEs governing Cartesian acoustic TTI wave propagation presented in Fletcher et al. (2009), and extending this formulation to generalized 3D coordinate systems through a change-of-variables transformation of the underlying partial differential operators. I then specify 2D/3D topographic coordinates and develop the corresponding 2D/3D TTI equations appropriate for this geometry. I present several examples for 2D and 3D elliptic and anelliptic

TTI media and compare impulse responses with those simulated through Cartesian meshes. Finally, I conclude with a discussion on several computational issues involving the implementation of acoustic TTI propagation operators.

3D CARTESIAN TTI PROPAGATION An efficient way of modeling 3D Cartesian acoustic wave propagation in TTI media is through a set of coupled second-order PDEs (Fletcher et al., 2009). I begin by rewriting equation 2 of Fletcher et al. (2009) and likewise setting the parameter α ¼ 1 to yield the following system of equations:

∂2 P x ¼ V 2Px H 2 ½P x  þ V 2Pz H 1 ½Qx  þ V 2Sz H 1 ½P x − Qx ; (1) ∂t2 ∂2 Qx ¼ V 2Pn H 2 ½P x  þ V 2Pz H 1 ½Qx  − V 2Sz H 2 ½P x − Qx ; (2) ∂t2 where V Pz is the vertical P-wave velocity in the direction normal to the symmetry plane, V Px is the horizontal P-wave velocity in the symmetry plane, V Pn is the P-wave NMO velocity in the symmetry plane, V Sz is the vertical S-wave velocity in the symmetry plane, P x is the acoustic pressure field, and Qx is an auxiliary pressure field. Subscripts x indicate fields specified in a Cartesian coordinate system, x ¼ ½x1 ; x2 ; x3 . The various P-wave velocity fields are related to the more common Thomsen’s (1986) anisotropy parameters ϵ and δ through

pffiffiffiffiffiffiffiffiffiffiffiffiffi V Pn ¼ V Pz 1 þ 2δ

pffiffiffiffiffiffiffiffiffiffiffiffiffi and V Px ¼ V Pz 1 þ 2ϵ:

(3)

Partial differential operators H 1 ½U x  and H 2 ½U x  act on a (generic) acoustic pressure field U x according to

H 1 ½U x  ¼ aij

∂2 U x ; ∂xi ∂xj

H 2 ½U x  ¼ ½δij − aij 

∂2 U x ; ∂xi ∂xj

i; j ¼ 1; 2; 3;

(4)

where summation is implied over repeated indices, δij is a Kronecker delta function, and aij are trigonometric coefficients given by

2

2 2 1 4 2 cos ϕ sin2 θ ij ½a  ¼ sin 2ϕ sin θ 2 cos ϕ sin 2θ

sin 2ϕ sin2 θ 2 sin2 ϕ sin2 θ sin ϕ sin 2θ

3 cos ϕ sin 2θ sin ϕ sin 2θ 5: 2 cos2 θ (5)

Angle θ is the dip measured between the symmetry and vertical axes, and ϕ is the azimuth of the symmetry axis with respect to the horizontal plane. To facilitate summation notation, the coefficients are symmetrized with respect to those specified in Fletcher et al. (2009). Substituting for the H 1 and H2 operators in equations 1 and 2 yields the desired form of the acoustic TTI wave equations in 3D Cartesian coordinates:

Topographic 3D TTI propagation

∂2 P x ∂2 P x ¼ ½V 2Px ðδij − aij Þ þ V 2Sz aij  2 ∂xi ∂xj ∂t

Downloaded 09/19/16 to 106.68.109.41. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

þ ½ðV 2Pz − V 2Sz Þaij 

∂ 2 Qx ; ∂xi ∂xj

(6)

∂2 Qx : ∂xi ∂xj

(7)

Fletcher et al. (2009) discuss a methodology for solving this system of equations in a stable manner for scenarios, where ϵ < δ (by including nonzero V Sz ), and present wavefield propagation and RTM examples that demonstrate the utility of this modeling approach in 3D TTI media.

3D GENERALIZED TTI PROPAGATION This section demonstrates how to extend the above system of coupled PDEs describing 3D Cartesian TTI wave propagation to generalized 3D coordinate meshes. Here, I denote a generalized coordinate mesh using variables ξ ¼ ½ξ1 ; ξ2 ; ξ3  and assume that the mapping relationship to Cartesian coordinates, ξðxÞ ¼ f, is known and one-to-one. Specifying the corresponding system of coupled PDEs on a generalized 3D mesh requires transforming the coordinate basis of the second-order tensor formed by the partial derivatives acting on wavefield U x (e.g., ∂2 U x ∕∂xi ∂xj ). A general curvilinear transformation (GCT) of the second-order partial derivatives from x- to ξ-coordinates is given by (Arnold, 1989)

∂2 U x ∂ξl ∂ξm ∂2 U ξ ∂2 ξl ∂U ξ ¼ þ : ∂xi ∂xj ∂xi ∂xj ∂ξl ∂ξm ∂xi ∂xj ∂ξl

∂2 P ξ ∂2 P ξ ∂P ξ ∂ 2 Qξ ∂Qξ ¼ slm þ f lp1 þ slm þ f lq1 ; p1 q1 2 ∂ξl ∂ξm ∂ξl ∂ξl ∂ξm ∂ξl ∂t (11) ∂ 2 Qξ ∂2 P ξ ∂P ξ ∂ 2 Qξ ∂Qξ lm l lm ¼ s þ f þ s þ f lq2 ; p2 p2 q2 2 ∂ξl ∂ξm ∂ξl ∂ξl ∂ξm ∂ξl ∂t (12)

∂2 Qx ∂2 P x ¼ ½ðV 2Pn − V 2Sx Þðδij − aij Þ 2 ∂xi ∂xj ∂t þ ½V 2Sz ðδij − aij Þ þ V 2Pz aij 

C267

(8)

I remark that the GCT transformation from Cartesian to a generalized coordinate system introduces the first-order partial differential operators into the generalized acoustic TTI wave equations. As rank-zero tensors, scalar fields (e.g., V Pz ) are invariant under coordinate transformation, and thus have the same values at point ξ0 , as they do at the corresponding Cartesian point x0 , and thereby require no additional geometric weighting factors. Applying the GCT in equation 8 to equations 6 and 7 yields

 l m 2  ∂2 P ξ ∂ Pξ ∂2 ξl ∂P ξ 2 ij ij 2 ij ∂ξ ∂ξ ¼ ½V Px ðδ − a Þ þ V Sz a  þ ∂xi ∂xj ∂ξl ∂ξm ∂xi ∂xj ∂ξl ∂t2  l m 2  ∂ Qξ ∂2 ξl ∂Qξ 2 2 ij ∂ξ ∂ξ þ þ ½ðV Pz − V Sz Þa  ; (9) ∂xi ∂xj ∂ξl ∂ξm ∂xi ∂xj ∂ξl  l m 2  ∂ 2 Qξ ∂ξ ∂ξ ∂ P ξ ∂2 ξl ∂P ξ 2 2 ij ij ¼½ðV −V Þðδ −a Þ þ Pn Sz ∂xi ∂xj ∂ξl ∂ξm ∂xi ∂xj ∂ξl ∂t2  l m 2  ∂ξ ∂ξ ∂ Qξ ∂2 ξl ∂Qξ þ þ½V 2Sz ðδij −aij ÞþV 2Pz aij  ; ∂xi ∂xj ∂ξl ∂ξm ∂xi ∂xj ∂ξl (10) where summation occurs over all i, j, l, and m indices. These expressions may be written compactly using the following notation:

where I use a convention in which coefficients f and s represent the first- and second-order derivative terms, subscripts p or q indicate the wavefield being acted upon (e.g., Pξ or Qξ ), and subscripts 1 or 2 provide the equation number. Note that f and s absorb the trigonometric aij coefficients and the corresponding velocity fields. The two main differences between acoustic TTI propagation in a generalized 3D coordinate system relative to a 3D Cartesian mesh are the introduction of: (1) first-order partial derivative terms that largely use the same wavefield values as the second-order derivatives but with different stencil coefficients and (2) more complex coefficients that premultiply the second-order partial differential operators. This suggests additional computational overhead when simulating acoustic TTI wavefields on 3D generalized coordinate meshes. However, for scenarios where one defines a generalized 3D coordinate system analytically, all partial derivatives found in the f and s coefficients likewise may be precomputed analytically, save for derivatives of the topographic surface, which must be evaluated numerically. Computing analytical derivatives leads to simplified expressions and forestalls the introduction of errors from imprecise numerical differentiation of the various chain-rule quantities in equation 8. This also prevents significant losses in computational efficiency because one need not repeatedly compute or store/ access 3D numerical representations of geometric coefficients or partial derivative fields held in global memory (see the “Discussion” section). I illustrate these points below for analytic 2D and 3D topographic coordinate systems.

TOPOGRAPHIC COORDINATES The generalized 3D acoustic TTI wave equations specified above are directly applicable to the problem of performing two-way acoustic propagation through a TTI medium underlying an irregular topographic surface. A topographic coordinate system may be specified using an analytic mapping based on a flattening transformation (Tessmer et al., 1992; Hestholm, 1999; Shragge, 2014b):

3 3 2 x 1 ξ1 7 4 ξ2 5 ¼ 6 4  x 2  5; x3 −T a a−T ξ3 2

(13)

where T ¼ Tðx1 ; x2 Þ is the topographic profile, the vertical coordinate ranges over ξ3 ¼ ½0; a, and a is the maximum depth of the coordinate system (here assumed to be a flat plane). I note that a more complex lower surface (i.e., a ¼ aðx1 ; x2 Þ) could be incorporated into the mapping equations; however, this approach has not been implemented in this work. The 3D topographic coordinate system in equation 13 represents a linear vertical stretch of the mesh and leads to a nonorthogonal basis that spans the 3D computational domain. The equations for acoustic TTI propagation in 3D topographic coordinates can be

Shragge

Downloaded 09/19/16 to 106.68.109.41. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

C268

found by substituting the analytic expression from equation 13 into equations 11 and 12. After somewhat involved computations, one can generate analytic expressions for the coefficient fields, which, for brevity, are reproduced in Appendix A. For 2D scenarios, one may assume a symmetry axis in the ξ1 − ξ3 plane (i.e., azimuth ϕ ¼ 0° ) and thereby eliminate the ξ2 -coordinate dependence. The corresponding 2D mapping equations, based on topographic profile Tðx1 Þ, are given by



ξ1 ξ3





  x1  ¼ : 3 −T a xa−T

(14)

Appendix A also presents the coefficients for generalized 2D acoustic TTI wave propagation.

4 ∂2 U 1 X ≈ sl ½U niþl;j;k þ U ni−l;j;k ; ∂ξ21 Δξ21 l¼0

(17)

where the sl coefficient vector is



205 ; ½sl  ¼ − 72

8 ; 5

1 − ; 5

8 ; 315

 1 − : 560

The remaining mixed second-order partial derivative terms are given by, e.g. 4 X 4 ∂2 U 1 X ≈ m ½U n − U niþl;j−p;k ∂ξ1 ∂ξ2 Δξ1 Δξ2 l¼1 p¼1 lp iþl;jþp;k

− U ni−l;jþp;k þ U ni−l;j−p;k ;

(19)

with the mlp coefficient matrix given by

NUMERICAL IMPLEMENTATION Numerical implementation of the 2D and 3D TTI propagators described above is based on the finite-difference time-domain (FDTD) operators found in the sfttifd2d and sfttifd3d Madagascar programs. I have extended these codes to incorporate finite-difference (FD) operators of OðΔx8 ; Δt2 Þ accuracy to approximate the continuous partial derivatives of equations 11 and 12. I replace continuous (generic) wavefield U with a discrete version, U ni;j;k , where superscript n represents the current time step and subscripts i, j, and k indicate spatial indices for the ξ1 -, ξ2 -, and ξ3 -axes, respectively. The FD expressions used here are based on familiar Taylor-series expansions for FD operators. First-order partial derivative operators are given by, e.g. 4 ∂U 1 X ≈ f ½U n − U ni;j;k−l ; ∂ξ3 Δξ3 l¼1 l i;j;kþl

(18)

(15)

2 ½mlp  ¼

16 25 6−4 6 1625 4 525 1 − 260

4 − 25

1 144 4 − 525 1 1040

16 525 4 − 525 16 11025 1 − 7350

1 − 260

3

1 7 1040 7: 1 5 − 7350 1 78400

(20)

The temporal partial derivative is discretized using a standard OðΔt2 Þ approximation

∂2 U 1 n n−1 ≈ 2 ½U nþ1 i;j;k − 2U i;j;k þ U i;j;k : 2 ∂t Δt

(21)

Substituting the above FD approximations into equations 11 and 12 and solving for the U nþ1 i;j;k wavefield yields a FDTD algorithm for propagating a 3D acoustic TTI wavefield.

Remarks on implementation where the f l coefficient vector is

 ½f l  ¼

4 ; 5

1 − ; 5

4 ; 105

 1 − : 280

(16)

Second-order partial derivative operators are given by, e.g.

Figure 1. Foothills coordinate system formed by the topographic mapping in equation 14.

I treat boundary reflections through a successive application of absorbing boundary condition (Clayton and Engquist, 1977) and exponential sponge operators. I smoothly extend the boundaries outward in the ξ1 and ξ2 dimensions to ensure that the spatial derivatives of the topographic surface found in the coefficients of Appendix A (e.g., ∂T∕∂ξ1 ) are identically zero in the boundary regions. I apply a zero-velocity layer free-surface boundary condition by fixing four fictitious zero-valued wavefield points above the boundary to enable evaluation of the FD stencils within the boundary halo region. Although the example below demonstrates that this approach is sufficient for simulating wavefields in topographic coordinates with an accuracy roughly equivalent to that simulated with the corresponding Cartesian implementation, I point out that this approach only approximately satisfies the free-surface condition at a lower order of accuracy than used in the interior region and will thus introduce numerical errors related to the mismatch in order of accuracy. In addition, the zero-velocity layer approach has an increased potential for long-term instability issues over many wavefield grid crossings due to the presence of variable geometric coefficients within the computational boundary halo. However, in the author’s experience, the zero-velocity layer approximation has proved to be computationally efficient and sufficiently accurate in several RTM applications.

Topographic 3D TTI propagation

Downloaded 09/19/16 to 106.68.109.41. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

NUMERICAL EXAMPLES This section presents numerical tests of the developed acoustic TTI propagators in 2D and 3D topographic coordinates for elliptic and anelliptic TTI media, as well as a demonstration of the accuracy of the approach in handling free-surface interactions.

2D examples The first example is based on the 2D Foothills synthetic model, which is a merger of common geologic features found in the

C269

Canadian Foothills region of northeastern British Columbia, Canada (Gray and Marfurt, 1995). Here, I use a 10 km section of the model’s relief profile to form the topographic coordinate mesh using equation 14 above. This mesh, shown in Figure 1, exhibits approximately 1.2 km of topographic variability and has a lower grid boundary set at a ¼ 4 km depth. The topographic mesh has N ξ1 × N ξ3 ¼ 1400 × 533 grid points regularly sampled in the computational domain, which includes an aforementioned padding of 40 grid points to either side of the topographic surface. The spatial sampling intervals are set at Δξ1 ¼ 7.5 m and range vertically between 8.0 and 10.3 m when viewed in the physical Cartesian

a)

b)

c)

d)

e)

f)

Figure 2. Example of 2D wave propagation through a homogeneous elliptical TTI medium defined by model parameters V Pz ¼ 3.0 km∕s, anisotropic parameters ϵ ¼ δ ¼ 0.1, and tilt angle θ ¼ 45° . Wavefield snapshots in the left panels are computed in a Cartesian domain with a horizontal free surface, whereas those shown in the right panels are simulated in topographic coordinates with an irregular free surface shown by the thick black line and then interpolated back to Cartesian. Elliptical TTI traveltimes are indicated by the blue lines. (a and b) t ¼ 0.45, (c and d) 0.90, and (e and f) 1.35 s.

Downloaded 09/19/16 to 106.68.109.41. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

C270

Shragge

domain. The time-stepping interval is set at Δt ¼ 0.75 ms, which globally satisfies the deformed mesh stability criterion discussed in Shragge (2014b). The 2D Cartesian meshes are sampled at Δx1 ¼ Δx3 ¼ 7.5 m. Visualizing the resulting wavefields requires accounting for the spatially variable Δx3 by interpolating between the topographic computational and Cartesian domains. I achieve this using 2D (or 3D as appropriate) sinc interpolators weighted by the Jacobian of the transformation (Shragge, 2014a). The first 2D wave-propagation example is for a homogeneous elliptical TTI medium with vertical velocity defined by V Pz ¼ 3.0 km∕s, anisotropy parameters ϵ ¼ δ ¼ 0.1 (i.e., V Px ¼ V Pn ¼ 3.28 km∕s), zero V Sz , and a tilt angle of θ ¼ 45° . Figure 2 presents

wavefield snapshots for a point source (using a 45 Hz Ricker wavelet) excited at 1.4 km depth and centered horizontally in the computational grid. The left panels present propagation snapshots for a wavefield computed directly in the Cartesian mesh with a free surface at x3 ¼ 0.0 m. The right panels show wavefield snapshots computed in topographic coordinates and then interpolated to the Cartesian domain. I indicate the topographic surfaces using a thick black line. The wavefield snapshots in the upper panels are taken at t ¼ 0.45 s after source excitation, the one in the middle panels is taken at t ¼ 0.90 s, and in the lower panels is taken at t ¼ 1.35 s. To demonstrate the accuracy of the above impulse responses, I overlay analytical elliptical TTI traveltime surfaces (Alkhalifah,

a)

b)

c)

d)

e)

f)

Figure 3. Snapshots of 2D propagation through an anelliptical TTI model with ϵ ¼ 0.4, δ ¼ 0.1, and θ ¼ 45° . (a) t ¼ 0.23, (b) 0.45, (c) 0.68, (d) 0.90, (e) 1.13, and (f) 1.35 s.

Downloaded 09/19/16 to 106.68.109.41. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

Topographic 3D TTI propagation 2011) computed with the corresponding model parameters used to simulate the wavefield snapshots. The wavefield sequence computed in the Cartesian model (left panels) shows the impulse propagating upward, interacting with the free surface, and generating the single free-surface reflected event with the required negative reflection coefficient. The wavefield sequence computed in topographic coordinates similarly propagates upward and interacts with the free surface; however, due to the topographic complexity, additional multiple-scattered events are now present in the downgoing wavefield. Wavefield components in all panels that have not interacted with the free-surface agree very well with the analytic elliptical traveltime surfaces, indicating very good kinematic accuracy of the numerical propagation in this TTI medium — even in the presence of spatially varying geometric coefficients. Figure 3 presents 2D wave-propagation examples for an anelliptic TTI medium that is identical to the previous example, save for an increase of ϵ ¼ 0.4 (i.e., V Px ¼ 4.02 km∕s) and the introduction pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi of nonzero S-wave velocity, V Sz ¼ V Pz 0.75ðϵ − δÞ, to avoid numerical instability (Fletcher et al., 2009). The source location and 45 Hz wavelet remain unaltered. Figure 3a–3f presents propagation snapshots starting at t ¼ 0.23 s and subsequent increments of Δt ¼ 0.23 s. The anellipticity of the impulse response and presence of an S-wave artifact are evident in these panels (Fletcher et al., 2009).

3D examples I next present an example of 3D TTI wave propagation for a topographic profile derived from relief in the San Francisco Bay Area. Figure 4 shows a topographic surface with elevations ranging from sea level at x3 ¼ 0.0 km to a maximum of x3 ¼ 0.55 km. The topographic mesh has N ξ1 × N ξ2 × N ξ3 ¼ 401 × 401 × 151 grid points that are regularly sampled in the computational domain with a spacing of Δξ1 ¼ Δξ2 ¼ 40 m and range vertically between 40.0 and 43.66 m when viewed in the physical Cartesian domain. Interpolated Cartesian wavefields are sampled at Δx1 ¼ Δx2 ¼ Δx3 ¼ 40.0 m. The first 3D wave-propagation example is for a homogeneous elliptical TTI medium defined by vertical velocity V Pz ¼ 2.0 km∕s, anisotropy parameters ϵ ¼ δ ¼ 0.1 (i.e., V Px ¼ V Pn ¼ 2.53 km∕s), and a θ ¼ 45° tilt oriented at azimuth ϕ ¼ 45° (i.e., with a northwest–southeast fast-axis orientation). I simulate an impulsive source with 10 Hz Ricker wavelet at 2.5 km depth horizontally centered in the computational grid. Figure 5 presents wavefield propagation snapshots, with the left and right panels showing propagation computed in the Cartesian and topographic coordinate domains, respectively. The snapshots in the upper, middle, and lower panels (taken at t ¼ 1.28, 2.24, and 3.2 s, respectively) show the elliptical impulse responses after free-surface reflection. The Cartesian impulse responses again reflect from the free surface in a straightforward manner; however, those generated in topographic coordinates show additional complexity due to topographic scattering (i.e., at points interior to the outer elliptical waveform). Figure 6 presents wavefield snapshots in an anelliptic medium defined by vertical velocity V Pz ¼2.0km∕s, anisotropy parameters ϵ¼0.4 and δ ¼ 0.1 (V Px ¼ 2.68 and 2.19 km∕s, respectively), tilt angle θ ¼ 45° , and azimuth ϕ ¼ 45° . The snapshots in Figure 6a–6d shows the increasing complexity of the near-surface reflection response, as well as the expected anelliptical moveout.

C271

Accuracy of boundary interactions Although the above examples demonstrate the accuracy of the acoustic TTI propagators within the interiors of 2D and 3D topographic meshes, these tests have not yet explicitly demonstrated an ability to accurately model free-surface interactions. One way to examine the kinematic and dynamic accuracy of the FDTD implementation is by comparing numerical and analytic solution for some canonical nonflat Cartesian surface. For example, de la Puente et al. (2014) compare numerical simulation results from a 2D mimetic elastic FDTD approach to analytic expressions from a tilted 2D Garvin’s (1956) test based on the response of a semi-infinite isotropic elastic medium. However because the acoustic TTI wave equations discussed above have no known exact solution, setting up a similar analytic comparison is not possible. Inspired by the tilted 2D Garvin test, I instead investigate a relative comparison to show the equivalence of 2D TTI wavefield solutions simulated in Cartesian and topographic coordinates. Figure 7 highlights the geometry and wavefield results for the comparative 2D tilted Garvin test using an elliptical 2D TTI pffiffiffi earth model defined by parameters V P ¼ 3 km∕s, V S ¼ 3∕ 3km∕s, ϵ ¼ δ ¼ 0.4, and θ ¼ 45° . The black line in Figure 7a shows a free-surface boundary tilting at a 45° angle (i.e., with a normal parallel to the TTI tilt axis), which is used in the 2D mapping transformation (equation 14) to construct a topographic coordinate system similar to that shown in Figure 1. The first step in the comparative test is to simulate a 20 Hz Ricker wavefield impulse response in the topographic coordinate system from selected point x1 ¼ x3 ¼ 1 km. After injecting the source wavefield and simulating 1.0 s of wavefield propagation on the topographic mesh, I interpolate the resulting wavefields back to a Cartesian reference frame for comparison purposes. Figure 7a shows the superposition of two wavefield snapshots extracted at t ¼ 0.05 and 0.625 s. I include the former snapshot to show the location and the polarity of the modeled source. By t ¼ 0.625 s, the impulse has propagated to, and reflected back from, the freesurface boundary with the expected polarity reversal. To test kinematic and dynamic accuracy of the above approach in handling boundary reflections, I compare the results in

Figure 4. Topographic model of the San Francisco Bay Area with 0.55 km relief.

Downloaded 09/19/16 to 106.68.109.41. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

C272

Shragge

a)

b)

c)

d)

e)

f)

Figure 5. Example of 3D wave propagation through a homogeneous elliptical TTI medium defined by V Pz ¼ 3.0 km∕s, ϵ ¼ δ ¼ 0.1, θ ¼ 45° , and ϕ ¼ 45° . Wavefield snapshots in the left panels are computed in a Cartesian domain with a horizontal free surface, whereas those shown in the right panels are simulated through a 3D topographic coordinate mesh and interpolated back to Cartesian. Velocity models have been superimposed to illustrate the free surface. (a and b) t ¼ 1.28, (c and d) 2.24, and (e and f) 3.2 s.

Downloaded 09/19/16 to 106.68.109.41. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

Topographic 3D TTI propagation Figure 7a with a mirror source of opposing polarity from point at x1 ¼ x3 ¼ −1 km. I simulate this mirror wavefield on a Cartesian mesh that is free of the topographic boundary. Importantly, the simulated Cartesian wavefield can be used to test of the equivalence of Cartesian and topographic coordinate modeling approaches by seeing whether the propagated Cartesian mirror source solution exactly negates the free-surface reflection observed in the topographic wavefield solution. Figure 7b shows the superposition and reverse polarity of two wavefield snapshots extracted at the same times as Figure 7a. Figure 7c shows the direct subtraction of the mirrorsource wavefield in Figure 7b from wavefield shown in Figure 7a. Importantly, no noticeable free-surface reflection energy remains indicating that the Cartesian mirror source successfully cancels the topographic boundary reflection as required. To more comprehensively illustrate the efficacy of wavefield cancellation, Figure 7d and 7e presents comparative traces extracted at horizontal location x3 ¼ −0.5 km and at depth x1 ¼ −0.25 km, respectively. In each figure part, the blue and yellow lines show the topographic and mirror Cartesian wavefield traces, whereas the black

C273

line presents the wavefield difference. Overall, a very good match is observed with only minor variations that are attributable to a combination of factors: (1) source wavefield interpolation into the computational grid, (2) different numerical dispersion due to inclusion of different FD terms with different coefficients, (3) errors potentially introduced through the zero-velocity layer approximation, and (4) interpolation of wavefield propagation results back to a Cartesian mesh. Thus, these results suggest that the full-wavefield TTI propagation operators and their respective adjoints would be useful as computational kernel functions for implementing RTM or FWI directly from a free-surface boundary influenced by topography.

DISCUSSION The computational efficiency of the above implementation is largely controlled by the memory access of mixed second-order partial derivative stencils which, for these OðΔx8 Þ operators, require memory access of 64 neighboring wavefield values. The firstand second-order partial-derivative calculations involve accessing

a)

b)

c)

d)

Figure 6. Snapshots of 3D wave propagation through an anelliptical TTI model defined by V Pz ¼ 2.0 km∕s, ϵ ¼ 0.4, δ ¼ 0.1, θ ¼ 45° , and ϕ ¼ 45° . (a) t ¼ 2.24, (b) 2.56, (c) 2.88, and (d) 3.20 s.

Shragge

Downloaded 09/19/16 to 106.68.109.41. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

C274

eight and nine neighboring wavefield values, respectively. Thus, for a single time step, the additional computational overhead in topographic coordinates relative to a Cartesian implementation in 2D is 10% (i.e., 90 versus 82 points) and in 3D is 4% (i.e., 227 versus 219 points). Because the ξ ¼ ½ξ1 ; ξ2  coordinates of the topographic mesh are independent of the topographic profile, it is computationally

a)

d)

advantageous to have the ξ3 dimension as the “fast axis” of memory array alignment. The array access of the topographic values and their derivatives (e.g., T and T xx in Appendix A) occurs only once per N ξ3 stencil operations, whereas accessing the ξ3 -dependent aij and velocity values is required at each Δξ3 location — even for Cartesian implementations. Thus, these topographic contributions represent fairly trivial overheads of factors 1 þ ð3∕7N ξ3 Þ and

b)

c)

e)

Figure 7. Snapshots of 2D propagation at t ¼ 0.625 s through an elliptical TTI model with ϵ ¼ δ ¼ 0.4 and θ ¼ 45° . (a) Simulation in a topographic coordinate systems formed by a linear ramp (shown in black) using a point source at x1 ¼ x3 ¼ 1 km and interpolated back to Cartesian. (b) Simulation in a Cartesian coordinate system of a mirror point source at x1 ¼ x3 ¼ −1 km. (c) Panels (a and b): wavefield difference showing the cancellation of free-surface reflected wave by the Cartesian mirror source. (d) Comparison for traces extracted at horizontal position x3 ¼ −0.5 km, where yellow, blue, and green lines indicate topographic, Cartesian, and wavefield difference from (a to c), respectively. (e) Comparison for traces extracted from wavefields at depth x1 ¼ −0.25 km, with the same color scheme as (d).

Downloaded 09/19/16 to 106.68.109.41. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

Topographic 3D TTI propagation 1 þ ð3∕5N ξ3 Þ for 2D and 3D applications, respectively. Evaluating the various derivative coefficients for the topographic profile introduces additional floating-point computations; however, these are largely hidden by the latency of memory access of wavefield values for the aforementioned stencil operations. Although the FD stencils described above represent a straightforward numerical implementation and generate accurate interior point and free-surface boundary interactions for the provided examples, they do not represent a panacea for all topographic scenarios. For example, the 3D coordinate mapping equations preclude scenarios, in which topography is effectively discontinuous (e.g., vertical cliff faces) or has overhanging strata. Additional tests using the 3D topographic surfaces described above indicate that stretching the elevations by a factor of two can lead to numerical instability within 5 s of numerical simulation using the current discretization intervals. However, apart from extreme scenarios, and because land seismic preprocessing efforts commonly remap the seismic data volume to a floating datum that is substantially smoother than the true topographic surface, the FD operators discussed here should efficiently generate wavefield solutions of sufficient accuracy for most RTM or FWI applications of interest. A second topic worth discussing is that the above topographic coordinate approach will inherit all of the drawbacks noted in the Cartesian implementation. For example, this approach does not explicitly attenuate S-wave artifacts, which would contaminate any TTI-RTM images without additional wavefield operations. However because Guan et al. (2011) explicitly use a general (i.e., coordinate independent) form of the H 1 and H2 differential operators in a preimaging wavefield filtering operation, this procedure should be equally applicable when using the generalized H1 and H 2 operators described above. In addition, I note that the partial differential operators in Bin Waheed and Alkhalifah (2015) lend themselves to reformulation in a topographic coordinate system, and following this approach would preclude S-wave artifacts by design. Both of these assertions, though, remain open research questions at this time. Another interesting question is how does the smoothness of coordinate transformations affect stability and numerical dispersion? The spatial variability of the free surface in a topographic coordinate system introduces three additional consequences. First, there are instability concerns in any FD method whenever the effective grid velocity Δx∕Δt drops below the physical velocity. The variable computational grid size, when viewed in a Cartesian coordinate system, effectively modifies the local grid velocity, such that local instabilities can arise in voxels with smaller Δx sampling. One approach is to use a Δt that satisfies the global stability requirements (see Shragge, 2014b). A second related issue is that the variable Δxi ðξj Þ sampling generates voxels of greater volume, which reduces the effective points per wavelength wavefield sampling and thereby increases numerical dispersion (see discussion in Shragge, 2014b). A more fundamental issue is that the FD operators used in this work (such as all FD operators based solely on Taylor-series expansions) do not explicitly satisfy underlying conservation equations, and are thus subject to nonphysical fluxes entering into through the boundaries of the computational domain (Castillo et al., 2001). Although these effects also are present in Cartesian FD approaches and can cause long-term instabilities, they are more acute in scenarios involving irregular boundaries because of an increased likelihood of introducing nonphysical fluxes and associated insta-

C275

bilities at earlier stages of propagation than in Cartesian scenarios. Following a mimetic FD approach for acoustic-/elastic-wave propagation (Castillo and Miranda, 2013; de la Puente et al., 2014; Solano-Feo et al., 2016) that uses fully staggered grids and a modified boundary layer representation to formally satisfy underlying conservation laws, should offer a more robust approach for handling boundary fluxes. However, the effectiveness of implementing mimetic FD operators for the TTI model scenarios described above remains an open research topic and would not be a straightforward technical undertaking. Finally, topographic meshes represent one class of coordinate system, for which the generalized acoustic TTI media propagation theory is applicable. For example, one could apply this approach in an “internal boundary” mesh (Shragge, 2014b) that partitions a volume of interest into isotropic and anisotropic regions (e.g., at the location of a geologic unconformity); however, this is beyond the scope of this work.

CONCLUSIONS This paper demonstrates that 2D/3D wave propagation in acoustic TTI media can be extended to generalized 3D computational geometries. This extension can be achieved through a GCT of the partial derivative operators forming a new standard governing system of coupled second-order PDEs. The GCT leads to a TTI wave equation with the first- and second-order partial derivative operators, which is different from Cartesian implementations and requires an altered FDTD scheme. I show that it is straightforward to extend acoustic wave propagation to 2D/3D TTI media using a topographic coordinate system formed by linear interpolation between a complex topographic surface and a uniform level at depth. Specifying an analytic topographic mesh affords allows for partial-differential coefficients to be evaluated analytically. These expressions can be precomputed and incorporated directly into FDTD solution schemes without introducing significant computational overheads relative to an equivalent Cartesian implementation. The 2D and 3D examples illustrate the ability to accurately model the kinematics of (an)elliptical TTI wave propagation in topographic and Cartesian coordinate meshes. Simulated impulse responses show very good agreement when compared with analytic traveltime expressions for homogeneous elliptical TTI media. A 2D boundary accuracy test demonstrates the accuracy of wavefields simulated on representative topographic mesh relative to ones computed on Cartesian grids, save for errors potentially introduced by approximate boundary conditions, interpolation noise, and difference in numerical dispersion. Accordingly, these generalized acoustic TTI propagators (and their adjoints) would be useful for undertaking most RTM or FWI applications using computational domains conforming to free-surface topography.

ACKNOWLEDGMENTS P. Sava and colleagues at Colorado School of Mines are thanked for making available the sfttifd2d and sfttifd3d codes, on which the sfttifd2d_topo and sfttifd3d_topo programs are based. I acknowledge the support of the UWA:RM Consortium sponsors. This work was supported by resources provided by the Pawsey Supercomputing Centre with funding from the Australian Government and the Government of Western Australia. Analytic coordinate system

Shragge

C276

Downloaded 09/19/16 to 106.68.109.41. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

geometry results were verified using the Mathematica software package (http://www.wolfram.com/mathematica/). Reproducible numerical examples were generated using the Madagascar package (http://www.ahay.org).

Coefficients f 3q1 and slm q1 are given by

ða − TÞ2 3 f ¼ a22 ð1 − ξ3 Þð2T 22 þ ða − TÞT 2;2 Þ V 2Pz − V 2Sz q1 þ a11 ð1 − ξ3 Þð2T 21 þ ða − TÞT 1;1 Þ − a13 T 1 − a23 T 2 þ a12 ð1 − ξ3 Þ

APPENDIX A ACOUSTIC TTI PROPAGATION COEFFICIENTS FOR TOPOGRAPHIC COORDINATE SYSTEMS This appendix presents the nonzero coefficients for TTI propagation in 2D and 3D topographic coordinates. I have computed these coefficients by evaluating the GCT approach described above. Note that: (1) these expressions explicitly depend on the topographic surface T ¼ Tðx1 ; x2 Þ and spatial derivatives thereof as well as the ξ3 coordinate and (2) the following shorthand notation is used: ð∂T∕∂xi Þ ¼ T i and ð∂2 T∕∂xi ∂xj Þ ¼ T i;j .

Coefficients for TTI propagation in 3D topographic coordinates Coefficients

f 3p1

and

slm p1

are given by

× ð1 − ξ3 Þð2T 22 þ ða − TÞT 2;2 Þ − ða11 ðV 2Sz − V 2Px Þ þ V 2Px Þð1 − ξ3 Þ þ ðV 2Px − V 2Sz Þ½a12 ð1 − ξ3 Þð2T 2 T 1 þ ða − TÞT 1;2 Þ

ða −

(A-1)

11 2 2 2 c11 p1 ¼ a ðV Sz − V Px Þ þ V Px ;

(A-2)

12 2 2 c12 p1 ¼ a ðV Sz − V Px Þ;

(A-3)

(A-9)

1 s12 ¼ a12 ; − V 2Sz q1

(A-10)

V 2Pz

a−T s13 ¼ a13 − ð1 − ξ3 Þða12 T 2 þ 2a11 T 1 Þ; (A-11) − V 2Sz q1

V 2Pz

1 s22 ¼ a22 ; V 2Pz − V 2Sz q1

− 2½ða11 ðV 2Sz − V 2Px Þ þ V 2Px Þð1 − ξ3 ÞT 1 ;(A-4)

þ ð1 − ξ3 Þ2 ða22 T 22 þ a11 T 21 þ a12 T 1 T 2 Þ: (A-14) Coefficients f 3p2 and slm p2 are given by

ða − TÞ2 3 f ¼ −a13 T 1 − a23 T 2 V 2Pn − V 2Sz p2 þ a12 ð1 − ξ3 Þð2T 1 T 2 þ ða − TÞT 1;2 Þ

(A-5)

þ ð1 − a11 Þð1 − ξ3 Þð2T 21 þ ða − TÞT 1;1 Þ; (A-15)

V 2Sz

1 s11 ¼ a11 − 1; − V 2Pn p2

(A-16)

1 s12 ¼ a12 ; − V 2Pn p2

(A-17)

2 2 23 12 ða − TÞc23 p1 ¼ ðV Sz − V Pz Þ½a − a ð1 − ξ3 ÞT 1 

− 2½ða22 ðV 2Sz − V 2Px Þ þ V 2Px Þð1 − ξ3 ÞT 2 ;(A-6) 33 2 2 2 ða − TÞ2 c33 p1 ¼ ða ðV Sz − V Px Þ þ V Px Þ

þ ða22 ðV 2Sz − V 2Px Þ þ V 2Px Þð1 − ξ3 Þ2 T 22 þ ða11 ðV 2Sz − V 2Px Þ þ V 2Px Þð1 − ξ3 Þ2 T 21

V 2Sz

a−T s13 ¼ a13 − ð1 − ξ3 Þða12 T 2 þ 2ða11 − 1ÞT 1 Þ; V 2Sz − V 2Pn p2 (A-18)

þ ðV 2Px − V 2Sz Þð1 − ξ3 Þ½a13 T 1 þ a23 T 2 −

a12 ð1

− ξ3 ÞT 1 T 2 :

(A-12)

þ ð1 − a22 Þð1 − ξ3 Þð2T 22 þ ða − TÞT 2;2 Þ

¼ ðV 2Sz − V 2Pz Þ½a13 − a12 ð1 − ξ3 ÞT 2 

22 2 2 2 c22 p1 ¼ a ðV Sz − V Px Þ þ V Px ;

1 s11 ¼ a11 ; V 2Pz − V 2Sz q1

ða − TÞ2 33 s ¼ a33 − ð1 − ξ3 Þða13 T 1 þ a23 T 2 Þ V 2Pz − V 2Sz q1

× ð2T 21 þ ða − TÞT 1;1 Þ

TÞc13 p1

(A-8)

a−T s23 ¼ a23 − ð1 − ξ3 Þða12 T 1 þ 2a22 T 2 Þ; (A-13) V 2Pz − V 2Sz q1

ða − TÞ2 f 3p1 ¼ −ða22 ðV 2Sz − V 2Px Þ þ V 2Px Þ

− a13 T 1 − a23 T 2 ;

× ð2T 1 T 2 þ ða − TÞT 1;2 Þ;

(A-7)

V 2Sz

1 s22 ¼ a22 − 1; − V 2Pn p2

(A-19)

Topographic 3D TTI propagation

a−T s23 ¼ a23 − ð1 − ξ3 Þða12 T 1 þ 2ða22 − 1ÞT 2 Þ; − V 2Pn p2 (A-20)

11 2 2 2 c11 p1 ¼ a ðV Sz − V Px Þ þ V Px ;

Downloaded 09/19/16 to 106.68.109.41. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

V 2Sz

ða − TÞ3 33 s ¼ ða33 − 1Þða − TÞ − ð1 − ξ3 Þða13 T 1 þ a23 T 2 Þ V 2Sz − V 2Pn p2 þ ð1 − ξ3 Þ2 ðða22 − 1ÞT 22 þ ða11 − 1ÞT 21 þ a12 ða − TÞT 1 T 2 Þ:

C277

13 2 2 ða − TÞc13 p1 ¼ a ðV Sz − V Pz Þ − 2ð1 − ξ3 ÞT 1

× ða11 ðV 2Sz − V 2Px Þ þ V 2Px Þ;

× ð1 − ξ3 ÞT 1 þ ða11 ðV 2Sz − V 2Px Þ þ V 2Px Þ × ð1 − ξ3 Þ2 T 21 :

þ ða − TÞT 1;2 Þ − ða22 ðV 2Pz − V 2Sz Þþ V 2Sz Þð1− ξ3 Þ

ða − TÞ2 3 f ¼ a13 T 1 − a11 ð1 − ξ3 Þ V 2Pz − V 2Sz q1

×ð2T 22 þ ða − TÞT 2;2 Þ − ða11 ðV 2Pz − V 2Sz Þ þ V 2Sz Þ ×ð1 − ξ3 Þð2T 21 þ ða − TÞT 1;1 Þ;

(A-22)

11 2 2 2 s11 q2 ¼ a ðV Pz − V Sz Þ þ V Sz ;

(A-23)

12 2 2 s12 q2 ¼ a ðV Pz − V Sz Þ;

(A-24)

13 12 2 2 ða − TÞs13 q2 ¼ ða − a ð1 − ξ3 ÞT 2 ÞðV Pz − V Sz Þ

− 2ða11 ðV 2Pz − V 2Sz Þ þ V 2Sz Þð1 − ξ3 ÞT 1 ; (A-25) 22 2 2 2 s22 q2 ¼ a ðV Pz − V Sz Þ þ V Sz ;

−

−

þ

V 2Sz Þð1

33 2 2 2 2 2 ða − TÞ2 s33 q2 ¼ a ðV Pz − V Sz Þ þ V Sz − ðV Pz − V Sz Þð1 − ξ3 Þ

× ða13 T 1 þ a23 T 2 − a12 ð1 − ξ3 ÞT 1 T 2 Þ þ ð1 − ξ3 Þ2 ½ða22 ðV 2Pz − V 2Sz Þ þ V 2Sz ÞT 22

1 s11 ¼ a11 ; − V 2Sz q1

(A-34)

a−T s13 ¼ a13 − 2a11 ð1 − ξ3 ÞT 1 ; − V 2Sz q1

V 2Pz

(A-35)

ða − TÞ2 33 s ¼ a33 − a13 ð1 − ξ3 ÞT 1 þ a11 ð1 − ξ3 Þ2 T 21 : V 2Pz − V 2Sz q1 (A-36)

ða − TÞ2 3 f ¼ a13 T 1 þ ð1 − a11 Þð1 − ξ3 Þ V 2Sz − V 2Pn p2 × ð2T 21 þ ða − TÞT 1;1 Þ; 1 s11 ¼ a11 − 1; V 2Sz − V 2Pn p2

(A-37) (A-38)

(A-28)

Note that in 3D Cartesian coordinates in which T and all its derivatives are equal to zero and ða − TÞ → 1, the coefficients simplify to the conventional expressions for 3D acoustic TTI wave propagation.

Coefficients for TTI propagation in 2D topographic coordinates Assuming that the azimuth of the fast axis lies in the ξ1 − ξ3 plane (i.e., ϕ ¼ 0° ), the nonzero coefficients for 2D TTI propagation in a topographic coordinate system become significantly less complex. Coefficients f 3p1 and slm p1 are given by

ða − TÞ2 f 3p1 ¼ a13 ðV 2Sz − V 2Px ÞT 1 − ða11 ðV 2Sz − V 2Px Þ þ V 2Px Þ × ð1 − ξ3 Þð2T 21 þ ða − TÞT 1;1 Þ;

(A-33)

Coefficients f 3p2 and slm p2 are given by

− ξ3 ÞT 2 ; (A-27)

þ ða11 ðV 2Pz − V 2Sz Þ þ V 2Sz ÞT 21 :

V 2Pz

× ð2T 21 þ ða − TÞT 1;1 Þ;

(A-26)

23 12 2 2 ða − TÞs23 q2 ¼ ða − a ð1 − ξ3 ÞT 1 ÞðV Pz − V Sz Þ

V 2Sz Þ

(A-32)

Coefficients f 3q1 and slm q1 are given by

ða− TÞ2 f 3q2 ¼ ðV 2Pz −V 2Sz Þ½a23 T 2 þ a13 T 1 þ a12 ð1− ξ3 Þð2T 1 T 2

2ða22 ðV 2Pz

(A-31)

33 2 2 2 13 2 2 ða − TÞ2 c33 p1 ¼ a ðV Sz − V Px Þ þ V Px − a ðV Sz − V Px Þ

(A-21)

Coefficients f 3q2 and slm q2 are given by

(A-30)

(A-29)

a−T s13 ¼ a13 þ 2ð1 − a11 Þð1 − ξ3 ÞT 1 ; V 2Sz − V 2Pn p2

(A-39)

ða − TÞ2 33 s ¼ ða33 − 1Þ − a13 ð1 − ξ3 ÞT 1 V 2Sz − V 2Pn p2 þ ða11 − 1Þð1 − ξ3 Þ2 T 21 :

(A-40)

Coefficients f 3q2 and slm q2 are given by

ða − TÞ2 f 3q2 ¼ a13 ðV 2Pz − V 2Sz ÞT 1 − ða11 ðV 2Pz − V 2Sz Þ þ V 2Sz Þ × ð1 − ξ3 Þð2T 21 þ ða − TÞT 1;1 Þ;

(A-41)

11 2 2 2 s11 q2 ¼ a ðV Pz − V Sz Þ þ V Sz ;

(A-42)

Shragge

C278

13 2 2 11 2 2 2 ða − TÞs13 q2 ¼ a ðV Pz − V Sz Þ − 2ða ðV Pz − V Sz Þ þ V Sz Þ

× ð1 − ξ3 ÞT 1 ;

(A-43)

Downloaded 09/19/16 to 106.68.109.41. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

33 2 2 2 13 2 2 ða − TÞ2 s33 q2 ¼ a ðV Pz − V Sz Þ þ V Sz − a ðV Pz − V Sz Þ

× ð1 − ξ3 ÞT 1 þ ða11 ðV 2Pz − V 2Sz Þ þ V 2Sz Þ × ð1 − ξ3 Þ2 T 21 :

(A-44)

As above, for 2D topographic coordinates in which T and all its partial derivatives are equal to zero and ða − TÞ → 1, the coefficients simplify to the conventional Cartesian expressions for 2D acoustic TTI wave propagation.

REFERENCES Alkhalifah, T., 2000, An acoustic wave equation for anisotropic media: Geophysics, 65, 1239–1250, doi: 10.1190/1.1444815. Alkhalifah, T., 2011, Traveltime approximations for transversely isotropic media with an inhomogeneous background: Geophysics, 76, no. 3, WA31–WA42, doi: 10.1190/1.3555040. Arnold, V., 1989, Mathematical methods of classical mechanics, graduate texts in mathematics: Springer-Verlag. Baysal, E., D. Kosloff, and J. Sherwood, 1983, Reverse-time migration: Geophysics, 48, 1514–1524, doi: 10.1190/1.1441434. Bin Waheed, U., and T. Alkhalifah, 2015, An efficient wave extrapolation method for anisotropic media with tilt: Geophysical Prospecting, 63, 1126–1141, doi: 10.1111/1365-2478.12233. Castillo, J., J. Hyman, M. Shashkov, and S. Steinberg, 2001, Fourth- and sixth-order conservative finite difference approximations of the divergence and gradient: Applied Numerical Mathematics, 37, 171–187, doi: 10.1016/S0168-9274(00)00033-7. Castillo, J. E., and G. F. Miranda, 2013, Mimetic discretization methods: CRC Press. Clayton, R., and B. Engquist, 1977, Absorbing boundary conditions for acoustic and elastic wave equations: Bulletin of the Seismological Society of America, 67, 1529–1540. de la Puente, J., M. Ferrer, M. Hanzich, J. Castillo, and J. Cela, 2014, Mimetic seismic wave modeling including topography on deformed staggered grids: Geophysics, 79, no. 3, T125–T141, doi: 10.1190/geo2013-0371.1. Fletcher, R., X. Du, and P. Fowler, 2009, Reverse time migration in tilted transversely isotropic (TTI) media: Geophysics, 74, no. 6, WCA179– WCA187, doi: 10.1190/1.3269902. Garvin, W., 1956, Exact transient solution of the buried line source problem: Proceedings of the Royal Society of London, Series A: Mathematical, Physical and Engineering Sciences, 234, 528–541. Gittins, J., and R. Vestrum, 2004, Overcoming thrust-belt imaging problems in Magdalena Valley, Colombia: 74th Annual International Meeting, SEG, Expanded Abstracts, 2021–2023.

Gray, S. H., and K. J. Marfurt, 1995, Migration from topography: Improving the near-surface image: Journal of the Canadian Society of Exploration Geophysicists, 31, 18–24. Guan, H., E. Dussaud, B. Denel, and P. Williamson, 2011, Techniques for an efficient implementation of RTM in TTI media: 81st Annual International Meeting, SEG, Expanded Abstracts, 3393–3397. Hestholm, S., 1999, Three-dimensional finite difference viscoelastic wave modeling including surface topography: Geophysical Journal International, 139, 852–878, doi: 10.1046/j.1365-246x.1999.00994.x. Lan, H., Z. Zhang, J. Chen, and Y. Liu, 2014, Reverse time migration from irregular surface by flattening surface topography: Tectonophysics, 627, 26–37, doi: 10.1016/j.tecto.2014.04.015. McMechan, G. A., 1983, Migration by extrapolation of time-dependent boundary values: Geophysical Prospecting, 31, 413–420, doi: 10.1111/ j.1365-2478.1983.tb01060.x. Pratt, R., C. Shin, and G. Hicks, 1998, Gauss-Newton and full Newton methods in frequency-space seismic waveform inversions: Geophysical Journal International, 133, 341–362, doi: 10.1046/j.1365-246X.1998 .00498.x. Shragge, J., 2014a, Reverse time migration from topography: Geophysics, 79, no. 4, S141–S152, doi: 10.1190/geo2013-0405.1. Shragge, J., 2014b, Solving the 3D acoustic wave equation on generalized structured meshes: A finite-difference time-domain approach: Geophysics, 79, no. 6, T363–T378, doi: 10.1190/geo2014-0172.1. Solano-Feo, F., J. Guevara-Jordan, O. Rojas, B. Otero, and R. Rodriguez, 2016, A new mimetic scheme for the acoustic wave equation: Journal of Computational and Applied Mathematics, 295, 2–12, doi: 10.1016/j.cam .2015.09.037. Tarantola, A., 1984, Inversion of seismic reflection data in the acoustic approximation: Geophysics, 49, 1259–1266, doi: 10.1190/1.1441754. Tessmer, E., D. Kosloff, and A. Behle, 1992, Elastic wave propagation simulation in the presence of surface topography: Geophysical Journal International, 108, 621–632, doi: 10.1111/j.1365-246X.1992.tb04641.x. Thomsen, L., 1986, Weak elastic anisotropy: Geophysics, 51, 1954–1966, doi: 10.1190/1.1442051. Tsvankin, I., 2001, Seismic signatures and analysis of reflection data in anisotropic media: Elsevier. Ursenbach, C., and R. Bale, 2009, TTI wave-equation migration for Canadian Foothills depth imaging: The Leading Edge, 28, 1344–1351, doi: 10 .1190/1.3259613. Vestrum, R., 2002, 2D and 3D anisotropic depth migration case histories: Recorder, 27, 31–32. Vestrum, R., V. Dolgov, G. Wittman, L. Csontos, and J. Gittins, 2011, 3D seismic imaging over two structurally complex surveys in the foothills of Pakistan: First Break, 29, 61–70. Whitmore, N., 1983, Iterative depth imaging by backward time propagation: 53rd Annual International Meeting, SEG, Expanded Abstracts, 382–384. Xu, S., and H. Zhou, 2014, Accurate simulations of pure quasi-P-waves in complex anisotropic media: Geophysics, 79, no. 6, T341–T348, doi: 10 .1190/geo2014-0242.1. Zhang, H., G. Zhang, and Y. Zhang, 2009, Removing S-wave noise in TTI reverse time migration: 79th Annual International Meeting, SEG, Expanded Abstracts, 2849–2853. Zhang, M., and Z. C. Li, 2013, True-amplitude angle domain RTM Using time-shift imaging condition under complex topographical conditions: 75th Annual International Conference and Exhibition, EAGE, Extended Abstracts, P0904.