# A hybrid eikonal solver for accurate first-arrival traveltime computation in
anisotropic media with strong contrasts

Kai Gao Corresponding Author; Email<EMAIL_ADDRESS>Los Alamos National
Laboratory, Geophysics Group, Los Alamos, NM 87545, USA  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 tomography, Kirchhoff migration, etc. There exist two major issues
in conventional eikonal solvers: 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 in medium
properties. We develop a stable and accurate hybrid eikonal solver for 2D and
3D transversely isotropic media with a tilted symmetry axis (TTI, or tilted
transversely isotropic media). Our new eikonal solver combines the traveltime
field factorization technique, the third-order Lax-Friedrichs update scheme,
and a new method for computing the base traveltime field. The solver has the
following three advantages. First, there is no need to assign exact traveltime
values in the near-source region, and the computed traveltime field near the
source location is accurate even for TTI media with strong anisotropy. Second,
the computed traveltime field is high-order accurate in space. Third, the
solver is numerically stable for 2D and 3D TTI media with strong anisotropy,
complex structures, and strong contrasts in medium properties. We verify the
stability and accuracy of our hybrid eikonal solver using several 2D and 3D
TTI medium examples. The results show that our solver is stable and accurate
in 2D and 3D complex TTI media, producing first-arrival traveltime fields that
are consistent with full-wavefield solutions.

Keywords: anisotropic media, first-arrival traveltime, eikonal equation,
strong medium property contrasts

##  1 Introduction

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 and May, (1994);
Buske, (1999) and first-arrival traveltime tomography Lin et al., (2009);
Taillandier et al., (2009) for reconstructing subsurface structures and medium
properties, rely on accurate and efficient 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 equations approximated from wave equations. A ray-
tracing system is then solved in the framework of either one-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éndez et al., (2015). These approaches are generally very efficient for
sparse source and receivers, but the computational costs may increase
dramatically as the number of source-receiver pairs increases. In addition,
ray tracing cannot trace rays to cover the entire model, and cannot handle
complex media with strong contrasts. A more flexible ray tracing approach is
the so-called the wavefront construction (WFC) method Vinje et al., (1993);
Lambaré et al., (1996); Gibson et al., (2005); Chambers and Kendall, (2008),
which computes traveltime and amplitude fields in a more intuitive manner by
dynamically inserting new rays where necessary. WFC generally requires more
programming efforts and is more computationally demanding compared with
conventional ray-based approaches.

Vidale, (1988) first developed an eikonal-equation-based approach to
generating first-arrival traveltimes 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. Currently, the most
widely used approaches include expanding wavefront methods Podvin and Lecomte,
(1991); Qin et al., (1992), fast marching methods Sethian and Popovici,
(1999); Rawlinson and 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. Eikonal solvers on
triangular or unstructured mesh (e.g., Qian et al.,, 2007; Le Bouteiller et
al.,, 2019) can handle complex interfaces and domain boundaries. Various high-
order and non-oscillatory numerical schemes Kim and Cook, (1999); Kim, (2002);
Luo and Qian, (2011); Luo et al., (2012) can improve the numerical accuracy
and stability of eikonal solvers in complex media with strong contrasts.
Different eikonal solvers have various computational complexities and
numerical accuracy. A comparison on several popular eikonal solvers can be
found in Gómez et al., (2019). One of the most distinct advantage of eikonal
solvers compared with ray-based methods is that the output of an eikonal
solver is a first-arrival traveltime field in the entire computational domain,
as opposed to ray tracing methods that compute only traveltimes on ray paths.
In addition, eikonal solvers generally allow models to be arbitrarily
heterogeneous and complex, whereas ray-based methods usually require simple or
smooth media.

Finite-difference eikonal solvers were first developed for isotropic media.
There is an intensive need for eikonal solvers in anisotropic media. Eikonal
equations for various kinds of anisotropic media are often significantly more
complex than that in isotropic media. Most of the eikonal solvers in isotropic
media require substantial modifications for anisotropic media if possible. In
anisotropic media, the phase and group velocity directions are generally not
the same Carcione, (2015). Therefore, eikonal solvers for anisotropic media
require sophisticated numerical schemes for updating traveltime fields. Eaton,
(1993) developed a high-order expanding-wavefront method on 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 the qP-wave traveltime field in TI media. Wang et al.,
(2006) developed an unconditionally-stable expanding wavefront method for the
eikonal equation in TI media. Their method explicitly tracks group 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. There exist several other methods for solving eikonal equations in
anisotropic media by assuming weak, elliptical anisotropy Ettrich and
Gajewski, (1998), by perturbing from elliptical reference medium Soukina et
al., (2003), or with only low-order numerical accuracy Lecomte, (1993).

The source singularity is a major problem in various eikonal solvers.
Conventional numerical schemes for solving eikonal equations are based on
local plane wave assumption for traveltime field update, and therefore cannot
accurately handle large curvatures of the traveltime field around the source
point. The numerical error in the near-source region can eventually
deteriorate the overall numerical accuracy in the entire computational domain.
The source singularity issue is generally solved with the traveltime
factorization method. The total traveltime field is factorized into an
addition or multiplication of a base traveltime and an additional or
multiplicative traveltime field, i.e., $t=t_{0}+\tau$ or $t=t_{0}\times\tau$,
where $t$ is the total traveltime field, $t_{0}$ is the base traveltime field,
and $\tau$ is the additional or multiplicative traveltime field Luo and Qian,
(2011, 2012); Luo et al., (2012). The base traveltime field is solved in a
homogeneous isotropic or elliptically anisotropic medium, while the traveltime
field $\tau$ is computed in the heterogeneous part of the medium. A notable
approach based 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 first decomposes the left-hand
side of the eikonal equation into a tilted elliptically anisotropic term and
an 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 TTI eikonal equation.

We develop a hybrid numerical scheme based on both the monotonic Godunov
scheme and the high-order weighted essentially non-oscillatory (WENO) scheme
to solve the eikonal equation in 2D and 3D anisotropic media. Our hybrid
eikonal solver has three advantages. First, our solver is free of source
singularity issue, and there is no need to assign traveltime around the source
by using multiplicative factorization of the traveltime field. Second, the
computed traveltime field of our 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, highly
complex structures, and strong medium property contrasts, by using the
weighted non-oscillatory scheme. Our eikonal solver employs both the
conventional first-order Godunov scheme and the third-order Lax-Friedrichs
scheme to achieve numerical stability and high-order accuracy. Therefore, we
call our solver a hybrid approach. To our knowledge, our hybrid eikonal solver
is the first fast-sweeping-based method to date that simultaneously holds
these three advantages.

Our paper is organized as follows. In the Methodology section, we describe the
three computational steps of our anisotropic eikonal solver, including the
first-order Godunov locking-sweeping step, the base traveltime field
computation step in the arbitrary TTI medium, and the third-order Lax-
Friedrichs fast sweeping step. We then use several numerical examples to
verify the stability and accuracy of our method. In the Conclusions section,
we summarize the most important features of our new anisotropic eikonal
solver.

##  2 Methodology

We derive the formulation for our hybrid eikonal solver in 2D anisotropic
media without loss of generality. We give the 3D formulation in Appendix A.

We adopt the following eikonal equation in TTI media Waheed et al., 2015b :

$v_{x}^{2}\left(\frac{\partial
t}{\partial\hat{x}}\right)^{2}+v_{z}^{2}\left(\frac{\partial
t}{\partial\hat{z}}\right)^{2}\left[1-2(\varepsilon-\delta)v_{z}^{2}\left(\frac{\partial
t}{\partial\hat{x}}\right)^{2}\right]=1,$ (1)

where $v_{x}(\mathbf{x})=V_{p}(\mathbf{x})\sqrt{1+2\varepsilon(\mathbf{x})}$
is the qP-wave velocity along the $x$-axis,
$v_{z}(\mathbf{x})=V_{p}(\mathbf{x})$ is the qP-wave velocity along the
$z$-axis (i.e., the anisotropy symmetry axis),
$\varepsilon=\varepsilon(\mathbf{x})$ and $\delta=\delta(\mathbf{x})$ are
Thomsen parameters describing anisotropy properties of a VTI medium. Equation
(1) is written in the rotated coordinates $\hat{x}-\hat{z}$, and the spatial
derivatives are combinations of the spatial derivatives in the unrotated
coordinates $x-z$:

$\displaystyle\frac{\partial t}{\partial\hat{x}}$
$\displaystyle=a_{x}\frac{\partial t}{\partial x}+c_{x}\frac{\partial
t}{\partial z},$ (2) $\displaystyle\frac{\partial t}{\partial\hat{z}}$
$\displaystyle=a_{z}\frac{\partial t}{\partial x}+c_{z}\frac{\partial
t}{\partial z},$ (3)

with the coordinate transformation matrix $\mathbf{R}$ written as

$\mathbf{R}=\begin{bmatrix}a_{x}&c_{x}\\\
a_{z}&c_{z}\end{bmatrix}=\begin{bmatrix}\cos\theta&\sin\theta\\\
-\sin\theta&\cos\theta\end{bmatrix},$ (4)

where $\theta=\theta(\mathbf{x})$ is the tilt angle of a VTI medium’s symmetry
axis (i.e., the counterclockwise rotation angle of the symmetry axis with
respect to the $y$-axis).

For notation clarity and derivation convenience, we further define

