Dispersion and Stability Condition of Seismic Wave Simulation in TTI Media

For seismic waveform simulation in tilted transversely isotropic (TTI) media, we derive explicitly the numerical dispersion relation and the stability condition for the computation of a 2D pseudo-acoustic wave equation. The numerical dispersion relation indicates that the number of sampling points per wavelength has the greatest influence on the dispersion, while the anisotropic parameters of the TTI media and the mesh rotation angle have little influence on the dispersion. Given an appropriate spatial sampling, the stability condition is for the selection of the time step for the implementation of the TTI wave equation. We partition a numerical model using quadrangle grids in Cartesian coordinates, and map it to a computing model in which any non-rectangular meshes in Cartesian coordinates become rectangular meshes. Then we reformulate the pseudo-acoustic wave equation for the TTI media accordingly in the computational space. We implement seismic waveform simulation using the second-order finite-difference method straightforwardly, and show examples with a desirable accuracy using a model with non-rectangular meshes in Cartesian coordinates along a curved surface and fluctuating interfaces in the TTI media.


Introduction
Seismic anisotropy commonly exists in the Earth's subsurface media (Tsvankin et al. 2010;Takanashi and Tsvankin 2012). Accurate waveform simulation in tilted transversely isotropic (TTI) media is of importance for seismic waveform inversion. The latter reconstructs the subsurface velocity model quantitatively based on seismic waveform data. Seismic waveform data routinely recorded by hydrocarbon explorations comprise of mainly P-wave reflection data. Therefore, we use an acoustic wave equation in this paper for waveform simulation through TTI media containing a curved surface and fluctuating interfaces. This waveform simulation scheme is a core engine employed by the iterative inversion of seismic P-wave reflection data (Wang 2003;Wang and Rao 2009).
Although the P-wave and S-wave are coupled in the elastic wave equation in TTI media, we can have a pseudo-acoustic wave equation if the S-wave velocity is fixed along the axis of symmetry (Alkhalifah 1998;Fletcher et al. 2009). This acoustic wave equation is defined by two anisotropic parameters e and d measuring the difference between two axes of the elliptic wavefront and the deviation from a perfect elliptical shape, respectively (Thomsen 1986). In the context of seismic waveform inversion, Pratt and Shipp (1999) and Rao and Wang (2011) use an acoustic wave equation defined by a single anisotropic parameter e. Rao et al. (2016) provide the derivation of this wave equation with a single anisotropic parameter. In this paper, we adopt the acoustic wave equation with two anisotropic parameters. That is the pseudo-acoustic wave equation.
We also consider the models with a curved surface and fluctuating interfaces in the subsurface. We use quadrangle grids to partition every individual layer, confined by two fluctuating interfaces, based on a body-fitting scheme. By solving Poisson's equation, these grids also satisfy the pseudo-orthogonal condition (Rao and Wang 2013) in which grids should have the acute angles [ 67°(90°for completely orthogonal grids). We transform non-rectangular meshes in Cartesian coordinates into rectangular ones through conformal mapping. These grids, as structured, keep the similar neighborhood relationships as they do in Cartesian coordinates. We reformulate the pseudo-acoustic wave equation accordingly in the computational space. For the reformulated TTI wave equation, we analytically derive the corresponding numerical dispersion relation and stability condition, which provide the basis of finite-difference parameter selection of the TTI media waveform simulation.

