Continuous theory

The eikonal equation, governing the traveltimes from a fixed source in an isotropic heterogeneous medium, has the form  

$$\left(\frac{\partial \tau}{\partial x}\right)^2 + \left(\f... ...{\partial \tau}{\partial z}\right)^2 = \frac{1}{v^2(x,y,z)}\;.$$ (1)

Here x, y, and z are spatial coordinates, $\tau$ is the traveltime (eikonal), and v stands for the velocity field. Constant-traveltime contours in the traveltime field $\tau (x,y,z)$,constrained by equation (1) and appropriate boundary conditions, correspond to wavefronts of the propagating wave. Additionally, each point on a wavefront can be parameterized by an arbitrarily chosen ray parameter $\gamma$. In three dimensions, $\gamma$includes a pair of independent parameters. For brevity, from now on we will restrict the analysis to two dimensions. One can easily generalize it to the 3-D case by considering $\gamma$ and x as vector quantities. Thus, we will refer to the following two-dimensional form of equation (1):  

$$\left(\frac{\partial \tau}{\partial x}\right)^2 + \left(\frac{\partial \tau}{\partial z}\right)^2 = \frac{1}{v^2(x,z)}\;.$$ (2)

For a point source, $\gamma$ can be chosen as the initial ray angle at the source. Zhang (1993) shows that $\gamma$ as a function of spatial coordinates satisfies the simple partial differential equation  

$$\frac{\partial \tau}{\partial x}\,\frac{\partial \gamma}{\p... ...al \tau}{\partial z}\,\frac{\partial \gamma}{\partial z} = 0\;.$$ (3)

Equation (3) merely expresses the fact that in an isotropic medium, rays are locally orthogonal to wavefronts. The field $\gamma (x,z)$ has not only theoretical interest as it provides one of the possible ways for evaluating propagation amplitudes. In particular, the geometrical spreading factor J (x,z) is connected to $\gamma$ by the simple relationship Zhang (1993)  

$$\left(\frac{\partial \gamma}{\partial x}\right)^2 + \left(... ...{\partial \gamma}{\partial z}\right)^2 = \frac{1}{J^2(x,z)}\;.$$ (4)

It is important to note that for complex velocity fields, both $\tau$and $\gamma$ as functions of x and z become multi-valued. In this case, the multi-valued character of the ray parameter $\gamma$corresponds to the situation, where more than one ray from the source passes through a particular point $\{x,z\}$ in the subsurface. This situation presents a very difficult problem when equations (2) and (3) are solved numerically. Typically, only the first-arrival branch of the traveltime is picked in the numerical calculation. The ray tracing method is free from that limitation because it operates in the ray coordinate system. Ray tracing computes the traveltime $\tau$ and the corresponding ray positions x and z for a fixed ray parameter $\gamma$.

Since $x (\tau,\gamma)$ and $z (\tau, \gamma)$ are uniquely defined for arbitrarily complex velocity fields, we can now make an important mathematical transformation. Considering equations (2) and (3) as a system and applying the general rules of calculus, we can transform this system by substituting the inverse functions $x (\tau,\gamma)$ and $z (\tau, \gamma)$ for the original fields $\tau (x,z)$ and $\gamma (x,z)$. The resultant expressions take the form  

$$\left(\frac{\partial x}{\partial \tau}\right)^2 + \left(\f... ...au}\right)^2 = v^2 \left(x(\tau,\gamma),z(\tau,\gamma)\right)$$ (5)

and  

$$\frac{\partial x}{\partial \tau}\,\frac{\partial x}{\partia... ...al z}{\partial \tau}\,\frac{\partial z}{\partial \gamma} = 0\;.$$ (6)

Comparing equations (5) and (6) with the original system (2-3) shows that equations (5) and (6) again represent the dependence of ray coordinates and Cartesian coordinates in the form of partial differential equations. However, the solutions of system (5-6) are better behaved and have a unique value for every $\tau$ and $\gamma$. These values can be computed with the conventional ray tracing. However, the ray-tracing approach is based on a system of ordinary differential equations, which represents a different mathematical model.

We use equations (5) and (6) as the basis of our wavefront tracing algorithm. The next section discusses the discretization of the differential equations and the physical interpretation we have given to the scheme.