%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ %+ + %+ NOTICE: This article is also available formatted with illustrations. + %+ + %+ SEE: http://sepwww.stanford.edu/oldreports/sep61/ + %+ + %+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ \def\be{\begin{equation}} \def\ee{\end{equation}} \def\bea{\begin{eqnarray}} \def\eea{\end{eqnarray}} \def\dsp{\displaystyle} \def\t{\tau} % traveltime \def\s{s} % slowness \def\u{u} % dt/dx \def\ubar{{\overline{\u}}} % u bar \def\f{F} % f \def\fplus{\f_+} % f+ \def\fmin{\f_-} % f- % derivatives \def\dtdx{\t_x} \def\dtdz{\t_z} \def\dtdzdx{\t_{zx}} \def\dudx{\u_x} \def\dudz{\u_z} \def\dsdx{\s_x} \def\dfdx{\f_x} \def\sqrtsu{\sqrt{\s^2 - \u^2}} % finite diff. \def\dplus{\Delta_+} \def\dmin{\Delta_-} \def\dx{\Delta x} \def\dz{\Delta z} \def\dzdx{{{\dz}\over{\dx}}} \def\unj{\u^n_j} \def\unpj{\u^{n+1}_j} \def\figs{/scr/jos/Figa} \title{Finite-difference calculation of traveltimes} \author{Jos van Trier} \lefthead{Van Trier} \righthead{Finite-difference traveltimes} \footer{SEP--61} \ABS{ The Eikonal equation can be transformed into a flux-conservative equation, which can then be solved with standard finite-difference techniques. A first-order upwind finite-difference scheme proves accurate enough for seismic applications. The method is limited to rays traveling in one direction, and only calculates single-valued traveltime functions. The algorithm efficiently calculates traveltimes on a regular grid, and is useful in Kirchhoff migration and modeling. } \INT Traveltime calculations play an important role in many methods in seismic data processing. For example, in the Kirchhoff implementation of migration and modeling of seismic data Green's functions need to be calculated that describe the traveltime between survey points on the surface and depth points in the velocity model. These traveltimes are often calculated by ray tracing, and, since 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 or not penetrate shadow zones, making the interpolation cumbersome and computationally expensive. Vidale~(1988) recently described a method that calculates traveltimes directly on a regular grid. His approach is as follows: given the initial traveltime value on a certain grid point, the traveltimes for neighbouring points are extrapolated using a planar or circular wavefront approximation. Although Vidale's method represents an elegant alternative to interpolating traveltimes between rays, it has some drawbacks. First, to ensure stability, extrapolation can only be initiated from certain points on the grid (the so-called ``relative minimum points''). Second, to achieve accuracy, at each grid point the curvature of the wavefront has to be determined, and, if necessary, an expensive circular wavefront extrapolation step has to be carried out. These drawbacks not only increase the computational cost, but also make it impossible to vectorize the computations. %Vectorization is normally an important advantage of finite-difference methods. W. Symes of Rice University introduced me to Vidale's method during a recent visit to Europe, and suggested an alternative scheme that overcomes the problems described above. Symes suggests writing the Eikonal equation as a flux equation, and then solving it with a standard upwind finite-difference scheme. This leads to a method that is computationally more efficient than Vidale's scheme. Just like Vidale's method, the finite-difference calculations correctly estimate traveltimes of first arrivals in caustics and of diffractions in shadow zones. One disadvantage however, is that traveltimes can only be calculated for rays traveling in one direction, and that the method does not calculate multi-valued traveltime curves. Although this may be a serious limitation for some applications, many Kirchhoff algorithms use single-valued traveltime functions anyway, and the method is readily applicable to those algorithms. In this paper I first describe the method and the finite-difference implementation. Then I demonstrate the calculations with some examples, where I compare the results both with analytical solutions and with Vidale's method. Finally, I discuss some of the limitations and applications of the method. \mhead{EIKONAL EQUATION} The Eikonal equation in two dimensions is as follows: \be \Bigl( \dtdx \Bigr) ^2 \ +\ \Bigl( \dtdz \Bigr) ^2 \ =\ \s^2(x,z), \label{eik} \ee where $s(x,z)$ is the 2-dimensional slowness model and $\t = \t(x,z)$ is the traveltime field. Subscripts $x$ and $z$ denote partial derivatives with respect to $x$ and $z$, respectively. \shead{Flux-conservative equation} Similar to the conversion of the full wave equation into a one-way wave equation, we can rewrite the eikonal equation as a one-way equation that can be solved with standard techniques. First, use $\u = \dtdx$ to rewrite equation~(\ref{eik}) as \be \dtdz \ =\ \sqrtsu, \label{eiku} \ee By choosing a positive sign in front of the square root, and by using $\u=\dtdx$ instead of $\u=\dtdz$ as the substitution variable, we limit ourselves to downward-traveling rays. I will come back to these choices in one of the later sections. Second, take the derivative of this equation with respect to $x$: \be \dtdzdx\ =\ \dudz \ =\ -\dfdx(\u). \label{flux} \ee where the function $\f(\u)$ is defined as \be \f(\u) \ =\ -\sqrtsu. \ee $\f(\u)$ is called the conserved flux; if $\f(\u) = 0$, the rays do not ``flow'' downward anymore, but travel horizontally. The (one-way) eikonal equation now has the form of a flux-conservative equation This is a well-known equation in computational fluid dynamics, and it can be solved in many ways (for example, see Roache,~1976). The next section discusses one particular solution. \mhead{FINITE-DIFFERENCE SCHEME} Although sophisticated higher-order finite-difference schemes exist for the solution of the flux-conservative equation (Centrella and Wilson,~1984; Hawley et al.,~1984), a first-order upwind finite-difference scheme is generally accurate enough for geophysical applications. I implemented a first-order method described by Engquist and Osher~(1980). Depending on the flow direction, upwind finite-difference schemes use a backward or forward finite-difference approximation to the derivative operator. The idea is that, although centered approximations are mathematically more accurate, the underlying physics may dictate otherwise. In other words, it is better to make a big error that is not transported, than to make a small one that blows up. The first-order forward and backward finite differences in space are, respectively: \be \dplus \u\ =\ \u_{j+1} - \u_j,\ \ \ \ \ \ \ \ \dmin \u =\ \u_j - \u_{j-1}, \ee where $\u_j$ is the discrete representation of $\u(x)$ on the grid $x_j = j\dx$. Using the above finite-difference scheme, equation~(\ref{flux}) is approximated by: \begin{equation} {{\unpj-\unj} \over {\dz}} \ =\ \left\{ \begin{array}{ll} \ -\dsp{{\dplus\ \f(\unj)} \over {\dx}}\ & {\rm if}\ \unj \leq 0; \\ \ \ \\ \ -\dsp{{\dmin \ \f(\unj)} \over {\dx}}\ & {\rm if}\ \unj > 0, \end{array} \right. \label{stupid} \end{equation} with $n$ the discrete-depth index of the grid: $z^n = n\dz$. However, the resulting scheme is unstable when $\u$ is discontinuous (when the field is going through a ``shock''). The discontinuity occurs at the source: going from left to right through the source, $\u$ changes sign from $-\s$ to $+\s$, i.e., rays change direction from left to right. Engquist and Osher solve the unstability problem in the shock region with the following scheme: \be \unpj = \unj\ -\ \dzdx\ \Bigl(\dplus \fmin(\unj)\ +\ \dmin \fplus(\unj) \Bigr), \ee where \be \fplus(u)\ =\ \f({\rm max}(\u,\ubar)),\ \ \ \ \ \ \ \fmin(u) \ =\ \f({\rm min}(\u,\ubar)), \label{eng} \ee and $\ubar$ is the stagnation point, i.e., $\f'(\ubar) = 0$. $\ubar = 0$ for the function $\f$ under consideration here. Engquist and Osher's scheme reduces to the standard first-order scheme~(\ref{stupid}) away from the source, and is nonlinearly stable near the source. \shead{Initial and boundary conditions} Equation~(\ref{eng}) is used to extrapolate $\u$ downward in $z$. Initial values for $\u$ have to be specified at the surface. In many cases the velocity model is constant at the surface, and the initial conditions can be calculated analytically for a given source position. For models with non-constant surface velocities, traveltimes can be computed along the surface with a simple 1-D finite-difference approximation. The derivative of the traveltime function with respect to $x$ then gives $\u$. When one assumes outgoing rays at the boundaries, one-sided finite differences can be used at the left and right side of the model. These boundary conditions are of the same order as the scheme inside the model, and generally do not cause any problems. After $\u$ has been computed, traveltimes are calculated by integrating $\u$ over $x$ using Simpson's rule. \plot{trmap}{3.75in}{\figs/fdtime}{Finite-difference traveltime function for a constant velocity model and a source on the surface in the middle of the model.} \plot{diffmap}{3.75in}{\figs/difftime}{Difference between the traveltime map of Figure~\protect\ref{trmap} and the analytical solution. Higher intensities in the plot represent larger errors. A quantitative measure of the errors is given in Figure~\protect\ref{bottom}.} \mhead{EXAMPLES} In the first example I calculate traveltimes for a constant-velocity medium with a velocity of 2.5 km/s. The grid is evenly sampled in depth and laterally, with a sample interval of 10 m. Figure~\ref{trmap} shows the traveltime function for a source on the surface at the middle of the model, and Figure~\ref{diffmap} displays the difference with the analytical solution. \plot{vid}{3.75in}{\figs/diffvidtime}{Difference between a traveltime function calculated with Vidale's plane wave extrapolation method and the analytical solution. The intensity scale is the same as in Figure~\protect\ref{diffmap}.} Errors accumulate as depth increases, and are largest on the outer sides of the model, where rays travel at angles closer to the horizontal. Since the finite-difference scheme is designed for downward traveling rays, this behavior is expected. The maximum error is about .6 ms, which is one order of magnitude smaller than the standard time sampling interval of 4 ms (see Figure~\ref{bottom}). Finite-difference traveltime calculations of one shot on a $100\times 100$ grid take about .1 s in CPU time on the Convex C-1. \plot{bottom}{3.in}{\figs/bottom}{Errors in the finite-difference traveltime calculations at the bottom of model. The solid line represents the error curve for the method described here, the dashed line denotes the error in Vidale's scheme.} Figure~\ref{vid} shows the difference between Vidale's scheme and the analytical solution. I have only implemented the plane wave extrapolation method, and errors can probably be reduced if a combination of plane and circular wave extrapolation is used. The errors are largest away from the vertical, diagonal and horizontal direction, where the plane wave approximation breaks down. The errors at the bottom of the model are of the same order of magnitude as the errors in the method described here (see again Figure~\ref{bottom}). However, Vidale's scheme does not vectorize, and on the Convex C-1 the finite-difference calculations are about 5-10 times faster than his method. \plot{strmod}{3.75in}{\figs/structmod}{Wedge model. The grid spacing is $10\times 10$ m. High intensities denote low velocities.} The next example illustrates the calculations for a more complicated model. The model is shown in Figure~\ref{strmod}; it consists of 3 layers and a wedge intrusion. The velocity in the top layer is 2 km/s, the middle layer has a velocity of 1.75 km/s, and the bottom velocity is 2.5 km/s. The velocity in the wedge that intrudes the middle layer from the right is 2.75 km/s. Figure~\ref{rays} shows the result of tracing rays from the surface downwards; because of the velocity contrasts, rays cross and shadow zones are apparent. \plot{rays}{3.75in}{\figs/raystr}{Rays traced through the model of Figure~\protect\ref{strmod}.} \plot{stmmap}{3.75in}{\figs/fdtimestr}{Finite-difference traveltimes calculated for the model of Figure~\protect\ref{strmod}.} As is obvious from Figure~\ref{rays}, interpolating traveltimes from the rays onto the grid is not easy for this model. However, the finite-difference calculation correctly fills in the problem areas as can be seen in Figure~\ref{stmmap}: the contour lines in the plot reveal the correct curvature of the wave fronts in the high- and low-velocity regions. Figure~\ref{sbottom} provides a more quantitative check; it compares interpolated traveltimes of the rays at the bottom of the model with the finite-difference traveltimes. Since the rays are calculated for a spline model that is fitted to the grid model (Van Trier, 1988), some discrepancies may be expected, especially if one realizes that a slight change in velocity contrast can drastically change the direction of the ray. Except for the region near the shadow zone, the traveltime curves match reasonably well. \plot{sbottom}{3.5in}{\figs/sbottom}{Comparison of finite-difference traveltimes with traveltimes interpolated from rays at the bottom of the model (fat line). The rays do not reach the rightmost part of the bottom, so interpolation can only be done up to $x=2400$ m.} \mhead{LIMITATIONS} The flux-conservative form of the Eikonal equation is essentially a one-way equation. This means that rays can only be propagated in one direction, and that only single-valued traveltime curves are computed. So far, I have been discussing downward traveling rays, but if $\u=\dtdz$ instead of $\u=\dtdx$ is used in equation~(\ref{eiku}), one can calculate traveltimes of rays that travel sideways. Initial conditions then have to be specified along a vertical line. In the downward-propagation finite-difference scheme, the numerical errors get large if \makebox{$\f \dz \ll \dx$}, which occurs when rays travel almost horizontally. This numerical noise can be reduced by decreasing the stepsize in $z$ when rays have large propagation angles at the current depth level, and increasing it when the propagation angles get smaller. If greater accuracy is still required, an other alternative is to use higher-order finite difference schemes. Since the method does not calculate ray paths, it is only of limited use in tomographic methods. It is possible, however, to calculate ray directions on the grid from the traveltime gradient. One component of the slowness vector, $\dtdx$, is already available as $\u$, the other one, $\dtdz$, is calculated using $\dtdz = \sqrtsu$. \mhead{APPLICATION IN KIRCHHOFF METHODS} My main use for the method has been in Kirchhoff migration and modeling. Although ray tracing in itself may be faster than finite-difference calculations, the interpolation of traveltimes along rays onto the grid can be quite cumbersome and expensive: special care is needed in handling crossing rays and shadow zones. The finite-difference method does not handle multi-valued traveltime functions, however, and ray tracing is still necessary if one wants to correctly include triplications in the migration. The Green's functions that are used in Kirchhoff methods describe traveltimes between depth points in the model and survey points on the surface. One can precompute these traveltime maps with the described method for densely spaced surface points, store them on disk (if enough disk space is available), and then read them when needed in the Kirchhoff summation. If a map for a certain surface point is not on disk, a simple linear interpolation of traveltimes in two neighboring maps is accurate enough for most purposes. The amplitudes terms in the Green's functions can be approximately estimated from the traveltime function: the geometrical spreading and obliquity terms are approximate functions of traveltime or its gradient. Ray propagation angles are readily available from the gradient vector of the traveltime field (see previous section), and can be used accordingly in amplitude calculations. \mhead{CONCLUSIONS} I have presented a method that greatly simplifies traveltime calculations on a regular grid. The algorithm has the advantage over Vidale's method that all points can be treated equally, i.e., no minimum traveltime points have to be found at each extrapolation step, and no choices based on local wavefront curvature have to be made. %advantage of Vidale: all directions, higher accuracy than first-order. The computations are vectorizable, and allow efficient calculation of Greens's functions in Kirchhoff methods. \ACK I would like to thank Bill Symes of Rice University for bringing this topic to my attention, and for the interesting discussions we had on velocity analysis. This paper is mostly due to his suggestions. \REF \reference{Centrella, J., and Wilson, J., 1984, Astrophysical Journal Supplement, {\bf 54}, 229-249.} \reference{Engquist, B., and Osher, S., 1980, Stable and entropy satisfying approximations for transonic flow calculations: Math. Comp., {\bf 34}, no. 149, 45-75.} \reference{Hawley, J., Smarr, L., and Wilson, J., 1984, Astrophysical Journal Supplement, {\bf 55}, 211-246.} \reference{Roache, P., 1976, Computational fluid dynamics: Hermosa, Albuquerque, New Mexico.} \reference{Van Trier, J., 1988, Vectorized initial-value ray tracing for arbitrary velocity models: SEP--{\bf 57}, 279-288.} \reference{Vidale, J., 1988, Finite-difference calculation of travel times: Bull., Seis. Soc. Am., {\bf 78}, no. 6, 2062-2076.}