Wave Equation
The pseudo-acoustic wave equation in TTI media is (Fletcher et al. 2009) where p is the P wavefield, q is an auxiliary wavefield, v pz is the P-wave velocity along the TTI symmetry axis, v px is the P-wave velocity perpendicular to the symmetry axis, v pn is the P-wave normal-moveout velocity, v sz is the SV-wave velocity along the TTI symmetry axis, and H x and H z are two 2D differential operators, H x o 2 =ox 2 and H z o 2 =oẑ 2 , with respect to the rotated coordinate system (x,ẑ). Denoting the anticlockwise rotation angle of the TTI symmetry axisẑ by / (Fig. 1), we can present the two 2D differential operators as where (x, z) are Cartesian coordinates. Note that we use H x and H z here to present differentials with respect tox andẑ, respectively. Fletcher et al. (2009) used H 1 and H 2 to present the two differentials with respect toẑ andx, respectively. For the completeness, we summarize the derivation of this pseudo-acoustic wave equation in the appendix.
In the pseudo-acoustic wave equation [Eq.
(1)], the two anisotropic velocities (v px , v pn ) are related to the TTI anisotropy parameters by , where e and d are two anisotropy parameters (Thomsen, 1986). Note that the SVwave velocity v sz in Eq. (1) does not have a significant effect on the P wavefront. It would have an effect on the SV-wave wavefront and in turn on the P-S converted wave imaging. However, from the P-P wave imaging/gradient calculation point of view, the SV-wave wavefront (traveling slowly) is simply an unwanted artefact and does not significantly affect the P-wave image.
We use body-fitted grids initially to partition the numerical model (Rao and Wang 2013). Because of the curved surface and fluctuating interfaces in the model, there must be non-rectangular meshes along the surface and the interfaces (Fig. 2a). Then, we map this initial Cartesian model into a computing model with a flat surface and flat interfaces (Fig. 2b). After this conformal mapping, the non-rectangular meshes in the Cartesian model become rectangular meshes in the computational space. Thus, we can adopt a simple finite-difference scheme straightforwardly for wave simulation.
Here, we transform the pseudo-acoustic wave equation from Cartesian coordinates (x, z) to the computational space (n; g). In the computational space, the first-order spatial derivatives become Figure 1 In homogeneous TTI media, the anticlockwise rotation angle of the symmetry axis of the wavefront is /: The ray vector points from the source to the receiver, and the corresponding wave (phase velocity) vector is normal to the wavefront. The wave angle with respect to the symmetry axisẑ is h: The angle between the wave vector and the positive direction of the g -axis in the computational space is c. The rotation angle of a local mesh from the conventional finite-difference mesh, which is parallel and perpendicular to the Cartesian coordinates, is u Y. Rao and Y. Wang Pure Appl. Geophys. where _ n x on=ox, _ g x og=ox, and so on. These derivatives are known parameters, because of the relationship between the model coordinates and the computational coordinates. The second-order spatial derivatives are Substituting these spatial derivatives into the 2D differential operators H x and H z , we rewrite wave equation [Eq. (1)] in the computational space (n; g).
We numerically solve this time-space domain wave equation in the computational space using a secondorder finite-difference method.

Numerical Dispersion
We analyze the numerical dispersion in this section, for the finite-difference implementation of the spatial derivatives in the pseudo-acoustic wave equation.
When e ¼ d; the P wavefield p is approximately equal to the auxiliary wavefield q; and we can use any of the two expressions in Eq. (1) for the dispersion analysis. We rewrite the first expression as for which we exploit the relation of v 2 px ¼ ð1 þ 2eÞv 2 pz . Assuming that the grids fully satisfy the orthogonality, and that the local mesh has a rotation of angle u from the conventional finite-difference mesh: We can express the spatial derivatives in Eq. (4) as Substituting these expressions into the differential operators H x and H z in Eq. (2), then we rewrite Eq. (5) as where A ¼ ð1 þ 2e cos 2 /Þ cos 2 u þ ð1 þ 2e sin 2 /Þ sin 2 u À e sin 2/ sin 2u ; B ¼ À2e sinð2/ þ 2uÞ; C ¼ ð1 þ 2e cos 2 /Þ sin 2 u þ ð1 þ 2e sin 2 /Þ cos 2 u þ e sin 2/ sin 2u : Figure 2 a Body-fitted grids in the physical space (x, z). b Non-rectangular meshes in the Cartesian coordinates are transformed into rectangular ones in the computational space (n; g) through conformal mapping

