# An accurate discontinuous Galerkin method for solving point-source Eikonal equation in 2-D heterogeneous anisotropic media

Ludovic Métivier, Philippe Le Bouteiller, Mondher Benjemaa, Jean Virieux

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Ludovic Métivier, Philippe Le Bouteiller, Mondher Benjemaa, Jean Virieux. An accurate discontinuous Galerkin method for solving point-source Eikonal equation in 2-D heterogeneous anisotropic media. Geophysical Journal International, 2018, 212 (3), pp.1498-1522. 10.1093/gji/ggx463. hal-02009739

HAL Id: hal-02009739

https://hal.science/hal-02009739

Submitted on 14 Oct 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# An accurate discontinuous Galerkin method for solving point-source Eikonal equation in 2-D heterogeneous anisotropic media

P. Le Bouteiller,<fnref n="1" /> M. Benjemaa,<fnref n="2" /> L. Métivier<fnref n="1" />,<fnref n="3" /> and J. Virieux<fnref n="1" />

Accepted 2017 October 20. Received 2017 October 17; in original form 2017 July 13

## SUMMARY

Accurate numerical computation of wave traveltimes in heterogeneous media is of major interest for a large range of applications in seismics, such as phase identification, data windowing, traveltime tomography and seismic imaging. A high level of precision is needed for traveltimes and their derivatives in applications which require quantities such as amplitude or take-off angle. Even more challenging is the anisotropic case, where the general Eikonal equation is a quartic in the derivatives of traveltimes. Despite their efficiency on Cartesian meshes, finite-difference solvers are inappropriate when dealing with unstructured meshes and irregular topographies. Moreover, reaching high orders of accuracy generally requires wide stencils and high additional computational load. To go beyond these limitations, we propose a discontinuous-finite-element-based strategy which has the following advantages: (1) the Hamiltonian formalism is general enough for handling the full anisotropic Eikonal equations; (2) the scheme is suitable for any desired high-order formulation or mixing of orders (*p-adaptivity*); (3) the solver is explicit whatever Hamiltonian is used (no need to find the roots of the quartic); (4) the use of unstructured meshes provides the flexibility for handling complex boundary geometries such as topographies (*h-adaptivity*) and radiation boundary conditions for mimicking an infinite medium. The point-source factorization principles are extended to this discontinuous Galerkin formulation. Extensive tests in smooth analytical media demonstrate the high accuracy of the method. Simulations in strongly heterogeneous media illustrate the solver robustness to realistic Earth-sciences-oriented applications.

**Key words:** Non-linear differential equations; Numerical approximations and analysis; Numerical modelling; Numerical solutions; Seismic anisotropy.

## 1 INTRODUCTION

Eikonal equation arises in many domains such as geometrical optics (Benamou 2003), optimality problems with shortest/geodesic path calculation (Moser 1991; Kimmel & Sethian 1998), computer vision problems with shape-from-shading (Rouy & Tourin 1992; Kimmel & Sethian 2001), as well as in geophysics and more particularly in seismic imaging for traveltime-based tomographic methods (Vidale 1988; Podvin & Lecomte 1991; Le Meur 1994; Hole & Zelt 1995; Le Meur et al. 1997; Billette & Lambaré 1998; Leung & Qian 2006; Virieux & Lambaré 2007; Taillandier et al. 2009). It establishes a mathematical link between the travel time of an oscillatory phenomenon from a source to any point of a given medium, and the speed at which this phenomenon propagates inside the medium. This equation is obtained in the high-frequency regime, where wavelengths are by far smaller than the characteristic lengths that describe the variations of the speed inside the medium. This asymptotic approach based on a ray *ansatz* allows for the smoothly varying quantities traveltime and amplitude to substitute for the highly oscillatory wavefield in the computation. This property is highly desirable for applications in Earth sciences where one generally needs to model the propagation of an oscillatory field along hundreds to thousands of wavelengths. Indeed, for such applications, standard volumetric methods, based on the discretization of the wave equation, yield large-scale problems occurring high computational costs.

A first method for computing traveltimes from a source point in the asymptotic approximation is the ray-tracing approach, which uses the method of characteristics in order to derive from the Eikonal equation a set of linear ordinary differential equations (ODEs). They are efficiently solved by integration along rays for given initial conditions following a Lagrangian approach (Courant & Hilbert 1966). A review of ray methods can be found in Virieux & Lambaré (2007). However, ray tracing suffers from inherent difficulties, such as non-uniform

<footnote n="1">¹ISTerre, Univ. Grenoble Alpes, F-38000 Grenoble, France. E-mail: philippe.le-bouteiller@univ-grenoble-alpes.fr</footnote>

<footnote n="2">²Laboratory of Stability and Control of Systems and Nonlinear PDE, University of Sfax, 3029 Sfax, Tunisia</footnote>

<footnote n="3">³LJK, Univ. Grenoble Alpes, CNRS, F-38000 Grenoble, France</footnote>
sampling of the medium, the complexity related to triplications, the existence of shadow zones and the difficulty to capture the ray arriving at a receiver from a given source.

On the other hand, solving directly the Eikonal equation in an Eulerian framework might be a good strategy when many source-receiver couples are considered. However, the Eikonal equation is nonlinear. Consequently, numerical grid-based schemes for solving Eikonal are more complex to build when compared to ray-tracing tools, and their convergence properties are generally not easy to establish. Fortunately, the robustness of this approach has been intensively studied inside a Hamilton–Jacobi framework. The Eikonal equation belongs to the general class of Hamilton–Jacobi equations for which existence and uniqueness of solutions have been established by Crandall & Lions (1983) and Lions (1982) by defining the mathematical concept of viscosity solution. In the case of the Eikonal equation, the viscosity solution corresponds to the shortest path, also known as the first-arrival traveltime. The first-arrival traveltime field $T(x)$ is the viscosity solution of an Eikonal equation expressed in a generic way by

$$ \mathcal{H}(x, \nabla_x u(x)) = 0. \qquad (1) $$

In this paper we restrain our study to these first-arrival traveltimes, leaving away multi-arrival considerations for future work, and in the following, for the sake of simplicity, we shall refer to first-arrival traveltimes as traveltimes.

Grid-based Eikonal solvers rely on two main ingredients, which allow for classifying them into general families. The first ingredient is the local discretization of the equation, for the computation of a traveltime value at a point, given the values of the traveltime at neighbouring points. The second one is the solution of the discretized system. The solving method is mainly ruled by the order in which the grid points are considered and the number of times each point is updated (single-pass methods/multipass methods). Grid-based solvers became popular in seismology since the work of Vidale (1988), which relies on a finite-difference discretization of the Eikonal equation (first ingredient) and an expanding square framework (second ingredient) which consists in updating traveltime values from the source point towards the boundaries of the domain along the edges of an expanding square. We refer the reader to the review by Rawlinson et al. (2007) for more details about the numerous recent developments on solving the Eikonal equation.

In presence of anisotropy, the most promising strategies that have been proposed so far rely on finite-difference discretization of the corresponding Eikonal equation and the Fast Sweeping Method (FSM) which is a multipass algorithm relying on global orderings of the nodes. All the nodes are updated at each Gauss–Seidel iteration (sweep), following alternating orderings (Boué & Dupuis 1999; Tsai et al. 2003; Kao et al. 2004; Zhao 2005). This approach has been extended to anisotropic media considering elliptical anisotropy at first (Tsai et al. 2003; Qian et al. 2007). Elliptical anisotropy is handled quite naturally since it can be understood as a dilation in one direction of the space. Therefore the design of a local solver relies on finding the roots of quadratic equation. However, the more general case of anelliptical tilted transversely isotropic (TTI) media is more challenging since spatial derivatives of the traveltime to the power of four appear in the Eikonal equation. Several approaches have been proposed, either by solving the quartic equation and selecting the appropriate root (Han et al. 2017), which may yield a high computational load, by treating the anellipticity as a perturbation to the elliptical case and solving for the corresponding time expansion (Waheed et al. 2015a), or by implementing a fixed-point iteration technique which solves an elliptical equation at each iteration with an appropriate right-hand side accounting for anellipticity (Tavakoli F. et al. 2015; Waheed et al. 2015b; Waheed & Alkhalifah 2017).

Most of implementations of all the Eikonal solvers are performed following finite-difference formalism. Reaching high-order accuracy generally involves a wide stencil, which can be sometimes in contradiction with singularities of the traveltime solution. Despite their efficiency on Cartesian meshes, another pitfall of finite differences is their difficulty in dealing with unstructured meshes and correctly handle complex geometries, such as in the case of an irregular topography. Therefore, in this work, we consider finite-element methods, and especially discontinuous Galerkin (DG) methods for their ability to handle unstructured meshes, and the so-called hp-adaptivity which allows for various orders and element sizes among a unique mesh. We shall consider numerical schemes using a general Hamiltonian formalism for TTI media in this paper by introducing a pseudo-time-dependent Hamilton–Jacobi evolution expressed by

$$ \partial_t u(x, t) + \mathcal{H}(x, \nabla_x u(x, t)) = 0, \qquad (2) $$

to be solved by a Runge–Kutta (RK) integration scheme (RKDG). A pseudo-time $t$ is introduced as well as a pseudo-time-dependent solution $u(x, t)$ which converges to a steady-state solution $u_\infty(x)$ we are looking for. The link between stationary and time-dependent Hamilton–Jacobi equations was suggested by Osher (1993) by using the level-set idea.

Initially developed for solving hyperbolic conservation laws (Cockburn & Shu 1998), the RKDG method was then applied to solve the conservation law verified by the derivatives of the solution of eq. (2), in order to reconstruct the solution itself afterward (Hu & Shu 1999). Efforts were made by Cheng & Shu (2007) to recover directly the solution of eq. (2) in order to simplify the scheme. However, the method would suffer from an entropy violation issue in some cells, which had to be corrected by a specific ad hoc procedure. Cheng & Wang (2014) achieved a new step for directly solving eq. (2) with a DG method. An entropy fix is embedded inside the weak formulation itself, which greatly simplifies the implementation with a compact scheme. For these reasons we shall consider the method proposed by Cheng & Wang (2014) in this work. Special attention will be devoted to boundary conditions: one is related to the source singularity where the solution is zero, the other one is the external boundary from where no information could come from outside. The first one is important as the solution accuracy depends on how it is handled (Pica 1997; Qian & Symes 2002a; Zhang et al. 2005a; Fomel et al. 2009) while the second one could be expressed quite naturally with the approach of Cheng & Wang (2014).

The main objective of this study is to assess convergence and accuracy of the solution and highlight the flexibility of this method, with specific focus on the steady-state solution (first-arrival traveltime) and its spatial derivatives (related to slowness vector) for isotropic and
anisotropic media. To develop such a method, we have designed and brought together the following novelties: (1) the design of outgoing flux conditions at the edges of a bounded domain; (2) the extension of the factorization technique to the DG discretization, and the resulting high-order accuracy regarding the point-source problem; (3) the consideration of anisotropy inside the numerical scheme with no causality-related complication; (4) the high accuracy reached in anelliptical anisotropic settings; (5) the use of unstructured DG meshes to account for variable topographies in a traveltime computation problem.

The remainder of the paper is organized as follows: in Section 2 we focus on the Hamiltonian formulation of the Eikonal equation for both isotropic and anisotropic cases. We then derive factored expressions for these equations in Section 3. In Section 4, we present the 2-D RKDG discretization adapted from Cheng & Wang (2014) for heterogeneous media with suitable boundary conditions and initial conditions. Section 5 illustrates convergence and flexibility properties of the method with a variety of numerical case studies. Concluding remarks are given in Section 6 along with a discussion about future work.

## 2 DYNAMIC HAMILTON-JACOBI FORMULATION FOR SOLVING EIKONAL EQUATIONS

### 2.1 Isotropic case

In an isotropic medium with a speed $c(x)$, the Hamiltonian can be written as

$$ \mathcal{H}(x, \nabla_x u) = \|\nabla_x u\| - \frac{1}{c(x)}, \qquad (3) $$

and the corresponding time-dependent Hamilton–Jacobi equation writes