$\displaystyle\alpha_{x}$ $\displaystyle=v_{x}a_{x}=v_{x}\cos\theta,$ (5)
$\displaystyle\gamma_{x}$ $\displaystyle=v_{x}c_{x}=v_{x}\sin\theta,$ (6)
$\displaystyle\alpha_{z}$ $\displaystyle=v_{z}a_{z}=-v_{z}\sin\theta,$ (7)
$\displaystyle\gamma_{z}$ $\displaystyle=v_{z}c_{z}=v_{z}\cos\theta,$ (8)

and

$\xi=\frac{2(\varepsilon-\delta)}{1+2\varepsilon},$ (9)

leading to

$\left(\frac{\partial t}{\partial\hat{x}}\right)^{2}+\left(\frac{\partial
t}{\partial\hat{z}}\right)^{2}\left[1-\xi\left(\frac{\partial
t}{\partial\hat{x}}\right)^{2}\right]=1,$ (10)

with

$\displaystyle\frac{\partial t}{\partial\hat{x}}$
$\displaystyle=\alpha_{x}\frac{\partial t}{\partial
x}+\gamma_{x}\frac{\partial t}{\partial z},$ (11) $\displaystyle\frac{\partial
t}{\partial\hat{z}}$ $\displaystyle=\alpha_{z}\frac{\partial t}{\partial
x}+\gamma_{z}\frac{\partial t}{\partial z}.$ (12)

Generally, $v_{x}\neq v_{z}$ and $\theta\neq 0$ in TTI media.

###  2.1 Step I: First-order Godunov fast locking-sweeping

The first step of our hybrid eikonal solver is to solve equation (10) using an
iterative first-order fast sweeping method. The iterative approach is based on
rewriting equation (10) in the form of

$\left[1-\frac{1}{2}\xi\left(\frac{\partial
t}{\partial\hat{z}}\right)^{2}\right]\left(\frac{\partial
t}{\partial\hat{x}}\right)^{2}+\left[1-\frac{1}{2}\xi\left(\frac{\partial
t}{\partial\hat{x}}\right)^{2}\right]\left(\frac{\partial
t}{\partial\hat{z}}\right)^{2}=1.$ (13)

By setting two coefficients

$\displaystyle c_{x}=1-\frac{1}{2}\xi\left(\frac{\partial
t}{\partial\hat{z}}\right)^{2},$ (14) $\displaystyle
c_{z}=1-\frac{1}{2}\xi\left(\frac{\partial t}{\partial\hat{x}}\right)^{2},$
(15)

and absorbing them into $\partial t/\partial\hat{x}$ and $\partial
t/\partial\hat{z}$, we have

$\left(\frac{\partial t}{\partial\tilde{x}}\right)^{2}+\left(\frac{\partial
t}{\partial\tilde{z}}\right)^{2}=1,$ (16)

with

$\displaystyle\frac{\partial t}{\partial\tilde{x}}$
$\displaystyle=\sqrt{|c_{x}|}\alpha_{x}\frac{\partial t}{\partial
x}+\sqrt{|c_{x}|}\gamma_{x}\frac{\partial t}{\partial z},$ (17)
$\displaystyle\frac{\partial t}{\partial\tilde{z}}$
$\displaystyle=\sqrt{|c_{z}|}\alpha_{z}\frac{\partial t}{\partial
x}+\sqrt{|c_{z}|}\gamma_{z}\frac{\partial t}{\partial z}.$ (18)

The left-hand side of equation (16) is the eikonal correspondence in an
elliptically transversely isotropic 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}(\mathbf{x})=c_{z}(\mathbf{x})=1$, and solve equation (16) using fast
sweeping; Then we update $c_{x}(\mathbf{x})$ and $c_{z}(\mathbf{x})$ using the
computed traveltime $t$ based on equations (14) and (15), and solve equation
(16) again. The reason for such a reformatting is that a Godunov scheme for
equation (10) can be fairly difficult to derive, and possibly leads to high
computational complexity for local solvers. Instead, a first-order Godunov
scheme 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 sweeping iteration, we simply lock the points where
$t(\mathbf{x})^{(m)}=t(\mathbf{x})^{(m-1)}$ where the superscript $(m)$
represents the $m$-th round of fast sweeping, and update traveltime values
only at unlocked points. The criterion
$t(\mathbf{x})^{(m)}=t(\mathbf{x})^{(m-1)}$ might not be as accurate as the
one given in Bak et al., (2010) and Gómez et al., (2019) based on checking the
changes of neighbor points of a certain spatial point, but is much more
efficient to compute. In practice, we find that even for complex media, this
criterion results in accurate traveltime fields. The locking-sweeping
procedure can significantly reduce the computational costs for the first step.

###  2.2 Step II: Base traveltime field computation

In Steps II and III, we express the traveltime field $t=t(\mathbf{x})$ in
heterogeneous TTI media using a multiplicative factorization Luo and Qian,
(2012); Waheed and Alkhalifah, (2017) as

$t(\mathbf{x})=t_{0}(\mathbf{x})\tau(\mathbf{x}),$ (19)

where $t_{0}(\mathbf{x})$ is the base traveltime field in the homogeneous TTI
media, and $\tau(\mathbf{x})$ is the multiplicative traveltime field that
accounts for heterogeneities of the model.

This multiplicative factorization results in

$\displaystyle\frac{\partial t}{\partial\hat{x}}$
$\displaystyle=\alpha_{x}\left(t_{0}\frac{\partial\tau}{\partial
x}+\frac{\partial t_{0}}{\partial
x}\tau\right)+\gamma_{x}\left(t_{0}\frac{\partial\tau}{\partial
z}+\frac{\partial t_{0}}{\partial z}\tau\right),$ (20)
$\displaystyle\frac{\partial t}{\partial\hat{z}}$
$\displaystyle=\alpha_{z}\left(t_{0}\frac{\partial\tau}{\partial
x}+\frac{\partial t_{0}}{\partial
x}\tau\right)+\gamma_{z}\left(t_{0}\frac{\partial\tau}{\partial
z}+\frac{\partial t_{0}}{\partial z}\tau\right),$ (21)

which transform equation (10) into a factorized eikonal equation for TTI
media.

For the factorized eikonal equation, the base traveltime field $t_{0}$ and its
first-order spatial derivatives, say, $t_{0x}=\partial t_{0}/\partial x$ and
$t_{0z}=\partial t_{0}/\partial z$, are computed analytically, and are fixed
during fast-sweeping iterations Luo and Qian, (2012); Waheed and Alkhalifah,
(2017). The base traveltime field is essential to avoid the source singularity
issue in eikonal-equation-based traveltime computation without specifying
exact traveltime values at the points near the point source location.

Existing techniques for solving the factorized isotropic or anisotropic
eikonal equation assume that the background media is either isotropic where
$\varepsilon=\delta=0$ Fomel et al., (2009); Luo and Qian, (2011); Luo et al.,
(2012) or elliptically anisotropic where $\varepsilon=\delta\neq 0$ Luo and
Qian, (2011, 2012); Waheed and Alkhalifah, (2017), because there exist closed-
form expressions for computing the traveltime and its spatial derivatives in
the case of $\varepsilon=\delta$.

However, the requirement of $\varepsilon=\delta$ also limits the application
of multiplicative factorization in complex or strong anisotropic media. For
instance, to apply multiplicative factorization to the case where
$\varepsilon\neq\delta$, Waheed and Alkhalifah, (2017) had to use an iterative
scheme to update the right-hand side of the eikonal equation, and update the
analytic $t_{0}$, $t_{0x}$ and $t_{0z}$ after several iterations. Even with
such an iterative scheme, in each iteration, the background anisotropic medium
is assumed to be elliptically anisotropic. Therefore, the base traveltime
field $t_{0}$ never truly approximates that in the anisotropic media where
$\varepsilon\neq\delta$. For strongly anisotropic media, the group velocity
curve can significantly deviate from an ellipse. The resulting base traveltime
field might significantly differ from the true solution, eventually leading to
suboptimal traveltime solutions, even though it helps avoid the source
singularity issue.

In contrast to existing techniques where $t_{0}$ is computed for an
elliptically anisotropic medium, we develop a semi-analytic approach to
directly computing $t_{0}$ for anelliptically anisotropic medium where
$\varepsilon\neq\delta$. Therefore, the background anisotropic medium in our
method is the non-degenerated TTI medium at the reference point. This is the
most important difference between our hybrid method and existing factorized
eikonal solvers.

In the following derivations, to distinguish the rotation angle $\theta$ of
the TTI symmetry axis in the following description, we use $\vartheta$ to
represent the phase velocity angle in a VTI medium, which measures the
deviation angle from the vertical symmetry axis of a VTI medium. The phase
velocity of the qP-wave in a VTI medium can be written as Tsvankin, (2012)

$v_{\text{phase}}(\vartheta)=V_{p}\sqrt{1+\varepsilon\sin^{2}\vartheta-\frac{f}{2}\left[1-\sqrt{\left(1+\frac{2\varepsilon\sin^{2}\vartheta}{f}\right)^{2}-\frac{2(\varepsilon-\delta)\sin^{2}2\vartheta}{f}}\right]},$
(22)

where $f=1-V_{s}^{2}/V_{p}^{2}$, 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_{\text{phase}}(\vartheta)=V_{p}\sqrt{\frac{1}{2}+\varepsilon\sin^{2}\vartheta+\sqrt{\left(1+2\varepsilon\sin^{2}\vartheta\right)^{2}-2(\varepsilon-\delta)\sin^{2}2\vartheta}}.$
(23)

Note that there is no approximation in equation (23), and therefore it is
accurate even for strong anisotropy. We also have the group velocity angle
$\psi$ in terms of the phase velocity angle $\vartheta$ as Tsvankin, (2012)

