Introduction

Geophysical methodology often requires computing numerical solutions of partial differential equations (PDEs). In many cases, computational efficiency and/or accuracy can be enhanced by posing the physical equations on generalized coordinate meshes rather than Cartesian grids. By following this approach, though, mesh generation becomes a necessary solution process step. An important consideration is ensuring that the developed mesh has no attributes that would cause numerical instability (e.g. grid clustering). This is not a straightforward task because mesh generation is controlled by non-linear geometric coupling. Hence, localized mesh regularization can be difficult to implement and developing new gridding techniques with built-in regularization remains an open and important topic.

Meshing techniques are less advanced in the geophysical community relative to other scientific disciplines (e.g. computer graphics and fluid flow). However, a number of different grid generation approaches have been reported. Alkhalifah (2003) performs anisotropic velocity analysis on structured, non-orthogonal $\tau$-coordinate meshes based on two-way travel times. Sava and Fomel (2005) developed Riemannian wavefield extrapolation that generalizes one-way wave-propagation to structured semi-orthogonal ray-based meshes. Shragge and Sava (2005) examine wave-equation migration from topography in acquisition coordinates using structured locally orthogonal meshes. Rüger and Hale (2006) use a non-structured gridding technique based on tesselation to break subsurface velocity models into logical units.

An important set of structured meshing techniques used in other scientific fields are based on differential methods Liseikin (2004). These approaches pose mesh generation as solving an elliptic boundary value problem (BVP) within a framework of differential geometry. Numerical solution of these generalized Laplacian systems is facilitated by recasting the elliptic problem into a set of parabolic equations for which well-developed and efficient solution techniques exist. The steady-state solutions of the parabolic equations are structured non-orthogonal meshes. Importantly, the gridding equations can be manipulated to ensure that meshes exhibit appropriate attributes for numerical solution of PDEs: i) piece-wise smoothness of interior cells; ii) non-propagation of boundary singularities; and iii) non-overlapping neighboring cells. In addition, problematic grid clustering can be controlled by introducing monitor functions that force the mesh to conform locally to minimum geometric standards.

In this paper, I generate structured non-orthogonal meshes appropriate for generalized Riemannian wavefield extrapolation Shragge (2006) using a differential method advocated by Liseikin (2004). The goal of this paper is not to simply re-develop Liseikin's method. Rather, it is to demonstrate its relevance to geophysical meshing problems by summarizing the general motivation behind the method and highlighting the essential theory and advocated numerical solution. I begin the paper by reviewing why differential methods form an elliptical BVP and describing how to pose the gridding problem within a N-D differential geometric framework. I then present gridding equations for 2-D geometry and provide meshing examples for two seismic imaging applications: wave-equation generation of Green's function estimates and wave-equation migration from topography. The attached appendix provides the numerical details required to implement the differential meshing technique.