A hybrid eikonal solver for accurate first-arrival traveltime computation in anisotropic media with strong contrasts
AA hybrid eikonal solver for accurate first-arrivaltraveltime computation in anisotropic media withstrong contrasts
Kai Gao ∗ ,1 and Lianjie Huang Los Alamos National Laboratory, Geophysics Group, Los Alamos, NM 87545,USA
Abstract
First-arrival traveltime computation is crucial for many applications such as traveltime to-mography, Kirchhoff migration, etc. There exist two major issues in conventional eikonalsolvers: the source singularity issue and insufficient numerical accuracy in complex media.Some existing eikonal solvers also exhibit the stability issue in media with strong contrasts inmedium properties. We develop a stable and accurate hybrid eikonal solver for 2D and 3Dtransversely isotropic media with a tilted symmetry axis (TTI, or tilted transversely isotropicmedia). Our new eikonal solver combines the traveltime field factorization technique, thethird-order Lax-Friedrichs update scheme, and a new method for computing the base travel-time field. The solver has the following three advantages. First, there is no need to assignexact traveltime values in the near-source region, and the computed traveltime field near thesource location is accurate even for TTI media with strong anisotropy. Second, the computedtraveltime field is high-order accurate in space. Third, the solver is numerically stable for 2Dand 3D TTI media with strong anisotropy, complex structures, and strong contrasts in mediumproperties. We verify the stability and accuracy of our hybrid eikonal solver using several 2D ∗ Corresponding Author; Email: [email protected] a r X i v : . [ phy s i c s . c o m p - ph ] A ug nd 3D TTI medium examples. The results show that our solver is stable and accurate in 2Dand 3D complex TTI media, producing first-arrival traveltime fields that are consistent withfull-wavefield solutions. Keywords: anisotropic media, first-arrival traveltime, eikonal equation, strong medium prop-erty contrasts
Traveltime computation is important for many applications, including underwater acoustics(Martinelli, 2012), geometrical optics (Qian and Leung, 2006), quantum mechanics (Jin et al.,2005), geophysics, etc. Many geophysical applications, such as Kirchhoff migration (Gray andMay, 1994; Buske, 1999) and first-arrival traveltime tomography (Lin et al., 2009; Taillandieret al., 2009) for reconstructing subsurface structures and medium properties, rely on accurate andefficient traveltime computation.There exist roughly two categories of numerical methods for traveltime computation: ray-based methods and eikonal-equation-based methods. Ray approaches are based on ray equationsapproximated from wave equations. A ray-tracing system is then solved in the framework of eitherone-point initial value problem or two-point boundary value problem using different techniques(Pereyra et al., 1980; Grechka and McMechan, 1996; Sadeghi et al., 1999; Mel´endez et al., 2015).These approaches are generally very efficient for sparse source and receivers, but the computationalcosts may increase dramatically as the number of source-receiver pairs increases. In addition, raytracing cannot trace rays to cover the entire model, and cannot handle complex media with strongcontrasts. A more flexible ray tracing approach is the so-called the wavefront construction (WFC)method (Vinje et al., 1993; Lambar´e et al., 1996; Gibson et al., 2005; Chambers and Kendall,2008), which computes traveltime and amplitude fields in a more intuitive manner by dynamicallyinserting new rays where necessary. WFC generally requires more programming efforts and ismore computationally demanding compared with conventional ray-based approaches.Vidale (1988) first developed an eikonal-equation-based approach to generating first-arrivaltraveltimes in heterogeneous media. His method is also known as the expanding-box method.Eikonal-equation-based approaches have since gained fast development and wide applications.2urrently, the most widely used approaches include expanding wavefront methods (Podvin andLecomte, 1991; Qin et al., 1992), fast marching methods (Sethian and Popovici, 1999; Rawlinsonand Sambridge, 2004; Zhang et al., 2006), fast sweeping methods (Tsai et al., 2003; Zhao, 2004;Kao et al., 2005; Fomel et al., 2009; Luo and Qian, 2011; Waheed et al., 2015b), etc. Eikonalsolvers on triangular or unstructured mesh (e.g., Qian et al., 2007; Le Bouteiller et al., 2019) canhandle complex interfaces and domain boundaries. Various high-order and non-oscillatory nu-merical schemes (Kim and Cook, 1999; Kim, 2002; Luo and Qian, 2011; Luo et al., 2012) canimprove the numerical accuracy and stability of eikonal solvers in complex media with strong con-trasts. Different eikonal solvers have various computational complexities and numerical accuracy.A comparison on several popular eikonal solvers can be found in G ´omez et al. (2019). One of themost distinct advantage of eikonal solvers compared with ray-based methods is that the output ofan eikonal solver is a first-arrival traveltime field in the entire computational domain, as opposedto ray tracing methods that compute only traveltimes on ray paths. In addition, eikonal solversgenerally allow models to be arbitrarily heterogeneous and complex, whereas ray-based methodsusually require simple or smooth media.Finite-difference eikonal solvers were first developed for isotropic media. There is an intensiveneed for eikonal solvers in anisotropic media. Eikonal equations for various kinds of anisotropicmedia are often significantly more complex than that in isotropic media. Most of the eikonalsolvers in isotropic media require substantial modifications for anisotropic media if possible. Inanisotropic media, the phase and group velocity directions are generally not the same (Carcione,2015). Therefore, eikonal solvers for anisotropic media require sophisticated numerical schemesfor updating traveltime fields. Eaton (1993) developed a high-order expanding-wavefront methodon a hexagonal mesh to compute qP-, qSV- or qSH-wave traveltime fields for transversely isotropic(TI) media. Qian and Symes (2002) developed a paraxial eikonal equation system to compute theqP-wave traveltime field in TI media. Wang et al. (2006) developed an unconditionally-stableexpanding wavefront method for the eikonal equation in TI media. Their method explicitly tracksgroup velocity propagation directions to ensure correct causal stencils. Waheed et al. (2015a)developed a method based on perturbation expansion to solve the TTI eikonal equation. Thereexist several other methods for solving eikonal equations in anisotropic media by assuming weak,elliptical anisotropy (Ettrich and Gajewski, 1998), by perturbing from elliptical reference medium3Soukina et al., 2003), or with only low-order numerical accuracy (Lecomte, 1993).The source singularity is a major problem in various eikonal solvers. Conventional numericalschemes for solving eikonal equations are based on local plane wave assumption for traveltime fieldupdate, and therefore cannot accurately handle large curvatures of the traveltime field around thesource point. The numerical error in the near-source region can eventually deteriorate the overallnumerical accuracy in the entire computational domain. The source singularity issue is generallysolved with the traveltime factorization method. The total traveltime field is factorized into anaddition or multiplication of a base traveltime and an additional or multiplicative traveltime field,i.e., t = t + τ or t = t × τ , where t is the total traveltime field, t is the base traveltime field, and τ is the additional or multiplicative traveltime field (Luo and Qian, 2011, 2012; Luo et al., 2012). Thebase traveltime field is solved in a homogeneous isotropic or elliptically anisotropic medium, whilethe traveltime field τ is computed in the heterogeneous part of the medium. A notable approachbased on the traveltime factorization is an iterative scheme developed by Waheed et al. (2015b)and Waheed and Alkhalifah (2017) to solve the eikonal equation in TTI media. Their method firstdecomposes the left-hand side of the eikonal equation into a tilted elliptically anisotropic term andan additional term, and moves the additional term to the right-hand side of the eikonal equation.During each iteration, their method updates the right-hand side term to approximate the true TTIeikonal equation.We develop a hybrid numerical scheme based on both the monotonic Godunov scheme and thehigh-order weighted essentially non-oscillatory (WENO) scheme to solve the eikonal equation in2D and 3D anisotropic media. Our hybrid eikonal solver has three advantages. First, our solveris free of source singularity issue, and there is no need to assign traveltime around the source byusing multiplicative factorization of the traveltime field. Second, the computed traveltime field ofour solver is high-order accurate in space by using the third-order Lax-Friedrichs WENO scheme.Third, our solver is numerically stable for 2D and 3D TTI media with strong anisotropy, highlycomplex structures, and strong medium property contrasts, by using the weighted non-oscillatoryscheme. Our eikonal solver employs both the conventional first-order Godunov scheme and thethird-order Lax-Friedrichs scheme to achieve numerical stability and high-order accuracy. There-fore, we call our solver a hybrid approach. To our knowledge, our hybrid eikonal solver is the firstfast-sweeping-based method to date that simultaneously holds these three advantages.4ur paper is organized as follows. In the Methodology section, we describe the three computa-tional steps of our anisotropic eikonal solver, including the first-order Godunov locking-sweepingstep, the base traveltime field computation step in the arbitrary TTI medium, and the third-orderLax-Friedrichs fast sweeping step. We then use several numerical examples to verify the stabil-ity and accuracy of our method. In the Conclusions section, we summarize the most importantfeatures of our new anisotropic eikonal solver. We derive the formulation for our hybrid eikonal solver in 2D anisotropic media without lossof generality. We give the 3D formulation in Appendix A.We adopt the following eikonal equation in TTI media (Waheed et al., 2015b): v x (cid:18) ∂t∂ ˆ x (cid:19) + v z (cid:18) ∂t∂ ˆ z (cid:19) (cid:34) − ε − δ ) v z (cid:18) ∂t∂ ˆ x (cid:19) (cid:35) = 1 , (1)where v x ( x ) = V p ( x ) (cid:112) ε ( x ) is the qP-wave velocity along the x -axis, v z ( x ) = V p ( x ) is theqP-wave velocity along the z -axis (i.e., the anisotropy symmetry axis), ε = ε ( x ) and δ = δ ( x ) areThomsen parameters describing anisotropy properties of a VTI medium. Equation (1) is written inthe rotated coordinates ˆ x − ˆ z , and the spatial derivatives are combinations of the spatial derivativesin the unrotated coordinates x − z : ∂t∂ ˆ x = a x ∂t∂x + c x ∂t∂z , (2) ∂t∂ ˆ z = a z ∂t∂x + c z ∂t∂z , (3)with the coordinate transformation matrix R written as R = a x c x a z c z = cos θ sin θ − sin θ cos θ , (4)where θ = θ ( x ) is the tilt angle of a VTI medium’s symmetry axis (i.e., the counterclockwiserotation angle of the symmetry axis with respect to the y -axis).For notation clarity and derivation convenience, we further define α x = v x a x = v x cos θ, (5)5 x = v x c x = v x sin θ, (6) α z = v z a z = − v z sin θ, (7) γ z = v z c z = v z cos θ, (8)and ξ = 2( ε − δ )1 + 2 ε , (9)leading to (cid:18) ∂t∂ ˆ x (cid:19) + (cid:18) ∂t∂ ˆ z (cid:19) (cid:34) − ξ (cid:18) ∂t∂ ˆ x (cid:19) (cid:35) = 1 , (10)with ∂t∂ ˆ x = α x ∂t∂x + γ x ∂t∂z , (11) ∂t∂ ˆ z = α z ∂t∂x + γ z ∂t∂z . (12)Generally, v x (cid:54) = v z and θ (cid:54) = 0 in TTI media. The first step of our hybrid eikonal solver is to solve equation (10) using an iterative first-orderfast sweeping method. The iterative approach is based on rewriting equation (10) in the form of (cid:34) − ξ (cid:18) ∂t∂ ˆ z (cid:19) (cid:35) (cid:18) ∂t∂ ˆ x (cid:19) + (cid:34) − ξ (cid:18) ∂t∂ ˆ x (cid:19) (cid:35) (cid:18) ∂t∂ ˆ z (cid:19) = 1 . (13)By setting two coefficients c x = 1 − ξ (cid:18) ∂t∂ ˆ z (cid:19) , (14) c z = 1 − ξ (cid:18) ∂t∂ ˆ x (cid:19) , (15)and absorbing them into ∂t/∂ ˆ x and ∂t/∂ ˆ z , we have (cid:18) ∂t∂ ˜ x (cid:19) + (cid:18) ∂t∂ ˜ z (cid:19) = 1 , (16)with ∂t∂ ˜ x = (cid:112) | c x | α x ∂t∂x + (cid:112) | c x | γ x ∂t∂z , (17)6 t∂ ˜ z = (cid:112) | c z | α z ∂t∂x + (cid:112) | c z | γ z ∂t∂z . (18)The left-hand side of equation (16) is the eikonal correspondence in an elliptically transverselyisotropic medium, and can be solved using the first-order Godunov scheme.We solve the eikonal equation (16) using the following procedure. We first set c x ( x ) = c z ( x ) =1 , and solve equation (16) using fast sweeping; Then we update c x ( x ) and c z ( x ) using the com-puted traveltime t based on equations (14) and (15), and solve equation (16) again. The reason forsuch a reformatting is that a Godunov scheme for equation (10) can be fairly difficult to derive, andpossibly leads to high computational complexity for local solvers. Instead, a first-order Godunovscheme for the left-hand side of equation (13) (or equation (16)) is usually simple to derive (e.g.,Tsai et al., 2003).We use a locking-sweeping procedure to reduce computational costs. At each fast sweepingiteration, we simply lock the points where t ( x ) ( m ) = t ( x ) ( m − where the superscript ( m ) repre-sents the m -th round of fast sweeping, and update traveltime values only at unlocked points. Thecriterion t ( x ) ( m ) = t ( x ) ( m − might not be as accurate as the one given in Bak et al. (2010) andG´omez et al. (2019) based on checking the changes of neighbor points of a certain spatial point, butis much more efficient to compute. In practice, we find that even for complex media, this criterionresults in accurate traveltime fields. The locking-sweeping procedure can significantly reduce thecomputational costs for the first step. In Steps II and III, we express the traveltime field t = t ( x ) in heterogeneous TTI media usinga multiplicative factorization (Luo and Qian, 2012; Waheed and Alkhalifah, 2017) as t ( x ) = t ( x ) τ ( x ) , (19)where t ( x ) is the base traveltime field in the homogeneous TTI media, and τ ( x ) is the multiplica-tive traveltime field that accounts for heterogeneities of the model.This multiplicative factorization results in ∂t∂ ˆ x = α x (cid:18) t ∂τ∂x + ∂t ∂x τ (cid:19) + γ x (cid:18) t ∂τ∂z + ∂t ∂z τ (cid:19) , (20)7 t∂ ˆ z = α z (cid:18) t ∂τ∂x + ∂t ∂x τ (cid:19) + γ z (cid:18) t ∂τ∂z + ∂t ∂z τ (cid:19) , (21)which transform equation (10) into a factorized eikonal equation for TTI media.For the factorized eikonal equation, the base traveltime field t and its first-order spatial deriva-tives, say, t x = ∂t /∂x and t z = ∂t /∂z , are computed analytically, and are fixed during fast-sweeping iterations (Luo and Qian, 2012; Waheed and Alkhalifah, 2017). The base traveltime fieldis essential to avoid the source singularity issue in eikonal-equation-based traveltime computationwithout specifying exact traveltime values at the points near the point source location.Existing techniques for solving the factorized isotropic or anisotropic eikonal equation assumethat the background media is either isotropic where ε = δ = 0 (Fomel et al., 2009; Luo andQian, 2011; Luo et al., 2012) or elliptically anisotropic where ε = δ (cid:54) = 0 (Luo and Qian, 2011,2012; Waheed and Alkhalifah, 2017), because there exist closed-form expressions for computingthe traveltime and its spatial derivatives in the case of ε = δ .However, the requirement of ε = δ also limits the application of multiplicative factorization incomplex or strong anisotropic media. For instance, to apply multiplicative factorization to the casewhere ε (cid:54) = δ , Waheed and Alkhalifah (2017) had to use an iterative scheme to update the right-handside of the eikonal equation, and update the analytic t , t x and t z after several iterations. Evenwith such an iterative scheme, in each iteration, the background anisotropic medium is assumed tobe elliptically anisotropic. Therefore, the base traveltime field t never truly approximates that inthe anisotropic media where ε (cid:54) = δ . For strongly anisotropic media, the group velocity curve cansignificantly deviate from an ellipse. The resulting base traveltime field might significantly differfrom the true solution, eventually leading to suboptimal traveltime solutions, even though it helpsavoid the source singularity issue.In contrast to existing techniques where t is computed for an elliptically anisotropic medium,we develop a semi-analytic approach to directly computing t for anelliptically anisotropic mediumwhere ε (cid:54) = δ . Therefore, the background anisotropic medium in our method is the non-degeneratedTTI medium at the reference point. This is the most important difference between our hybridmethod and existing factorized eikonal solvers.In the following derivations, to distinguish the rotation angle θ of the TTI symmetry axis inthe following description, we use ϑ to represent the phase velocity angle in a VTI medium, which8easures the deviation angle from the vertical symmetry axis of a VTI medium. The phase velocityof the qP-wave in a VTI medium can be written as (Tsvankin, 2012) v phase ( ϑ ) = V p (cid:118)(cid:117)(cid:117)(cid:117)(cid:116) ε sin ϑ − f − (cid:115)(cid:18) ε sin ϑf (cid:19) − ε − δ ) sin ϑf , (22)where f = 1 − V s /V p , and V p and V s are qP- and qS-wave velocities along the symmetry axis,respectively. By setting V s = 0 , we have v phase ( ϑ ) = V p (cid:114)
12 + ε sin ϑ + (cid:113)(cid:0) ε sin ϑ (cid:1) − ε − δ ) sin ϑ. (23)Note that there is no approximation in equation (23), and therefore it is accurate even for stronganisotropy. We also have the group velocity angle ψ in terms of the phase velocity angle ϑ as(Tsvankin, 2012) ψ ( ϑ ) = arctan tan ϑ + 1 v phase ( ϑ ) dv phase ( ϑ ) dϑ − tan ϑv phase ( ϑ ) dv phase ( ϑ ) dϑ , (24)and the magnitude of the group velocity in terms of the phase velocity angle ϑ as v group ( ϑ ) = v phase ( ϑ ) (cid:115) v phase ( ϑ ) dv phase ( ϑ ) dϑ , (25)where based on equation (23), we have dv phase ( ϑ ) dϑ = V p sin 2 ϑ v phase ( ϑ ) ε (cid:0) ε sin ϑ (cid:1) − ε − δ ) cos 2 ϑ (cid:113)(cid:0) ε sin ϑ (cid:1) − ε − δ ) sin ϑ + ε . (26)For an arbitrary spatial point x = ( x, z ) in the computational domain, the group velocity angleat this point can be evaluated as ψ = arctan x − x z − z , (27)where x = ( x , z ) is the position of the point source.Our goal is to compute the semi-analytic magnitude of the group velocity at this point, sothat we can find the exact first-arrival traveltime at x . This requires the determination of thecorresponding phase velocity angle ϑ at x , by which we can evaluate the magnitude of the groupvelocity at x using equation (25). Unfortunately, there is no closed-form expression to compute9 from ψ based on equation (24), because equation (24) is a complicated transcendental equationand is extremely difficult to solve analytically, if not impossible.We therefore adopt a numerical method to compute the magnitude of group velocity v group at x . We first compute a series of group velocity values v group ( ϑ ) , v group ( ϑ ) , · · · , v group ( ϑ n ) , where ϑ , ϑ , · · · , ϑ n is an equal division of the phase angle range [0 , π ] . Meanwhile, we compute thecorresponding group velocity angles ψ ( ϑ ) , ψ ( ϑ ) , · · · , ψ ( ϑ n ) based on these phase velocity an-gles. In any VTI medium, the group velocity angles ψ ( ϑ ) , ψ ( ϑ ) , · · · , ψ ( ϑ n ) range exactly from [0 , π ] , but are generally not equally distributed within this range. Therefore, we obtain a seriesof group angle-velocity pairs { ψ ( ϑ i ) , v group ( ϑ i ) } with i = 1 , , · · · , n . We then use cubic splineinterpolation to obtain the group velocity for the spatial location x , which corresponds to a groupvelocity angle ψ x based on equation (27). The interpolant function of this cubic spline interpola-tion is built from the group angle-velocity pairs { ψ ( ϑ i ) , v group ( ϑ i ) } . In practical computations, weuse a large n to divide the range [0 , π/ , leading to high accuracy for cubic spline interpolation.We repeat the process until all the spatial points in the model are covered.Therefore, there is no need to analytically compute the phase velocity angle ϑ for the spatialpoint x in our numerical scheme. The group velocity value at any spatial point is obtained throughan 1D interpolation process with smooth and continuous interpolants built from the analytic groupangle-velocity pairs. As a result, the computed group velocities are practically of analytic accuracy.The computational cost associated with this part is small in the entire eikonal equation solvingprocess.Because for any VTI medium, the phase or group velocity is symmetric with respect to bothaxes, it is sufficient to build a complete group velocity profile for the entire π range based onthe computed group velocity values in [0 , π/ . In addition, because any TTI medium is simplya coordinate rotation result of some VTI medium, the group velocity values for the TTI mediumcan be easily computed using the scheme described above. Assume that the tilt angle of a TTImedium is θ , then for a normalized spatial location x = ( x − x , z − z ) where ( x , z ) is thesource location, the corresponding directional vector in the unrotated coordinate is x (cid:48) = x (cid:48) z (cid:48) = cos θ sin θ − sin θ cos θ x − x z − z , (28)10hich indicates the group velocity angle corresponding to x should be ψ (cid:48) = arctan | x (cid:48) || z (cid:48) | . (29)We take absolute value in equation (29) to ensure that the angle ψ (cid:48) falls in [0 , π/ .We then compute the group velocity value v group | ψ (cid:48) for x by the aforementioned interpolationprocedure at ψ (cid:48) , and the traveltime at x is t ( x ) = (cid:112) ( x − x ) + ( z − z ) v group | ψ (cid:48) . (30)Finally, we use a high-order centered finite-difference scheme to compute the spatial derivativesof the base traveltime field t x and t z : t x ( i, j ) = 1∆ x M (cid:88) l =1 c l [ t ( i + l, j ) − t ( i − l, j )] , (31) t z ( i, j ) = 1∆ z M (cid:88) l =1 c l [ t ( i, j + l ) − t ( i, j − l )] , (32)where c l are finite-difference coefficients, M is the half length of the finite-difference operator(Fornberg, 1988), and ∆ x and ∆ z are the grid sizes in the x - and z -directions, respectively. Inour computation, we use M = 10 , and compute the associate coefficients c l using the proceduredescribed in Fornberg (1988). Again, these fields are practically of analytic accuracy because t ispractically of analytic accuracy and M is large.In Figure 1, we show three examples for base traveltime computation in anelliptically anisotropicTI media, including a VTI medium (Figure 1a), a HTI medium (Figure 1b), and a TTI media (Fig-ure 1c), all containing strong anisotropies. All group velocity curves significantly deviate from anellipse. The results indicate that our semi-analytic approach can accurately compute the base trav-eltime fields in complex anisotropic media. The numerical scheme to compute the base traveltimefield in arbitrary TTI media is also applicable to 3D scenario as shown in Figure 2.When a model contains multiple simultaneous sources, we need to compute the base traveltimefields for all the point sources, and compute a minimum base traveltime field by finding the minimalvalue from all the base traveltime field values at each point. That is, t ( x ) = min { t s ( x ) , t s ( x ) , · · · , t s n ( x ) } , (33)11here t s i represents the base traveltime field for the i -th point source. The medium propertiesat different point source locations can be different. Figure 1d shows a simple example of basetraveltime field computation with multiple simultaneous point sources in a 2D TTI medium.The above procedure also implies that, in homogeneous media, we can obtain the semi-analyticsolution to the eikonal equation regardless of TTI anisotropy type and the number of simultaneouspoint sources, without the need of the aforementioned Step I or the Step III described below. For heterogeneous media, once the initial traveltime filed t , the base traveltime field t andits derivatives are computed, we solve the factorized eikonal equation (10) along with equa-tions (20) and (21) in TTI media using the Lax-Friedrichs scheme based on a third-order WENOdiscretization (Zhang et al., 2006). The initial multiplicative traveltime field τ is computed as τ ( x ) = t ( x ) /t ( x ) , with τ ( x , z ) ≡ . The value of τ at location ( x , z ) is kept unchanged dur-ing iterations. For multiple simultaneous point source applications, the values of τ at all the pointsource locations are kept unchanged during iterations. In contrast to existing high-order schemeswhere neighbor points of the source point should be assigned and kept unchanged during iterations,our algorithm requires only the value at the source point fixed during iterations thanks to the useof the multiplicative factorization.To facilitate our description, we define two functionals associated with the multiplicative timefield τ and its spatial derivatives τ x and τ z : F x ( τ, τ x , τ z ) = α x ( t τ x + t x τ ) + γ x ( t τ z + t z τ ) , (34) F z ( τ, τ x , τ z ) = α z ( t τ x + t x τ ) + γ z ( t τ z + t z τ ) . (35)Then the Hamiltonian for equation (10) can be written as H ( τ, τ x , τ z ) = F x ( τ, τ x , τ z ) + F z ( τ, τ x , τ z ) − ξF x ( τ, τ x , τ z ) F z ( t, τ x , τ z ) . (36)In the Lax-Friedrichs update scheme, it is necessary to compute a set of artificial viscosities.To improve numerical stability, we use the following artificial viscosities: ω x = max Ω (cid:18)(cid:12)(cid:12)(cid:12)(cid:12) ∂H∂τ (cid:12)(cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12)(cid:12) ∂H∂τ x (cid:12)(cid:12)(cid:12)(cid:12)(cid:19) , (37)12 z = max Ω (cid:18)(cid:12)(cid:12)(cid:12)(cid:12) ∂H∂τ (cid:12)(cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12)(cid:12) ∂H∂τ z (cid:12)(cid:12)(cid:12)(cid:12)(cid:19) , (38)where Ω represents the entire computational domain. Based on equation (36), we have ∂H∂τ = 2 (cid:20) (1 − ξF z ) ∂F x ∂τ F x + (1 − ξF x ) ∂F z ∂τ F z (cid:21) = 2 (cid:2) (1 − ξF z )( α x t x + γ x t z ) F x + (1 − ξF x )( α z t x + γ z t z ) F z (cid:3) , (39) ∂H∂τ x = 2 (cid:20) (1 − ξF z ) ∂F x ∂τ x F x + (1 − ξF x ) ∂F z ∂τ x F z (cid:21) = 2 t (cid:2) (1 − ξF z ) α x F x + (1 − ξF x ) α z F z (cid:3) , (40) ∂H∂τ z = 2 (cid:20) (1 − ξF z ) ∂F x ∂τ z F x + (1 − ξF x ) ∂F z ∂τ z F z (cid:21) = 2 t (cid:2) (1 − ξF z ) γ x F x + (1 − ξF x ) γ z F z (cid:3) . (41)In our algorithm, we compute the quantities | ∂H/∂t | , | ∂H/∂τ x | and | ∂H/∂τ z | at every spatialpoints of the model, and find ω x and ω z based on the maximum values of these quantities. Notethat we do not adopt the artificial viscosities defined by Luo and Qian (2012) and Luo et al. (2012),which may lead to unstable results in complex anisotropic media with large medium propertycontrasts.We then obtain the following Lax-Friedrichs scheme to update the multiplicative traveltimefield τ at the spatial grid point ( i, j ) : τ ( m +1) i,j = 1 − H (cid:16) τ ( m ) i,j , τ ( m ) , † x ; i,j , τ ( m ) , † z ; i,j (cid:17) + ω x τ ( m ) , ∗ x ; i,j + ω z τ ( m ) , ∗ z ; i,j ω x / ∆ x + ω z / ∆ z + τ ( m ) i,j , (42)where the superscripts ( m ) and ( m +1) represent the values of τ i,j at the m th and ( m +1) th sweepingiterations, respectively, and according to Zhang et al. (2006), τ † x ; i,j = 12 (cid:0) τ +0 x ; i,j + τ − x ; i,j (cid:1) , (43) τ ∗ x ; i,j = 12 (cid:0) τ +0 x ; i,j − τ − x ; i,j (cid:1) , (44)with the third-order WENO discretizations τ +0 x ; i,j = (1 − w + x ) τ i +1 ,j − τ i − ,j x + w + x − τ i +2 ,j + 4 τ i +1 ,j − τ i,j x , (45) τ − x ; i,j = (1 − w − x ) τ i +1 ,j − τ i − ,j x − w − x − τ i − ,j + 4 τ i − ,j − τ i,j x , (46)13 + x = (cid:40) (cid:20) (cid:15) + ( τ i +2 ,j − τ i +1 ,j + τ i,j ) (cid:15) + ( τ i +1 ,j − τ i,j + τ i − ,j ) (cid:21) (cid:41) − , (47) w − x = (cid:40) (cid:20) (cid:15) + ( τ i − ,j − τ i − ,j + τ i,j ) (cid:15) + ( τ i +1 ,j − τ i,j + τ i − ,j ) (cid:21) (cid:41) − , (48)where (cid:15) is a small number to avoid singularity. The expressions for τ † z ; i,j and τ ∗ z ; i,j can be analo-gously derived.Because we actually use the first-order traveltime field computed at Step I as the initial solutionfor Step III, Step III in our algorithm requires much fewer iterations to achieve accurate resultscompared with the LF-3 method that directly solves the eikonal equation from a rough or constantinitial guess. We summarize the workflow of our hybrid eikonal solver as follows:1. Compute an initial, first-order accurate traveltime field t ( x ) using the first-order Godunovalgorithm based on the locking-sweeping procedure in the following order: i = 1 , · · · , N x , j = 1 , · · · , N z , (49) i = N x , · · · , , j = 1 , · · · , N z , (50) i = 1 , · · · , N x , j = N z , · · · , , (51) i = N x , · · · , , j = N z , · · · , , (52)where i and j are indices of the finite-difference grids, and N x and N z are the number ofgrids in the model in the x - and z -directions, respectively. The order of fast-sweeping is notimportant.2. Compute the base traveltime field t ( x ) and its spatial derivatives using the semi-analyticapproach described in Step II.3. Compute the multiplicative traveltime field τ ( x ) using the numerical scheme described inStep III. The sweeping follows the order listed in equations (49)-(52).14 Numerical Results
We use five numerical examples to verify the stability and accuracy of our hybrid eikonal solverfor 2D and 3D TTI anisotropic media. We compare the results from four different methods:1. Godunov: the first-order Godunov method based on the iterative scheme developed in Wa-heed et al. (2015b).2. Factorized Godunov: the first-order Godunov method based on the iterative scheme for thefactorized eikonal equation developed in Waheed and Alkhalifah (2017).3. LF-3: the third-order direct Lax-Friedrichs method without traveltime factorization devel-oped in Luo and Qian (2012).4. Hybrid: our hybrid eikonal solver in this paper.In our tests, we adapt and program all methods to solve the eikonal equation (10). The mean-ing of “iterative” in the first-two approaches is that we need to update the right-hand side of thedegenerated TTI eikonal equation (13) during sweeping iterations. The meaning of “factorized”in the second approach is that we use the traveltime field factorization scheme. The meaning of“direct” in the Lax-Friedrichs approach is that we directly discretize equation (10) based on theLax-Friedrichs update scheme and the third-order WENO finite-difference scheme, without anydegeneration or right-hand-side iteration as in the first two approaches.In all the implementations, we do not assign exact values for the points around the source tostudy the efficacy of these methods in realistic computational tasks. In practical applications, themedia around the source can be heterogeneous, where assigning exact traveltime values can bevery difficult, if not impossible, particularly for heterogeneous, anelliptically anisotropic media. Inall the implementations, we only fix the traveltime value at the source point (i.e., where t = 0 ) oversweeping iterations. Traveltime field values at all other spatial points can change during iterations.Note that without accurate and fixed values around the point source, it can be very difficult toachieve convergence for LF-3. Therefore, we use the result from the Godunov method as the initialguess for LF-3. 15 .1 Homogeneous model Since we are able to obtain traveltime field with analytic accuracy solely using the algorithm inStep II, traveltime computation in homogeneous TTI media is trivial for our hybrid method. How-ever, it can still be challenging for conventional eikonal solvers. In the first numerical example,we first compare the traveltime field computed using our method with those obtained using threeconventional methods.The model parameters for a homogeneous TTI medium are V p = 2000 m/s, Thomsen param-eters ε = 0 . and δ = 0 . , and TTI symmetry axis tilt angle θ = π/ . Figures 3a-d show thetraveltime fields in the homogeneous TTI model computed using Godunov, factorized Godunov,LF-3 and our hybrid method, respectively. The solution computed using our method is taken as thereference solution.Visually the solutions from different methods are close to one another, except the near-sourcetraveltime field contours in Figure 3a (Godunov method), which clearly deviate from the referencetraveltime contours shown in Figure 3d (our hybrid method). We compute the differences betweenthe traveltime fields obtained using the three conventional methods and the reference solution(Figures 4a-c). Different conventional methods have different error levels. The result of the LF-3method is the least accurate, partially because we do not assign exact values in the near-sourceregion for this method. The two Godunov methods have higher accuracy compared with the LF-3method. Nevertheless, obvious errors occur at the near-source region in both solutions, and theerrors become larger with the increased distance away from the source position. In the second example, we study the convergence of our hybrid method using an isotropicconstant gradient model. We choose this model because we can analytically compute the first-arrival traveltime in such a medium.The model is 4 km in both spatial directions, with a constant gradient of velocity s ( x ) = 1 s + G · ( x − x ) , (53)where s ( x ) is the spatially variant slowness, s is the slowness at the source point x , and G isthe constant gradient of the model. The analytical first-arrival traveltime for this medium is (Fomel16t al., 2009): t ( x ) = 1 | G | arccosh (cid:18) s ( x ) s | G | | x − x | (cid:19) , (54)with arccosh( x ) = ln (cid:16) x + √ x − (cid:17) . (55)We set a point source at x = (2 . , . km in the model, with a background constant slowness s = 1 / s/m. The constant gradient is G = ( G x , G z ) = (0 . , . s − . The velocityin the model varies from 1500 m/s to 4500 m/s as shown in Figure 5a. Figure 5b depicts thecorresponding analytical traveltime field computed using equation (54).We compute the traveltime fields with different model grid spacing using the Godunov andour hybrid method in this isotropic model. We compare the relative L -norm misfit between thenumerical solutions and the analytical solution in Figure 6, demonstrating that our method is moreaccurate than the Godunov method. The convergence order of our method is approximately 3.13while that of the Godunov method is approximately 0.79. Even at the largest grid size where thenumber of grids along each direction is 10, our hybrid method is almost two orders of magnitudemore accurate than the Godunov method. We use the third numerical example to verify the numerical stability of our method in ananisotropic medium with a strong contrast as depicted in Figure 7. The background medium isan isotropic homogeneous medium with V p = 5 , m/s and ε = δ = θ = 0 . The blue regionat the center of the model indicates the location of a low-velocity strong TTI anisotropic anomalywith V p = 1 , m/s, ε = 0 . , δ = − . and θ = π/ . The model size is 3.2 km in bothdimensions. The grid size is 10 m in both directions.Figures 7a-d show the traveltime fields computed using Godunov, factorized Godunov, LF-3and our hybrid method, respectively. The Godunov solution exhibits weak instabilities around theboundaries of the TTI ball. The factorized Godunov solution in Figure 7b shows evident numer-ical instabilities. The computed traveltime field inside the TTI ball indicates that the factorizedGodunov scheme becomes unstable for this TTI medium, and the traveltime field outside of theTTI ball is therefore mostly wrong with error propagating from the inside of the ball.17igure 7c is the LF-3 solution. The solution is more stable than those computed using theGodunov methods. Figure 7d displays our hybrid solution. Similar with the LF-3 solution, ourhybrid method produces a stable solution with the help of the weighted non-oscillatory scheme.It is important to compare the eikonal equation solution with the full-wavefield solution. Wecompute the full-wavefield solution using the fully staggered-grid finite-difference method (Lisitsaand Vishnevskiy, 2010) with a high-order stencil, plus an optimal multi-axial perfectly matchedlayers (Gao and Huang, 2018). Figures 8a-d show the full-wavefield solution at 0.55 s after thesource excitation and the corresponding eikonal equation solutions in black curves computed usingthe Godunov method, factorized Godunov method, LF-3 method and our hybrid method, respec-tively. We find that the Godunov and the LF-3 methods are stable, yet are not consistent with thefull-wavefield solution. The two solutions have an obvious delay compared with the full-wavefieldsolution wavefront. The factorized Godunov solution is completely inconsistent with the full-wavefield solution because it is not numerically stable. Only our hybrid method produces a stableand accurate solution that is highly consistent with the full-wavefield solution shown in Figure 8d. The fourth example in Figure 9 is a five-block anisotropic heterogeneous model. The modelis 20 km in the X direction and 5 km in the Z direction, with a uniform grid sampling of 25 m inboth directions. The model has uniform Thomsen parameters of ε = 0 . and δ = − . , but hasa strongly contrasted V p and TTI symmetry axis tilt angle θ . The V p contrast at each interface is2000 m/s, and that of the tilt angle θ is at least π/ .Figures 10a-d show the traveltime filed solutions computed using the Godunov method, factor-ized Godunov method, LF-3 method and our hybrid method, respectively. The Godunov, factorizedGodunov and LF-3 solutions have obvious spurious modes starting from the X position of 10 km.This artifact is generated at the interface between the second and the third block and propagates tothe far end of the model. In the factorized Godunov solution in Figure 10b, we also observe someinstabilities near the interface between the first and the second block. These instabilities propagatesfrom the first interface to the positive direction of X, eventually deteriorating the traveltime fieldin the entire computational domain. By contrast, the solution computed with our hybrid method18hown in Figure 10d is the only one of the four solutions that is stable in all the five TTI blocks.We further compare the accuracy of different solutions against the full-wavefield solution inFigures 11-14. Figure 11 shows the computed traveltime field overlying on the full-wavefield so-lution at 1 s after source excitation. Figures 11a-d are the solutions computed using the Godunovmethod, the factorized Godunov method, the LF-3 method, and our hybrid method, respectively.The factorized Godunov solution contains some weak instabilities near the 5 km interface. TheGodunov and LF-3 solutions give a slight delay along the direction perpendicular to the TTI sym-metry axis compared with the full-wavefield solution. In comparison, our hybrid method producesa solution in Figure 11d that is both stable and accurate, and is highly consistent with the full-wavefield solution wavefront along all propagation directions.At snapshot time 2 s depicted in Figure 12, the inaccuracy caused by the instability of the fac-torized Godunov solution (Figure 12b) becomes fairly apparent, while the Godunov (Figure 12a)and LF-3 solutions (Figure 12c) start to show inconsistency with the full-wavefield solution at theZ position of approximately 3 km. This inconsistency is in fact the artifact in Figures 10a andc. Only our hybrid method produces a stable and accurate solution (Figure 12d) that is highlyconsistent with the full-wavefield solution.The consistency check between the full-wavefield solution and the eikonal equation solutionat two other time steps shown in Figures 13 and 14 further verifies that our hybrid method is ableto produce stable and accurate solutions to the anisotropic eikonal equation with strong mediumproperty contrasts where conventional methods fail. We verify the accuracy and stability of our hybrid method using a 3D anisotropic model inFigure 15 modified from the SEG/EAGE salt model. The model dimension is 6.76 km in both theX and Y directions and 2 km in the Z direction. The P-wave velocity model shown in Figure 15ahas a value range from 1500 m/s to 4500 m/s. We create the models of Thomsen parameters ε and δ , with values varying from 0 to 0.4 and -0.3 to 0.3, respectively, from the original velocitymodel. We create the TTI symmetry axis tilt angle θ and φ models with a value range from 0 to ◦ and 0 to ◦ , respectively. The φ model has the same spatial pattern as the θ model shown in19igure 15d.Figures 16a and b compares the full-wavefield solution at 0.2 s with the Godunov and ourhybrid method solutions, respectively. There exist obvious inconsistency between the Godunovsolution and the full-wavefield solution in the X-Z slice of Figure 16a at a depth of approximately0.3 km. The Godunov solution is faster than the wavefront of the full-wavefield solution aroundthis depth. At a depth of approximately 1.7 km, the Godunov solution is slower than the wavefrontof the full-wavefield solution. By contrast, our hybrid method solution in Figure 16b shows goodconsistency with the full-wavefield solution in both shallow and deep regions.The full-wavefield and Godunov solution consistency check in Figure 17a for the snapshot of0.3 s shows that the Godunov solution is faster than the wavefront of the full-wavefield solutionat the depth around 1.1 km in the X-Z slice. There exist obvious inconsistency between the twosolutions in the X-Y slice in Figure 17a. By contrast, our hybrid method produces a solution inFigure 17b that is consistent with the full-wavefield solution in all three slices.Figure 18 depicts a full-wavefield snapshot at 0.4 s superimposed with the corresponding trav-eltime contours obtained using the Godunov method and our hybrid method. The results furtherverify that that our hybrid method is stable and accurate for 3D heterogeneous anisotropic mediawith strong contrasts. We have developed a hybrid eikonal solver for computing first-arrival traveltime in 2D and3D anisotropic media. The numerical scheme of our hybrid eikonal solver consists of three steps:the Godunov fast locking-sweeping step, the base traveltime computation step in anellipticallyanisotropic media, and the third-order Lax-Friedrichs fast sweeping step. There are three advan-tages in our hybrid eikonal solver compared with existing eikonal solvers for anisotropic media.(1) The solver avoids the source singularity issue by multiplicative traveltime factorization and re-quires no specification of near-source traveltime values. (2) It is high-order accurate in space. (3)It can produces stable and accurate solution in models with strong anelliptically anisotropy, strongmedium property contrasts, and complex structures. We have used five numerical examples, in-cluding four 2D examples and one 3D example, to verify the stability and high-order accuracy20f our hybrid eikonal solver. The results show that our new method is advantageous in terms ofstability and accuracy compared with conventional approaches. Future work aims at extendingour method to address media with more complex anisotropies such as orthorhombic and rotatedorthorhombic anisotropies.
This work was supported by the U.S. Department of Energy (DOE) Geothermal TechnologiesOffice through the Los Alamos National Laboratory (LANL). LANL is operated by Triad NationalSecurity, LLC, for the U.S. DOE National Nuclear Security Administration (NNSA) under Con-tract No. 89233218CNA000001. This research used resources provided by the LANL InstitutionalComputing Program supported by the U.S. DOE NNSA under Contract No. 89233218CNA000001.
References
Bak, S., J. McLaughlin, and D. Renzi, 2010, Some improvements for the fast sweeping method:SIAM Journal on Scientific Computing, , no. 5, 2853–2874, doi: 10.1137/090749645.Buske, S., 1999, Three-dimensional pre-stack Kirchhoff migration of deep seismic reflec-tion data: Geophysical Journal International, , no. 1, 243–260, doi: 10.1046/j.1365-246x.1999.00789.x.Carcione, J. M., 2015, Wave fields in real media: Wave propagation in anisotropic, anelastic,porous and electromagnetic media (third edition): Elsevier, Amsterdam, Netherlands.Chambers, K., and J.-M. Kendall, 2008, A practical implementation of wave front construc-tion for 3-D isotropic media: Geophysical Journal International, , no. 3, 1030–1038, doi:10.1111/j.1365-246X.2008.03790.x.Eaton, D. W. S., 1993, Finite difference traveltime calculation for anisotropic media: GeophysicalJournal International, , no. 2, 273–280, doi: 10.1111/j.1365-246X.1993.tb03915.x.Ettrich, N., and D. Gajewski, 1998, Traveltime computation by perturbation with fd-eikonalsolvers in isotropic and weakly anisotropic media: Geophysics, , no. 3, 1066–1078, doi:10.1190/1.1444385. 21omel, S., S. Luo, and H. Zhao, 2009, Fast sweeping method for the factored eikonal equation:Journal of Computational Physics, , no. 17, 6440–6455, doi: 10.1016/j.jcp.2009.05.029.Fornberg, B., 1988, Generation of finite difference formulas on arbitrarily spaced grids: Mathe-matics of Computation, , no. 184, 699–706, doi: 10.2307/2008770.Gao, K., and L. Huang, 2018, Optimal damping profile ratios for stabilization of perfectly matchedlayers in general anisotropic media: Geophysics, , no. 1, T15–T30, doi: 10.1190/geo2017-0430.1.Gibson, R. L., V. Durussel, and K.-J. Lee, 2005, Modeling and velocity analysis with awavefront-construction algorithm for anisotropic media: Geophysics, , no. 4, T63–T74, doi:10.1190/1.1988188.G´omez, J. V., D. ´Alvarez, S. Garrido, and L. Moreno, 2019, Fast methods for eikonal equations:An experimental survey: IEEE Access, , 39005–39029, doi: 10.1109/access.2019.2906782.Gray, S. H., and W. P. May, 1994, Kirchhoff migration using eikonal equation traveltimes: Geo-physics, , no. 5, 810–817, doi: 10.1190/1.1443639.Grechka, V. Y., and G. A. McMechan, 1996, 3-D twopoint ray tracing for heterogeneous, weaklytransversely isotropic media: Geophysics, , no. 6, 1883–1894, doi: 10.1190/1.1444103.Jin, S., H. Liu, S. Osher, and Y.-H. R. Tsai, 2005, Computing multivalued physical observables forthe semiclassical limit of the Schr¨odinger equation: Journal of Computational Physics, , no.1, 222 – 241, doi: 10.1016/j.jcp.2004.11.008.Kao, C.-Y., S. Osher, and Y.-H. Tsai, 2005, Fast Sweeping Methods for Static Hamilton–Jacobi Equations: SIAM Journal on Numerical Analysis, , no. 6, 2612–2632, doi:10.1137/S0036142902419600.Kim, S., 2002, 3D eikonal solvers: Firstarrival traveltimes: Geophysics, , no. 4, 1225–1231, doi:10.1190/1.1500384.Kim, S., and R. Cook, 1999, 3-d traveltime computation using secondorder eno scheme: Geo-physics, , no. 6, 1867–1876, doi: 10.1190/1.1444693.Lambar´e, G., P. S. Lucio, and A. Hanyga, 1996, Two-dimensional multivalued traveltime andamplitude maps by uniform sampling of a ray field: Geophysical Journal International, , no.2, 584–598, doi: 10.1111/j.1365-246X.1996.tb00021.x.Le Bouteiller, P., M. Benjemaa, L. M´etivier, and J. Virieux, 2019, A discontinuous Galerkin fast-22weeping eikonal solver for fast and accurate traveltime computation in 3D tilted anisotropicmedia: Geophysics, , no. 2, C107–C118, doi: 10.1190/geo2018-0555.1.Lecomte, I., 1993, Finite difference calculation of first traveltimes in anisotropic media1: Geo-physical Journal International, , no. 2, 318–342, doi: 10.1111/j.1365-246X.1993.tb00890.x.Lin, F.-C., M. H. Ritzwoller, and R. Snieder, 2009, Eikonal tomography: surface wave tomogra-phy by phase front tracking across a regional broad-band seismic array: Geophysical JournalInternational, , no. 3, 1091–1110, doi: 10.1111/j.1365-246X.2009.04105.x.Lisitsa, V., and D. Vishnevskiy, 2010, Lebedev scheme for the numerical simulation of wavepropagation in 3D anisotropic elasticity: Geophysical Prospecting, , no. 4, 619–635, doi:10.1111/j.1365-2478.2009.00862.x.Luo, S., and J. Qian, 2011, Factored singularities and high-order Lax-Friedrichs sweeping schemesfor point-source traveltimes and amplitudes: Journal of Computational Physics, , no. 12,4742–4755, doi: 10.1016/j.jcp.2011.02.043.——–, 2012, Fast sweeping methods for factored anisotropic eikonal equations: Multiplicativeand additive factors: Journal of Scientific Computing, , no. 2, 360–382, doi: 10.1007/s10915-011-9550-y.Luo, S., J. Qian, and H. Zhao, 2012, Higher-order schemes for 3D first-arrival traveltimes andamplitudes: Geophysics, , no. 2, T47–T56, doi: 10.1190/geo2010-0363.1.Martinelli, S. L., 2012, An application of the level set method to underwater acous-tic propagation: Communications in Computational Physics, , no. 5, 1359–1391, doi:10.1017/S1815240600003479.Mel´endez, A., J. Korenaga, V. Sallar`es, A. Miniussi, and C. Ranero, 2015, TOMO3D: 3-D jointrefraction and reflection traveltime tomography parallel code for active-source seismic data—synthetic test: Geophysical Journal International, , no. 1, 158–174, doi: 10.1093/gji/ggv292.Pereyra, V., W. H. K. Lee, and H. B. Keller, 1980, Solving two-point seismic-ray tracing problemsin a heterogeneous medium: Part 1. A general adaptive finite difference method: Bulletin of theSeismological Society of America, , no. 1, 79–99.Podvin, P., and I. Lecomte, 1991, Finite difference computation of traveltimes in very contrastedvelocity models: a massively parallel approach and its associated tools: Geophysical JournalInternational, , no. 1, 271–284, doi: 10.1111/j.1365-246X.1991.tb03461.x.23ian, J., and S. Leung, 2006, A local level set method for paraxial geometrical optics: SIAMJournal on Scientific Computing, , no. 1, 206–223, doi: 10.1137/030601673.Qian, J., and W. W. Symes, 2002, Finitedifference quasiP traveltimes for anisotropic media: Geo-physics, , no. 1, 147–155, doi: 10.1190/1.1451438.Qian, J., Y.-T. Zhang, and H.-K. Zhao, 2007, Fast sweeping methods for eikonal equations ontriangular meshes: SIAM J. Numer. Anal., , no. 1, 83–107, doi: 10.1137/050627083.Qin, F., Y. Luo, K. B. Olsen, W. Cai, and G. T. Schuster, 1992, Finite-difference solutionof the eikonal equation along expanding wavefronts: Geophysics, , no. 3, 478–487, doi:10.1190/1.1443263.Rawlinson, N., and M. Sambridge, 2004, Wave front evolution in strongly heterogeneous layeredmedia using the fast marching method: Geophysical Journal International, , no. 3, 631–647,doi: 10.1111/j.1365-246X.2004.02153.x.Sadeghi, H., S. Suzuki, and H. Takenaka, 1999, A two-point, three-dimensional seismic ray tracingusing genetic algorithms: Physics of the Earth and Planetary Interiors, , no. 1, 355 – 365,doi: 10.1016/S0031-9201(99)00011-4.Sethian, J. A., and A. M. Popovici, 1999, 3-D traveltime computation using the fast marchingmethod: Geophysics, , no. 2, 516–523, doi: 10.1190/1.1444558.Soukina, S., D. Gajewski, and B. Kashtan, 2003, Traveltime computation for 3D anisotropic mediaby a finite-difference perturbation method: Geophysical Prospecting, , no. 5, 431–441, doi:10.1046/j.1365-2478.2003.00385.x.Taillandier, C., M. Noble, H. Chauris, and H. Calandra, 2009, First-arrival traveltime to-mography based on the adjoint-state method: Geophysics, , no. 6, WCB1–WCB10, doi:10.1190/1.3250266.Tsai, Y., L. Cheng, S. Osher, and H. Zhao, 2003, Fast Sweeping Algorithms for a Class ofHamilton-Jacobi Equations: SIAM Journal on Numerical Analysis, , no. 2, 673–694, doi:10.1137/S0036142901396533.Tsvankin, I., 2012, Seismic Signatures and Analysis of Reflection Data in Anisotropic Media, 3rded.: Society of Exploration Geophyscists.Vidale, J., 1988, Finite-difference calculation of travel times: Bulletin of the Seismological Societyof America, , no. 6, 2062–2076. 24inje, V., E. Iversen, and H. Gjøystdal, 1993, Traveltime and amplitude estimation using wavefrontconstruction: Geophysics, , no. 8, 1157–1166, doi: 10.1190/1.1443499.Waheed, U. B., and T. Alkhalifah, 2017, A fast sweeping algorithm for accurate solution of thetilted transversely isotropic eikonal equation using factorization: Geophysics, , no. 6, WB1–WB8, doi: 10.1190/geo2016-0712.1.Waheed, U. B., T. Alkhalifah, and H. Wang, 2015a, Efficient traveltime solutions of theacoustic TI eikonal equation: Journal of Computational Physics, , 62 – 76, doi:10.1016/j.jcp.2014.11.006.Waheed, U. B., C. E. Yarman, and G. Flagg, 2015b, An iterative, fast-sweeping-based eikonalsolver for 3D tilted anisotropic media: Geophysics, , no. 3, C49–C58, doi: 10.1190/geo2014-0375.1.Wang, Y., T. Nemeth, and R. T. Langan, 2006, An expanding-wavefront method for solvingthe eikonal equations in general anisotropic media: Geophysics, , no. 5, T129–T135, doi:10.1190/1.2235563.Zhang, Y.-T., H.-K. Zhao, and J. Qian, 2006, High order fast sweeping methods for static Hamilton-Jacobi equations: Journal of Scientific Computing, , no. 1, 25–56, doi: 10.1007/s10915-005-9014-3.Zhao, H., 2004, A fast sweeping method for eikonal equations: Mathematics of Computation, ,no. 250, 603–627, doi: 10.1090/S0025-5718-04-01678-3. Appendix A: Hybrid eikonal solver for 3D TTI media
In 3D, our hybrid eikonal solver is based on the following eikonal equation in TTI media: v x (cid:18) ∂t∂ ˆ x (cid:19) + v y (cid:18) ∂t∂ ˆ y (cid:19) + v z (cid:18) ∂t∂ ˆ z (cid:19) (cid:40) − ε − δ ) v z (cid:34)(cid:18) ∂t∂ ˆ x (cid:19) + (cid:18) ∂t∂ ˆ y (cid:19) (cid:35)(cid:41) = 1 , (56)where v x ( x ) = v y ( x ) = V p ( x ) (cid:112) ε ( x ) is the qP-wave velocity along the x - and y -axes, v z ( x ) = V p ( x ) is the qP-wave velocity along the z -axis (i.e., the anisotropy symmetry axis), ε = ε ( x ) and δ = δ ( x ) are Thomsen parameters describing a VTI medium’s anisotropy properties.Equation (56) is written in the rotated coordinates, and the spatial derivatives are combinations25f the spatial derivatives in the original coordinates: ∂t∂ ˆ x = α x ∂t∂x + β x ∂t∂y + γ x ∂t∂z , (57) ∂t∂ ˆ y = α y ∂t∂x + β y ∂t∂y + γ y ∂t∂z , (58) ∂t∂ ˆ z = α z ∂t∂x + β z ∂t∂y + γ z ∂t∂z , (59)with the coordinate transformation matrix R written as R = α x β x γ x α y β y γ y α z β z γ z = cos θ cos φ cos θ sin φ sin θ − sin φ cos φ − sin θ cos φ − sin θ sin φ cos θ , (60)where θ = θ ( x ) is the tilt angle of a VTI medium’s symmetry axis (i.e., the counterclockwiserotation angle of the symmetry axis w.r.t. the y -axis), φ = φ ( x ) is the azimuth angle of a VTImedium’s symmetry axis (i.e., the counterclockwise rotation angle of the symmetry axis w.r.t. the z -axis).We then solve the factorized eikonal equation using the Lax-Friedrichs scheme based on athird-order WENO discretization. To facilitate our description, we define F x ( τ, τ x , τ y , τ z ) = α x ( t τ x + t x τ ) + β x ( t τ y + t y τ ) + γ x ( t τ z + t z τ ) , (61) F y ( τ, τ x , τ y , τ z ) = α y ( t τ x + t x τ ) + β y ( t τ y + t y τ ) + γ y ( t τ z + t z τ ) , (62) F z ( τ, τ x , τ y , τ z ) = α z ( t τ x + t x τ ) + β z ( t τ y + t y τ ) + γ z ( t τ z + t z τ ) . (63)Then the Hamiltonian for equation (1) can be written as H ( τ, τ x , τ y , τ z ) = F x ( τ, τ x , τ y , τ z ) + F y ( τ, τ x , τ y , τ z ) + F z ( τ, τ x , τ y , τ z ) − ξF x ( τ, τ x , τ y , τ z ) F z ( t, τ x , τ y , τ z ) − ξF y ( τ, τ x , τ y , τ z ) F z ( t, τ x , τ y , τ z ) , (64)which leads to ∂H∂τ = 2 (cid:20) (1 − ξF z ) ∂F x ∂τ F x + (1 − ξF y ) ∂F y ∂τ F y + (1 − ξF x − ξF y ) ∂F z ∂τ F z (cid:21) = 2 (cid:2) (1 − ξF z )( α x t x + β x t y + γ x t z ) F x + (1 − ξF z )( α y t x + β y t y + γ y t z ) F y +(1 − ξF x − ξF y )( α z t x + β y t y + γ z t z ) F z (cid:3) , (65)26 H∂τ x = 2 (cid:20) (1 − ξF z ) ∂F x ∂τ x F x + (1 − ξF z ) ∂F y ∂τ x F y + (1 − ξF x − ξF y ) ∂F z ∂τ x F z (cid:21) = 2 t (cid:2) (1 − ξF z ) α x F x + (1 − ξF z ) α y F y + (1 − ξF x − ξF y ) α z F z (cid:3) , (66) ∂H∂τ y = 2 (cid:20) (1 − ξF z ) ∂F x ∂τ y F x + (1 − ξF z ) ∂F y ∂τ y F y + (1 − ξF x − ξF y ) ∂F z ∂τ y F z (cid:21) = 2 t (cid:2) (1 − ξF z ) β x F x + (1 − ξF z ) β y F y + (1 − ξF x − ξF y ) β z F z (cid:3) , (67) ∂H∂τ z = 2 (cid:20) (1 − ξF z ) ∂F x ∂τ z F x + (1 − ξF z ) ∂F y ∂τ z F y + (1 − ξF x − ξF y ) ∂F z ∂τ z F z (cid:21) = 2 t (cid:2) (1 − ξF z ) γ x F x + (1 − ξF z ) γ y F y + (1 − ξF x − ξF y ) γ z F z (cid:3) . (68)We then compute the following artificial viscosities along the three spatial axes: ω x = max Ω (cid:18)(cid:12)(cid:12)(cid:12)(cid:12) ∂H∂τ (cid:12)(cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12)(cid:12) ∂H∂τ x (cid:12)(cid:12)(cid:12)(cid:12)(cid:19) , (69) ω y = max Ω (cid:18)(cid:12)(cid:12)(cid:12)(cid:12) ∂H∂τ (cid:12)(cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12)(cid:12) ∂H∂τ y (cid:12)(cid:12)(cid:12)(cid:12)(cid:19) , (70) ω z = max Ω (cid:18)(cid:12)(cid:12)(cid:12)(cid:12) ∂H∂τ (cid:12)(cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12)(cid:12) ∂H∂τ z (cid:12)(cid:12)(cid:12)(cid:12)(cid:19) , (71)where Ω is the entire computational domain.This results in the following Lax-Friedrichs update scheme for the traveltime field τ : τ ( m +1) i,j,k = 1 − H (cid:16) τ ( m ) i,j,k , τ ( m ) , † x ; i,j,k , τ ( m ) , † y ; i,j,k , τ ( m ) , † z ; i,j,k (cid:17) + ω x τ ( m ) , ∗ x ; i,j,k + ω y τ ( m ) , ∗ y ; i,j,k + ω z τ ( m ) , ∗ z ; i,j,k ω x / ∆ x + ω y / ∆ y + ω z / ∆ z + τ ( m ) i,j,k , (72)where ∆ x , ∆ y and ∆ z are the regular grid sample intervals along the x -, y - and z -axis, respectively.The fast sweepings in 3D consist of the following sweepings: i = 1 , · · · , N x , j = 1 , · · · , N y , k = 1 , · · · , N z , (73) i = N x , · · · , , j = 1 , · · · , N y , k = 1 , · · · , N z , (74) i = 1 , · · · , N x , j = N y , · · · , , k = 1 , · · · , N z , (75) i = N x , · · · , , j = N y , · · · , , k = 1 , · · · , N z , (76) i = 1 , · · · , N x , j = 1 , · · · , N y , k = N z , · · · , , (77) i = N x , · · · , , j = 1 , · · · , N y , k = N z , · · · , , (78) i = 1 , · · · , N x , j = N y , · · · , , k = N z , · · · , , (79)27 = N x , · · · , , j = N y , · · · , , k = N z , · · · , . (80)where, again, the sweeping order is not important as along as all the listed sweepings are imple-mented. 28 ist of Figures V p = 2000 m/s, ε = 0 . , δ = 0 . , θ = 0 , (b) a HTI mediumwith V p = 2000 m/s, ε = 0 . , δ = − . , θ = π/ , (c) a TTI medium with V p = 2000 m/s, ε = 0 . , δ = 0 . , θ = π/ for single point sources, and (d) V p = 2000 m/s, ε = 0 . , δ = 0 . , θ = π/ for three random point sources. In allthe panels, the interval between any two adjacent contours is 0.03 s. The modelsare homogeneous in all the panels. . . . . . . . . . . . . . . . . . . . . . . . . . . 312 Base traveltime fields computed using the semi-analytic approach for a 3D TTImedium with V p = 2000 m/s, ε = 0 . , δ = − . , θ = π/ and φ = 5 π/ for apoint source. The interval between any two adjacent contours or isosurfaces is 0.1 s. 323 Traveltime fields in the homogeneous TTI model computed using (a) Godunovmethod, (b) factorized Godunov method, (c) LF-3 method, and (d) our hybridmethod. The interval between any two adjacent contours in all panels is 0.1 s. . . . 334 Difference between the traveltime field computed using our hybrid method (thesemi-analytic solution) and the traveltime field computed using (a) Godunov method,(b) factorized Godunov method, and (c) LF-3 method. . . . . . . . . . . . . . . . . 345 (a) A constant-gradient velocity model. (b) The analytical traveltime field com-puted using equation (54). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 356 Comparison between the convergence curves for the Godunov method (blue) andour hybrid method (red). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 367 Traveltime fields computed using (a) Godunov method, (b) factorized Godunovmethod, (c) LF-3 method and (d) our hybrid method. The interval between anytwo contours in all panels is 0.02 s, with contours starting from zero. Center blue-colored ball region is a TTI anomaly. . . . . . . . . . . . . . . . . . . . . . . . . 378 Comparison between the wavefield snapshot taken at 0.55 s after source excita-tion and the corresponding traveltime contours (red curves) computed using (a)Godunov method, (b) factorized Godunov method, (c) LF-3 method and (d) ourhybrid method. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 389 An TTI anisotropic model with block heterogeneities. The Thomsen parametersare ε = 0 . and δ = − . for all the blocks. . . . . . . . . . . . . . . . . . . . . . 3910 Traveltime field solutions computed using (a) Godunov method, (b) factorized Go-dunov method, (c) LF-3 method and (d) our hybrid method. The contours startfrom 0 s and the interval between any two adjacent contours is 0.15 s. . . . . . . . 4011 Comparison between the full-wavefield snapshot at 1 s and the correspondingeikonal equation solution contour (red curve) computed using (a) Godunov method,(b) factorized Godunov method, (c) LF-3 method and (d) our hybrid method. . . . 4112 Comparison between the full-wavefield snapshot at 3 s and the correspondingeikonal equation solution contour (red curve) computed using (a) Godunov method,(b) factorized Godunov method, (c) LF-3 method and (d) our hybrid method. . . . 4213 Comparison between the full-wavefield snapshot at 5 s and the correspondingeikonal equation solution contour (red curve) computed using (a) Godunov method,(b) factorized Godunov method, (c) LF-3 method and (d) our hybrid method. . . . 43294 Comparison between the full-wavefield snapshot at 7 s and the correspondingeikonal equation solution contour (red curve) computed using (a) Godunov method,(b) factorized Godunov method, (c) LF-3 method and (d) our hybrid method. . . . 4415 The anisotropic salt model. Panels (a)-(d) show the P-wave velocity, Thomsenparameters ε and δ , and TTI symmetry axis tilt angle θ , respectively. The TTIsymmetry axis tilt angle φ has a same spatial variation pattern with θ , with a valuerange of [0 , π ] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4516 Comparison between the full-wavefield snapshot at 0.2 s and the correspondingtraveltime contour (red curves and isosurface) computed using (a) Godunov methodand (b) our hybrid method. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4617 Comparison between the full-wavefield snapshot at 0.3 s and the correspondingtraveltime contour (red curves and isosurface) computed using (a) Godunov methodand (b) our hybrid method. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4718 Comparison between the full-wavefield snapshot at 0.4 s and the correspondingtraveltime contour (red curves and isosurface) computed using (a) Godunov methodand (b) our hybrid method. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4830 X (km) Z ( k m ) (a) X (km) Z ( k m ) (b) X (km) Z ( k m ) (c) X (km) Z ( k m ) (d) Figure 1: Base traveltime fields computed using the semi-analytic approach for (a) a VTI mediumwith V p = 2000 m/s, ε = 0 . , δ = 0 . , θ = 0 , (b) a HTI medium with V p = 2000 m/s, ε = 0 . , δ = − . , θ = π/ , (c) a TTI medium with V p = 2000 m/s, ε = 0 . , δ = 0 . , θ = π/ forsingle point sources, and (d) V p = 2000 m/s, ε = 0 . , δ = 0 . , θ = π/ for three random pointsources. In all the panels, the interval between any two adjacent contours is 0.03 s. The models arehomogeneous in all the panels. 31 X (km) Z ( k m ) Y ( k m ) Y (km)
Figure 2: Base traveltime fields computed using the semi-analytic approach for a 3D TTI mediumwith V p = 2000 m/s, ε = 0 . , δ = − . , θ = π/ and φ = 5 π/ for a point source. The intervalbetween any two adjacent contours or isosurfaces is 0.1 s.32 Horizontal Position (km) D ep t h ( k m ) (a) Horizontal Position (km) D ep t h ( k m ) (b) Horizontal Position (km) D ep t h ( k m ) (c) Horizontal Position (km) D ep t h ( k m ) (d) Figure 3: Traveltime fields in the homogeneous TTI model computed using (a) Godunov method,(b) factorized Godunov method, (c) LF-3 method, and (d) our hybrid method. The interval betweenany two adjacent contours in all panels is 0.1 s. 33
Horizontal Position (km) D ep t h ( k m ) -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08 Relative Traveltime Error (a)
Horizontal Position (km) D ep t h ( k m ) -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08 Relative Traveltime Error (b)
Horizontal Position (km) D ep t h ( k m ) -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08 Relative Traveltime Error (c)
Figure 4: Difference between the traveltime field computed using our hybrid method (the semi-analytic solution) and the traveltime field computed using (a) Godunov method, (b) factorizedGodunov method, and (c) LF-3 method. 34
X (km) Z ( k m ) V e l o c i t y ( m / s ) (a) X (km) Z ( k m ) T r a v e l t i m e ( s ) (b) Figure 5: (a) A constant-gradient velocity model. (b) The analytical traveltime field computedusing equation (54). 35igure 6: Comparison between the convergence curves for the Godunov method (blue) and ourhybrid method (red). 36 a) (b)(c) (d)
Figure 7: Traveltime fields computed using (a) Godunov method, (b) factorized Godunov method,(c) LF-3 method and (d) our hybrid method. The interval between any two contours in all panelsis 0.02 s, with contours starting from zero. Center blue-colored ball region is a TTI anomaly.37
X (km) Z ( k m ) (a) X (km) Z ( k m ) (b) X (km) Z ( k m ) (c) X (km) Z ( k m ) (d) Figure 8: Comparison between the wavefield snapshot taken at 0.55 s after source excitation andthe corresponding traveltime contours (red curves) computed using (a) Godunov method, (b) fac-torized Godunov method, (c) LF-3 method and (d) our hybrid method.38
X (km) Z ( k m ) V p = 2000 m/s ‚ = 45 – V p = 4000 m/s ‚ = 0 – V p = 2000 m/s ‚ = 60 – V p = 4000 m/s ‚ = 90 – V p = 2000 m/s ‚ = 45 – Figure 9: An TTI anisotropic model with block heterogeneities. The Thomsen parameters are ε = 0 . and δ = − . for all the blocks. 39 X (km) Z ( k m ) (a) X (km) Z ( k m ) (b) X (km) Z ( k m ) (c) X (km) Z ( k m ) (d) Figure 10: Traveltime field solutions computed using (a) Godunov method, (b) factorized Godunovmethod, (c) LF-3 method and (d) our hybrid method. The contours start from 0 s and the intervalbetween any two adjacent contours is 0.15 s. 40
X (km) Z ( k m ) (a) X (km) Z ( k m ) (b) X (km) Z ( k m ) (c) X (km) Z ( k m ) (d) Figure 11: Comparison between the full-wavefield snapshot at 1 s and the corresponding eikonalequation solution contour (red curve) computed using (a) Godunov method, (b) factorized Go-dunov method, (c) LF-3 method and (d) our hybrid method.41
X (km) Z ( k m ) (a) X (km) Z ( k m ) (b) X (km) Z ( k m ) (c) X (km) Z ( k m ) (d) Figure 12: Comparison between the full-wavefield snapshot at 3 s and the corresponding eikonalequation solution contour (red curve) computed using (a) Godunov method, (b) factorized Go-dunov method, (c) LF-3 method and (d) our hybrid method.42
X (km) Z ( k m ) (a) X (km) Z ( k m ) (b) X (km) Z ( k m ) (c) X (km) Z ( k m ) (d) Figure 13: Comparison between the full-wavefield snapshot at 5 s and the corresponding eikonalequation solution contour (red curve) computed using (a) Godunov method, (b) factorized Go-dunov method, (c) LF-3 method and (d) our hybrid method.43
X (km) Z ( k m ) (a) X (km) Z ( k m ) (b) X (km) Z ( k m ) (c) X (km) Z ( k m ) (d) Figure 14: Comparison between the full-wavefield snapshot at 7 s and the corresponding eikonalequation solution contour (red curve) computed using (a) Godunov method, (b) factorized Go-dunov method, (c) LF-3 method and (d) our hybrid method.44 Z ( k m ) Y ( k m ) P-wave Velocity (m/s) (a) Z ( k m ) Y ( k m ) Thomsen Parameter (cid:181) (b) Z ( k m ) Y ( k m ) -0.2 -0.1 0 0.1 0.2 0.3 Thomsen Parameter · (c) Z ( k m ) Y ( k m ) TTI Symmetry Axis Tilt Angle ‚ (rad) (d) Figure 15: The anisotropic salt model. Panels (a)-(d) show the P-wave velocity, Thomsen parame-ters ε and δ , and TTI symmetry axis tilt angle θ , respectively. The TTI symmetry axis tilt angle φ has a same spatial variation pattern with θ , with a value range of [0 , π ] .45 X (km) Z ( k m ) Y ( k m ) Y (km) (a)
X (km) Z ( k m ) Y ( k m ) Y (km) (b)(b)
X (km) Z ( k m ) Y ( k m ) Y (km) (b)(b)
Figure 16: Comparison between the full-wavefield snapshot at 0.2 s and the corresponding travel-time contour (red curves and isosurface) computed using (a) Godunov method and (b) our hybridmethod. 46
X (km) Z ( k m ) Y ( k m ) Y (km) (a)
X (km) Z ( k m ) Y ( k m ) Y (km) (b)(b)
X (km) Z ( k m ) Y ( k m ) Y (km) (b)(b)
Figure 17: Comparison between the full-wavefield snapshot at 0.3 s and the corresponding travel-time contour (red curves and isosurface) computed using (a) Godunov method and (b) our hybridmethod. 47
X (km) Z ( k m ) Y ( k m ) Y (km) (a)
X (km) Z ( k m ) Y ( k m ) Y (km) (b)(b)