$\psi(\vartheta)=\arctan\dfrac{\tan\vartheta+\dfrac{1}{v_{\text{phase}}(\vartheta)}\dfrac{dv_{\text{phase}}(\vartheta)}{d\vartheta}}{1-\dfrac{\tan\vartheta}{v_{\text{phase}}(\vartheta)}\dfrac{dv_{\text{phase}}(\vartheta)}{d\vartheta}},$
(24)

and the magnitude of the group velocity in terms of the phase velocity angle
$\vartheta$ as

$v_{\text{group}}(\vartheta)=v_{\text{phase}}(\vartheta)\sqrt{1+\frac{1}{v_{\text{phase}}(\vartheta)}\frac{dv_{\text{phase}}(\vartheta)}{d\vartheta}},$
(25)

where based on equation (23), we have

$\frac{dv_{\text{phase}}(\vartheta)}{d\vartheta}=\frac{V_{p}^{2}\sin
2\vartheta}{2v_{\text{phase}}(\vartheta)}\left[\frac{\varepsilon\left(1+2\varepsilon\sin^{2}\vartheta\right)-2(\varepsilon-\delta)\cos
2\vartheta}{\sqrt{\left(1+2\varepsilon\sin^{2}\vartheta\right)^{2}-2(\varepsilon-\delta)\sin^{2}2\vartheta}}+\varepsilon\right].$
(26)

For an arbitrary spatial point $\mathbf{x}=(x,z)$ in the computational domain,
the group velocity angle at this point can be evaluated as

$\psi=\arctan\frac{x-x_{0}}{z-z_{0}},$ (27)

where $\mathbf{x}_{0}=(x_{0},z_{0})$ is the position of the point source.

Our goal is to compute the semi-analytic magnitude of the group velocity at
this point, so that we can find the exact first-arrival traveltime at
$\mathbf{x}$. This requires the determination of the corresponding phase
velocity angle $\vartheta$ at $\mathbf{x}$, by which we can evaluate the
magnitude of the group velocity at $\mathbf{x}$ using equation (25).
Unfortunately, there is no closed-form expression to compute $\vartheta$ from
$\psi$ based on equation (24), because equation (24) is a complicated
transcendental equation and is extremely difficult to solve analytically, if
not impossible.

We therefore adopt a numerical method to compute the magnitude of group
velocity $v_{\text{group}}$ at $\mathbf{x}$. We first compute a series of
group velocity values
$v_{\text{group}}(\vartheta_{1}),v_{\text{group}}(\vartheta_{2}),\cdots,v_{\text{group}}(\vartheta_{n})$,
where $\vartheta_{1},\vartheta_{2},\cdots,\vartheta_{n}$ is an equal division
of the phase angle range $[0,\frac{\pi}{2}]$. Meanwhile, we compute the
corresponding group velocity angles
$\psi(\vartheta_{1}),\psi(\vartheta_{2}),\cdots,\psi(\vartheta_{n})$ based on
these phase velocity angles. In any VTI medium, the group velocity angles
$\psi(\vartheta_{1}),\psi(\vartheta_{2}),\cdots,\psi(\vartheta_{n})$ range
exactly from $[0,\frac{\pi}{2}]$, but are generally not equally distributed
within this range. Therefore, we obtain a series of group angle-velocity pairs
$\left\\{\psi(\vartheta_{i}),v_{\text{group}}(\vartheta_{i})\right\\}$ with
$i=1,2,\cdots,n$. We then use cubic spline interpolation to obtain the group
velocity for the spatial location $\mathbf{x}$, which corresponds to a group
velocity angle $\psi_{\mathbf{x}}$ based on equation (27). The interpolant
function of this cubic spline interpolation is built from the group angle-
velocity pairs
$\left\\{\psi(\vartheta_{i}),v_{\text{group}}(\vartheta_{i})\right\\}$. In
practical computations, we use a large $n$ to divide the range $[0,\pi/2]$,
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
$\vartheta$ for the spatial point $\mathbf{x}$ in our numerical scheme. The
group velocity value at any spatial point is obtained through an 1D
interpolation process with smooth and continuous interpolants built from the
analytic group angle-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 solving
process.

Because for any VTI medium, the phase or group velocity is symmetric with
respect to both axes, it is sufficient to build a complete group velocity
profile for the entire $2\pi$ range based on the computed group velocity
values in $[0,\pi/2]$. In addition, because any TTI medium is simply a
coordinate rotation result of some VTI medium, the group velocity values for
the TTI medium can be easily computed using the scheme described above. Assume
that the tilt angle of a TTI medium is $\theta$, then for a normalized spatial
location $\mathbf{x}=(x-x_{0},z-z_{0})$ where $(x_{0},z_{0})$ is the source
location, the corresponding directional vector in the unrotated coordinate is

$\mathbf{x}^{\prime}=\begin{bmatrix}x^{\prime}\\\
z^{\prime}\end{bmatrix}=\begin{bmatrix}\cos\theta&\sin\theta\\\
-\sin\theta&\cos\theta\end{bmatrix}\begin{bmatrix}x-x_{0}\\\
z-z_{0},\end{bmatrix}$ (28)

which indicates the group velocity angle corresponding to $\mathbf{x}$ should
be

$\psi^{\prime}=\arctan\frac{\left|x^{\prime}\right|}{\left|z^{\prime}\right|}.$
(29)

We take absolute value in equation (29) to ensure that the angle
$\psi^{\prime}$ falls in $[0,\pi/2]$.

We then compute the group velocity value $v_{\text{group}}|_{\psi^{\prime}}$
for $\mathbf{x}$ by the aforementioned interpolation procedure at
$\psi^{\prime}$, and the traveltime at $\mathbf{x}$ is

$t_{0}(\mathbf{x})=\frac{\sqrt{(x-x_{0})^{2}+(z-z_{0})^{2}}}{v_{\text{group}}|_{\psi^{\prime}}}.$
(30)

Finally, we use a high-order centered finite-difference scheme to compute the
spatial derivatives of the base traveltime field $t_{0x}$ and $t_{0z}$:

$\displaystyle t_{0x}(i,j)$ $\displaystyle=\frac{1}{\Delta
x}\sum_{l=1}^{M}c_{l}[t_{0}(i+l,j)-t_{0}(i-l,j)],$ (31) $\displaystyle
t_{0z}(i,j)$ $\displaystyle=\frac{1}{\Delta
z}\sum_{l=1}^{M}c_{l}[t_{0}(i,j+l)-t_{0}(i,j-l)],$ (32)

where $c_{l}$ are finite-difference coefficients, $M$ is the half length of
the finite-difference operator Fornberg, (1988), and $\Delta x$ and $\Delta z$
are the grid sizes in the $x$\- and $z$-directions, respectively. In our
computation, we use $M=10$, and compute the associate coefficients $c_{l}$
using the procedure described in Fornberg, (1988). Again, these fields are
practically of analytic accuracy because $t_{0}$ is practically of analytic
accuracy and $M$ is large.

In Figure 1, we show three examples for base traveltime computation in
anelliptically anisotropic TI media, including a VTI medium (Figure 1a), a HTI
medium (Figure 1b), and a TTI media (Figure 1c), all containing strong
anisotropies. All group velocity curves significantly deviate from an ellipse.
The results indicate that our semi-analytic approach can accurately compute
the base traveltime fields in complex anisotropic media. The numerical scheme
to compute the base traveltime field in arbitrary TTI media is also applicable
to 3D scenario as shown in Figure 2.

(a)  (b)

(c)  (d)

Figure 1: Base traveltime fields computed using the semi-analytic approach for
(a) a VTI medium with $V_{p}=2000$ m/s, $\varepsilon=0.25$, $\delta=0.1$,
$\theta=0$, (b) a HTI medium with $V_{p}=2000$ m/s, $\varepsilon=0.3$,
$\delta=-0.3$, $\theta=\pi/2$, (c) a TTI medium with $V_{p}=2000$ m/s,
$\varepsilon=0.4$, $\delta=0.1$, $\theta=\pi/6$ for single point sources, and
(d) $V_{p}=2000$ m/s, $\varepsilon=0.35$, $\delta=0.05$, $\theta=\pi/4$ for
three random point sources. In all the panels, the interval between any two
adjacent contours is 0.03 s. The models are homogeneous in all the panels.
Figure 2: Base traveltime fields computed using the semi-analytic approach for
a 3D TTI medium with $V_{p}=2000$ m/s, $\varepsilon=0.3$, $\delta=-0.3$,
$\theta=\pi/9$ and $\phi=5\pi/12$ for a point source. The interval between any
two adjacent contours or isosurfaces is 0.1 s.

When a model contains multiple simultaneous sources, we need to compute the
base traveltime fields for all the point sources, and compute a minimum base
traveltime field by finding the minimal value from all the base traveltime
field values at each point. That is,

$t_{0}(\mathbf{x})=\min\left\\{t_{0}^{s_{1}}(\mathbf{x}),t_{0}^{s_{2}}(\mathbf{x}),\cdots,t_{0}^{s_{n}}(\mathbf{x})\right\\},$
(33)

where $t_{0}^{s_{i}}$ represents the base traveltime field for the $i$-th
point source. The medium properties at different point source locations can be
different. Figure 1d shows a simple example of base traveltime 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-analytic solution to the eikonal equation regardless of TTI anisotropy
type and the number of simultaneous point sources, without the need of the
aforementioned Step I or the Step III described below.

###  2.3 Step III: Third-order Lax-Friedrichs fast sweeping

