Wave-equation imaging is computed almost exclusively on Cartesian meshes both for computational simplicity and because imaged subsurface volumes are usually rectangular parallelepipeds. However, situations exist where performing wave-equation imaging on more generalized coordinate system meshes is warranted. For example, extrapolating wavefields directly from an undulating topographic surface Shragge and Sava (2004a), orienting the extrapolation axis of lower-order operators in the direction of wavefield propagation to improve higher angle accuracy Sava and Fomel (2004), or imaging with overturning waves not currently modeled by one-way extrapolation operators Shan and Biondi (2004). Employing non-Cartesian meshes, though, necessitates resolving these three issues: which coordinate system should be chosen? how is the preferred coordinate system generated? and what are the appropriate governing wavefield extrapolation equations for this coordinate system choice?
Non-Cartesian wavefield extrapolation theory has advanced in recent years in the context of both global and exploration seismology. All of these methods locally transform the coordinate system and the corresponding governing propagation equations to a more appropriate reference frame linked to an underlying Cartesian grid through one-to-one mappings. For example, Bostock et al. (1993) formulate an orthogonal, 2-D, plane-wave-centric coordinate system appropriate for propagating overturning teleseismic wavefields. Etgen (2002) and Shan and Biondi (2004) discuss the use of tilted coordinate systems that enable propagation of overturning waves with one-way extrapolation operators. Brandsberg-Dahl and Etgen (2003) and Sava and Fomel (2004) present generalized curvilinear coordinate and Riemannian metric formulations of the governing wavefield propagation equations, respectively.
Emerging with these developments are two classes of coordinate systems: those conformal with the propagating wavefield direction, and those conformal with pre-existing topology. Examples of the former include coordinate systems constructed by plane-wave decomposition Bostock et al. (1993), and those generated using Eikonal equation ray-tracing Sava and Fomel (2004), whereas examples of the latter include a conformal mapping transform to incorporate surface topography Shragge and Sava (2004a). These methods, though, neither are universally applicable, nor always practical for 3D wavefield extrapolation. For example, ray-based methods can generate semi-orthogonal, triplicating coordinate systems that lead to numerical instability when propagating wavefields through coordinate system caustics. Another example is that the complex mathematics of conformal mapping surface topography does not extend easily to 3D. Hence, formulating a general approach for constructing 3D orthogonal coordinate systems appropriate for non-Cartesian wavefield extrapolation remains an open research topic.
This paper examines how potential function (PF) solutions of Laplace's equation, coupled with phase-ray tracing Shragge and Sava (2004b), can be used to generate orthogonal coordinate system meshes. Starting with appropriate PF boundary conditions, I generate a solution to Laplace's equation for any simply connected domain using conjugate gradient solvers Claerbout (1999). Importantly, the PF equipotentials define surfaces equivalent to wavefield extrapolation steps. The PF gradient field, by definition orthogonal to equipotentials, similarly forms geometric rays originating from the acquisition surface that collective form a ray coordinate system. The most obvious example is the Cartesian coordinate system, where each depth step surface represents an equipotential surface, and each vertical line projecting downward defines a coordinate system ray. Mathematical properties of a PF also guarantee that the 3D orthogonality of equipotential coordinate systems. Hence, this approach constructs computational meshes conformal to 3D geometric boundaries, including undulating 2D topographic surfaces or deviated wells found in vertical seismic profiling.
I begin this paper by reviewing some important properties of PF solutions of Laplace's equation, and then outline an approach for defining a Laplace's equation appropriate for incorporating generalized topology of a solution domain. Subsequently, I show how modified phase-ray tracing generates a suite of geometric rays that collectively form a coordinate system mesh. I then demonstrate the utility of the approach by developing a 2D coordinate system for the Canadian Rocky Mountain Foothills model discussed in Shragge and Sava (2004a) and a 3D coordinate system conformal to the 2D topography of the San Francisco Bay area.