Dispersion and Stability Condition of Seismic Wave Simulation in TTI Media
We can approximate the temporal and spatial derivatives in Eq. (8) by finite differencing. Following the procedure of von Neumann stability analysis, we insert a plane wave pðn; g; tÞ ¼ p 0 e Àiðk n nþk g gÞ e ixt into the finite-difference operators (Nilsson et al. 2007). For instance, Then, we approximate the four derivatives in Eq. (8) as harmonic functions: Now, assuming Dn ¼ Dg; we express Eq. (8) as A Dn 2 sin 2 k n Dn 2 þ B 4Dn 2 sinðk n DnÞ sinðk g DgÞ where x ¼ 2pv h s=Dn is the angular frequency with velocity of wave propagation v h ; following De Basabe and Sen (2007), s ¼ Dn=k represents the reciprocal number of sampling points per wavelength k; k n ¼ 2p sin c=k and k g ¼ 2p cos c=k are wavenumbers, and c is the angle between the wave vector and the positive direction of the g-axis in the computational space (Fig. 1). Defining the Courant number as we rewrite Eq. (12) as where K ¼ A sin 2 ðps sin cÞ þ B 4 sinð2ps sin cÞ sinð2ps cos cÞ þ C sin 2 ðps cos cÞ: ð15Þ In the appendix, we have derived the P-wave phase velocity vðhÞ as where L ¼ 1 À v 2 sz =v 2 pz , and h is the angle between the wave vector and the TTI symmetry axis (Fig. 1).
As we assume e ¼ d, and assuming L ¼ 1, the maximum possible value of L, we can simplify Eq. (17) where Therefore, we obtain the following dispersion relation: Figure 3 shows that the dispersion increases with the increase of the sampling rate s, which is the reciprocal number of sampling points per wavelength. In this test, the Courant number Q ¼ 0:1, the anisotropy parameter e ¼ 0:1, the angle between the symmetry axis and the vertical direction of the TTI medium is / ¼ 30 ; and the meshing rotation angle u ¼ 30 : Numerical tests indicate that the dispersion is weakest when the propagation direction is perpendicular to the axis of TTI media (h ¼ 90 ). The dispersion is slightly increased when the propagation direction is close to the axis of symmetry (h ! 0 ).
As it is well known, a general requirement for forward modeling in the isotropic media with the second-order finite-difference simulations is 10 grids per wavelength (Kreiss and Oliger 1972;Wu et al. 1996). Figure 3 indicates that satisfying the condition s 0:1, that is ! 10 grids per wavelength, means that the dispersion is less than 1% in wavefield simulation through anisotropic media. The smallest dispersion ratio is 0.99, for the anisotropic parameters e and d between À0:25 and 0:25 and for the mesh rotation angle u between 0 and 90 o . These numerical tests confirm that, as long as s ¼ 0:1, the dispersion requirement of forward modeling is still satisfied.
Note that the points per wavelength quoted above is for wave propagation within one wavelength. Because of the error accumulation, the number of points required for a given tolerance depends on the number of wavelengths the wave is to propagate. According to Kreiss and Petersson (2012), a typical scaling is ðk=lÞ 1=n , where (k, l) here are the Lame parameters, and n is the order of accuracy of the method.