$$ \partial_t u + \|\nabla_x u\| - \frac{1}{c(x)} = 0. \qquad (4) $$

The stationary state of (4) verifies Eikonal equation $\mathcal{H} = 0$. Since we are only interested in the stationary state of (2), we are free to consider other Hamiltonians which yield the same final state, such as

$$ \mathcal{H}(x, \nabla_x u) = c(x) \|\nabla_x u\| - 1, \qquad (5) $$

and the corresponding time-dependent Hamilton–Jacobi equation

$$ \partial_t u + c(x) \|\nabla_x u\| - 1 = 0. \qquad (6) $$

Eq. (6) describes the propagation of a front with the local speed $v(x) = c(x)$, whereas in (4) the front propagates with a uniform speed $v = 1$. This might impact upon the computational efficiency, although the steady state would be the same in both cases. The best choice is discussed further in Section 4.

### 2.2 Anisotropic case

From Christoffel's dispersion relation in an elastic medium and considering only the coupled P-SV propagation mode, we may write the corresponding Eikonal equation for 2-D vertical transversely isotropic (VTI) medium as

$$ ap_x^4 + bp_z^4 + cp_x^2 p_z^2 + dp_x^2 + ep_z^2 - 1 = 0, \qquad (7) $$

where

$$ \left\{ \begin{array}{l} a = -(1 + 2\epsilon)V_P^2 V_S^2, \\ b = -V_P^2 V_S^2, \\ c = -(1 + 2\epsilon)V_P^4 - V_S^4 + (V_P^2 - V_S^2)[V_P^2(1 + 2\delta) - V_S^2], \\ d = V_S^2 + (1 + 2\epsilon)V_P^2, \\ e = V_P^2 + V_S^2, \end{array} \right. \qquad (8) $$

and with the slowness vector $\mathbf{p} = (p_x, p_z)$ defined by

$$ \begin{aligned} n_x^2 &= v^2 p_x, \\ n_z^2 &= v^2 p_z. \end{aligned} \qquad (9) $$

The quantities $V_P$ and $V_S$ denote the P- and SV-wave velocities along the rotation-symmetry axis, while $\epsilon$ and $\delta$ are known as the Thomsen's parameters (Thomsen 1986). Their expressions in terms of elastic parameters are given in Appendix A, as well as the derivation of eq. (7).

The identification of the components of $\mathbf{p}$ as spatial derivatives of the solution of Hamilton–Jacobi eq. (2) leads to the VTI Hamiltonian

$$ \mathcal{H}_{\text{VTI}} = a(u_{,x})^4 + b(u_{,z})^4 + c(u_{,x})^2(u_{,z})^2 + d(u_{,x})^2 + e(u_{,z})^2 - 1, \qquad (10) $$
where the derivatives of $u(x, z, t)$ with respect to $x$ and $z$ are respectively denoted by $u_{,x}$ and $u_{,z}$. Eq. (10) is equivalent to the one proposed in (Postma 1955; Payton 1983; Carcione et al. 1988; Červený 2001; Slawinski 2003). Setting $\epsilon = \delta = 0$, we retrieve the isotropic case for two roots corresponding to P waves and S waves. When $\epsilon = \delta \neq 0$, the coefficient $c$ in front of the cross-term $u_{,x}^2 u_{,z}^2$ cancels out, so that we obtain the so-called elliptical anisotropy. This particular case is equivalent to a simple dilation applied to the isotropic case along the axis orthogonal to the rotation-symmetry axis. Such a particular case is of little physical interest as discussed by Levin (1979). In a typical VTI medium, we have $\epsilon > \delta$, and corresponding anelliptic effects are important to account for.

As considered by Alkhalifah (2000), eq. (10) can be simplified for the acoustic case which does not perturb the numerical solution since we only consider first-arrival traveltimes always related to P-wave propagation (we assume that compressional velocity is higher than shear velocity). Hence, for the acoustic case, $V_S = 0$, and we have

$$ \left\{ \begin{array}{l} a = 0, \\ b = 0, \\ c = -2(\epsilon - \delta)V_p^4, \\ d = (1 + 2\epsilon)V_p^2, \\ e = V_p^2. \end{array} \right. \qquad (11) $$

The VTI Hamiltonian becomes

$$ \mathcal{H}_{\text{VTI}} = d(u_{,x})^2 + e(u_{,z})^2 + c(u_{,x})^2(u_{,z})^2 - 1. \qquad (12) $$

Introducing a TTI medium implies an additional local tilt angle $\theta(\mathbf{x})$, which is the angle between the local rotation-symmetry axis and the vertical axis. Hamiltonian (12) becomes

$$ \mathcal{H}_{\text{TTI}} = d(u_{,x} \cos \theta + u_{,z} \sin \theta)^2 + e(u_{,z} \cos \theta - u_{,x} \sin \theta)^2 + c(u_{,x} \cos \theta + u_{,z} \sin \theta)^2 (u_{,z} \cos \theta - u_{,x} \sin \theta)^2 - 1. \qquad (13) $$

Other types of anisotropy might be derived the same way in a Hamiltonian formalism, such as tilted orthorhombic (TOR) for 3-D anisotropic structures (Waheed et al. 2015b).

# 3 POINT-SOURCE FACTORIZATION

## 3.1 Derivation of the equations

The point-source condition is known to impair accuracy and convergence orders of numerical solutions of Eikonal equation. More precisely, high-order numerical schemes exhibit large errors and are only first-order convergent if no special treatment is applied, due to the singularity of the solution close to the source (see e.g. the analysis by Qian & Symes (2002a)). Several studies have been performed in order to mitigate this effect in finite-difference schemes. The celerity transform, initially described by Pica (1997), was then further developed in Zhang et al. (2005a) and promoted under the word factorization in Fomel et al. (2009); Luo & Qian (2011). We propose to extend these studies to DG discretizations. By the use of a reference solution $u_0(\mathbf{x})$, which does not depend on time, the factorization consists in computing a multiplicative numerical solution $\tau(\mathbf{x}, t)$ such that

$$ u(\mathbf{x}, t) = u_0(\mathbf{x})\tau(\mathbf{x}, t), \qquad (14) $$

or an additive numerical solution, as proposed in Luo & Qian (2012), such that

$$ u(\mathbf{x}, t) = u_0(\mathbf{x}) + \tau(\mathbf{x}, t). \qquad (15) $$

Plugging expression (14) or (15) into eq. (2) yields a new equation to be solved for the new solution field $\tau(\mathbf{x}, t)$. With the additive formulation (15), the new equation in $\tau(\mathbf{x}, t)$ is of the same Hamilton-Jacobi type as the original one (2) in $u(\mathbf{x}, t)$, while the multiplicative formulation yields additional complexities out of the frame of eq. (2). For this reason, we shall only consider the additive formulation (15) in what follows. Rewriting (2) in terms of $u_0$ and $\tau$ for the isotropic Hamiltonian of eq. (3) yields

$$ \partial_t \tau + \|\nabla u_0 + \nabla_x \tau\| - \frac{1}{c} = 0. \qquad (16) $$

Applying the same strategy to the TTI Hamiltonian of eq. (13) gives

$$ \partial_t \tau + d [(u_{0,x} + \tau_{,x}) \cos \theta + (u_{0,z} + \tau_{,z}) \sin \theta]^2 + e [(u_{0,z} + \tau_{,z}) \cos \theta - (u_{0,x} + \tau_{,x}) \sin \theta]^2 \\ + c [(u_{0,x} + \tau_{,x}) \cos \theta + (u_{0,z} + \tau_{,z}) \sin \theta]^2 [(u_{0,z} + \tau_{,z}) \cos \theta - (u_{0,x} + \tau_{,x}) \sin \theta]^2 - 1 = 0. \qquad (17) $$

Using eqs (16) and (17), $\tau(\mathbf{x}, t)$ can be computed for a given reference solution $u_0(\mathbf{x})$.
## 3.2 Choice of the reference solution $u_0(x)$

### 3.2.1 Isotropic case

In an arbitrarily heterogeneous isotropic medium, the simplest assumption one can make is that of a homogeneous speed in a small vicinity of the source. This seems to be a reasonable approximation in the high-frequency regime where the velocity model should be smooth. Therefore the reference solution $u_0(x)$ to be used in eq. (16) is chosen to be the exact solution in a homogeneous medium taking the speed value at the source. Such solution is the distance function to the source divided by the source speed:

$$u_{0\text{ISO}}(\mathbf{x}) = \frac{\text{dist}(\mathbf{x}, \mathbf{x}_s)}{c(\mathbf{x}_s)} \qquad (18)$$

This reference solution captures the source singularity and allows the gain of one order of convergence and high level of accuracy. This property is illustrated in our numerical case studies.

### 3.2.2 TI case

For the anisotropic case, Luo & Qian (2012) have considered a factored anisotropic Eikonal equation, but only for elliptical media ($\epsilon = \delta$). Following their approach, the reference solution $u_0(x)$ in eq. (17) might be taken as the analytical solution in a homogeneous elliptical TI medium of velocity $V_p(x_s)$, thus accounting for $V_p(x_s)$, $\epsilon(x_s)$ and $\theta(x_s)$, but not $\delta(x_s)$. This yields

$$u_{0\text{VTI}}(\mathbf{x}) = \sqrt{\frac{(x - x_s)^2}{V_p(\mathbf{x}_s)\sqrt{1 + 2\epsilon}} + \frac{(z - z_s)^2}{V_p(\mathbf{x}_s)}} \qquad (19)$$

for the VTI case ($\theta(\mathbf{x}_s) = 0$). The TTI case is retrieved by the local rotation of coordinates

$$\begin{aligned} x - x_s &\rightarrow (x - x_s) \cos \theta(\mathbf{x}_s) + (z - z_s) \sin \theta(\mathbf{x}_s), \\ z - z_s &\rightarrow (z - z_s) \cos \theta(\mathbf{x}_s) - (x - x_s) \sin \theta(\mathbf{x}_s). \end{aligned} \qquad (20)$$

The reference solution would thus capture a part of the source singularity, although the remaining source effect due to the anellipticity would not be well accounted for near the source. The efficiency of such a factorization for the TTI case is illustrated in our numerical results.

In order to increase the accuracy in an anelliptical TI medium, we can also use exact solution of the anelliptical equation as the reference solution $u_0(x)$. However such a solution is not known analytically but only parametrically (Carcione et al. 1988). As detailed in Case study 3 in Section 5, the pre-computation of the anelliptical reference solution at all discrete points where values are needed allows for a significant improvement of the results.

# 4 RUNGE-KUTTA DISCONTINUOUS GALERKIN METHOD FOR HAMILTON-JACOBI EQUATIONS

## 4.1 Space discretization: the discontinuous Galerkin solver

We present the DG solver designed for time-dependent Hamilton–Jacobi eq. (2). The RKDG method relies on the method of lines (Schliesser 1991), where we first discretize the problem in space, yielding a system of ODEs in time that we solve with an RK solver. The first RKDG solver able to directly compute the solution of this equation was proposed by Cheng & Shu (2007). However, it exhibits several numerical complications such as the need for an $L^2$ reconstruction of the derivatives of the solution at cell interfaces and a posterior entropy fix. These complications are avoided in a more recent scheme proposed by Cheng & Wang (2014) so that implementation is made simpler. This new scheme is shown to work either on structured or unstructured meshes, convex and non-convex Hamiltonians. This is the one we implement.

In the following, the formulation of this scheme for a 2-D unstructured triangular mesh is described, and the meaning of the contributions of each term of the scheme is detailed. In the general case, the two-dimensional spatial domain $\Omega$ is partitioned into $n$ triangular elements denoted by $K_i$, $i = 1, \dots, n$. A local approximation space $P_i$ of dimension $d_i$ is chosen for each element $K_i$ together with a basis of shape functions $\phi_i^j(x, z)$, $j \in \{1, \dots, d_i\}$ spanning this space. The choice of the approximation space is local: it is attached to one element. This property is referred to as p-adaptivity in the literature, an interesting feature of the DG strategy allowing to adjust the numerical accuracy locally. In our numerical tests, we use standard polynomial approximation spaces $P^k$ containing all polynomials of degree at most $k$ with $k \in \{1, 2, 3\}$ (Cockburn & Shu 1998).

Following Cheng & Wang (2014), we define $n_{K_i}$ to be the outward unit normal to the $K_i$ cell boundary and $t_{K_i}$ the unit tangential vector. At cell interfaces, traces $v_h^\pm$, jumps $[v_h]$ and means $\overline{v_h}$, of any numerical quantity $v_h$ defined inside two neighbouring cells are given respectively by

$$v_h^\pm(\mathbf{x}) = \lim_{\epsilon \downarrow 0} v_h(\mathbf{x} \pm \epsilon n_{K_i}),$$
$$[v_h](\mathbf{x}) = v_h^+(\mathbf{x}) - v_h^-(\mathbf{x}), \\ \bar{v}_h(\mathbf{x}) = \frac{1}{2}(v_h^+(\mathbf{x}) + \bar{v}_h^-(\mathbf{x})). \tag{21}$$

At cell interface, a two-component vector is defined by the expression

$$\nabla_x u_{hK_i}^{\pm} = \left( \frac{(\nabla_x u_h \cdot \mathbf{n}_{K_i})^{\pm}}{\overline{\nabla_x u_h} \cdot \mathbf{t}_{K_i}} \right). \tag{22}$$

The first component $(\nabla_x u_h \cdot \mathbf{n}_{K_i})^{\pm}$ is the projection onto the normal $\mathbf{n}_{K_i}$, of the gradient of the numerical solution computed inside the $K_i$ cell (-), or inside its neighbour (+). The second component $\overline{\nabla_x u_h} \cdot \mathbf{t}_{K_i}$ is the mean of the projections onto the tangential vector $\mathbf{t}_{K_i}$ of the gradient of the numerical solution computed inside the $K_i$ cell and inside its corresponding neighbour.

The weak formulation of eq. (2) can be stated as follows:

Find $u_h(., t) \in \{v: v|_{K_i} \in \mathcal{P}_i, \forall i \in \{1, ..., n\}\} \forall t \ge 0$ such that

$$\int_{K_i} (\partial_t u_h(\mathbf{x}, t) + \mathcal{H}(\mathbf{x}, \nabla_x u_h(\mathbf{x}, t))) v_i(\mathbf{x}) d\mathbf{x} + \int_{\partial K_i} \min(\tilde{\mathcal{H}}_{n_{K_i}}(\mathbf{x}, t), 0) [u_h](\mathbf{x}, t) v_i^-(\mathbf{x}) d\mathbf{s} \\ - C \Delta K_i \sum_{S_i' \in \partial K_i} \frac{1}{\Delta S_i'} \int_{S_i'} (\chi_{n_{K_i}}(\mathbf{x}, t) - |\tilde{\mathcal{H}}_{n_{K_i}}(\mathbf{x}, t)|) [\nabla_x u_h \cdot \mathbf{n}_{K_i}](\mathbf{x}, t) v_i^-(\mathbf{x}) d\mathbf{s} \\ - 2C \Delta K_i \sum_{\tilde{S}_i' \in \bar{\partial} K_i} \frac{1}{\Delta \tilde{S}_i'} \int_{\tilde{S}_i'} \min(\tilde{\mathcal{H}}_{n_{K_i}}^-(\mathbf{x}, t), 0) (\nabla_x u_h^- \cdot \mathbf{n}_{K_i})(\mathbf{x}, t) v_i^-(\mathbf{x}) d\mathbf{s} = 0, \tag{23}$$

for each $i \in \{1, ..., n\}$ and for any test function $v_i \in \mathcal{P}_i$,

where $\Delta K_i$ (respectively $\Delta S_i'$) is the size of the element $K_i$ (respectively the length of the edge $j$ of element $K_i$). The set $\partial K_i$ denotes the internal edges of element $K_i$ which are shared with other cells. The set $\bar{\partial} K_i$ denotes the external edges of element $K_i$ which are part of the domain boundary $\bar{\partial} \Omega$. The test functions $v_i$ are shape functions as usual for Galerkin approaches (Zienkewicz & Morgan 1983).

In Scheme (23), the following quantities are introduced:

$$\begin{align*}
\mathcal{H}_{K_i}^{\pm} &= \mathcal{H}(\mathbf{x}^{\pm}, \nabla_{\mathbf{x}} u_{hK_i}^{\pm}), \\
\mathcal{H}_{n_{K_i}} &= \nabla_{\nabla_{\mathbf{x}}} \mathcal{H} \cdot \mathbf{n}_{K_i}, \\
\mathcal{H}_{n_{K_i}}^{\pm} &= \mathcal{H}_{n_{K_i}}(\mathbf{x}^{\pm}, \nabla_{\mathbf{x}} u_{hK_i}^{\pm}), \\
\tilde{\mathcal{H}}_{n_{K_i}}(\mathbf{x}) &= \begin{cases} \frac{\mathcal{H}_{K_i}^{\pm} - \tilde{\mathcal{H}}_{K_i}^{-}}{|\nabla_{\mathbf{x}} u_h \cdot \mathbf{n}_{K_i}|}(\mathbf{x}), & \text{if } [\nabla_{\mathbf{x}} u_h \cdot \mathbf{n}_{K_i}](\mathbf{x}) \neq 0, \\ \frac{1}{2}(\mathcal{H}_{n_{K_i}}^{+} + \tilde{\mathcal{H}}_{n_{K_i}}^{-}), & \text{otherwise,} \end{cases} \\
\delta_{n_{K_i}}(\mathbf{x}) &= \max(0, \tilde{\mathcal{H}}_{n_{K_i}}(\mathbf{x}) - \mathcal{H}_{n_{K_i}}^{-}, \mathcal{H}_{n_{K_i}}^{+} - \tilde{\mathcal{H}}_{n_{K_i}}(\mathbf{x})), \\
\chi_{n_{K_i}}(\mathbf{x}) &= \max(\delta_{n_{K_i}}(\mathbf{x}), |\tilde{\mathcal{H}}_{n_{K_i}}(\mathbf{x})|). \tag{24}
\end{align*}$$

The first term of Scheme (23) ensures consistency. It embeds a weak formulation of the Hamilton-Jacobi partial differential equation inside each element. The scheme accounts for causality thanks to the quantity $\tilde{\mathcal{H}}_{n_{K_i}}$, referred to as the Roe speed, estimated across interfaces between elements. Its sign determines the information flow direction at an interface between two cells. At a point located at an interface between two cells $i$ and $j$, we have

$$\tilde{\mathcal{H}}_{n_{K_i}} = -\tilde{\mathcal{H}}_{n_{K_j}}. \tag{25}$$

The second term of Scheme (23) penalizes the jump of the solution at the interface, $[u_h]$, only in the cells where the Roe speed is negative

$$\tilde{\mathcal{H}}_{n_K} < 0. \tag{26}$$

In other words, the downwind cell receives information from the upwind cell, and the upwind cell is not influenced by the downwind cell. When flows from both cells, determined by $\mathcal{H}_{n_K}$, are oriented towards the other cell, the situation is equivalent to a shock in the case of hyperbolic conservation laws. If both flows are equivalent then the shock is captured at the interface. In this case, we have

$$\tilde{\mathcal{H}}_{n_{K_i}} = \tilde{\mathcal{H}}_{n_{K_j}} = 0, \tag{27}$$

therefore the second term of Scheme (23) is equal to zero in both cells. In the case where flows from both cells are inward, this is similar to what is called a rarefaction (see Qian & Symes (2001, appendix)). The entropy condition is violated in such cells.
Only in such a case, the third term of Scheme (23), referred to as the viscosity term, is non-zero. This yields a penalization on the jump of the normal component of spatial derivatives of the solution at the interface, $[\nabla_x u_h \cdot n_{K_i}]$, so as to correct the entropy violation. Thanks to this third term, the entropy fix is directly embedded in the scheme. This viscosity term is weighted by a ratio related to the geometry of the element, and balanced by an empirical constant $C$. Like Cheng & Wang (2014), We have verified that the choice $C = 0.25$ yields satisfactory results in our numerical experiments.

We have designed the fourth term as an additional term to the original scheme, acting on external edges only, in order to enforce suitable boundary conditions (see Section 4.2).

For numerical implementation, we use efficient quadrature rules with enough accuracy together with associated Gauss points the same way as in Cockburn & Shu (1998) to evaluate edge and volume integrals. We emphasize that Scheme (23) is compact in the sense that, when considering an element $K_i$, the only values needed from other elements are the values of $u_h$ and its derivatives at edges, and more precisely at the quadrature points of the edge integrals. For example, the mean values of the solution in the neighbouring elements are not required.

## 4.2 Boundary conditions and initial conditions

Two types of boundary conditions are imposed and should take benefit from the discontinuous finite-element formulation.

The first one is related to the outside limit of the domain $\partial\Omega$, expressed through the fourth term that we add into Scheme (23). When considering a propagation problem in a unbounded physical domain, we need to control that the outside of the computational domain has no influence on the inside. In other words, information must not enter the computational domain. In terms of traveltimes, all the information we put into the system comes from the source, and it propagates towards the outside. The fourth term, acting on external edges only, is in fact a combination of the second and third terms of the scheme. Heuristically, it aims at finding the flow direction, and penalizing the derivative of the solution at the boundary if the flow is inward. In other words, the Roe speed $H_{\text{RK}}^{-}$ is computed and a penalization term is applied if it is negative. As far as we know, this is the first implementation of such boundary conditions in a RKDG scheme for the Eikonal equation.

A second boundary condition needs to be imposed at the source point where the traveltime value is zero. The source point falls into one element of the mesh, referred to as the source element. Values of all the degrees of freedom of this source element are set to a pre-computed value and these values are kept unchanged during the computation process: they do contribute to flux estimations and, therefore, they are boundary conditions, while being set at the initial time. These values correspond to a projection of a specified solution around the source into the approximation space of the source element. When we use the factorization, then the solution $\tau$ we compute is a perturbation to the reference solution $u_0$ which is supposed to be a suitable approximation around the source. Consequently, the solution $\tau(\mathbf{x}, t)$ is set to zero everywhere inside the source element when the factorization is applied.

Finally, we use the reference solution $u_0(\mathbf{x})$ as the initial guess, so that we set the initial condition $\tau(\mathbf{x}, 0) = 0$.

## 4.3 Time discretization

Decomposing the numerical solution over each element $K_i$ in terms of degrees of freedom $u_i^j$ we get

$$u_{h|K_i}(x, z, t) = \sum_{j=1}^{d_i} u_i^j(t) \phi_i^j(x, z). \qquad (28)$$

Replacing $u_h$ by expression (28) in Scheme (23) yields a system of ODE in $\partial_i u_i^j(t)$, namely the time derivatives of the degrees of freedom. Time integration is then performed with a standard explicit RK scheme. Since we are looking for the steady-state, we do not need a high level of accuracy on the intermediate states, so that a second-order RK scheme is enough. In our numerical experiments, we observe that higher-order RK schemes do not modify the steady state while increasing the computational cost. In practice, we detect the steady state by comparing the relative evolution of the solution between two successive time steps.

## 4.4 Courant–Friedrichs–Lewy condition, best Hamiltonian choice

### 4.4.1 Formulation of the CFL condition

The stability of the time integration is constrained by a CFL condition, which can be a severe constraint in terms of computational cost as we are only interested in reaching the steady-state solution. This CFL condition establishes a proportional relationship between the maximum size of time steps and the characteristic length $\lambda$ of the mesh. For a regular triangular mesh the characteristic length $\lambda$ is generally taken as the radius of the circumcircle of a cell. For a rectangular Cartesian mesh, it is defined by half the length of the longest edge of a cell. To ensure stability when using a non-regular mesh, we must consider the cell where the CFL constraint is the strongest since the RK integration we perform is global. Therefore the characteristic length $\lambda$ to consider for the CFL condition is the smallest one among all cells of the mesh. The CFL constraint also depends on the polynomial degree $k$ we use for numerical approximation. We may write the CFL constraint into general form

$$\Delta t \le \frac{1}{2k+1} \lambda Q. \qquad (29)$$

The $\frac{1}{2k+1}$ factor comes from the analysis performed in the case of hyperbolic conservation laws (Cockburn & Shu 1989).
The $Q$ factor depends on the Hamiltonian we use and connects the constraint on the time steps to the propagation velocity of the solution. In other words, we need to ensure that the numerical solution is able to propagate as fast as the exact solution. The stability condition for RKDG schemes can be established in the case of hyperbolic conservation laws (see e.g., Cockburn & Shu (1989)) which write in 1-D

$$u_t + (f(u))_{,x} = 0. \qquad (30)$$

The $Q$ factor is given by

$$Q = \frac{\alpha}{\max |f'(u)|}, \qquad (31)$$

where $\alpha$ is a constant.

This can be extended to the case of Hamilton-Jacobi equations in one dimension with a convex Hamiltonian (Cheng & Shu 2007; Cheng & Wang 2014). We must consider the quantities $\max |\tilde{H}_{K_i}|$ at cell boundaries and $\max |\partial H(x, u_{,x})/\partial u_{,x}|$ inside cells, yielding

$$Q_{1D} = \frac{\alpha}{\max (\max |\partial H(x, u_{,x})/\partial u_{,x}|, \max |\tilde{H}_{K_i}|)} \qquad (32)$$

For 2-D and/or non-convex Hamiltonians, we have not found general proofs but we rely on the behaviours observed in our numerical tests. In a 2-D case, the vector $H$ must be considered, the components of which are the derivatives of the Hamiltonian with respect to both components of the gradient of $u$, namely,

$$\begin{aligned} H_1(x, z, u_{,x}, u_{,z}) &= \frac{\partial H(x, z, u_{,x}, u_{,z})}{\partial u_{,x}}, \\ H_2(x, z, u_{,x}, u_{,z}) &= \frac{\partial H(x, z, u_{,x}, u_{,z})}{\partial u_{,z}}. \end{aligned} \qquad (33)$$

A criterion for defining $Q$ in the 2-D case is given by Zhang et al. (2005b) in a finite-difference framework, using $H_1$ and $H_2$ as well as characteristic lengths of the grid along both $x$ and $z$ directions. However, since we want to use unstructured meshes, we need to keep a general criterion and consider the norm of the vector $H$, instead of considering $x$ and $z$ directions as well as $H_1$ and $H_2$ separately as in Zhang et al. (2005b). This yields

$$Q_{2D} = \frac{\alpha}{\max (\max ||H||, \max |\tilde{H}_{K_i}|)} \qquad (34)$$

Such a norm is considered for instance by Qian & Symes (2002b) for the CFL derivation in a paraxial approach, from where we find $\alpha = \sqrt{2}/2$ in the general 2-D case.

It is interesting to note that for some choices of the Hamiltonian, the vector $H$ is equivalent to the group velocity vector, which gives a physical meaning to the propagation of the solution inside the computational domain (Červený 2001, eq. 4.2.8). In our DG approach, we might also be careful about the quantity $\max |\tilde{H}_{n_{K_i}}|$ at cell boundaries, which appears in the scheme. In the isotropic case, it can be shown from the definitions (24) that values of $|\tilde{H}_{n_{K_i}}|$ at cell boundaries are bounded by those of $H_1$ and $H_2$ in neighbouring cells. We make the assumption that this also holds in the TTI case. The factor $Q$ then simplifies:

$$Q_{2D} = \frac{\alpha}{\max ||H||} \qquad (35)$$

### 4.4.2 Isotropic case

Isotropic Hamiltonian (3) gives the expression

$$\|H\| = 1 \quad \forall \mathbf{x}, \qquad (36)$$

whereas Hamiltonian (5) yields

$$\|H\| = c(\mathbf{x}). \qquad (37)$$

The factor $Q$ in eq. (29) is inversely proportional to the maximum of these quantities. We see that, in both cases, $Q$ does not depend on the solution $u$ nor its derivatives $u_j$. If we consider a heterogeneous medium, the global constraint on time steps in (37) is imposed by the highest speed value since

$$Q = \frac{\alpha}{\max_{\Omega} c(\mathbf{x})}. \qquad (38)$$

This maximum speed value occurs only in a subdomain of $\Omega$, so that the computation is not optimal in terms of number of iterations. For this reason, although (5) has the meaning of mimicking the wave-front evolution in the physical medium from the source, we prefer to use (3) for computational efficiency. The CFL constraint is uniform in space in such a formulation, due to the uniform propagation velocity of the solution.
### 4.4.3 VTI/TTI case

Performing the same analysis is less straightforward in the anisotropic case. For the sake of clarity, we first consider the VTI case. The differentiation of Hamiltonian (12) yields

$$ \mathcal{H}_1 = 2du_{,x}^2 + 2cu_{,x}^2 u_{,z}^2, \\ \mathcal{H}_2 = 2eu_{,x}^2 + 2cu_{,x}^2 u_{,z}^2, \qquad (39) $$

$$ \|\mathcal{H}\| = 2\sqrt{d^2 u_{,x}^2 + e^2 u_{,z}^2 + (2c(d+e) + c^2(u_{,x}^2 + u_{,z}^2))u_{,x}^2 u_{,z}^2}. \qquad (40) $$

The VTI Hamiltonian is not Lipschitz continuous. The value of $\|\mathcal{H}\|$ depends on the derivatives of the solution $u_{,x}$ and $u_{,z}$ in an unbounded way. This is not desirable because we cannot assign a value to $\Delta t$ in (29) since the Q factor is virtually equal to zero (eq. 35). For this reason we switch to another VTI Hamiltonian, which yields the same steady state as (12). The new VTI Hamiltonian is similar to the one given by Zhang et al. (2006, eq. 3.12) and writes

$$ \mathcal{H}_{\text{VTI}} = \frac{1}{\sqrt{d}} \left( \sqrt{\frac{1}{4}(du_{,x}^2 + eu_{,z}^2)^2 + cu_{,x}^2 u_{,z}^2 + \frac{1}{2}(du_{,x}^2 + eu_{,z}^2) - 1} \right). \qquad (41) $$

The equivalence of Hamiltonians (12) and (41) at steady state is demonstrated in Appendix B. As in the case of isotropic Hamiltonian (3), this new VTI Hamiltonian (41) is Lipschitz continuous of Lipschitz constant 1, allowing a predictable stable $\Delta t$ for integration. This property is demonstrated in Appendix C. The TTI case is retrieved by introducing the tilt angle $\theta$ as in (13). Finally, the factorization strategy for Hamiltonian (41) is the same as for Hamiltonian (12). It is now possible to estimate a time step $\Delta t$ for performing stable computation in arbitrary heterogeneous TTI media.

We have shown how to optimize the CFL constraint in both isotropic and anisotropic cases. We emphasize that the use of the additive factorization technique does not change these conclusions.

## 5 NUMERICAL RESULTS

In the first case study, we illustrate the convergence properties of the RKDG scheme in a smooth isotropic medium where the exact solution is known. We focus on $L^2$ error computed on the solution and its derivatives, convergence orders, and point-source treatment (no treatment, enforcement of the exact solution on a fixed area around the source, point-source factorization). An analysis is performed for several polynomial orders in both Cartesian and triangular configurations. We also illustrate the impact of the choice of the Hamiltonian, and a comparison is performed with a finite-difference code from Noble et al. (2014).

In the second case study, we test a higher-order point-source factorization where we use an exact solution known in a virtual medium the speed of which is similar at the source to two first terms of the power series expansion of the velocity model of the real medium at the source. In other words, the new reference solution accounts for the value of the speed at the source but also for its gradient. Doing so, we show that we gain one more order of convergence.

The third case study illustrates the accuracy of the method in a homogeneous TTI case. We compare our results with a finite-difference code from Waheed et al. (2015b) and Tavakoli F. et al. (2015) with respect to an exact parametric solution. We exhibit the efficiency of the point-source factorization in an aneliptical case.

With a volcano model, the fourth case study illustrates the flexibility of the method for handling complex topographies. We test the scheme on an unstructured mesh to highlight the ability of our method to deal with complex topographies.

The fifth case study is an application inside a realistic TTI medium (BP TTI).

For all these test cases, the parameter C appearing in Scheme (23) is set to the value 0.25.

### 5.1 Case study 1: smooth isotropic velocity

In this first case study, the computation is performed inside a 4 × 4 km square. Inside this domain, a constant vertical gradient is defined for the speed $c(z)$ (km s$^{-1}$), such that

$$ c(z) = 1 + 0.5z. \qquad (42) $$

The point source is located in the centre ($x_s = 2$ km, $z_s = 2$ km). The analytical solution of the Eikonal equation as well as its spatial derivatives are known in this case (Fomel et al. 2009). Isochrones of the solution are shown in Fig. 1 together with the velocity model. The $L^2$ error is
<figure><img src="image_4.png" /><figcaption>Figure 1. Constant vertical gradient of velocity model of Case study 1. Left: velocity model. Right: isochrones of the analytical solution at different times (s).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 2. L² error of the solution (left) and L² error of the x-derivative of the solution (right) with respect to the number of degrees of freedom (#dof) in the Cartesian case for P¹, P² and P³ polynomial approximations, for Case study 1 with a source area of radius 0.4 km: optimal k + 1 convergence of the solution and k convergence of its x-derivative in a setting with no point source.</figcaption></figure>

computed as

$$L^2 \text{ error} = \sqrt{\frac{\sum_{i=1}^{K} \sum_{j=1}^{G} (u_h(\mathbf{x}_i^j) - u_a(\mathbf{x}_i^j))^2}{\sum_{i=1}^{K} \sum_{j=1}^{G} u_a^2(\mathbf{x}_i^j)}}, \quad (43)$$

where the values of the numerical solution $u_h$ are compared to those of the analytical solution $u_a$ at each of the G Gauss points of each of the K elements of the discretized domain. The L² errors of the x- and z-derivatives of the solution are also computed in the same way, replacing $u_h$ and $u_a$ by their x- and z-derivatives respectively in (43). The domain is first discretized in a rectangular Cartesian frame with $N_x = N_z = N$ elements along x- and z-axes so that the total number of elements is $N^2$.

In order to avoid the point-source numerical issue and first illustrate the convergence behaviour of the solver, a first set of simulations is performed with an extended source area. This area contains all the elements distant of less than 0.4 km from the source point. These source elements are excluded from the computational domain, and the analytical solution is used for computing their contribution in the integrals at the edges shared with the elements within the computational domain. Hamiltonian (3) without factorization is considered. Fig. 2 illustrates the optimal k + 1 convergence orders reached for P¹, P² and P³ approximation spaces regarding the L² error of the numerical solution. The optimal k convergence orders for the x-derivative of the solution can also be observed. Here and after, results concerning the z-derivative of
<table><caption>**Table 1.** Number of degrees of freedom (#dof), L² error of the solution and convergence orders for several values of N in the Cartesian case for P¹, P² and P³ polynomial approximations, for Case study 1 without factorization: non-optimal first-order-only convergence due to the source singularity.</caption><thead><tr><th rowspan="2">N</th><th colspan="3">p¹</th><th colspan="3">p²</th><th colspan="3">p³</th></tr><tr><th>#dof</th><th>L² Error</th><th>Order</th><th>#dof</th><th>L² Error</th><th>Order</th><th>#dof</th><th>L² Error</th><th>Order</th></tr></thead><tbody><tr><td>20</td><td>1200</td><td>4.84E-03</td><td></td><td>2400</td><td>1.77E-03</td><td></td><td>4000</td><td>7.18E-04</td><td></td></tr><tr><td>40</td><td>4800</td><td>2.24E-03</td><td>1.11</td><td>9600</td><td>8.60E-04</td><td>1.04</td><td>16000</td><td>2.18E-04</td><td>1.72</td></tr><tr><td>80</td><td>19200</td><td>1.09E-03</td><td>1.03</td><td>38400</td><td>4.29E-04</td><td>1.00</td><td>64000</td><td>8.17E-05</td><td>1.41</td></tr><tr><td>160</td><td>76800</td><td>5.45E-04</td><td>1.01</td><td>153600</td><td>2.15E-04</td><td>1.00</td><td>256000</td><td>3.69E-05</td><td>1.15</td></tr><tr><td>320</td><td>307200</td><td>2.72E-04</td><td>1.00</td><td>614400</td><td>1.08E-04</td><td>1.00</td><td>1024000</td><td>1.80E-05</td><td>1.04</td></tr></tbody></table>

<table><caption>**Table 2.** Number of degrees of freedom (#dof), L² error of the solution and convergence orders for several values of N in the Cartesian case for P¹, P² and P³ polynomial approximations, for Case study 1 with factorization: second-order convergence is achieved, which is optimal in the P¹ case.</caption><thead><tr><th rowspan="2">N</th><th colspan="3">p¹</th><th colspan="3">p²</th><th colspan="3">p³</th></tr><tr><th>#dof</th><th>L² Error</th><th>Order</th><th>#dof</th><th>L² Error</th><th>Order</th><th>#dof</th><th>L² Error</th><th>Order</th></tr></thead><tbody><tr><td>20</td><td>1200</td><td>5.24E-04</td><td></td><td>2400</td><td>1.30E-04</td><td></td><td>4000</td><td>2.21E-05</td><td></td></tr><tr><td>40</td><td>4800</td><td>1.27E-04</td><td>2.05</td><td>9600</td><td>3.30E-05</td><td>1.98</td><td>16000</td><td>5.59E-06</td><td>1.98</td></tr><tr><td>80</td><td>19200</td><td>3.11E-05</td><td>2.03</td><td>38400</td><td>8.29E-06</td><td>1.99</td><td>64000</td><td>1.41E-06</td><td>1.99</td></tr><tr><td>160</td><td>76800</td><td>7.71E-06</td><td>2.01</td><td>153600</td><td>2.08E-06</td><td>2.00</td><td>256000</td><td>3.52E-07</td><td>2.00</td></tr><tr><td>320</td><td>307200</td><td>1.92E-06</td><td>2.01</td><td>614400</td><td>5.20E-07</td><td>2.00</td><td>1024000</td><td>8.81E-08</td><td>2.00</td></tr></tbody></table>

the solution are not presented, however the z-derivative always yields the same convergence orders and thus the same conclusions as for the x-derivative.

As discussed in Section 3, in a realistic application, the exact solution around the source might not be analytically calculated, so that the previous treatment is not applicable. The above convergence behaviour is expected to collapse when the point-source singularity is introduced. If no special treatment is performed at the source point, when using Hamiltonian (3), Table 1 exhibits the expected non-optimal first-order-only convergence of the solver whatever the degree of polynomials used. However, the error magnitude decreases when the polynomial order increases for a given N, as well as for a given number of degrees of freedom.

When the factorization technique is applied with the use of Hamiltonian (16), a gain of one order of convergence is observed, as shown in Table 2. For P¹, P² and P³ polynomial approximations a second-order accuracy is achieved, which is optimal for P¹. Here again, even if the convergence orders are the same, the error magnitude decreases when the polynomial order increases.

When looking at the spatial derivatives of the solution, with no factorization, our experiments exhibit a degenerated first-order convergence which is controlled by the first-order-only convergence of the solution itself. When the factorization technique is used, an optimal first-order convergence for the x-derivative is observed in the P¹ case. Degenerated second-order convergences are observed for P² and higher-order polynomial spaces, dominated by the second-order convergence of the solution itself (Fig. 3).

In conclusion, the factorization technique allows for an optimal P¹ second-order solver as well as non-optimal second-order P² and P³ solvers. The P² solver nearly reaches second-order optimality in terms of derivatives. In all the cases, the factorization yields a significant decrease of the error magnitude.

Similar results are obtained when a structured triangulation of the domain is used. The Union-Jack (UJ) pattern is obtained from the Cartesian grid by cutting each rectangular element into two triangles in an alternating diagonal direction, as shown in Fig. 4. With $N_x = N_z = N$,

<figure><img src="image_6.png" /><figcaption>**Figure 3.** L² error of the solution (left) and L² error of the x-derivative of the solution (right) with respect to the number of degrees of freedom (#dof) in the Cartesian case for P¹ and P² polynomial approximations, for Case study 1 with factorization (fact.) and without factorization (no fact.). Left: non-optimal first-order convergence of the solver without factorization; second-order convergence with factorization (optimal for P¹). Right: degenerated first-order convergence of the x-derivative without factorization; optimal fully first-order convergence (P¹) and nearly second-order convergence (P²) with factorization.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 4. The Union–Jack triangulation shown for $N_x = N_z = 10$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 5. L² error of the numerical solution of Case study 1 with P² polynomial approximations, in both Cartesian and UJ triangular cases, and both without factorization (Cart. and UJ) and with factorization (Cart. fact. and UJ fact.). The first-order convergence of the standard case is improved to a second-order convergence when the factorization is applied. Both Cartesian and UJ discretizations achieve the same orders, although the magnitude of error is higher in the UJ case.</figcaption></figure>

the number of elements of such a triangulation is now $2N^2$. Results are the same as in the Cartesian case in terms of convergence orders. However, the magnitude of error with respect to the number of degrees of freedom is higher, which is illustrated in the P² case in Fig. 5.

Finally, our results in the Cartesian P¹ discretization with the factorization technique are compared to those obtained with the 2-D version of the Finite-Difference (FD) solver from Noble et al. (2014). The FD solver is based on local operators which also make use of the factorization principles. Fig. 6 highlights the much lower error level of the DG solver compared to the one of the FD solver with respect to the number of degrees of freedom. The different slopes of the two lines exhibit the first-order convergence of the FD solver compared to the second-order convergence of the DG solver. The huge difference between the two lines has to be balanced by the fact that the FD solver is a faster solver, thanks to the use of a Fast Sweeping strategy.

## 5.2 Case study 2: higher-order factorization in a constant gradient of isotropic squared slowness field

The purpose of this case study is to exhibit a way to gain one more order of convergence thanks to a higher-order source factorization. The principle holds in choosing a suitable reference solution $u_0(x)$ accounting for higher-order terms of the power series expansion of the velocity model at the source point. Luo et al. (2014b) proposed this kind of approach for the squared slowness and the squared Eikonal.

In this simple example, the same domain of computation and source location as in the first case study are considered, but with a different speed. Slowness $s(x)$ is defined by

$$s(x) = 1/c(x). \qquad (44)$$
<figure><img src="image_9.png" /><figcaption>**Figure 6.** L² error of the numerical solution of Case study 1 with the DG solver, P¹ approximation, Cartesian discretization and factorization (DG), compared to the FD solver from Noble *et al.* (2014).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** L² error of the solution (left) and L² error of the x-derivative of the solution (right) with respect to the number of degrees of freedom (#dof) in the Cartesian case for P¹, P² and P³ polynomial approximations, for Case study 2 with a second-order factorization. Left: optimal second-order convergence (P¹); optimal third-order convergence (P²); non-optimal third-order convergence (P³). Right: optimal first-, second- and third-order convergences of the x-derivative for P¹, P² and P³, respectively.</figcaption></figure>

The medium is defined by a constant vertical gradient of the squared slowness (slowness unit s.km⁻¹):

$$s^2(z) = 0.25 - 0.1(z - z_s), \qquad (45)$$

where the depth is expressed in km. The slowness is indeed strictly positive inside the domain. The analytical solution of the problem is known (see e.g. Fomel *et al.* (2009)) as well as its spatial derivatives. We can write the speed as

$$c(z) = \frac{1}{\sqrt{0.25 - 0.1(z - z_s)}}, \qquad (46)$$

which can be expanded around the source point as

$$c(z) = 2 + 0.4(z - z_s) + O((z - z_s)^2). \qquad (47)$$

The standard factorization technique would account for the first term of the above expansion, using the solution in a constant velocity model of value $c_0 = 2$ km s⁻¹ as the reference solution. Here we propose to account for both terms of the expansion using the exact solution in a constant gradient of velocity model as the reference solution. Therefore in this example $u_0(x)$ is the analytical solution for the point-source problem in a velocity model

$$c_0(z) = 2 + 0.4(z - z_s), \qquad (48)$$

with units of the eq. (45). Computations are performed with P¹, P² and P³ in Cartesian meshes and convergence behaviours are compatible with the order of selected polynomials (Fig. 7). Again, results concerning the z-derivative of the solution are not presented for the sake of concision, but we mention that it yields the same convergence orders and thus the same conclusions as the x-derivative.
<figure><img src="image_11.png" /><figcaption>**Figure 8.** At any point $(x, z)$ of a homogeneous TTI medium, the phase (ray) angle $\psi$ defined by $x, z$ and the source point differs from the group angle $\varphi$ which is defined by the normal to the wave front, except on the axes where they are equal. The parametric relationship between point coordinates and traveltime allows for explicitly computing a phase angle from a given group angle. We solve the inverse problem iteratively.</figcaption></figure>

This case study illustrates that this second-order factorization allows for reaching third-order convergence. In a general case, for a given velocity model which can be expanded in power series at the source point, we are able to design a second-order reference solution based on the analytical solution in a constant gradient of speed oriented in the direction of the gradient of the true velocity model. The resulting third-order scheme with a $P^2$ approximation is particularly advantageous for applications requiring quantities related to first- and second-order derivatives of the traveltime, such as amplitudes and take-off angles (Luo *et al.* 2012, 2014a).

## 5.3 Case study 3: homogeneous TTI velocity model

In this case study, the ability of our solver to compute traveltimes in a TTI medium using the tilted formulation of Hamiltonian (41) as well as the corresponding factorized formulation is validated. The medium is defined by the homogeneous fields of vertical speed, Thomsen's parameters and tilt angle, respectively $V_p = 2$ km s$^{-1}$, $\epsilon = 0.4$, $\delta = 0.2$ and $\theta = 40$ deg. The computational domain is a rectangle of 32 km length and 8 km depth, and the point source is located at $x = 2.025$ km, $z = 2.025$ km. For this anelliptic homogeneous medium, the exact solution is not known directly but a parametric formulation allows for computing the position of the wave front at a given time for a given phase angle. This formulation was established by Payton (1983, eqs 2.8.8 and 2.8.9) and by Carcione *et al.* (1988, eq. 5.9). Hence we can build the wave front at a given time $t$ with a dense sampling of phase velocities, and visually compare the isochrones of the traveltime maps computed by the solver with these wave fronts.

The medium is discretized in a Cartesian frame with $N_x = 640$ and $N_z = 160$ and a $P^1$ approximation, yielding 307 200 degrees of freedom. The first computation is performed without any factorization. The second computation uses the factorization technique with the elliptical reference solution from (19) and (20). Since $\epsilon \neq \delta$, the medium is anelliptic, so that the elliptical reference solution does not account for the whole source singularity. However, we observe an improvement of the solution compared to the first computation.

In the third computation, we pre-compute values for the exact anelliptical solution and its spatial derivatives at the points where they are required, and we use them as the reference solution. Since the reference solution and its derivatives are involved inside the integrals of Scheme (23), we must obtain their values at each Gauss point necessary for integral estimation by quadrature rule. We pre-compute these values by retrieving the group angle from the phase angle at a given point iteratively in a similar way as in Qian & Symes (2001, alg. 1), using the parametric formulation that we mentioned above (see Fig. 8). We then approximate $x$- and $z$-derivatives using central differences. Once the required values are pre-computed and properly stored, we are able to proceed with the DG solver which calls these values when needed.

Finally, a finite-difference solution is computed inside the same medium using the iterative fast-sweeping factored TTI Eikonal solver detailed in Waheed *et al.* (2015b) and Tavakoli F. *et al.* (2015). For that purpose, a finite-difference grid composed of 277 × 1105 points is considered, so that the number of degrees of freedom is equivalent to the DG discretization: 306 085.

Analytical wave fronts in the whole medium are plotted in Fig. 9. A zoom on the isochrone corresponding to the time $t = 2$ s is shown in Fig. 10. There is an obvious improvement between the first computation with no factorization and the second computation which uses the factorization with an elliptical reference solution. The third plot is nearly mingled with the exact solution. This illustrates the great advantage of using pre-computed anelliptical values as the reference solution. Finally, the FD computation, with a comparable number of degrees of freedom, exhibits a solution which is more or less equivalent to the DG computation with the elliptical solution as reference for the factorization. These results illustrate the good behaviour of our solver regarding the TTI configuration. Factorization yields a great
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Isochrones of the solutions for the homogeneous TTI medium of Case study 3. Analytical and numerical solutions are superimposed. Isochrones are plotted every second.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Isochrones at *t* = 2 s of the solutions for the homogeneous TTI medium of Case study 3. 1: DG computation with no factorization. 2: DG computation with elliptical reference solution. 3: DG computation with anelliptical reference solution. 4: Analytical solution. 5: FD computation.</figcaption></figure>

improvement of the solution, and we have shown that it is possible to use a highly accurate approximation of the anelliptical exact solution as the reference solution for the factorization.

## 5.4 Case study 4: volcano with variable topography

The flexibility of the DG approach and its ability to handle topographies is illustrated with a Gaussian topography simulating a volcanic dome. The TTI model is shown in Fig. 11. A magma chamber is represented by a high-speed zone located at the vertical under the dome. The speed reaches 4 km s⁻¹ at its maximum, while it is 2 km s⁻¹ away from the high-speed structure. Thomsen's parameters ϵ and δ take values respectively 0.4 and 0.2 away from the high-speed structure, while their values decrease inside the magma chamber, reaching 0 at their minimum. The tilt angle θ follows the topography, its absolute value decreases with depth until 0 at z = 4 km. To compute traveltimes in this medium, a triangular mesh composed of 105442 P¹ triangular elements is built; a subset of the mesh is shown in Fig. 12. Computations are performed for various source locations, at the surface as well as in depth. Isochrones are superimposed over the velocity model in Fig. 13. In Fig. 14, isochrones are superimposed over the maps of the spatial derivatives of the traveltime, for a source located at x_s = 3.04 km, z_s = 0.63 km. For this source location, a line of discontinuity of both derivatives is visible on these maps. This discontinuity matches with a singular part of the traveltime field, namely the kink at the junction between the direct wave on the right part and the refracted wave passing through the high-speed structure.

## 5.5 Case study 5: realistic TTI model

In this last case study, we test our solver on the 2-D BP TTI benchmark model from Shah (2007), which is used in the geophysics community for testing and validating modelling tools. The model is described by a highly contrasted P-wave velocity over a distance of 79 km and a depth of 11 km and corresponding Thomsen's parameters ϵ, δ and tilt angle θ shown in Fig. 15. In order to mitigate the impact of the discretization, the original model has been smoothed with Gaussian characteristic lengths of 200 m in both x and z directions. The DG solver proceeds over a Cartesian mesh of 179 200 P¹ elements with Nₓ = 1120 and Nᵧ = 160. The corresponding number of degrees of freedom is 537 600. The same finite-difference solver as in the previous case studies is used for comparison purpose, which proceeds over a 1938 × 278 grid, so that the number of degrees of freedom is similar (538 764). We proceed with various source locations, at the surface as well as in depth. Isochrones obtained with the two methods are shown in Fig. 16. In all the cases, they are very similar. The complexity of the medium yields complex propagation phenomena illustrated by the tortuous shape of the isochrones in some parts of the medium. Among many phenomena, we can mention that a refracted wave is visible on top-left part of the first panel (source located at x_s = 0 km, z_s = 0 km) due to the high-speed salt body at x = 7 km.
<figure><img src="image_14.png" /><figcaption>**Figure 11.** The volcano TTI model of Case study 4: (a) vertical velocity, (b) Thomsen's epsilon parameter, (c) Thomsen's delta parameter, (d) tilt angle.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Detail of the triangular mesh built for the volcanic structure of Case study 4. The topography is sampled by edges of triangles.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Isochrones computed in the volcano model of Case study 4 for several source positions in km: (a) $x_s = 0, z_s = 1.3$, (b) $x_s = 1.7, z_s = 2.2$, (c) $x_s = 2, z_s = 0$, (d) $x_s = 3.04, z_s = 0.63$, (e) $x_s = 3.3, z_s = 3$, (f) $x_s = 1, z_s = 4$, (g) $x_s = 2, z_s = 3.7$. Isochrones are plotted every 0.1 second.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 14.** Maps of spatial derivatives of the traveltime and isochrones superimposed (blue line), computed in the volcano model of Case study 4 for a source point located at $x_s = 3.04, z_s = 0.63$ km. Top: x-derivative. Bottom: z-derivative. Isochrones are plotted every 0.1 second.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure 15.** Smoothed 2-D BP TTI model of Case study 5. From top to bottom: (a) vertical velocity, (b) Thomsen's epsilon parameter, (c) Thomsen's delta parameter, (d) tilt angle.</figcaption></figure>

Isochrones are superimposed over the maps of the spatial derivatives of the traveltime in Fig. 17, for a source located at $x_s = 33$ km, $z_s = 3$ km. A detailed view in a smaller area is shown in Fig. 18. Although the traveltime field is continuous, discontinuities of the spatial derivatives of the traveltime are prominent on these maps. They match with angles observed in the isochrones (singularities). The singularities of the solution are due to the viscosity solution which selects the lowest traveltime value where different phases compete. A line of discontinuity occurs where two branches of the solution meet (shock). Since the directions of propagation from both sides differ, the traveltime derivatives are discontinuous at such shock. The resolution of a shock is related to the size of the element inside which it occurs. Outside of this element, the solution is not affected. We have an illustration of the viscosity solution.

The same computation is performed using $P^2$ elements, with $N_x = 791$ and $N_z = 113$, yielding 536 298 degrees of freedom. Profiles of the traveltime and its derivatives along $x = 47$ km and $z = 6.7$ km are shown in Fig. 19. In Fig. 20, the profiles along $x = 47$ km of the traveltime derivatives computed with $P^1$ and $P^2$ are compared near the shock occurring at 6.8 km. $P^1$ yields a piecewise constant approximation of the derivatives, whereas $P^2$ yields a piecewise linear one. Regarding the derivatives, the shock is poorly approximated inside the element where it occurs but it does not affect the solution elsewhere. For practical applications, one has to be careful if the solution inside such an element is needed. A criterion based on the variations of the derivatives of the solution could be designed in order to refine the mesh and recompute locally the solution with a better resolution of the shock.

In a finite-difference approach, the traveltime field, as well as its derivatives, are computed at grid points, then interpolated if needed somewhere else, whereas in the DG approach, these quantities are computed everywhere inside each element to a given order of accuracy. We emphasize that no interpolation is required to obtain the maps of Fig. 17 nor the profiles of Figs 19 and 20.

# 6 CONCLUDING REMARKS

We have presented an RK DG method for solving time-dependent Hamilton–Jacobi equations with a point-source condition. The steady state is the solution to the corresponding static Eikonal equation. This high-order accurate, compact and flexible method computes traveltimes and its spatial derivatives in heterogeneous anisotropic media. We emphasize that Scheme (23) is written in a general Hamiltonian formulation, thus it may hold for a large variety of Hamiltonians, opening doors to other types of anisotropy. The extension to 3-D does not introduce new specific methodological issues and is part of a forthcoming work.

The method we have presented is computationally intensive. At each time iteration, every cell of the mesh is updated, while only a small subset of cells is concerned by the front propagation. To give a comparison, the DG computation for one source in Case study 5 took about
<figure><img src="image_19.png" /><figcaption>**Figure 16.** Isochrones computed in the smoothed BP TTI model of Case study 5 for several source positions in km: (a) $x_s = 0, z_s = 0$, (b) $x_s = 0, z_s = 10$, (c) $x_s = 33, z_s = 3$, (d) $x_s = 33, z_s = 11$, (e) $x_s = 48, z_s = 0$, (f) $x_s = 48, z_s = 11$, (g) $x_s = 72, z_s = 11$, (h) $x_s = 78, z_s = 0$. Blue plain line: DG P¹ computation. Red dashed line: FD computation. Isochrones are plotted every second.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 17.** Maps of spatial derivatives of the traveltime and isochrones superimposed (blue line), computed in the smoothed BP TTI model of Case study 5 for a source point located at $x_s = 33, z_s = 3$ km. P¹ computation with $N_x = 1120$ and $N_z = 160$. Top: x-derivative. Bottom: z-derivative. Isochrones are plotted every second.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**Figure 18.** Maps of spatial derivatives of the traveltime and isochrones superimposed (blue line), computed in the smoothed BP TTI model of Case study 5 for a source point located at $x_s = 33$, $z_s = 3$ km. P¹ computation with $N_x = 1120$ and $N_z = 160$. Zoom in the [28, 55] x [0, 11] rectangle. Top: x-derivative. Bottom: z-derivative. Isochrones are plotted every second.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Figure 19.** Profiles computed in the smoothed BP TTI model of Case study 5 for a source point located at $x_s = 33$, $z_s = 3$ km. P² computation with $N_x = 791$ and $N_z = 113$. Profile along $x = 47$ km: (a) Traveltime field; (b) x-derivative (plain line) and z-derivative (dashed line) of the traveltime. Profile along $z = 6.7$ km: (c) Traveltime field; (d) x-derivative (plain line) and z-derivative (dashed line) of the traveltime. These plots highlight the smoothness of the traveltime field compared to its discontinuous derivatives.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>**Figure 20.** Profile of the traveltime derivatives along *x* = 47 km computed in the smoothed BP TTI model of Case study 5 for a source point located at *x*<sub>s</sub> = 33, *z*<sub>s</sub> = 3 km. P<sup>1</sup> computation (blue line) and P<sup>2</sup> computation (red line) with similar numbers of degrees of freedom. Top lines: x-derivative; bottom lines: z-derivative. P<sup>1</sup> yields a piecewise constant approximation of the derivatives, whereas P<sup>2</sup> yields a piecewise linear one. Note the local variation of the viscous solution quite sensitive to the element size and the polynomial interpolation while the solution accuracy is not impacted elsewhere.</figcaption></figure>

fifty minutes on a 2.6 GHz machine with 8 GB of RAM, while the FD computation took 20 seconds with the same number of degrees of freedom. However, in order to provide a more meaningful comparison between methods, we will investigate on the following computational improvements in the future.

We first highlight that, in (23), the evolution of the degrees of freedom in time $\partial_t u_i^f(t)$ of element $K_i$ can be computed without knowledge of any $\partial_t u_k^f(t)$, $k \neq i$. Therefore, the DG scheme is easily parallelizable in space as a block-Jacobi method: one can compute $\partial_t u_i^f(t)$ independently for each element $K_i$ before updating the solution everywhere.

Another way to deal with the computational cost will be through a different solver for the DG discretization. The new solver would work as a block-Gauss-Seidel method, exploiting the directions of propagation of the solution by a suitable ordering of the elements. The elements of the mesh would then be updated sequentially instead of all together. Some related works in Zhang *et al.* (2005b); Li *et al.* (2008); Zhang *et al.* (2011) draw a path for the integration of a RKDG local solver into a global sweeping procedure in order to reach more rapidly the steady state. Zhang *et al.* (2005b) embedded a finite-difference fixed-point iterative method inside a Gauss-Seidel sweeping method in order to accelerate the convergence. Extensions of fast-sweeping methods to DG discretizations are performed in Li *et al.* (2008); Zhang *et al.* (2011), but only for piecewise-linear Cartesian (implicit) isotropic cases. Applying such a fast-sweeping approach to our DG scheme, we would quickly reach the steady state and hence reduce drastically the computation time, which could eventually attain the same order of magnitude as the FD computation. These expectations prompt us to further investigate this topic in the future.

## ACKNOWLEDGEMENTS

This study was partially funded by the SEISCPEO consortium (http://seiscope2.osug.fr), sponsored by AKERBP, CGG, CHEVRON, EXXON-MOBIL, JGI, SHELL, SINOPEC, STATOIL, TOTAL and WOODSIDE. This study was granted access to the HPC resources of the Froggy platform of the CIMENT infrastructure (https://ciment.ujf-grenoble.fr), which is supported by the Rhône-Alpes region (GRANT CPER07_13 CIRA), the OSUG@2020 labex (reference ANR10 LABX56) and the Equip@Meso project (reference ANR-10-EQPX-29-01) of the programme Investissements d'Avenir supervised by the Agence Nationale pour la Recherche, and the HPC resources of CINES/IDRIS/TGCC under the allocation 046091 made by GENCI.

## REFERENCES

Alkhalifah, T., 2000. An acoustic wave equation for anisotropic media, *Geophysics*, **65**, 1239–1250.

Benamou, J.D., 2003. An introduction to Eulerian geometrical optics (1992–2002), *J. Sci. Comput.*, **19**(1), 63–93.

Billette, F. & Lambaré, G., 1998. Velocity macro-model estimation from seismic reflection data by stereotomography, *Geophys. J. Int.*, **135**(2), 671–680.

Boué, M. & Dupuis, P., 1999. Markov chain approximations for deterministic control problems with affine dynamics and quadratic cost in the control, *SIAM J. Numer. Anal.*, **36**(3), 667–695.

Carcione, J., Kosloff, D. & Kosloff, R., 1988. Wave-propagation simulation in an elastic anisotropic (transversely isotropic) solid, *Q. J. Mech. Appl. Math.*, **41**(3), 319–345.

Červený, V., 2001. *Seismic Ray Theory*, Cambridge Univ. Press.

Cheng, Y. & Shu, C.-W., 2007. A discontinuous Galerkin finite element method for directly solving the Hamilton–Jacobi equations, *J. Comput. Phys.*, **223**(1), 398–415.

Cheng, Y. & Wang, Z., 2014. A new discontinuous Galerkin finite element method for directly solving the Hamilton–Jacobi equations, *J. Comput. Phys.*, **268**, 134–153.
Cockburn, B. & Shu, C.-W., 1989. TVB Runge-Kutta local projection dis-continuous Galerkin finite element method for conservation laws II: Gen-eral framework, *Math. Comput.*, **52**, 411–435.

Cockburn, B. & Shu, C.-W., 1998. The Runge-Kutta discontinuous Galerkin method for conservation laws V, *J. Comput. Phys.*, **141**(2), 199–224.

Courant, R. & Hilbert, D., 1966. *Methods of Mathematical Physics*, John Wiley.

Crandall, M.G. & Lions, P.L., 1983. Viscosity solutions of Hamilton-Jacobi equations, *Trans. Am. Math. Soc.*, **277**(1), 1–42.

Fomel, S., Luo, S. & Zhao, H.-K., 2009. Fast sweeping method for the factored Eikonal equation, *J. Comput. Phys.*, **228**, 6440–6455.

Han, S., Zhang, W. & Zhang, J., 2017. Calculating *qp*-wave travel times in 2D TTI media by high-order fast sweeping methods with a numerical quartic equation solver, *Geophys. J. Int.*, **210**(3), 1560–1569.

Hole, D. & Zelt, B., 1995. 3-D finite difference reflection traveltimes, *Geophys. J. Int.*, **121**, 427–434.

Hu, C. & Shu, C.-W., 1999. A discontinuous Galerkin finite element method for Hamilton-Jacobi equations, *SIAM J. Sci. Comput.*, **21**, 666–690.

Kao, C.Y., Osher, S. & Qian, J., 2004. Lax-friedrichs sweeping schemes for static hamilton-jacobi equations, *J. Comput. Phys.*, **196**, 367–391.

Kimmel, R. & Sethian, J.A., 1998. Computing geodesic paths on manifolds, *Proc. Natl. Acad. Sci. USA*, **95**, 8431–8435.

Kimmel, R. & Sethian, J.A., 2001. Optimal algorithm for shape from shading and path planning, *J. Math. Imaging Vis.*, **14**(3), 237–244.

Le Meur, H., 1994. Tomographie tridimensionnelle à partir des temps des premières arrivées des ondes P et S, *PhD thesis*, Université Paris VII.

Le Meur, H., Virieux, J. & Podvin, P., 1997. Seismic tomography of the gulf of Corinth: a comparison of methods, *Ann. Geofis.*, **XL**(1), 1–24.

Leung, S. & Qian, J., 2006. An adjoint state method for three-dimensional transmission traveltime tomography using first-arrivals, *Commun. Math. Sci.*, **4**(1), 249–266.

Levin, F.K., 1979. Reply to K. Helbig by F. K. Levin, *Geophysics*, **44**(5), 990–990.

Li, F., Shu, C.W., Zhang, Y.T. & Zhao, H., 2008. A second order discontinuous Galerkin fast sweeping method for Eikonal equations, *J. Comput. Phys.*, **227**, 8191–8208.

Lions, P.-L., 1982. *Generalized Solutions of Hamilton-Jacobi Equations*, Pitman.

Luo, S. & Qian, J., 2011. Factored singularities and high-order Lax-Friedrichs sweeping schemes for point-source traveltimes and amplitudes, *J. Comput. Phys.*, **230**, 4742–4755.

Luo, S. & Qian, J., 2012. Fast sweeping method for factored anisotropic Eikonal equations: multiplicative and additive factors, *J. Sci. Comput.*, **52**, 360–382.

Luo, S., Qian, J. & Zhao, H.-K., 2012. Higher-order schemes for 3D first-arrival traveltimes and amplitudes, *Geophysics*, **77**, T47–T56.

Luo, S., Qian, J. & Burridge, R., 2014a. Fast Huygens sweeping methods for Helmholtz equations in inhomogeneous media in the high frequency regime, *J. Comput. Phys.*, **270**, 378–401.

Luo, S., Qian, J. & Burridge, R., 2014b. High-order factorization based high-order hybrid fast sweeping methods for point source Eikonal equations, *SIAM J. Numer. Anal.*, **52**, 23–44.

Moser, T.J., 1991. Shortest path calculation of seismic rays, *Geophysics*, **56**(1), 59–67.

Noble, M., Gesret, A. & Belayouni, N., 2014. Accurate 3-D finite difference computation of travel time in strongly heterogeneous media, *Geophys. J. Int.*, **199**, 1572–1585.

Osher, S., 1993. A level set formulation for the solution of the Dirichlet problem for Hamilton-Jacobi equations, *SIAM J. Math. Anal.*, **24**, 1145–1152.

Payton, R., 1983. *Elastic Wave Propagation in Transversely Isotropic Media*, Springer.

Pica, A., 1997. Fast and accurate finite-difference solutions of the 3-D Eikonal equation parameterized in celerity, in *Expanded Abstracts*, pp. 1774–1777, Society of Exploration Geophysics.

Podvin, P. & Lecomte, I., 1991. Finite difference computation of traveltimes in very contrasted velocity model: a massively parallel approach and its associated tools, *Geophys. J. Int.*, **105**, 271–284.

Postma, G., 1955. Wave propagation in a stratified medium, *Geophysics*, **20**(4), 780–806.

Qian, J. & Symes, W.W., 2001. Paraxial Eikonal solvers for anisotropic quasi-P travel times, *J. Comput. Phys.*, **173**, 256–278.

Qian, J. & Symes, W., 2002a. An adaptive finite-difference method for traveltimes and amplitudes, *Geophysics*, **67**, 167–176.

Qian, J. & Symes, W.W., 2002b. Finite-difference quasi-P traveltimes for anisotropic media, *Geophysics*, **67**, 147–155.

Qian, J., Zhang, Y.-T. & Zhao, H.-K., 2007. A fast sweeping method for static convex Hamilton-Jacobi equations, *J. Sci. Comput.*, **31**, 237–271.

Rawlinson, N., Hauser, J. & Sambridge, M., 2007. Seismic ray tracing and wavefront tracking in laterally heterogeneous media, *Adv. Geophys.*, **49**, 203–267.

Rouy, E. & Tourin, A., 1992. A viscosity solution approach to shape-from-shading, *SIAM J. Numer. Anal.*, **29**(3), 867–884.

Schiesser, W.E., 1991. *The numerical Method of Lines: Integration of Partial Differential Equations*, Academic Press.

Shah, H., 2007. The 2007 bp anisotropic velocity-analysis benchmark, in *Expanded Abstracts*, EAGE workshop.

Slawinski, M., 2003. *Seismic Waves and Rays in Elastic Media*, Elsevier Science.

Taillandier, C., Noble, M., Chauris, H. & Calandra, H., 2009. First-arrival travel time tomography based on the adjoint state method, *Geophysics*, **74**(6), WCB1-WCB10.

Tavakoli F., B., Ribodetti, A., Virieux, J. & Operto, S., 2015. An iterative fac-tored Eikonal solver for TTI media, in *SEG Technical Program Expanded Abstracts 2015*, Vol. 687, pp. 3576–3581.

Thomsen, L.A., 1986. Weak elastic anisotropy, *Geophysics*, **51**, 1954–1966.

Tsai, Y.-H.R., Chen, L.-T., Osher, S. & Zhao, H.-K., 2003. Fast sweeping algorithms for a class of Hamilton-Jacobi equations, *SIAM J. Numer. Anal.*, **41**(2), 673–694.

Vidale, D., 1988. Finite-difference calculation of travel time, *Bull. seism. Soc. Am.*, **78**, 2062–2076.

Virieux, J. & Lambaré, G., 2007. Theory and observations - body waves: ray methods and finite frequency effects, in *Treatise of Geophysics, Volume 1: Seismology and Structure of the Earth*, pp. 127–155, eds Romanowicz, B. & Diewonski, A., Elsevier.

Waheed, U. & Alkhalifah, T., 2017. A fast sweeping algorithm for accurate solution of the tilted transversely isotropic Eikonal equation using factorization, *Geophysics*, **82**(6), WB1-WB8.

Waheed, U., Alkhalifah, T. & Wang, H., 2015a. Efficient traveltime solutions of the acoustic TI Eikonal equation, *J. Comput. Phys.*, **282** (Supplement C), 62–76.

Waheed, U.B., Yarman, C.E. & Flagg, G., 2015b. An iterative, fast-sweeping-based Eikonal solver for 3D tilted anisotropic media, *Geophysics*, **80**, C49-C58.

Zhang, L., Rector III, J.W. & Hoversten, G.M., 2005a. Eikonal solver in the celerity domain, *Geophys. J. Int.*, **162**, 1–8.

Zhang, Y., Zhao, H. & Chen, S., 2005b. Fixed-point iterative sweeping methods for static Hamilton-Jacobi equations, *Methods Appl. Anal.*, **13**(3), 299–320.

Zhang, Y.-T., Zhao, H.-K. & Qian, J., 2006. High order fast sweeping methods for static Hamilton-Jacobi equations, *J. Sci. Comput.*, **29**, 25–56.

Zhang, Y.-T., Chen, S., Li, F., Zhao, H. & Shu, C.-W., 2011. Uniformly accurate discontinuous Galerkin fast sweeping methods for Eikonal equations, *SIAM J. Sci. Comput.*, **33**(4), 1873–1896.

Zhao, H., 2005. A fast sweeping method for Eikonal equations, *Math. Comput.*, **74**, 603–627.

Zienkewicz, O. & Morgan, K., 1983. *Finite Elements and Approximation*, J. Wiley and Sons.
APPENDIX A: EIKONAL DERIVATION IN THE VTI CASE

For a VTI medium, let us find out a first expression of the Hamiltonian. Transverse isotropy involves five independent parameters in the elastic stiffness tensor, which writes

$$ C_{\text{VTI}} = \begin{bmatrix} C_{11} & C_{12} & C_{13} & 0 & 0 & 0 \\ C_{12} & C_{11} & C_{13} & 0 & 0 & 0 \\ C_{13} & C_{13} & C_{33} & 0 & 0 & 0 \\ 0 & 0 & 0 & C_{44} & 0 & 0 \\ 0 & 0 & 0 & 0 & C_{44} & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{C_{11}-C_{12}}{2} \end{bmatrix}, \qquad (A1) $$

where the elastic stiffness components $C_{ij}$ are expressed with the Voigt notation. Christoffel's dispersion relation for an arbitrary direction of propagation for the three wave modes yields

$$ \begin{aligned} & (C_{66}(n_x^2 + n_y^2) + C_{44}n_z^2 - \rho v^2) \\ & (-C_{13}^2(n_x^2 + n_y^2)n_z^2 - 2C_{13}C_{44}(n_x^2 + n_y^2)n_z^2 + C_{33}C_{44}n_z^4 - C_{44}(n_x^2 + n_y^2)\rho v^2 - C_{33}n_z^2\rho v^2 - C_{44}n_z^2\rho v^2 \\ & + C_{11}(n_x^2 + n_y^2)(C_{44}(n_x^2 + n_y^2) + C_{33}n_z^2 - \rho v^2) + \rho^2 v^4) = 0. \end{aligned} \qquad (A2) $$

In eq. (A2), $\mathbf{n} = (n_x, n_y, n_z)$ is the unit vector oriented towards the direction of propagation, and the z-axis is the rotation-symmetry axis in the VTI case (Slawinski 2003, p. 229). We now assume that the elastic parameters are invariant along the y direction and we consider a 2-D propagation inside the xz-plane, letting $n_y = 0$. Equation A2 exhibits two factors: the first one corresponding to SH propagation (shear wave propagating inside the xz-plane with a purely orthogonal displacement direction parallel to the y-axis), and the second one coupling both qP (quasi-pure compressional wave with a displacement direction contained in the xz-plane and quasi-parallel to the direction of propagation) and qSV (quasi-pure shear wave with a displacement direction contained in the xz-plane and quasi-orthogonal to the direction of propagation) modes. These two factors can be solved independently. We consider here the coupled P-SV mode which contains the compressional mode. The slowness vector $\mathbf{p} = (p_x, p_z)$ is defined by

$$ \begin{aligned} n_x^2 &= v^2 p_x, \\ n_z^2 &= v^2 p_z. \end{aligned} \qquad (A3) $$

We also define

$$ \begin{aligned} V_P &= \sqrt{\frac{C_{33}}{\rho}}, & V_S &= V_{SV} = \sqrt{\frac{C_{44}}{\rho}}, \\ \epsilon &= \frac{C_{11} - C_{33}}{2C_{33}}, & \delta &= \frac{(C_{13} + C_{44})^2 - (C_{33} - C_{44})^2}{2C_{33}(C_{33} - C_{44})}. \end{aligned} \qquad (A4) $$

The quantities $V_P$ and $V_S$ denote the P- and SV-wave velocities along the rotation-symmetry axis, while $\epsilon$ and $\delta$ are known as the Thomsen's parameters (Thomsen 1986). The P-SV Eikonal for 2-D VTI medium can be written as

$$ ap_x^4 + bp_z^4 + cp_x^2p_z^2 + dp_x^2 + ep_z^2 - 1 = 0, \qquad (A5) $$

where

$$ \begin{cases} a = -(1 + 2\epsilon)V_P^2V_S^2, \\ b = -V_P^2V_S^2, \\ c = -(1 + 2\epsilon)V_P^4 - V_S^4 + (V_P^2 - V_S^2)[V_P^2(1 + 2\delta) - V_S^2], \\ d = V_S^2 + (1 + 2\epsilon)V_P^2, \\ e = V_P^2 + V_S^2. \end{cases} \qquad (A6) $$

APPENDIX B: EQUIVALENCE OF VTI HAMILTONIANS

The two VTI Hamiltonians (12) and (41) are expected to yield the same steady state. This is demonstrated in this appendix. At steady state, with the notations $X = u_{,x}$, $z = u_{,z}$, inserting (12) inside (2) writes

$$ dX^2 + eZ^2 + cX^2Z^2 - 1 = 0. \qquad (B1) $$
Adding a term on both sides and rearranging (B1) yields

$$1 + \frac{1}{4} (dX^2 + eZ^2)^2 - dX^2 - eZ^2 = \frac{1}{4} (dX^2 + eZ^2)^2 + cX^2 Z^2, \qquad (B2)$$

$$\left[ 1 - \frac{1}{2} (dX^2 + eZ^2) \right]^2 = \frac{d^2}{4} X^4 + \frac{e^2}{4} Z^4 + \left( \frac{de}{2} + c \right) X^2 Z^2. \qquad (B3)$$

From (11) we get

$$\frac{de}{2} + c = \frac{V_p^4}{2} (1 - 2\epsilon + 4\delta), \qquad (B4)$$

which is positive in practice for all realistic applications in geophysics. Therefore the square root of (B3) can be taken without loss of generality, yielding

$$\sqrt{\frac{1}{4}(dX^2 + eZ^2)^2 + cX^2 Z^2} + \frac{1}{2}(dX^2 + eZ^2) = 1. \qquad (B5)$$

Since $d > 0$ and $e > 0$, the square root can be taken again, and dividing by $\sqrt{d}$, the VTI Hamiltonian (41) is obtained:

$$\mathcal{H}_{\text{VTI}} = \frac{1}{\sqrt{d}} \left( \sqrt{\frac{1}{4} (du_{,x}^2 + eu_{,z}^2)^2 + cu_{,x}^2 u_{,z}^2} + \frac{1}{2} (du_{,x}^2 + eu_{,z}^2) - 1 \right). \qquad (B6)$$

## APPENDIX C: LIPSCHITZ CONTINUITY OF THE VTI HAMILTONIAN

With the notations $X = u_{,x}$, $z = u_{,z}$, for $(X, z) \neq (0, 0)$, the derivatives of VTI Hamiltonian (41) write

$$\mathcal{H}_1(X, Z) = \frac{dX + \frac{2cXZ^2 + dX(dX^2 + eZ^2)}{\sqrt{(dX^2 + eZ^2)^2 + 4cX^2Z^2}}}{\sqrt{2\sqrt{(dX^2 + eZ^2)^2 + 4cX^2Z^2} + 2(dX^2 + eZ^2)}}, \qquad (C1)$$

$$\mathcal{H}_2(X, Z) = \frac{eZ + \frac{2cX^2Z + eZ(dX^2 + eZ^2)}{\sqrt{(dX^2 + eZ^2)^2 + 4cX^2Z^2}}}{\sqrt{2\sqrt{(dX^2 + eZ^2)^2 + 4cX^2Z^2} + 2(dX^2 + eZ^2)}}. \qquad (C2)$$

This yields

$$\|\mathcal{H}\|^2 = \frac{1}{2d} \frac{\left( dX + \frac{2cX^2 + dX(dX^2 + eZ^2)}{\sqrt{(dX^2 + eZ^2)^2 + 4cX^2Z^2}} \right)^2 + \left( eZ + \frac{2cX^2Z + eZ(dX^2 + eZ^2)}{\sqrt{(dX^2 + eZ^2)^2 + 4cX^2Z^2}} \right)^2}{\sqrt{(dX^2 + eZ^2)^2 + 4cX^2Z^2 + dX^2 + eZ^2}}. \qquad (C3)$$

In the elliptical case, $c = 0$ thus (C3) simplifies:

$$\|\mathcal{H}\|_{\text{ELL}}^2 = \frac{d^2 X^2 + e^2 Z^2}{d^2 X^2 + d e Z^2}. \qquad (C4)$$

Since $\epsilon \ge 0$, then $d > e$ and we have

$$\max \|\mathcal{H}\|_{\text{ELL}} = 1, \qquad (C5)$$

which is obtained for $Z = 0$.

In the general anelliptical case, we use the polar coordinate system

$$\begin{aligned} \sqrt{d}X &= r \cos \gamma, \\ \sqrt{e}Z &= r \sin \gamma. \end{aligned} \qquad (C6)$$

Changing variables in (C3) yields

$$\|\mathcal{H}\|^2 = \frac{A + B + C}{2d \left( 1 + \sqrt{1 + \frac{4c}{de} \cos^2 \gamma \sin^2 \gamma} \right)}, \qquad (C7)$$
with

$$A = d \cos^2 \gamma + e \sin^2 \gamma,$$

$$B = \frac{d \cos^2 \gamma + e \sin^2 \gamma + 4c \cos^2 \gamma \sin^2 \gamma \left(\frac{1}{d} + \frac{1}{e}\right) + \frac{4c^2}{de} \cos^2 \gamma \sin^2 \gamma \left(\frac{\cos^2 \gamma}{d} + \frac{\sin^2 \gamma}{e}\right)}{1 + \frac{4c}{de} \cos^2 \gamma \sin^2 \gamma},$$

$$C = \frac{2(d \cos^2 \gamma + e \sin^2 \gamma) + 4c \cos^2 \gamma \sin^2 \gamma \left(\frac{1}{d} + \frac{1}{e}\right)}{\sqrt{1 + \frac{4c}{de} \cos^2 \gamma \sin^2 \gamma}} \quad (C8)$$

Expression (C7) holds for $r \neq 0$. The variable $r$ simplifies so that $||\mathcal{H}||^2$ is a function of one variable $\gamma$. This function is $\pi$-periodic. Getting back to the Thomsen's parameters $\epsilon$ and $\delta$ using (11), some calculus that we do not reproduce here gives

$$\|\mathcal{H}\|^2 = \frac{D + E}{(2 + 4\epsilon) \left(1 - \frac{2(\epsilon - \delta)}{1+2\epsilon} \sin^2 2\gamma\right) \left(1 + \sqrt{1 - \frac{2(\epsilon - \delta)}{1+2\epsilon} \sin^2 2\gamma}\right)} + 1, \quad (C9)$$

with

$$D = \frac{4\epsilon(\epsilon - \delta)(1 + 4\epsilon - 2\delta)}{(1 + 2\epsilon)^2} \sin^2 \gamma \sin^2 2\gamma + \frac{4(\epsilon - \delta)(\epsilon(1 + 2\epsilon) + \epsilon - \delta)}{(1 + 2\epsilon)^2} \sin^2 2\gamma - (6\epsilon - 2\delta) \sin^2 \gamma,$$

$$E = 4\epsilon \sqrt{1 - \frac{2(\epsilon - \delta)}{1 + 2\epsilon} \sin^2 2\gamma} \left( \frac{\epsilon - \delta}{1 + 2\epsilon} \sin^2 2\gamma - \sin^2 \gamma \right). \quad (C10)$$

Since $\epsilon \ge \delta \ge 0$ and assuming (B4) is positive, we have

$$\frac{\epsilon - \delta}{1 + 2\epsilon} \sin^2 2\gamma - \sin^2 \gamma = \sin^2 \gamma \left( \frac{4(\epsilon - \delta)}{1 + 2\epsilon} (1 - \sin^2 \gamma) - 1 \right) \\ \le \frac{4(\epsilon - \delta)}{1 + 2\epsilon} - 1 \le 0. \quad (C11)$$

Therefore $E \le 0$.

To state $D \le 0$, we define $y = \sin^2 2\gamma$. Thus

$$D = y(4A'y(1-y) + 4B'(1-y) - C') \\ = y(-4A'y^2 + (4A' - 4B')y + 4B' - C'), \quad (C12)$$

with the notations

$$A' = \frac{4\epsilon(\epsilon - \delta)(1 + 4\epsilon - 2\delta)}{(1 + 2\epsilon)^2},$$

$$B' = \frac{4(\epsilon - \delta)(\epsilon(1 + 2\epsilon) + \epsilon - \delta)}{(1 + 2\epsilon)^2},$$

$$C' = 6\epsilon - 2\delta. \quad (C13)$$

Since $y \ge 0$, $D$ is of the sign of the second-order polynomial in $y$ in (C12). Its discriminant writes

$$\Delta = 16((A' - B')^2 + A'(4B' - C')) = 16((A' + B')^2 - A'C'). \quad (C14)$$

Using (C13), we obtain

$$(A' + B')^2 - A'C' = \frac{4(\epsilon - \delta)(3\epsilon - \delta)}{(1 + 2\epsilon)^2} (4(\epsilon - \delta)(3\epsilon - \delta) - 2\epsilon(1 + 4\epsilon - 2\delta)) \\ = \frac{4(\epsilon - \delta)(3\epsilon - \delta)}{(1 + 2\epsilon)^2} (2\epsilon(2\epsilon - 4\delta - 1) - 4\delta(\epsilon - \delta)) \\ \le 0, \quad (C15)$$

since (B4) is positive. It follows that $D \le 0$.