For heterogeneous media, once the initial traveltime filed $t$, the base
traveltime field $t_{0}$ and its derivatives are computed, we solve the
factorized eikonal equation (10) along with equations (20) and (21) in TTI
media using the Lax-Friedrichs scheme based on a third-order WENO
discretization Zhang et al., (2006). The initial multiplicative traveltime
field $\tau$ is computed as
$\tau(\mathbf{x})=t(\mathbf{x})/t_{0}(\mathbf{x})$, with
$\tau(x_{0},z_{0})\equiv 1$. The value of $\tau$ at location $(x_{0},z_{0})$
is kept unchanged during iterations. For multiple simultaneous point source
applications, the values of $\tau$ at all the point source locations are kept
unchanged during iterations. In contrast to existing high-order schemes where
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 use of the multiplicative factorization.

To facilitate our description, we define two functionals associated with the
multiplicative time field $\tau$ and its spatial derivatives $\tau_{x}$ and
$\tau_{z}$:

$\displaystyle F_{x}(\tau,\tau_{x},\tau_{z})$
$\displaystyle=\alpha_{x}\left(t_{0}\tau_{x}+t_{0x}\tau\right)+\gamma_{x}\left(t_{0}\tau_{z}+t_{0z}\tau\right),$
(34) $\displaystyle F_{z}(\tau,\tau_{x},\tau_{z})$
$\displaystyle=\alpha_{z}\left(t_{0}\tau_{x}+t_{0x}\tau\right)+\gamma_{z}\left(t_{0}\tau_{z}+t_{0z}\tau\right).$
(35)

Then the Hamiltonian for equation (10) can be written as