Numerical Stability
We derive the stability condition in this section, for the numerical implementation of the time derivatives in the wave equation, given appropriate sampling in the spatial domain. We investigate this numerical stability when the spatial grids are discretized in Cartesian coordinates.
Assuming e ¼ d; a similar procedure to that in the previous section, we approximate the auxiliary wavefield q by the P wavefield p, and obtain the condition of numerical stability from the implementation of Eq. (5). Using the second-order finitedifferencing operator to present the time-domain derivative, we express Eq. (5) as where n is the time step. After Fourier transformation, it becomes whereH x andH z are the wavenumber-domain 2D differential operators, andp is the wavenumber-domain P wavefield. Now, we rewrite Eq. (22) in the form of a growth matrix as Denoting the 2 Â 2 growth matrix by A, the stability condition of this differential equation is that the maximum absolute eigenvalue of matrix A is not greater than one. This maximum value defines the spectral radius of the matrix. Thus, the stability condition is The eigen-equation of matrix A is which has two eigenvalues: with D ¼ f v 2 pz Dt 2 ½ð1 þ 2eÞH x þH z À 2 g 2 À 4: ð27Þ When D [ 0; both eigenvalues are real valued, and RðAÞ is larger than one. When D 0, both solutions k 1;2 are complex valued. The moduli of both the complex eigenvalues k 1;2 are identically equal to one. Therefore, to guarantee a stable iteration, the condition of D 0 must be satisfied. That is, According to Eq. (2), the wavenumber-domain 2D differential operatorsH x andH z may be expressed as The wavenumbers in Eq. (29) can be approximated as harmonic functions in the wavenumber domains, similarly to Eq. (11). Denoting the wavenumber-domain wavefield as pð'; m; tÞ !p 0 ðk x ; k z ; tÞ; pð' þ 1; m; tÞ !p 0 ðk x ; k z ; tÞe Àik x Dx ; pð'; m þ 1; tÞ !p 0 ðk x ; k z ; tÞe Àik z Dz ; ð30Þ the wavenumber factors can be approximated by second-order finite-differencing as Then, the 2D differential operatorsH x andH z are rewritten as where sin 2/ [ 0 is assumed, and k h is either k x or k z , whichever Dx or Dz is the smallest ð¼ hÞ. The denominator will be maximal only when cosðk h hÞ ¼ À1: Hence, setting k h ¼ AEp=h; the Nyquist wave number, and replacing v pz with v max , we express the stability condition finally as where v max is the maximum value of the P-wave velocity along the symmetry axis. Note that the condition of e ! d commonly exists and expression (34) is also applicable to parameter d: Expression (34) is the stability condition for the OðDt 2 ; h 2 Þ scheme, the second-order finite-differencing in both time and space. For selecting the time step using this expression, we set h as the smallest cell size of non-rectangular grids generated by body fitting. It is worth mentioning that, setting e ¼ 0, the stability condition will reduce to the well-known formula of an OðDt 2 ; h 2 Þ scheme in isotropic media (Lines et al. 1999). Figure 4 demonstrates the waveform simulation in TTI media. This is a homogeneous model (Fig. 4a) with a constant P-wave velocity and constant anisotropy parameters. The objective here is to show that the wave simulation method is able to generate an accurate wavefield, even for a model discretized with non-rectangular grids. The P-wave velocity is v pz ¼ 2000 m/s, the S-wave velocity is v sz ¼ v pz = ffiffi ffi 3 p m/s, the dip angle of the symmetry axis is / ¼ 45 ; and the two anisotropy parameters are e ¼ 0:24 and d ¼ 0:1. Figure 4a shows that the body-fitted grids coincide well with the curved surface and a fluctuating interface at the middle of the model, and, meanwhile, there are non-rectangular meshes unavoidably existing on either side of the interface. As long as the media are discretized by body-fitted grids, there is no explicit enforcement of the normal stress and displacement continuity conditions at the interface. In order to avoid any instability caused by the low meshing precision, we adapt a summation-by-parts (SBP) finite-difference method (Kreiss and Scherer 1974;Nilsson et al. 2007) to the case here with the cell size variation of body-fitted grids. We use the finite-difference operators with the second-order accuracy in both temporal and spatial directions. According to Sjögreen and Petersson (2012), the SBP method is automatically stable for partial differential equations with the second-order spatial derivatives. The pseudo-acoustic wave equation is a system of partial differential equations with an initial condition which defines the source signature. For the absorbing boundary, we use the perfectly matched layer method, presented in the computational space, and discretized by finite differencing with a second-order accuracy (Rao and Wang 2013).

Waveform Simulation in TTI Media
The snapshots (Fig. 4b, c) indicate that the nonrectangularity in the mesh does not have visible effect deteriorating wavefield simulation, once the numerical dispersion relation and the stability condition are satisfied. This example demonstrates that this numerical simulation method for a model with such irregular grids is able to produce an accurate wavefield, without any artificial reflections from the interface. In contrast, if using a standard finite-difference method, there must be some artificial reflections caused by strong variation in the cell sizes of body-fitted grids, and even the two-layer velocities are assumed to be constant.
We measure the amplitudes at a number of time samples, and compare them to the theoretical amplitudes as shown in Fig. 5. Measuring along the directions parallel and perpendicular to the TTI symmetry axis (/ ¼ 45 ), the theoretical amplitude is ffiffiffiffiffiffiffi ffi t 0 =t p , where t 0 ¼ 100 ms is used in this example. Normalized amplitudes measured at time steps t ¼ 100, 150, 200, 250, 300 ms show minor errors, which probably include picking errors. Figure 6 demonstrates wave simulation in a twolayer TTI media. The P-wave velocities of two layers are 2500 and 3000 m/s. The S-wave velocity is v sz ¼ v pz = ffiffi ffi 3 p m/s. The rest of the parameters are the same as that used in Fig. 4. The two anisotropy parameters are constant for two layers, e ¼ 0:24 and d ¼ 0:1. The dip angle of the symmetry axis is / ¼ 45 . The body-fitted grid, plotted by each set of 10 grids, is a well-coincided irregular surface and a Wave simulation in TTI media. a A homogeneous model with constant velocity and anisotropy parameters. The body-fitted grid, plotted by each set of 10 grids, is a well-coincided irregular surface and a curved interface at the middle of model (red curves). b, c The snapshots of the wavefield at 200 ms and 450 ms, respectively, demonstrate that the numerical simulation method with nonrectangular grids is able to produce an accurate wavefield, without any artificial reflections, by strong variation in the cell sizes of body-fitted grids curved interface at the middle of the model. The snapshots of the wavefield at 350 and 480 ms clearly show reflections from the fluctuating interface. Note that the geometrical configuration of Figs. 4a and 6a is the same. In waveform simulation, we always implement two steps. In the first step, we use a model with a constant velocity for all layers, such as Fig. 4a, to check the effectiveness of grids. In the second step, we use the proper model with layered velocities, such as Fig. 6a, to generate the desired waveform.

Conclusions
For seismic waveform simulation with the TTI wave equation, we have derived a numerical dispersion relation, and demonstrated that the number of sampling points per wavelength has the greatest influence on the dispersion, and both the anisotropic parameters and the mesh rotation angle have little influence on the dispersion. Further, we have derived the condition for the numerical stability for the selection of time step, given the spatial sampling is appropriately selected in Cartesian coordinates. In the derivation, we have used the second-order finite-differencing operators in both time and space and assumed an elliptical anisotropy. Appendix A: P-wave Phase Velocity The wave equation for the homogeneous anisotropic media, following Newton's second law, can be represented as (Tsvankin 2001 where s ij are the elements of the stress tensor, x j are Cartesian coordinates, u i are the components of the displacement vector, t is time, and q is density. According to the generalized Hooke's law, the stress tensor is linearly related to the strain tensor by s ij ¼ c ijm' e m' , where e m' are the elements of the strain tensor, and c ijm' are the elastic constants in the stiffness tensor. The elements of the strain tensor are defined by e m' ¼ 1 2 ðou m =ox ' þ ou ' =ox m Þ. Assuming the elastic constants c ijm' are constant in the space, equation (A1) can be written as In this wave equation, the media anisotropy property is accounted for by the stiffness tensor c ijm' .