Introduction

Traveltime calculations play an important role in many methods of seismic data processing. For example, the Kirchhoff methods of migrating and modeling seismic data require the calculation of Green's functions that describe the traveltimes between survey points on the surface and depth points in the velocity model. These traveltimes are often calculated by ray tracing (Cervený, 1987), and, because depth and surface points are normally distributed on a regular grid, the traveltimes along the rays are then interpolated onto the grid. For complicated velocity models, rays may cross each other or not penetrate shadow zones; the interpolation is thus cumbersome and computationally expensive.

Recently, several methods have been introduced to calculate traveltimes directly on a regular grid. In calculating the depth gradient of traveltimes, Reshef and Kosloff (1986) use a finite-difference approximation to the eikonal equation, which they then integrate with a Runge-Kutta method. Vidale (1988) also approximates the eikonal equation by finite differences, but he solves directly for traveltime using a planar or circular wavefront extrapolation.

Although these methods are elegant alternatives to interpolating traveltimes between rays, and far more efficient computationally, they have some drawbacks. Reshef and Kosloff's method is unstable for models with large velocity contrasts; in regions with velocity discontinuities, traveltimes have to be calculated by alternative methods, such as ray tracing. The principal drawback in Vidale's method is that the points on the computational front must be ordered, with calculations starting at the so-called ``relative minimum points''. This ordering step complicates the algorithm, and, because it is not vectorizable, makes the method computationally expensive.

Here we suggest an alternative scheme that overcomes the problems described above. We observe that the components of the traveltime gradient satisfy so-called hyperbolic conservation laws. Conservation laws are extremely important in fluid mechanics -- the Euler equations and the gas dynamics equations are examples -- and a great deal of effort has been devoted to their numerical solution. Accordingly we have used an ``off-the-shelf,'' upwind finite-difference scheme as the basis for our method. The resulting scheme is stable, involves none of the sorting implicit in Vidale's approach, and vectorizes fully.

Just like Vidale's method (which is in fact also an upwind finite-difference scheme), our scheme correctly estimates the traveltimes of first arrivals in caustics and of diffractions in shadow zones. A disadvantage is that neither method calculates multi-valued traveltime curves. This restriction is intrinsic to the approaches pursued here and by Vidale. Although this limitation can be serious for some applications, many Kirchhoff algorithms use single-valued traveltime functions anyway, and the method is readily applicable to these algorithms. Also, many tomography codes assume that picked times are first arrivals. To the extent that this assumption is correct, finite-difference traveltime calculation might serve as the basis for very efficient tomography codes.

In this paper we first describe the notion of viscosity solutions, which underlies the finite-difference calculations, and explain the need for upwind differencing, which plays a key role in both our scheme and Vidale's. Next we state the eikonal equation and describe the related conservation laws and upwind finite-difference schemes. Then we demonstrate the calculations with some examples, where we compare the results both with analytical solutions and with those obtained by Vidale's method. Finally, we discuss some of the limitations and applications of the method.