$H(\tau,\tau_{x},\tau_{z})=F_{x}^{2}(\tau,\tau_{x},\tau_{z})+F_{z}^{2}(\tau,\tau_{x},\tau_{z})-\xi
F_{x}^{2}(\tau,\tau_{x},\tau_{z})F_{z}^{2}(t,\tau_{x},\tau_{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:

$\displaystyle\omega_{x}$
$\displaystyle=\max\limits_{\Omega}\left(\left|\frac{\partial
H}{\partial\tau}\right|+\left|\frac{\partial
H}{\partial\tau_{x}}\right|\right),$ (37) $\displaystyle\omega_{z}$
$\displaystyle=\max\limits_{\Omega}\left(\left|\frac{\partial
H}{\partial\tau}\right|+\left|\frac{\partial
H}{\partial\tau_{z}}\right|\right),$ (38)

where $\Omega$ represents the entire computational domain. Based on equation
(36), we have

$\displaystyle\frac{\partial H}{\partial\tau}$ $\displaystyle=2\left[(1-\xi
F_{z}^{2})\frac{\partial F_{x}}{\partial\tau}F_{x}+(1-\xi
F_{x}^{2})\frac{\partial F_{z}}{\partial\tau}F_{z}\right]$
$\displaystyle=2\left[(1-\xi
F_{z}^{2})(\alpha_{x}t_{0x}+\gamma_{x}t_{0z})F_{x}+(1-\xi
F_{x}^{2})(\alpha_{z}t_{0x}+\gamma_{z}t_{0z})F_{z}\right],$ (39)
$\displaystyle\frac{\partial H}{\partial\tau_{x}}$
$\displaystyle=2\left[(1-\xi F_{z}^{2})\frac{\partial
F_{x}}{\partial\tau_{x}}F_{x}+(1-\xi F_{x}^{2})\frac{\partial
F_{z}}{\partial\tau_{x}}F_{z}\right]$ $\displaystyle=2t_{0}\left[(1-\xi
F_{z}^{2})\alpha_{x}F_{x}+(1-\xi F_{x}^{2})\alpha_{z}F_{z}\right],$ (40)
$\displaystyle\frac{\partial H}{\partial\tau_{z}}$
$\displaystyle=2\left[(1-\xi F_{z}^{2})\frac{\partial
F_{x}}{\partial\tau_{z}}F_{x}+(1-\xi F_{x}^{2})\frac{\partial
F_{z}}{\partial\tau_{z}}F_{z}\right]$ $\displaystyle=2t_{0}\left[(1-\xi
F_{z}^{2})\gamma_{x}F_{x}+(1-\xi F_{x}^{2})\gamma_{z}F_{z}\right].$ (41)

In our algorithm, we compute the quantities $\left|\partial H/\partial
t\right|$, $\left|\partial H/\partial\tau_{x}\right|$ and $\left|\partial
H/\partial\tau_{z}\right|$ at every spatial points of the model, and find
$\omega_{x}$ and $\omega_{z}$ based on the maximum values of these quantities.
Note that 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 property contrasts.

We then obtain the following Lax-Friedrichs scheme to update the
multiplicative traveltime field $\tau$ at the spatial grid point $(i,j)$:

$\displaystyle\tau_{i,j}^{(m+1)}=\dfrac{1-H\left(\tau_{i,j}^{(m)},\tau_{0x;i,j}^{(m),\dagger},\tau_{0z;i,j}^{(m),\dagger}\right)+\omega_{x}\tau_{0x;i,j}^{(m),*}+\omega_{z}\tau_{0z;i,j}^{(m),*}}{\omega_{x}/\Delta
x+\omega_{z}/\Delta z}+\tau_{i,j}^{(m)},$ (42)

where the superscripts $(m)$ and $(m+1)$ represent the values of $\tau_{i,j}$
at the $m^{\text{th}}$ and $(m+1)^{\text{th}}$ sweeping iterations,
respectively, and according to Zhang et al., (2006),

$\displaystyle\tau_{0x;i,j}^{\dagger}$
$\displaystyle=\frac{1}{2}\left(\tau_{0x;i,j}^{+}+\tau_{0x;i,j}^{-}\right),$
(43) $\displaystyle\tau_{0x;i,j}^{*}$
$\displaystyle=\frac{1}{2}\left(\tau_{0x;i,j}^{+}-\tau_{0x;i,j}^{-}\right),$
(44)

with the third-order WENO discretizations

$\displaystyle\tau_{0x;i,j}^{+}$
$\displaystyle=(1-w_{x}^{+})\frac{\tau_{i+1,j}-\tau_{i-1,j}}{2\Delta
x}+w_{x}^{+}\frac{-\tau_{i+2,j}+4\tau_{i+1,j}-3\tau_{i,j}}{2\Delta x},$ (45)
$\displaystyle\tau_{0x;i,j}^{-}$
$\displaystyle=(1-w_{x}^{-})\frac{\tau_{i+1,j}-\tau_{i-1,j}}{2\Delta
x}-w_{x}^{-}\frac{-\tau_{i-2,j}+4\tau_{i-1,j}-3\tau_{i,j}}{2\Delta x},$ (46)
$\displaystyle w_{x}^{+}$
$\displaystyle=\left\\{1+2\left[\frac{\epsilon+(\tau_{i+2,j}-2\tau_{i+1,j}+\tau_{i,j})^{2}}{\epsilon+(\tau_{i+1,j}-2\tau_{i,j}+\tau_{i-1,j})^{2}}\right]^{2}\right\\}^{-1},$
(47) $\displaystyle w_{x}^{-}$
$\displaystyle=\left\\{1+2\left[\frac{\epsilon+(\tau_{i-2,j}-2\tau_{i-1,j}+\tau_{i,j})^{2}}{\epsilon+(\tau_{i+1,j}-2\tau_{i,j}+\tau_{i-1,j})^{2}}\right]^{2}\right\\}^{-1},$
(48)

where $\epsilon$ is a small number to avoid singularity. The expressions for
$\tau_{0z;i,j}^{\dagger}$ and $\tau_{0z;i,j}^{*}$ can be analogously derived.

Because we actually use the first-order traveltime field computed at Step I as
the initial solution for Step III, Step III in our algorithm requires much
fewer iterations to achieve accurate results compared with the LF-3 method
that directly solves the eikonal equation from a rough or constant initial
guess.

###  2.4 Workflow of our hybrid eikonal solver

We summarize the workflow of our hybrid eikonal solver as follows:

  1. 1.

Compute an initial, first-order accurate traveltime field $t(\mathbf{x})$
using the first-order Godunov algorithm based on the locking-sweeping
procedure in the following order:

$\displaystyle i=1,\cdots,N_{x},\quad j=1,\cdots,N_{z},$ (49) $\displaystyle
i=N_{x},\cdots,1,\quad j=1,\cdots,N_{z},$ (50) $\displaystyle
i=1,\cdots,N_{x},\quad j=N_{z},\cdots,1,$ (51) $\displaystyle
i=N_{x},\cdots,1,\quad j=N_{z},\cdots,1,$ (52)

where $i$ and $j$ are indices of the finite-difference grids, and $N_{x}$ and
$N_{z}$ are the number of grids in the model in the $x$\- and $z$-directions,
respectively. The order of fast-sweeping is not important.

  2. 2.

Compute the base traveltime field $t_{0}(\mathbf{x})$ and its spatial
derivatives using the semi-analytic approach described in Step II.

  3. 3.

Compute the multiplicative traveltime field $\tau(\mathbf{x})$ using the
numerical scheme described in Step III. The sweeping follows the order listed
in equations (49)-(52).

##  3 Numerical Results

We use five numerical examples to verify the stability and accuracy of our
hybrid eikonal solver for 2D and 3D TTI anisotropic media. We compare the
results from four different methods:

  1. 1.

Godunov: the first-order Godunov method based on the iterative scheme
developed in Waheed et al., 2015b .

  2. 2.

Factorized Godunov: the first-order Godunov method based on the iterative
scheme for the factorized eikonal equation developed in Waheed and Alkhalifah,
(2017).

  3. 3.

LF-3: the third-order direct Lax-Friedrichs method without traveltime
factorization developed in Luo and Qian, (2012).

  4. 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 meaning of “iterative” in the first-two approaches is that we need
to update the right-hand side of the degenerated 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 the Lax-Friedrichs update scheme and the third-order
WENO finite-difference scheme, without any degeneration 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 to study the efficacy of these methods in realistic
computational tasks. In practical applications, the media around the source
can be heterogeneous, where assigning exact traveltime values can be very
difficult, if not impossible, particularly for heterogeneous, anelliptically
anisotropic media. In all the implementations, we only fix the traveltime
value at the source point (i.e., where $t=0$) over sweeping 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 to achieve convergence for LF-3. Therefore, we use the result
from the Godunov method as the initial guess for LF-3.

###  3.1 Homogeneous model

Since we are able to obtain traveltime field with analytic accuracy solely
using the algorithm in Step II, traveltime computation in homogeneous TTI
media is trivial for our hybrid method. However, 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 three conventional methods.

The model parameters for a homogeneous TTI medium are $V_{p}=2000$ m/s,
Thomsen parameters $\varepsilon=0.25$ and $\delta=0.05$, and TTI symmetry axis
tilt angle $\theta=\pi/4$. Figures 3a-d show the traveltime 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 the reference solution.

Visually the solutions from different methods are close to one another, except
the near-source traveltime field contours in Figure 3a (Godunov method), which
clearly deviate from the reference traveltime contours shown in Figure 3d (our
hybrid method). We compute the differences between the 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-3 method is the least accurate, partially because we do
not assign exact values in the near-source region for this method. The two
Godunov methods have higher accuracy compared with the LF-3 method.
Nevertheless, obvious errors occur at the near-source region in both
solutions, and the errors become larger with the increased distance away from
the source position.

(a)

(b)

(c)

(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 between any two adjacent contours in all panels is
0.1 s.

(a)

(b)

(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) factorized Godunov method, and (c) LF-3 method.

###  3.2 Gradient model

In the second example, we study the convergence of our hybrid method using an
isotropic constant 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

$\frac{1}{s(\mathbf{x})}=\frac{1}{s_{0}}+\mathbf{G}_{0}\cdot(\mathbf{x}-\mathbf{x}_{0}),$
(53)

where $s(\mathbf{x})$ is the spatially variant slowness, $s_{0}$ is the
slowness at the source point $\mathbf{x}_{0}$, and $\mathbf{G}_{0}$ is the
constant gradient of the model. The analytical first-arrival traveltime for
this medium is Fomel et al., (2009):

$t(\mathbf{x})=\frac{1}{|\mathbf{G}_{0}|}\mathrm{arccosh}\left(1+\frac{1}{2}s(\mathbf{x})s_{0}|\mathbf{G}_{0}|^{2}|\mathbf{x}-\mathbf{x}_{0}|^{2}\right),$
(54)

with

$\mathrm{arccosh}(x)=\ln\left(x+\sqrt{x^{2}-1}\right).$ (55)

We set a point source at $\mathbf{x}_{0}=(2.5,2.5)$ km in the model, with a
background constant slowness $s_{0}=1/3000$ s/m. The constant gradient is
$\mathbf{G}_{0}=(G_{x},G_{z})=(0.1,0.5)$ s-1. The velocity in the model varies
from 1500 m/s to 4500 m/s as shown in Figure 5a. Figure 5b depicts the
corresponding analytical traveltime field computed using equation (54).

We compute the traveltime fields with different model grid spacing using the
Godunov and our hybrid method in this isotropic model. We compare the relative
$L_{2}$-norm misfit between the numerical solutions and the analytical
solution in Figure 6, demonstrating that our method is more accurate than the
Godunov method. The convergence order of our method is approximately 3.13
while that of the Godunov method is approximately 0.79. Even at the largest
grid size where the number of grids along each direction is 10, our hybrid
method is almost two orders of magnitude more accurate than the Godunov
method.

(a)  (b)

Figure 5: (a) A constant-gradient velocity model. (b) The analytical
traveltime field computed using equation (54). Figure 6: Comparison between
the convergence curves for the Godunov method (blue) and our hybrid method
(red).

###  3.3 Ball model

We use the third numerical example to verify the numerical stability of our
method in an anisotropic medium with a strong contrast as depicted in Figure
7. The background medium is an isotropic homogeneous medium with $V_{p}=5,000$
m/s and $\varepsilon=\delta=\theta=0$. The blue region at the center of the
model indicates the location of a low-velocity strong TTI anisotropic anomaly
with $V_{p}=1,800$ m/s, $\varepsilon=0.3$, $\delta=-0.3$ and $\theta=\pi/4$.
The model size is 3.2 km in both dimensions. The grid size is 10 m in both
directions.

Figures 7a-d show the traveltime fields computed using Godunov, factorized
Godunov, LF-3 and our hybrid method, respectively. The Godunov solution
exhibits weak instabilities around the boundaries of the TTI ball. The
factorized Godunov solution in Figure 7b shows evident numerical
instabilities. The computed traveltime field inside the TTI ball indicates
that the factorized Godunov scheme becomes unstable for this TTI medium, and
the traveltime field outside of the TTI ball is therefore mostly wrong with
error propagating from the inside of the ball.

Figure 7c is the LF-3 solution. The solution is more stable than those
computed using the Godunov methods. Figure 7d displays our hybrid solution.
Similar with the LF-3 solution, our hybrid 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. We compute the full-wavefield solution using the fully
staggered-grid finite-difference method Lisitsa and Vishnevskiy, (2010) with a
high-order stencil, plus an optimal multi-axial perfectly matched layers Gao
and Huang, (2018). Figures 8a-d show the full-wavefield solution at 0.55 s
after the source excitation and the corresponding eikonal equation solutions
in black curves computed using the Godunov method, factorized Godunov method,
LF-3 method and our hybrid method, respectively. We find that the Godunov and
the LF-3 methods are stable, yet are not consistent with the full-wavefield
solution. The two solutions have an obvious delay compared with the full-
wavefield solution 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 stable and accurate solution that is
highly consistent with the full-wavefield solution shown in Figure 8d.

(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 panels is 0.02 s, with contours starting from
zero. Center blue-colored ball region is a TTI anomaly.

(a)  (b)

(c)  (d)

Figure 8: Comparison between the wavefield snapshot taken at 0.55 s after
source excitation and the corresponding traveltime contours (red curves)
computed using (a) Godunov method, (b) factorized Godunov method, (c) LF-3
method and (d) our hybrid method.

###  3.4 Block model

The fourth example in Figure 9 is a five-block anisotropic heterogeneous
model. The model is 20 km in the X direction and 5 km in the Z direction, with
a uniform grid sampling of 25 m in both directions. The model has uniform
Thomsen parameters of $\varepsilon=0.4$ and $\delta=-0.2$, but has a strongly
contrasted $V_{p}$ and TTI symmetry axis tilt angle $\theta$. The $V_{p}$
contrast at each interface is 2000 m/s, and that of the tilt angle $\theta$ is
at least $\pi/6$.

Figures 10a-d show the traveltime filed solutions computed using the Godunov
method, factorized Godunov method, LF-3 method and our hybrid method,
respectively. The Godunov, factorized Godunov 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 to the far end of the model. In the factorized Godunov solution in
Figure 10b, we also observe some instabilities near the interface between the
first and the second block. These instabilities propagates from the first
interface to the positive direction of X, eventually deteriorating the
traveltime field in the entire computational domain. By contrast, the solution
computed with our hybrid method shown 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 in Figures 11-14. Figure 11 shows the computed traveltime
field overlying on the full-wavefield solution at 1 s after source excitation.
Figures 11a-d are the solutions computed using the Godunov method, 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. The Godunov and LF-3 solutions give a slight delay
along the direction perpendicular to the TTI symmetry axis compared with the
full-wavefield solution. In comparison, our hybrid method produces a 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 factorized 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 the Z position of
approximately 3 km. This inconsistency is in fact the artifact in Figures 10a
and c. Only our hybrid method produces a stable and accurate solution (Figure
12d) that is highly consistent with the full-wavefield solution.

The consistency check between the full-wavefield solution and the eikonal
equation solution at two other time steps shown in Figures 13 and 14 further
verifies that our hybrid method is able to produce stable and accurate
solutions to the anisotropic eikonal equation with strong medium property
contrasts where conventional methods fail.

Figure 9: An TTI anisotropic model with block heterogeneities. The Thomsen
parameters are $\varepsilon=0.4$ and $\delta=-0.2$ for all the blocks.

(a)

(b)

(c)

(d)

Figure 10: Traveltime field solutions computed using (a) Godunov method, (b)
factorized Godunov method, (c) LF-3 method and (d) our hybrid method. The
contours start from 0 s and the interval between any two adjacent contours is
0.15 s.

(a)

(b)

(c)

(d)

Figure 11: Comparison between the full-wavefield snapshot at 1 s and the
corresponding eikonal equation solution contour (red curve) computed using (a)
Godunov method, (b) factorized Godunov method, (c) LF-3 method and (d) our
hybrid method.

(a)

(b)

(c)

(d)

Figure 12: Comparison between the full-wavefield snapshot at 3 s and the
corresponding eikonal equation solution contour (red curve) computed using (a)
Godunov method, (b) factorized Godunov method, (c) LF-3 method and (d) our
hybrid method.

(a)

(b)

(c)

(d)

Figure 13: Comparison between the full-wavefield snapshot at 5 s and the
corresponding eikonal equation solution contour (red curve) computed using (a)
Godunov method, (b) factorized Godunov method, (c) LF-3 method and (d) our
hybrid method.

(a)

(b)

(c)

(d)

Figure 14: Comparison between the full-wavefield snapshot at 7 s and the
corresponding eikonal equation solution contour (red curve) computed using (a)
Godunov method, (b) factorized Godunov method, (c) LF-3 method and (d) our
hybrid method.

###  3.5 Salt model

We verify the accuracy and stability of our hybrid method using a 3D
anisotropic model in Figure 15 modified from the SEG/EAGE salt model. The
model dimension is 6.76 km in both the X and Y directions and 2 km in the Z
direction. The P-wave velocity model shown in Figure 15a has a value range
from 1500 m/s to 4500 m/s. We create the models of Thomsen parameters
$\varepsilon$ and $\delta$, with values varying from 0 to 0.4 and -0.3 to 0.3,
respectively, from the original velocity model. We create the TTI symmetry
axis tilt angle $\theta$ and $\phi$ models with a value range from 0 to
$90^{\circ}$ and 0 to $180^{\circ}$, respectively. The $\phi$ model has the
same spatial pattern as the $\theta$ model shown in Figure 15d.

Figures 16a and b compares the full-wavefield solution at 0.2 s with the
Godunov and our hybrid method solutions, respectively. There exist obvious
inconsistency between the Godunov solution and the full-wavefield solution in
the X-Z slice of Figure 16a at a depth of approximately 0.3 km. The Godunov
solution is faster than the wavefront of the full-wavefield solution around
this depth. At a depth of approximately 1.7 km, the Godunov solution is slower
than the wavefront of the full-wavefield solution. By contrast, our hybrid
method solution in Figure 16b shows good consistency 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 of 0.3 s shows that the Godunov solution is faster than the
wavefront of the full-wavefield solution at the depth around 1.1 km in the X-Z
slice. There exist obvious inconsistency between the two solutions in the X-Y
slice in Figure 17a. By contrast, our hybrid method produces a solution in
Figure 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 traveltime contours obtained using the Godunov method and our
hybrid method. The results further verify that that our hybrid method is
stable and accurate for 3D heterogeneous anisotropic media with strong
contrasts.

(a)  (b)

(c)  (d)

Figure 15: The anisotropic salt model. Panels (a)-(d) show the P-wave
velocity, Thomsen parameters $\varepsilon$ and $\delta$, and TTI symmetry axis
tilt angle $\theta$, respectively. The TTI symmetry axis tilt angle $\phi$ has
a same spatial variation pattern with $\theta$, with a value range of
$[0,\pi]$.

(a)

(b)

Figure 16: Comparison between the full-wavefield snapshot at 0.2 s and the
corresponding traveltime contour (red curves and isosurface) computed using
(a) Godunov method and (b) our hybrid method.

(a)

(b)

Figure 17: Comparison between the full-wavefield snapshot at 0.3 s and the
corresponding traveltime contour (red curves and isosurface) computed using
(a) Godunov method and (b) our hybrid method.

(a)

(b)

Figure 18: Comparison between the full-wavefield snapshot at 0.4 s and the
corresponding traveltime contour (red curves and isosurface) computed using
(a) Godunov method and (b) our hybrid method.

##  4 Conclusions

We have developed a hybrid eikonal solver for computing first-arrival
traveltime in 2D and 3D 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 anelliptically anisotropic
media, and the third-order Lax-Friedrichs fast sweeping step. There are three
advantages 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 requires 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, strong medium property contrasts, and complex structures. We have
used five numerical examples, including four 2D examples and one 3D example,
to verify the stability and high-order accuracy of our hybrid eikonal solver.
The results show that our new method is advantageous in terms of stability and
accuracy compared with conventional approaches. Future work aims at extending
our method to address media with more complex anisotropies such as
orthorhombic and rotated orthorhombic anisotropies.

##  5 Acknowledgments

This work was supported by the U.S. Department of Energy (DOE) Geothermal
Technologies Office through the Los Alamos National Laboratory (LANL). LANL is
operated by Triad National Security, LLC, for the U.S. DOE National Nuclear
Security Administration (NNSA) under Contract No. 89233218CNA000001. This
research used resources provided by the LANL Institutional Computing Program
supported by the U.S. DOE NNSA under Contract No. 89233218CNA000001.

## References

  * Bak et al., (2010) Bak, S., J. McLaughlin, and D. Renzi, 2010, Some improvements for the fast sweeping method: SIAM Journal on Scientific Computing, 32, no. 5, 2853–2874, doi: 10.1137/090749645.
  * Buske, (1999) Buske, S., 1999, Three-dimensional pre-stack Kirchhoff migration of deep seismic reflection data: Geophysical Journal International, 137, no. 1, 243–260, doi: 10.1046/j.1365-246x.1999.00789.x.
  * Carcione, (2015) Carcione, J. M., 2015, Wave fields in real media: Wave propagation in anisotropic, anelastic, porous and electromagnetic media (third edition): Elsevier, Amsterdam, Netherlands.
  * Chambers and Kendall, (2008) Chambers, K., and J.-M. Kendall, 2008, A practical implementation of wave front construction for 3-D isotropic media: Geophysical Journal International, 173, no. 3, 1030–1038, doi: 10.1111/j.1365-246X.2008.03790.x.
  * Eaton, (1993) Eaton, D. W. S., 1993, Finite difference traveltime calculation for anisotropic media: Geophysical Journal International, 114, no. 2, 273–280, doi: 10.1111/j.1365-246X.1993.tb03915.x.
  * Ettrich and Gajewski, (1998) Ettrich, N., and D. Gajewski, 1998, Traveltime computation by perturbation with fd-eikonal solvers in isotropic and weakly anisotropic media: Geophysics, 63, no. 3, 1066–1078, doi: 10.1190/1.1444385.
  * Fomel et al., (2009) Fomel, S., S. Luo, and H. Zhao, 2009, Fast sweeping method for the factored eikonal equation: Journal of Computational Physics, 228, no. 17, 6440–6455, doi: 10.1016/j.jcp.2009.05.029.
  * Fornberg, (1988) Fornberg, B., 1988, Generation of finite difference formulas on arbitrarily spaced grids: Mathematics of Computation, 51, no. 184, 699–706, doi: 10.2307/2008770.
  * Gao and Huang, (2018) Gao, K., and L. Huang, 2018, Optimal damping profile ratios for stabilization of perfectly matched layers in general anisotropic media: Geophysics, 83, no. 1, T15–T30, doi: 10.1190/geo2017-0430.1.
  * Gibson et al., (2005) Gibson, R. L., V. Durussel, and K.-J. Lee, 2005, Modeling and velocity analysis with a wavefront-construction algorithm for anisotropic media: Geophysics, 70, no. 4, T63–T74, doi: 10.1190/1.1988188.
  * Gómez et al., (2019) Gómez, J. V., D. Álvarez, S. Garrido, and L. Moreno, 2019, Fast methods for eikonal equations: An experimental survey: IEEE Access, 7, 39005–39029, doi: 10.1109/access.2019.2906782.
  * Gray and May, (1994) Gray, S. H., and W. P. May, 1994, Kirchhoff migration using eikonal equation traveltimes: Geophysics, 59, no. 5, 810–817, doi: 10.1190/1.1443639.
  * Grechka and McMechan, (1996) Grechka, V. Y., and G. A. McMechan, 1996, 3-D two‐point ray tracing for heterogeneous, weakly transversely isotropic media: Geophysics, 61, no. 6, 1883–1894, doi: 10.1190/1.1444103.
  * Jin et al., (2005) Jin, S., H. Liu, S. Osher, and Y.-H. R. Tsai, 2005, Computing multivalued physical observables for the semiclassical limit of the Schrödinger equation: Journal of Computational Physics, 205, no. 1, 222 – 241, doi: 10.1016/j.jcp.2004.11.008.
  * Kao et al., (2005) Kao, C.-Y., S. Osher, and Y.-H. Tsai, 2005, Fast Sweeping Methods for Static Hamilton–Jacobi Equations: SIAM Journal on Numerical Analysis, 42, no. 6, 2612–2632, doi: 10.1137/S0036142902419600.
  * Kim, (2002) Kim, S., 2002, 3‐D eikonal solvers: First‐arrival traveltimes: Geophysics, 67, no. 4, 1225–1231, doi: 10.1190/1.1500384.
  * Kim and Cook, (1999) Kim, S., and R. Cook, 1999, 3-d traveltime computation using second‐order eno scheme: Geophysics, 64, no. 6, 1867–1876, doi: 10.1190/1.1444693.
  * Lambaré et al., (1996) Lambaré, G., P. S. Lucio, and A. Hanyga, 1996, Two-dimensional multivalued traveltime and amplitude maps by uniform sampling of a ray field: Geophysical Journal International, 125, no. 2, 584–598, doi: 10.1111/j.1365-246X.1996.tb00021.x.
  * Le Bouteiller et al., (2019) Le Bouteiller, P., M. Benjemaa, L. Métivier, and J. Virieux, 2019, A discontinuous Galerkin fast-sweeping eikonal solver for fast and accurate traveltime computation in 3D tilted anisotropic media: Geophysics, 84, no. 2, C107–C118, doi: 10.1190/geo2018-0555.1.
  * Lecomte, (1993) Lecomte, I., 1993, Finite difference calculation of first traveltimes in anisotropic media1: Geophysical Journal International, 113, no. 2, 318–342, doi: 10.1111/j.1365-246X.1993.tb00890.x.
  * Lin et al., (2009) Lin, F.-C., M. H. Ritzwoller, and R. Snieder, 2009, Eikonal tomography: surface wave tomography by phase front tracking across a regional broad-band seismic array: Geophysical Journal International, 177, no. 3, 1091–1110, doi: 10.1111/j.1365-246X.2009.04105.x.
  * Lisitsa and Vishnevskiy, (2010) Lisitsa, V., and D. Vishnevskiy, 2010, Lebedev scheme for the numerical simulation of wave propagation in 3D anisotropic elasticity‡: Geophysical Prospecting, 58, no. 4, 619–635, doi: 10.1111/j.1365-2478.2009.00862.x.
  * Luo and Qian, (2011) Luo, S., and J. Qian, 2011, Factored singularities and high-order Lax-Friedrichs sweeping schemes for point-source traveltimes and amplitudes: Journal of Computational Physics, 230, no. 12, 4742–4755, doi: 10.1016/j.jcp.2011.02.043.
  * Luo and Qian, (2012) ——–, 2012, Fast sweeping methods for factored anisotropic eikonal equations: Multiplicative and additive factors: Journal of Scientific Computing, 52, no. 2, 360–382, doi: 10.1007/s10915-011-9550-y.
  * Luo et al., (2012) Luo, S., J. Qian, and H. Zhao, 2012, Higher-order schemes for 3D first-arrival traveltimes and amplitudes: Geophysics, 77, no. 2, T47–T56, doi: 10.1190/geo2010-0363.1.
  * Martinelli, (2012) Martinelli, S. L., 2012, An application of the level set method to underwater acoustic propagation: Communications in Computational Physics, 12, no. 5, 1359–1391, doi: 10.1017/S1815240600003479.
  * Meléndez et al., (2015) Meléndez, A., J. Korenaga, V. Sallarès, A. Miniussi, and C. Ranero, 2015, TOMO3D: 3-D joint refraction and reflection traveltime tomography parallel code for active-source seismic data—synthetic test: Geophysical Journal International, 203, no. 1, 158–174, doi: 10.1093/gji/ggv292.
  * Pereyra et al., (1980) Pereyra, V., W. H. K. Lee, and H. B. Keller, 1980, Solving two-point seismic-ray tracing problems in a heterogeneous medium: Part 1. A general adaptive finite difference method: Bulletin of the Seismological Society of America, 70, no. 1, 79–99.
  * Podvin and Lecomte, (1991) Podvin, P., and I. Lecomte, 1991, Finite difference computation of traveltimes in very contrasted velocity models: a massively parallel approach and its associated tools: Geophysical Journal International, 105, no. 1, 271–284, doi: 10.1111/j.1365-246X.1991.tb03461.x.
  * Qian and Leung, (2006) Qian, J., and S. Leung, 2006, A local level set method for paraxial geometrical optics: SIAM Journal on Scientific Computing, 28, no. 1, 206–223, doi: 10.1137/030601673.
  * Qian and Symes, (2002) Qian, J., and W. W. Symes, 2002, Finite‐difference quasi‐P traveltimes for anisotropic media: Geophysics, 67, no. 1, 147–155, doi: 10.1190/1.1451438.
  * Qian et al., (2007) Qian, J., Y.-T. Zhang, and H.-K. Zhao, 2007, Fast sweeping methods for eikonal equations on triangular meshes: SIAM J. Numer. Anal., 45, no. 1, 83–107, doi: 10.1137/050627083.
  * Qin et al., (1992) Qin, F., Y. Luo, K. B. Olsen, W. Cai, and G. T. Schuster, 1992, Finite-difference solution of the eikonal equation along expanding wavefronts: Geophysics, 57, no. 3, 478–487, doi: 10.1190/1.1443263.
  * Rawlinson and Sambridge, (2004) Rawlinson, N., and M. Sambridge, 2004, Wave front evolution in strongly heterogeneous layered media using the fast marching method: Geophysical Journal International, 156, no. 3, 631–647, doi: 10.1111/j.1365-246X.2004.02153.x.
  * Sadeghi et al., (1999) Sadeghi, H., S. Suzuki, and H. Takenaka, 1999, A two-point, three-dimensional seismic ray tracing using genetic algorithms: Physics of the Earth and Planetary Interiors, 113, no. 1, 355 – 365, doi: 10.1016/S0031-9201(99)00011-4.
  * Sethian and Popovici, (1999) Sethian, J. A., and A. M. Popovici, 1999, 3-D traveltime computation using the fast marching method: Geophysics, 64, no. 2, 516–523, doi: 10.1190/1.1444558.
  * Soukina et al., (2003) Soukina, S., D. Gajewski, and B. Kashtan, 2003, Traveltime computation for 3D anisotropic media by a finite-difference perturbation method: Geophysical Prospecting, 51, no. 5, 431–441, doi: 10.1046/j.1365-2478.2003.00385.x.
  * Taillandier et al., (2009) Taillandier, C., M. Noble, H. Chauris, and H. Calandra, 2009, First-arrival traveltime tomography based on the adjoint-state method: Geophysics, 74, no. 6, WCB1–WCB10, doi: 10.1190/1.3250266.
  * Tsai et al., (2003) Tsai, Y., L. Cheng, S. Osher, and H. Zhao, 2003, Fast Sweeping Algorithms for a Class of Hamilton-Jacobi Equations: SIAM Journal on Numerical Analysis, 41, no. 2, 673–694, doi: 10.1137/S0036142901396533.
  * Tsvankin, (2012) Tsvankin, I., 2012, Seismic Signatures and Analysis of Reflection Data in Anisotropic Media, 3rd ed.: Society of Exploration Geophyscists.
  * Vidale, (1988) Vidale, J., 1988, Finite-difference calculation of travel times: Bulletin of the Seismological Society of America, 78, no. 6, 2062–2076.
  * Vinje et al., (1993) Vinje, V., E. Iversen, and H. Gjøystdal, 1993, Traveltime and amplitude estimation using wavefront construction: Geophysics, 58, no. 8, 1157–1166, doi: 10.1190/1.1443499.
  * Waheed and Alkhalifah, (2017) Waheed, U. B., and T. Alkhalifah, 2017, A fast sweeping algorithm for accurate solution of the tilted transversely isotropic eikonal equation using factorization: Geophysics, 82, no. 6, WB1–WB8, doi: 10.1190/geo2016-0712.1.
  * (44) Waheed, U. B., T. Alkhalifah, and H. Wang, 2015a, Efficient traveltime solutions of the acoustic TI eikonal equation: Journal of Computational Physics, 282, 62 – 76, doi: 10.1016/j.jcp.2014.11.006.
  * (45) Waheed, U. B., C. E. Yarman, and G. Flagg, 2015b, An iterative, fast-sweeping-based eikonal solver for 3D tilted anisotropic media: Geophysics, 80, no. 3, C49–C58, doi: 10.1190/geo2014-0375.1.
  * Wang et al., (2006) Wang, Y., T. Nemeth, and R. T. Langan, 2006, An expanding-wavefront method for solving the eikonal equations in general anisotropic media: Geophysics, 71, no. 5, T129–T135, doi: 10.1190/1.2235563.
  * Zhang et al., (2006) Zhang, Y.-T., H.-K. Zhao, and J. Qian, 2006, High order fast sweeping methods for static Hamilton-Jacobi equations: Journal of Scientific Computing, 29, no. 1, 25–56, doi: 10.1007/s10915-005-9014-3.
  * Zhao, (2004) Zhao, H., 2004, A fast sweeping method for eikonal equations: Mathematics of Computation, 74, 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}^{2}\left(\frac{\partial
t}{\partial\hat{x}}\right)^{2}+v_{y}^{2}\left(\frac{\partial
t}{\partial\hat{y}}\right)^{2}+v_{z}^{2}\left(\frac{\partial
t}{\partial\hat{z}}\right)^{2}\left\\{1-2(\varepsilon-\delta)v_{z}^{2}\left[\left(\frac{\partial
t}{\partial\hat{x}}\right)^{2}+\left(\frac{\partial
t}{\partial\hat{y}}\right)^{2}\right]\right\\}=1,$ (56)

where
$v_{x}(\mathbf{x})=v_{y}(\mathbf{x})=V_{p}(\mathbf{x})\sqrt{1+2\varepsilon(\mathbf{x})}$
is the qP-wave velocity along the $x$\- and $y$-axes,
$v_{z}(\mathbf{x})=V_{p}(\mathbf{x})$ is the qP-wave velocity along the
$z$-axis (i.e., the anisotropy symmetry axis),
$\varepsilon=\varepsilon(\mathbf{x})$ and $\delta=\delta(\mathbf{x})$ are
Thomsen parameters describing a VTI medium’s anisotropy properties.

Equation (56) is written in the rotated coordinates, and the spatial
derivatives are combinations of the spatial derivatives in the original
coordinates:

$\displaystyle\frac{\partial t}{\partial\hat{x}}$
$\displaystyle=\alpha_{x}\frac{\partial t}{\partial x}+\beta_{x}\frac{\partial
t}{\partial y}+\gamma_{x}\frac{\partial t}{\partial z},$ (57)
$\displaystyle\frac{\partial t}{\partial\hat{y}}$
$\displaystyle=\alpha_{y}\frac{\partial t}{\partial x}+\beta_{y}\frac{\partial
t}{\partial y}+\gamma_{y}\frac{\partial t}{\partial z},$ (58)
$\displaystyle\frac{\partial t}{\partial\hat{z}}$
$\displaystyle=\alpha_{z}\frac{\partial t}{\partial x}+\beta_{z}\frac{\partial
t}{\partial y}+\gamma_{z}\frac{\partial t}{\partial z},$ (59)

with the coordinate transformation matrix $\mathbf{R}$ written as

$\mathbf{R}=\begin{bmatrix}\alpha_{x}&\beta_{x}&\gamma_{x}\\\
\alpha_{y}&\beta_{y}&\gamma_{y}\\\
\alpha_{z}&\beta_{z}&\gamma_{z}\end{bmatrix}=\begin{bmatrix}\cos\theta\cos\phi&\cos\theta\sin\phi&\sin\theta\\\
-\sin\phi&\cos\phi&0\\\
-\sin\theta\cos\phi&-\sin\theta\sin\phi&\cos\theta\end{bmatrix},$ (60)

where $\theta=\theta(\mathbf{x})$ is the tilt angle of a VTI medium’s symmetry
axis (i.e., the counterclockwise rotation angle of the symmetry axis w.r.t.
the $y$-axis), $\phi=\phi(\mathbf{x})$ is the azimuth angle of a VTI medium’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 a third-order WENO discretization. To facilitate our description, we
define

$\displaystyle F_{x}(\tau,\tau_{x},\tau_{y},\tau_{z})$
$\displaystyle=\alpha_{x}\left(t_{0}\tau_{x}+t_{0x}\tau\right)+\beta_{x}\left(t_{0}\tau_{y}+t_{0y}\tau\right)+\gamma_{x}\left(t_{0}\tau_{z}+t_{0z}\tau\right),$
(61) $\displaystyle F_{y}(\tau,\tau_{x},\tau_{y},\tau_{z})$
$\displaystyle=\alpha_{y}\left(t_{0}\tau_{x}+t_{0x}\tau\right)+\beta_{y}\left(t_{0}\tau_{y}+t_{0y}\tau\right)+\gamma_{y}\left(t_{0}\tau_{z}+t_{0z}\tau\right),$
(62) $\displaystyle F_{z}(\tau,\tau_{x},\tau_{y},\tau_{z})$
$\displaystyle=\alpha_{z}\left(t_{0}\tau_{x}+t_{0x}\tau\right)+\beta_{z}\left(t_{0}\tau_{y}+t_{0y}\tau\right)+\gamma_{z}\left(t_{0}\tau_{z}+t_{0z}\tau\right).$
(63)

Then the Hamiltonian for equation (1) can be written as

$\displaystyle H(\tau,\tau_{x},\tau_{y},\tau_{z})$
$\displaystyle=F_{x}^{2}(\tau,\tau_{x},\tau_{y},\tau_{z})+F_{y}^{2}(\tau,\tau_{x},\tau_{y},\tau_{z})+F_{z}^{2}(\tau,\tau_{x},\tau_{y},\tau_{z})$
$\displaystyle-\xi
F_{x}^{2}(\tau,\tau_{x},\tau_{y},\tau_{z})F_{z}^{2}(t,\tau_{x},\tau_{y},\tau_{z})-\xi
F_{y}^{2}(\tau,\tau_{x},\tau_{y},\tau_{z})F_{z}^{2}(t,\tau_{x},\tau_{y},\tau_{z}),$
(64)

which leads to

$\displaystyle\frac{\partial H}{\partial\tau}$ $\displaystyle=2\left[(1-\xi
F_{z}^{2})\frac{\partial F_{x}}{\partial\tau}F_{x}+(1-\xi
F_{y}^{2})\frac{\partial F_{y}}{\partial\tau}F_{y}+(1-\xi F_{x}^{2}-\xi
F_{y}^{2})\frac{\partial F_{z}}{\partial\tau}F_{z}\right]$
$\displaystyle=2\left[(1-\xi
F_{z}^{2})(\alpha_{x}t_{0x}+\beta_{x}t_{0y}+\gamma_{x}t_{0z})F_{x}+(1-\xi
F_{z}^{2})(\alpha_{y}t_{0x}+\beta_{y}t_{0y}+\gamma_{y}t_{0z})F_{y}\right.$
$\displaystyle\quad\left.+(1-\xi F_{x}^{2}-\xi
F_{y}^{2})(\alpha_{z}t_{0x}+\beta_{y}t_{0y}+\gamma_{z}t_{0z})F_{z}\right],$
(65) $\displaystyle\frac{\partial H}{\partial\tau_{x}}$
$\displaystyle=2\left[(1-\xi F_{z}^{2})\frac{\partial
F_{x}}{\partial\tau_{x}}F_{x}+(1-\xi F_{z}^{2})\frac{\partial
F_{y}}{\partial\tau_{x}}F_{y}+(1-\xi F_{x}^{2}-\xi F_{y}^{2})\frac{\partial
F_{z}}{\partial\tau_{x}}F_{z}\right]$ $\displaystyle=2t_{0}\left[(1-\xi
F_{z}^{2})\alpha_{x}F_{x}+(1-\xi F_{z}^{2})\alpha_{y}F_{y}+(1-\xi
F_{x}^{2}-\xi F_{y}^{2})\alpha_{z}F_{z}\right],$ (66)
$\displaystyle\frac{\partial H}{\partial\tau_{y}}$
$\displaystyle=2\left[(1-\xi F_{z}^{2})\frac{\partial
F_{x}}{\partial\tau_{y}}F_{x}+(1-\xi F_{z}^{2})\frac{\partial
F_{y}}{\partial\tau_{y}}F_{y}+(1-\xi F_{x}^{2}-\xi F_{y}^{2})\frac{\partial
F_{z}}{\partial\tau_{y}}F_{z}\right]$ $\displaystyle=2t_{0}\left[(1-\xi
F_{z}^{2})\beta_{x}F_{x}+(1-\xi F_{z}^{2})\beta_{y}F_{y}+(1-\xi F_{x}^{2}-\xi
F_{y}^{2})\beta_{z}F_{z}\right],$ (67) $\displaystyle\frac{\partial
H}{\partial\tau_{z}}$ $\displaystyle=2\left[(1-\xi F_{z}^{2})\frac{\partial
F_{x}}{\partial\tau_{z}}F_{x}+(1-\xi F_{z}^{2})\frac{\partial
F_{y}}{\partial\tau_{z}}F_{y}+(1-\xi F_{x}^{2}-\xi F_{y}^{2})\frac{\partial
F_{z}}{\partial\tau_{z}}F_{z}\right]$ $\displaystyle=2t_{0}\left[(1-\xi
F_{z}^{2})\gamma_{x}F_{x}+(1-\xi F_{z}^{2})\gamma_{y}F_{y}+(1-\xi
F_{x}^{2}-\xi F_{y}^{2})\gamma_{z}F_{z}\right].$ (68)

We then compute the following artificial viscosities along the three spatial
axes:

$\displaystyle\omega_{x}$
$\displaystyle=\max\limits_{\Omega}\left(\left|\frac{\partial
H}{\partial\tau}\right|+\left|\frac{\partial
H}{\partial\tau_{x}}\right|\right),$ (69) $\displaystyle\omega_{y}$
$\displaystyle=\max\limits_{\Omega}\left(\left|\frac{\partial
H}{\partial\tau}\right|+\left|\frac{\partial
H}{\partial\tau_{y}}\right|\right),$ (70) $\displaystyle\omega_{z}$
$\displaystyle=\max\limits_{\Omega}\left(\left|\frac{\partial
H}{\partial\tau}\right|+\left|\frac{\partial
H}{\partial\tau_{z}}\right|\right),$ (71)

where $\Omega$ is the entire computational domain.

This results in the following Lax-Friedrichs update scheme for the traveltime
field $\tau$:

$\displaystyle\tau_{i,j,k}^{(m+1)}=\dfrac{1-H\left(\tau_{i,j,k}^{(m)},\tau_{0x;i,j,k}^{(m),\dagger},\tau_{0y;i,j,k}^{(m),\dagger},\tau_{0z;i,j,k}^{(m),\dagger}\right)+\omega_{x}\tau_{0x;i,j,k}^{(m),*}+\omega_{y}\tau_{0y;i,j,k}^{(m),*}+\omega_{z}\tau_{0z;i,j,k}^{(m),*}}{\omega_{x}/\Delta
x+\omega_{y}/\Delta y+\omega_{z}/\Delta z}+\tau_{i,j,k}^{(m)},$ (72)

where $\Delta x$, $\Delta y$ and $\Delta 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:

$\displaystyle i=1,\cdots,N_{x},\quad j=1,\cdots,N_{y},\quad
k=1,\cdots,N_{z},$ (73) $\displaystyle i=N_{x},\cdots,1,\quad
j=1,\cdots,N_{y},\quad k=1,\cdots,N_{z},$ (74) $\displaystyle
i=1,\cdots,N_{x},\quad j=N_{y},\cdots,1,\quad k=1,\cdots,N_{z},$ (75)
$\displaystyle i=N_{x},\cdots,1,\quad j=N_{y},\cdots,1,\quad
k=1,\cdots,N_{z},$ (76) $\displaystyle i=1,\cdots,N_{x},\quad
j=1,\cdots,N_{y},\quad k=N_{z},\cdots,1,$ (77) $\displaystyle
i=N_{x},\cdots,1,\quad j=1,\cdots,N_{y},\quad k=N_{z},\cdots,1,$ (78)
$\displaystyle i=1,\cdots,N_{x},\quad j=N_{y},\cdots,1,\quad
k=N_{z},\cdots,1,$ (79) $\displaystyle i=N_{x},\cdots,1,\quad
j=N_{y},\cdots,1,\quad k=N_{z},\cdots,1.$ (80)

where, again, the sweeping order is not important as along as all the listed
sweepings are implemented.