\lefthead{Shragge}
\righthead{RWE: NOC}
\footer{SEP--124}
\usepackage{color}
\title{Riemannian Wavefield Extrapolation: Non-orthogonal coordinate systems}
\email{jeff@sep.stanford.edu}
\keywords{}
\author{Jeff Shragge}

\sepinput{mainmacros}
\sepinput{macros}

\ABS{Riemannian wavefield extrapolation (RWE) is a method for
  performing one-way wave propagation on generalized coordinate
  meshes.  Previous RWE implementations assume that  coordinate
  systems are either orthogonal or semi-orthogonal, which can lead to
  situations where meshes suffer from problematic bunching,
  singularities, or multivaluedness.  This paper develops a procedure
  for avoiding these problems that is based on wavefield extrapolation
  on smoother, but generally non-orthogonal, coordinate system meshes.
  The corresponding  extrapolation wavenumber is comprised of a number
  of mixed-domain fields, in addition to velocity, that encode
  coordinate system geometry.  An extended split-step Fourier
  approximation of the extrapolation operator is developed that
  provides accurate results when multiple reference parameter sets are
  used.  Two analytic coordinate system examples are presented to
  validate the theory.  A procedure for generating triplication-free
  coordinate systems in complex media is developed, which is followed
  by calculations of 2D and 3D Green's \textcolor{red}{function} estimates in
  cylindrical and near-spherical geometries.   An important
  consequence is that greater emphasis can be placed on generating
  smooth computational meshes for RWE rather than satisfying more
  restrictive semi-orthogonal criteria.}        
%
\section{Introduction}
%
A persistent goal of wave-equation migration research is to improve
seismic imaging capabilities in complex geologic settings.  Although
ubiquitous velocity model uncertainty and uneven illumination can
contribute greatly to image interpretation ambiguity in these
contexts, extrapolation operator inaccuracy remains a significant
problem.  The central issues with conventional one-way wave-equation
extrapolation operators are well documented: while naturally handling
wavefield multipathing in the presence of lateral velocity variation,
they are of limited accuracy at steep propagation angles and cannot
propagate overturning waves by design (though extensions including
phase-shift table lookup \cite{Hale} or two-pass migration
\cite{Zhang} obviate some of these concerns).  Propagation errors are
subsequently manifest in migration images as defocused or misplaced
reflectors or even by a complete absence of interpretable
reflectivity. Accordingly, minimizing these deleterious effects should
improve image quality and any subsequent interpretation based thereon.   
%
\par
%
One strategy for reducing extrapolation operator inaccuracy is to
decompose the complete computational grid into subdomains oriented in
the wave propagation direction (see, e.g., Gaussian beams
\cite{GEO66-04-12401250}, beam-waves 
\cite{BDEtgen-SEG03}, coherent states \cite{Albertin-SEG01} or tilted
Cartesian meshes \cite{Etgen-SEG02}).  The key concept in each of
these approaches is that a judicious choice of reference frame lowers
the effective propagation angle, reducing the need for expensive
extrapolation operators and enabling imaging of overturned waves.
\longcite{Sava.geo.rwepub} followed this approach in developing
Riemannian wavefield extrapolation (RWE), a theory of one-way
wavefield propagation through a generalized Riemannian space.  This
formulation specifies wave-equation operators appropriate for
wavefield extrapolation on generalized computational meshes.  One
important ramification is that the user is free to specify the degree
to which wave-propagation effects are incorporated into the
computational mesh.  However, finding the optimal trade-off between
computational mesh simplicity, how well the mesh conforms to the
wavefield propagation direction, and the computational cost is not a
straightforward task. 
%
\par
%
RWE was initially implemented to model high-quality Green's functions.
This process involved extrapolating wavefields on a point-source
coordinate mesh comprised of a suite of rays traced beforehand through
a smoothed version of the migration velocity model.  Hence, RWE
computational meshes explicitly were asserted to exhibit ray-field
characteristics: semi-orthogonal geometry with an extrapolation
direction (i.e. travel-time along a ray) orthogonal to the two other
axes (i.e. shooting angles) that are not necessarily mutually
orthogonal.  This supposition results in a wave-equation dispersion
relationship for semi-orthogonal meshes that contains a number of
mixed-domain fields, in addition to velocity, that encoded coordinate
system geometry.  In most examples, modeled Green's function estimates
interpolated into the Cartesian domain are highly accurate at steep
propagation angles; however, accuracy is compromised in certain
situations exhibiting certain unfavorable mesh characteristics. 
%
\par
%
In general, the semi-orthogonality assertion is overly restrictive
because it precludes using non-orthogonal computational meshes
designed to be smoothly varying, with singularities and single-valued.
One problematic example occurs in coordinate system triplications (see
upper panel of figure~\ref{fig:Problem}), that arise wherever a mesh
is developed from a triplicating field of rays.  This generates
spatial singularities that lead to zero-division of the wavefield
during extrapolation and cause infinite amplitudes.  Although
ray-coordinate triplications can be avoided by iterative velocity
model smoothing, this less-than-ideal solution counters the goal of
having a coordinate system conformal to the wavefield propagation
direction.   A second example of restrictive semi-orthogonal geometry
is illustrated in \longcite{Shragge.seg.2005}, who formulate a
wave-equation migration from topography strategy that poses wavefield
extrapolation directly in locally orthogonal meshes conformal to the
acquisition surface.  This approach successfully generates subsurface
images beneath areas exhibiting longer wavelength and lower amplitude
relief; however, imaging results in situations involving more rugged
acquisition topography degrade due to the grid compression/extension
demanded by semi-orthogonality geometry (see lower panel of
figure~\ref{fig:Problem}b).  
%
\par
%
In this paper, I argue that the RWE approach can be made more broadly
applicable by resolving the issues related to semi-orthogonal geometry
using a combination of localized computational mesh smoothing and the
judicious removal of all coordinate triplication branches.  In doing
so, though,  the asserted semi-orthogonality mesh geometry is lost,
and one must propagate wavefields on non-orthogonal grids.
Fortunately, a non-orthogonal extrapolation wavenumber can be isolated
from the general RWE dispersion relationship.  This wavenumber is used
in forming a basis for a one-way extrapolation operator that is used
in the usual sense of  wave-equation imaging. 
% In comparison to
%semi-orthogonal geometry, propagation on non-orthogonal meshes
%introduces two additional dispersion relationship terms, and makes the
%existing coefficients slightly more involved.  These spatially varying
%coefficients, though, can be handled through an extended split-step
%Fourier approach \cite{GEO55-04-04100421}.   
%
\par
%
The goal of this paper is develop and implement a one-way
wave-equation extrapolation operator appropriate for RWE in 3D
non-orthogonal coordinates.  A second goal is to specify a procedure
for generating unconditionally triplication-free computational meshes.
The development presented herein naturally follows that of
\longcite{Sava.geo.rwepub}; however, I recast the theory in a more
compact notation that leads to a closer analytic connection of the
generalized computation geometry with the underlying Cartesian grid.
The paper begins with the formulation of the 3D Riemannian acoustic
wave-equation and the corresponding non-orthogonal one-way wavefield
extrapolation wavenumber.  Appendix~A presents an overview of the
required differential geometry theory, while the split-step Fourier
approximation corresponding to the extrapolation operator developed
herein is presented in Appendix~B.  Two analytic 2D non-orthogonal
coordinate system examples are then provided to validate the theory. 
The final sections detail the procedure for generating
triplication-free coordinate systems, and present 2D and 3D Green's
functions estimates modeled in cylindrical and near-spherical
coordinates, respectively. 
%
\plot{Problem}{width=7in}{Illustration of problems with RWE
  computational grids.  Top panel: Ray-coordinate system
  triplications. Bottom panel: Grid bunching for topographically
  conformal coordinate system.} 
%
\section{Acoustic wave-equation in 3D generalized Riemannian spaces}
%
Specifying the acoustic wave-equation in a 3D Riemannian space
requires formulating the physics of wave-propagation in a generalized
coordinate system framework.  By definition, generalized Riemannian
coordinates are related to the underlying Cartesian mesh by unique
transformations (i.e. non-triplicating and one-to-one).   In this
paper, I use a notation where a generalized coordinate system $\bxi =
\{ \xi_1,\xi_2,\xi_3 \}$ to mapped to a Cartesian grid  $\mathbf{x}=\{
x_1,x_2,x_3 \}$ through transformation $\xi_i (x_j) = f_i$.  Provided
these conditions are met, the acoustic wave-equation for a wavefield,
$\W$, in a generalized Riemannian space is,  
\beq \label{eqn:helm}
\dd_\bxi \W= - \ww^2\ss^2 \lp \mathbf{\bxi} \rp \W,
\eeq
where $\dd_\bxi$ is the Laplacian operator applied in coordinates
$\bxi$, $\ww$ is frequency, and $\ss$ is the propagation slowness. 
%
\par
%
Correctly formulating the wave-equation in coordinate system $\bxi$
requires that Laplacian operator $\dd_\bxi$ be specified by
differential geometry relationships.  (An overview of necessary
differential geometry theory is provided in Appendix~A.) The Laplacian
operator in generalized coordinates is, 
\beq \label{eqn:laplac2}
  \dd_\bxi \W = \frac{1}{\sqrt{|\mathbf{g}|}}\,
    \frac{\partial}{\partial \xi_i}\,
    \left(\,\sqrt{|\mathbf{g}|} g^{ij}\,\frac{\partial \W}{\partial \xi_j}\right)
    =
  \frac{1}{\sqrt{|\mathbf{g}|}}\,
    \frac{\partial}{\partial \xi_i}\,
    \left(\,m^{ij}\,\frac{\partial \W}{\partial \xi_j}\right), \quad i,j=1,2,3,
\eeq
where $g^{ij}$ is an element of the metric tensor $\mathbf{g}$,
$|\mathbf{g}|$ is the metric tensor discriminant, and
$m^{ij}=\sqrt{|\mathbf{g}|} g^{ij}$ is weighted metric tensor element
that enables a more compact notation.  Unless otherwise stated,
summation over all repeated indices (i.e. $i,j=1,2,3$) is assumed
throughout this paper.  Substituting equation~\ref{eqn:laplac2}\ into
\ref{eqn:helm} generates a Helmholtz equation appropriate for
propagating waves through a 3D Riemannian space,  
\beq \label{eqn:we1}
\frac{1}{\gdetf} \de\qi \lp \mmm\ii\jj \done\W\qj \rp = - \ww^2\ss^2 \W.
\eeq
%
\par
%
The first step in developing a generalized RWE wave-equation
dispersion relationship is to expand the derivative terms in
equation~\ref{eqn:we1} and multiply through by $\gdetf$ to obtain,   
\beq \label{eqn:helm1}
\done{\mmm\ii\jj}\qi \done\W\qj + \mmm\ii\jj \mtwo\W\qi\qj
= - \gdetf \ww^2 \ss^2 \W.
\eeq
The spatial derivative of the weighted metric tensor in the first term
of equation~\ref{eqn:helm1} is written concisely using the following
substitution, 
\beq
n^j=\done{\mmm\ii\jj}\qi = 
\done{m^{1j}}\qx + 
\done{m^{2j}}\qy +
\done{m^{3j}}\qz.
\eeq
Fields $n^j$ are interpreted as measures of the rates by which space
expands, compresses and/or shears in the $j^{th}$ direction and can be
non-zero even for orthogonal coordinate systems.  Using this
substitution, equation~\ref{eqn:helm1} is rewritten, 
\beq \label{eqn:helm2}
n^j \done\W\qj + \mmm\ii\jj \mtwo\W\qi\qj
= -\gdetf \ww^2 \ss^2 \W.
\eeq
A wave-equation dispersion relation is developed by replacing the
partial differential operators acting on wavefield $\W$ with their
Fourier domain wavenumber duals \cite{Claerbout85}, 
\beq \label{eqn:disp}
\lp m^{ij} \ki - {\rm i} n^j \rp \kj = \gdetf \ww^2\ss^2,
\eeq
where ${\rm i} \ki$ is the Fourier domain dual of differential operator
$\frac{\partial}{\partial \xi_i}$.  
%
\par
%
Equation~\ref{eqn:disp} represents the dispersion relationship
required to propagate a wavefield through a generalized 3D Riemannian
space.   Quantity $m^{ij}$ in the first term, $m^{ij}\ki\kj$, is a
measure of the dot product between wavenumber vectors in the $\ki$ and
$\kj$ directions (i.e. orthogonal wavenumbers will have a zero
$m^{ij}$ coefficient).  Fields $n^j$ in the second term, ${\rm i} n^j
\kj$, represent a scaling of wavenumber $\kj$ caused by local
expansion, contraction and/or shearing of the coordinate system in the
j$^{th}$ direction. 
%
\par
%
Note that the expression in equation~\ref{eqn:disp} reduces to the more
familiar Cartesian expression when introducing $n^{j}=0$ and
$m^{ij}=\delta^{ij}$, which leaves
\beq
\ki \ki= \kx^2 + \ky^2 + \kz^2 = \ww^2 s^2.
\eeq
%
\subsection{Extrapolation wavenumber isolation}
%
Specifying a one-way extrapolation operator requires isolating one of
the wavenumbers in equation~\ref{eqn:disp}.  In this paper, I
associate the extrapolation direction with coordinate $\xi_3$.
Expanding equation~\ref{eqn:disp} by evaluating indices and
rearranging the results yields,    
\begin{small}
\beq \label{eqn:disp2}
m^{33}\kz^2 + \lp 2m^{13}\kx + 2m^{23}\ky - {\rm i} \,n^3 \rp \kz =
\gdetf \ww^2 \ss^2 + {\rm i} \lp n^1\kx + n^2\ky \rp - m^{11}\kx^2 -
m^{22}\ky^2 - 2m^{12}\kx\ky.
\eeq
\end{small}
An expression for wavenumber $\kz$ is obtained through a
complete-the-square transform,  
\begin{small}
\beqa
\,& m^{33} \lp \kz + 
\frac{ \lp 2 m^{13}\kx + 2m^{23}\ky- {\rm i} \,n^3 \rp }{ 2m^{33} } 
\rp^2 =  \,
\gdetf\ww^2\ss^2 -
\kx^2 \lp m^{11}-\frac{ \lp m^{13} \rp^2 }{ m^{33}} \rp - \\
\,& \ky^2 \lp m^{22}-\frac{ \lp m^{23} \rp^2 }{ m^{33}} \rp - 
\kx\ky \lp 2 m^{12}-\frac{ 2 m^{13}m^{23}}{ m^{33}} \rp +
{\rm i} \, \kx \lp n^1 - \frac{m^{13}\,n^3}{ m^{33} } \rp + {\rm i} \, \ky \lp n^2
- \frac{m^{23}\,n^3}{ m^{33} }\rp - \frac{(n^3)^2 }{ m^{33}}.
\,\nonumber
\eeqa
\end{small}
Isolating the contributions of wavenumber $\kz$ yields,
\beq \label{eqn:kzfinal}
\kz =
- a_1 \kx
- a_2 \ky 
+ {\rm i} a_3  
\pm
\lb
  a_4^2 \ww^2 
- a_5^2 \kx^2 
- a_6^2 \ky^2
- a_7 \kx\ky 
+ a_8 {\rm i}\,\kx 
+ a_9 {\rm i}\,\ky  
- a_{10}^2  
\rb^{\frac{1}{2}},
\eeq
where $a_i$ are non-stationary coefficients given by,
\begin{tiny}
\beq \label{eqn:ametric2}
\mathbf{a}
=
\left[
\frac{g^{13}}{ g^{33}} \;\;\;\;\;
\frac{g^{23}}{ g^{33} }\;\;\;\;\;
\frac{ n^3 } {2 m^{33} } \;\;\;\;\;
\frac{ \ss }{ \sqrt{g^{33}} } \;\;\;\;\;\;
\sqrt{ \frac{g^{11}}{ g^{33}} - \lp \frac{g^{13}}{g^{33}}  \rp^2 } \;\;\;\;\;
\sqrt{\frac{g^{22}}{ g^{33} } - \lp \frac{ g^{23}}{g^{33}}  \rp^2 }
\;\;\;\;\; 
\lb \frac{2\,g^{12}}{g^{33}} - \frac{2\, g^{13}\,g^{23}}{ \lp g^{33}
  \rp^2} \rb \;\;\;\;\; 
\lb \frac{n^1}{m^{33}} - \frac{m^{13}\, n^3}{ \lp m^{33}\rp^2}\rb \;\;\;\;\;
\lb \frac{n^2}{m^{33}} - \frac{m^{23}\, n^3}{ \lp m^{33}\rp^2}\rb \;\;\;\;\;
\frac{ n^3 }{ m^{33} }
\right]^{\mathbf{T}}.
\eeq
\end{tiny}
Note that the coefficients contain mixed $m^{ij}$ and $g^{ij}$
terms and positive definite terms $a_4, a_5, a_6$ and $a_{10}$ are
squared.  
%
\par
%
The special Cartesian case is again recovered from the two equations
above by substitution $n^j=0$ and $g^{ij}=\delta^{ij}$ (for all values
$i$ and $j$ indicies) in the coefficients in
equation~\ref{eqn:ametric2}
\beq
\kz=\left[ s^2\ww^2-\kx^2-\ky^2  \right]^{\frac{1}{2}}.
\eeq  
%
\par
%
The dispersion relationship specified by equations~\ref{eqn:kzfinal}
and \ref{eqn:ametric2} contains ten coefficients that represent
mixed-domain fields.  Similar to Cartesian-based wavefield
extrapolation through homogeneous media, a constant-coefficient
Fourier-domain ($\omega-\boldsymbol{k_\xi}$) phase-shift extrapolation
scheme can be developed to recursively advance a wavefield from level
$\xi_3$ to level $\xi_3+\Delta \xi_3$, 
%
\beq \label{eqn:bulk}
\W(\xi_3+\Delta \xi_3, \kx,\ky ; \omega) = 
\W(\xi_3, \kx, \ky ; \omega) {\rm e}^{ {\rm i} \kz \Delta \xi_3}. 
\eeq
%
If $\W$ represents a post-stack wavefield, an image $\mathcal{I}$ can
be produced from the propagated wavefield through an imaging condition
\cite{Claerbout85}, 
%
\beq \label{eqn:ZOI}
\mathcal{I}(\xi_3,\xi_1,\xi_2) = \sum_{\omega} \W(\xi_3,\xi_1,\xi_2; \omega).
\eeq
%
\par
%
Situations where coefficients vary across an extrapolation step,
though, required additional  approximations.  One straightforward
approach is to use an extended split-step Fourier (SSF) method
\cite{GEO55-04-04100421}, which is detailed in Appendix~B.  This
method uses Taylor expansions of the dispersion relation about a set
of reference parameters to form a bulk phase-shift operator in the
Fourier domain ($\omega-\boldsymbol{k_\bxi}$).  Differences between
the reference and true parameters then form a correction term applied
in the mixed $\omega-\bxi$ domain.  In addition, the phase-shift
extrapolation scheme can be extended to incorporate multiple reference
media followed by interpolation similar to the phase-shift plus
interpolation (PSPI) technique of \longcite{Gazdag.1984}. 
% 
\subsection{Discussion}
%
The accuracy of the extended SSF approach is directly related to the
degree to which coefficients in equation~\ref{eqn:ametric2} vary
across a propagation step.  At a 
first glance, one might expect that far too many expansions are
required to make a PSPI approach practical.  (For example, three
reference expansions for each of the ten terms would seemingly require
$3^{10}=22599$ separate wavefield extrapolations.)  However, three
factors combine to greatly reduce the total number of required
reference coefficient sets.  
%
\par
%
First, the $a_i$ coefficients in equation~\ref{eqn:kzfinal} tend to be
highly correlated because they are composed of similar metric tensor
elements $g^{ij}$.  Thus, the central issue is how accurately can
we characterize vector coefficient fields.  Fortunately, this problem
is similar to the quantization problem in computer graphics: What is
the fewest number of colors by which an image can be represented given
a maximum allowable error?  In this paper, I calculate reference
coefficients using a multi-dimensional Lloyd's algorithm developed by
\longcite{Clapp.Lloyds}.  This iterative procedure allows for
variable rate quantization that represents the multi-dimensional
histogram of the coefficients with the sparsest number of points
within a specified error tolerance. 
% 
\par
%
Second, numerous situations exists where some coefficients are zero or
negligible.  For example, the mixed-domain fields for a 3D weakly
non-orthogonal coordinate system (i.e. ${\rm max}
|\{g^{12},g^{13},g^{23}\}| << {\rm max} \{g^{11},g^{22},g^{33} \}$)
within a kinematic approximation (i.e. no imaginary terms) can be
accurately approximated with only four coefficients.  (Appendix~C
details situations where additional approximations are appropriate.)
Third, one may apply algorithms that locally smooth the coordinate
system mesh, which reduces the spatial variability of the coefficients
and allows a more reliable representation of wavenumber $\kz$.   
%
\section{Numerical Modeling Examples}
%
\subsection{2D sheared Cartesian coordinate system}
%
An instructive analytic coordinate system to examine is a 2D sheared
Cartesian grid formed by a uniform shearing action on a 2D Cartesian
mesh (see left panel of figure~\ref{fig:NewRT}).  This coordinate
system is uniquely specified by one additional degree of freedom and
is related to an underlying Cartesian mesh through the following
transformation,  
%
\beq \label{eqn:2Dred}
  \left[\begin{array}{c}
  	\xone \\
  	\xthr
  	\end{array}\right] 
  	= 
  \left[\begin{array}{cc}
      1 & \rm{cos}\, \theta      \\
      0 & \rm{sin}\, \theta
  \end{array}\right]
  \left[\begin{array}{c}
    \qx \\
    \qz \\
  \end{array}\right],
\eeq
where $\theta$ is the shear angle of the coordinate system
($\theta=90^{\circ}$ is Cartesian).  The metric tensor of this
transformation is, 
\beq \label{eqn:2Dmet}
\left[g_{ij}\right]
=
\left[\begin{array}{cc}
\hxxdef & \hxzdef \\
\hxzdef & \hzzdef
\end{array}\right]
=
\left[\begin{array}{cc}
    \EE & \HH \\
    \HH & \AA
\end{array}\right] 
= 
\left[\begin{array}{cc}
    1                & \rm{cos}\, \theta \\
    \rm{cos}\,\theta & 1  
\end{array}\right],
\eeq
and has a discriminant $\gdet = \rm{sin}^2\,\theta$ and an associated
metric tensor $g^{ij}$ given by, 
\beq \label{eqn:2Dmet6}
\left[g^{ij}\right]
=
\frac{1}{\rm{sin}^2\,\theta}
\left[\begin{array}{cc}
    1                 & - \rm{cos}\,\theta \\
    -\rm{cos}\,\theta & 1
\end{array}\right].
\eeq
%
Note that because the tensor in equation~\ref{eqn:2Dmet6} is
coordinate invariant, equation~\ref{eqn:helm2} simplifies to,   
\beq \label{eqn:laplac}
  \dd_\bxi \W = g^{ij} \,\frac{\partial^2 \W}{\partial \qi \partial \qj}
  = - \ww^2\ss^2\W,
\eeq
which generates the following dispersion relation,
\beq
g^{ij}\ki\kj = \ww^2\ss^2.
\eeq
Expanding out these terms leads to an expression for wavenumber $\kz$, 
\beq \label{eqn:kz2Dex}
\kz =
-\frac{ g^{13} }{ g^{33} } \kx
\pm
\sqrt{\frac{\ss^2\ww^2}{g^{33}}- \lp\frac{g^{11}}{g^{33}}-\lp
  \frac{g^{13}}{g^{33}} \rp^2\rp \kx^2}. 
\eeq
Substituting the values of the associated metric tensor in
equation~\ref{eqn:2Dmet6} into equation~\ref{eqn:kz2Dex} yields, 
\beq
\kz =
\rm{cos} \, \theta \;\kx
\pm \rm{sin} \,\theta
\sqrt{\ss^2\ww^2 - \kx^2}, 
\eeq 
which is appropriate for performing RWE on the 2D sheared Cartesian
coordinate system shown in left panel in figure~\ref{fig:NewRT}. 
%
\par
%
The right panel of figure~\ref{fig:NewRT} shows the results of
extrapolating plane waves in a Cartesian coordinate system sheared at
25$^{\circ}$ from vertical.  The background velocity model is 1500~m/s
and the zero-offset data consist of 4 flat plane-waves at times
$t=0.2, 0.4, 0.6$ and $0.8$ s.  Zero-offset migration results
generated by equation~\ref{eqn:ZOI} show migrated reflectors at the
expected depths of $z$=300, 600, 900, and 1200~m.  Note that the
propagation creates explainable boundary artifacts: reflections on the
left are due to a truncated coordinate system and hyperbolic
diffractions on the right are caused by truncated plane-waves.   
%
\sideplot{NewRT}{width=7in}{Sheared Cartesian coordinate system test.
  Coordinate system (left panel) shear angle and velocity are
  $\theta=25^{\circ}$ and 1500 ms$^{-1}$, respectively.  Zero-offset
  data consist of 4 flat plane-wave impulses at $t=0.2, 0.4, 0.6$ and
  $0.8$ s that are correctly imaged (right panel) at depths $z=300,
  600, 900,$ and $1200~m$.}    
%
\subsection{Polar Ellipsoidal Coordinates}
%
A second example is an ellipsoidal polar coordinate system (see upper
left panel in figure~\ref{fig:Polar}) appropriate for migrating
overturning plane-waves.  A polar ellipsoidal coordinate system is
specified by,  
\beq \label{eqn:2DPE1}
  \left[\begin{array}{c}
  	\xone \\
  	\xthr
  	\end{array}\right] 
  	= 
  \left[\begin{array}{c}
    a(\xi_3)\,\xi_1\, \rm{cos}\,\xi_3 \\
    a(\xi_3)\,\xi_1\, \rm{sin}\,\xi_3 \\
  \end{array}\right],
\eeq
where parameter coordinate $\xi_1$ is the radius from the center
focus, $\xi_3$ is polar angle, and is $a=a(\xi_3)$ is a smooth
function that controls coordinate system ellipticity and has curvature
parameters $b=\frac{\partial a}{\partial \xi_3}$ and
$c=\frac{\partial^2 a}{\partial \xi_3^2}$.  The metric tensor $g_{ij}$
for the polar ellipsoidal coordinate system defined in
equation~\ref{eqn:2DPE1} is,
\beq \label{eqn:2DPE2}
\left[g_{ij}\right]
=
\left[\begin{array}{cc}
    a^2 & \xi_1\, a\, b \\
    \xi_1 \, a \, b & \xi_1^2(b^2+a^2) 
\end{array}\right],
\eeq
and has a metric discriminant given by $\gdet = a^4 \xi_1^2$. The
associated and weighted associated metric tensors are given by,  
\beq 
\left[g^{ij}\right]
=
\left[\begin{array}{cc}
    \frac{b^2+a^2}{a^4}   & -\frac{b}{a^3\,\xi_1} \\
    -\frac{b}{a^3\,\xi_1} & \frac{1}{a^2\, \xi_1^2}
\end{array}\right]
\quad
{\rm and} 
\quad \label{eqn:2DPE4}
\left[m^{ij}\right]
=
\left[\begin{array}{cc}
    \frac{\xi_1(b^2+a^2)}{a^2} & -\frac{b}{a} \\
    -\frac{b}{a} & \frac{1}{\xi_1}
\end{array}\right].
\eeq
%	
Tensors $g^{ij}$ and $m^{ij}$ are used to form the extrapolation
wavenumber appropriate for one-way wavefield propagation on a 2D polar
ellipsoidal mesh.  However, because the computational mesh is
non-stationary, we must also compute the $n^i$ fields:
$n^1=\frac{a^2+2b^2-ac}{a^2}$ and $n^3=0$.  Inserting these values
leads to the following extrapolation wavenumber expression (see
equation~\ref{eqn:2DNO}),     
\beq \label{eqn:2Dex5}
\kz = \frac{\xi_1 b}{a} \kx \pm \sqrt{ a^2 \xi_1^2 s^2 \omega^2 -
  \xi^2_1 \kx^2 - i \kx \xi_1 \lp \frac{a^2+2b^2-ac}{a^2} \rp }.
\eeq
The kinematic approximation of equation~\ref{eqn:2Dex5} is given by,
\beq \label{eqn:2Dex6}
\kz = \xi_1 \lb \frac{b}{a} \kx \pm \sqrt{ a^2 s^2 \omega^2 -
  \kx^2 }\rb, 
\eeq
and the orthogonal polar case (i.e. $a=1$ and $b=0$) yields, 
\beq \label{eqn:2Dex7}
\kz = \pm \xi_1 \sqrt{ s^2 \omega^2 -\kx^2 },
\eeq
which is examined in \longcite{Nichols}.
%
\par
%
Figure~\ref{fig:Polar} shows a wavefield extrapolation example for an
ellipsoidal polar coordinate system in equation~\ref{eqn:2DPE1}
defined by ellipticity parameter
$a(\xi_3)=1+0.2\,\xi_3-0.05\,\xi_3^2$.  The upper and lower panels of
figure~\ref{fig:Polar} correspond to velocity/coordinate and wavefield
domains, respectively.  Similarly, the left and right panels represent
the Cartesian and Riemannian domains.  Note that wavefield
interpolation between the latter two domains is possible because of
the established mapping relationships.  The upper left panel shows the
polar ellipsoidal coordinate system mesh overlying a linear
$v(z)=1500+0.35\,z$ m/s velocity function.  The upper right panel
presents the velocity model mapped into the RWE domain under the
transformations defined in equation~\ref{eqn:2DPE1}.    
%
\plot{Polar}{width=6.5in}{Ellipsoidal polar coordinate system test
  example.  Upper left: $v(z)=1500+0.35\,z$ velocity  function
  overlain by a polar ellipsoidal coordinate system defined  by
  parameter $a=1+0.2\,\xi_3-0.05\,\xi_3^2$.  Upper right: Velocity
  model mapped in the RWE domain.  Bottom right: Imaged reflectors in
  RWE domain. Bottom left: the RWE domain image mapped to a Cartesian 
  mesh.}  
%
The data used in this test consist of ten plane-waves defined by a ray
parameter $p_x=-0.5$ s/km.  At greater depths the waves propagate
faster and go through a turning point before propagating upward to the
left. The tops of the waves, though, still travel
through slower material and have not yet overturned.  One interesting
observation is that if a propagating plane-wave wavepath is well
represented by a single ellipticity parameter $a=a(\xi_3)$, then a
polar elliptical mesh can form the coordinate basis for a plane-wave
migration strategy that is potentially more dynamic than tilted
Cartesian coordinates. 
%
\section{Generating Triplication-free Coordinate Meshes}
%
A RWE computational mesh design challenge is finding a coordinate
system that is fairly conformal to the wavefield propagation
direction, but is unconditionally triplication-free.   This problem is
illustrated by the panels in figure~\ref{fig:GaussExample}.  The top
panel shows a v(z) velocity model with three Gaussian anomaly
inclusions overlain by a ray-coordinate system calculated by the
Huygens rayfront tracing method \cite{Sava.hwt}.  Note that the
anomalies cause both mesh triplications to the left and right of the
model, as well as a grid rarefaction directly beneath the shot-point.  
%
\par
%
The center panel shows a first-arrival Eikonal equation solution for
the same shot-point presented in the top panel.  Plotted over-top are
single-valued Eikonal solution isochrons.  Note that isochron
generally conform to the propagation direction and can be used to
construct the extrapolation steps of a RWE computational mesh.  The
first mesh generation step is to extract initial and final isochron
surfaces from the Eikonal equation solution to form the inner and
outer mesh boundaries.  The mesh domain is then enclosed by
interpolating between the edges of the inner and outer bounding
surfaces.  The interior mesh can then be formed through bi-linear
interpolation methods, such as blending functions \cite{Liseikin,ShraggeMesh}. 
%
\par
%
A corresponding triplication-free mesh is presented in the bottom
panel of figure~\ref{fig:GaussExample}.  The grid is regularly-spaced
on the outer isochron and has a couple of dimples at the locations of
removed triplications.   These discontinuities have been reduced by
applying a smoothing operator to the Eikonal equation solution before
calculating the mesh.  Importantly, smoothing does not greatly affect
propagation accuracy because the coordinate system mesh forms only the
skeleton on which wavefield extrapolation occurs.   
%
\plot{GaussExample}{width=6.5in}{Example of triplication-free mesh
  generation. Top panel: Velocity model with three Gaussian velocity
  perturbations.  Overlain is a coordinate mesh generated from
  ray-tracing.  Note the triplication to either side of the
  shot-point, as well as the spreading beneath the shot point.  Middle
  panel:  Velocity model overlain by isochrons of an Eikonal equation
  solution for same shot-point.  Bottom panel:  Triplication-free
  computational mesh generated by through the approach specified in
  this section. }  
%
\subsection{2D Green's Function Generation}
%
The third test uses RWE to model 2D Green's functions on coordinate
systems constructed by the aforementioned approach.
Figure~\ref{fig:SPexamp} presents the slice of the SEG-EAGE salt
velocity model chosen for the test.  Importantly, the contrast between
the salt body and sediment velocities causes complicated wavefield
propagation including triplication and multipathing.  The upper left
panel shows the velocity model with an overlain coordinate system
generated by meshing procedure discussed in the section above.  The
velocity model in the RWE domain is illustrated in the upper right
panel.    
%
\par
%
The lower right panel shows the impulse response tests in the RWE
domain.  The impulses conform fairly well to the travel-time steps,
except where they enter the salt body in the lower left of the image.
The migration results mapped back to Cartesian space are shown in the
lower left panel.  The complicated wavefield left of the shot point
advances through the salt body and is subsequently refracted upward.
Note that wide-angle reflections off the top-salt/sediment boundary
are also strongly present. 
%
\plot{SPexamp}{width=6.5in}{Example of wave-equation generated Green's
  Functions on structured non-orthogonal mesh for a slice through the
  SEG-EAGE salt velocity model.  Top left: Salt model  in physical
  space with an overlain ray-based mesh. Top right: Velocity model in
  the RWE domain.  Bottom right: Wavefield propagated in
  ray-coordinates through velocity model shown in the upper right.
  Bottom left: Wavefield in bottom right interpolated back to
  Cartesian space.}  
%
Figure~\ref{fig:SPcompare} presents a comparison test between RWE and
Cartesian extrapolation.  Beneath and to the right of the shot-point
the wavefields are fairly similar except for a $90^{\circ}$
phase-change associated with differences between modeling a Cartesian
point-source and a RWE plane-wave.  However, the wavefields are
significantly different to the left of the shot-point because
Cartesian-based extrapolation cannot propagate wavefields laterally
with the same accuracy and upward at all. Hence, this energy is absent
from the propagating wavefield in the lower panel.  
%
\plot{SPcompare}{width=6.5in}{Results of a comparison test between
  Riemannian wavefield extrapolation (top panel) and Cartesian-based
  extrapolation (bottom panel).}  

\subsection{3D Point Source Coordinates}
%
The procedure discussed above can be extended to 3D point-source
coordinate systems.  However, any wavefield extrapolation procedure in
near-spherical coordinates must account for the apparent singularity
at the poles.  One approach is to rotate the poles such that the small
arc-length regions are at the Earth's surface \cite{Schneider}.  A
second approach is to avoid sampling at the existing polar singularity
\cite{Fowler}.  Both solutions, though, have highly variable spatial
sampling because regular spherical angle discritization leads to
highly variable sampling in Cartesian space. 
%
\par
%
Another solution is to choose a discritization of the spherical
surface that has a fairly uniform distribution in Cartesian domain.
Herein, I use a 3D point source gridding technique based on the
Winkel-Tripel cartographic projection \cite{Winkel}.  Note that he
mapping of a spherical surface illustrated in
figure~\ref{fig:Starting} does not eliminate distortions in area,
direction or distance; rather, it tries to minimize the sum of all
three.  The gridding equations for a Winkel-Tripel surface are the
following,
\beqa
\,&\,& x_1 = \xi_3 (2 \xi_1 / \pi + 2 \, w\, \alpha {\rm cos} \, \xi_2
\, {\rm sin}(\xi_1 /2)) /2 \nonumber \\ 
\,&\,& x_2 = \xi_3 ( \xi_2 + 2\, \alpha \, {\rm sin} \, \xi_2 ) /2 \\
\,&\,& x_3 = \sqrt{\xi_3^2 - x_1^2 - x_2^2} \nonumber
\eeqa
where $\xi_3$ is in the radial direction, $\xi_1$ and $\xi_2$ are a
latitude parallel and longitude meridian, respectively, and, 
\beqa
\,&\,& \alpha = {\rm cos}^{-1} \, ({\rm cos} \, \xi_2 \, {\rm cos}
(\xi_1/2)) \nonumber \\ 
\,&\,& w = \left\{ \begin{array}{l}  
			0, \quad \quad \quad {\rm sin} \,\alpha = 0 \\
			({\rm sin} \alpha)^{-1}, \quad {\rm otherwise} 
			\end{array} \right. 
\eeqa
Note that the 3D grid is created by stepping outward in the radial
direction $\xi_3$ (or equally by the Eikonal equation solution method
described above). 
%
\plot{Starting}{width=5in}{Starting surface for Winkel-Tripel
  projection of unit hemisphere.  Note that the cells on the unit
  hemisphere surface have nearly consistent area indicating fairly
  uniform sampling, and that not all parts of the unit hemisphere are
  sampled with this approach.} 
%
\par
%
Figure~\ref{fig:ElfVel} presents an example of a wavefield
extrapolated on 3D point source coordinate mesh.   
%
\plot{ElfVel}{width=7in}{Wavefield extrapolation example through
  Elf-IFP-CGG synthetic velocity model.} 
%
The Elf-IFP-CGG synthetic velocity model used in this example is
inspired from field data observations from  North Sea Block L7d.  The
sediments located beneath the shot point form a mostly v(z) profile,
while the salt body to the right generates a more complex geologic
setting.  Overlain are a number of impulse responses associated with
different travel times.  Note that the wavefield  does not exhibit
azimuthal anisotropy found in 3D impulse response tests based on
finite-difference splitting methods.  Also, the coordinate system
boundaries are restricted at steep angles, which generates boundary
reflection artifacts.  Overall, though, this extrapolation results
demonstrates the stability of wavefield propagation on a  3D grid
defined by the Winkel-Tripel spatial discritization. 
%
\section{Concluding Remarks}
%
This paper extends the implementation of Riemannian wavefield
extrapolation theory to 3D non-orthogonal coordinate systems.  The
extrapolation wavenumber decouples from the other wavenumbers enabling
the use an extended split-step Fourier approximate solution of the
one-way wavefield extrapolation operator.  Also documented is a
methodology for generating computational meshes that are
unconditionally triplication-free.  The provided examples illustrate
that wavefields can be extrapolated on non-orthogonal coordinate
meshes, which opens up a range of possibilities including Riemannian
plane-wave migration in elliptical polar coordinates. 
%
\section{Acknowledgments}
%
I acknowledge the contributions of Paul Sava and Sergey Fomel in
laying the groundwork for the current theory and for ongoing RWE
conversations, and I thank Biondo Biondi, Bob Clapp, Brad Artman, Paul
Fowler, Tom Dickens, and Peter Traynor for enlightening discussions.
The sponsors of the SEP consortium provided funding for this work. 
%
\bibliographystyle{sep}
\bibliography{SEP,MISC,GEOTLE,refs,refs2}

\newpage

\APPENDIX{A}
%
Geometry in a generalized 3D Riemannian space is described by a symmetric
metric tensor, $g_{ij}=g_{ji}$, that relates the geometry in a general
non-orthogonal coordinate system, $\{\xi_1,\xi_2,\xi_3\}$ , to an underlying
Cartesian mesh, $\{x_1,x_2,x_3\}$ \cite{Guggenheimer}.  In
matrix form, the metric tensor is written,    
\beq \label{eqn:metric}
\left[g_{ij}\right] 
= 
\left[\begin{array}{ccc}
g_{11} & g_{12} & g_{13} \\
g_{21} & g_{22} & g_{23} \\
g_{31} & g_{32} & g_{33} 
\end{array}\right] 
= 
\left[\begin{array}{ccc}
g_{11} & g_{12} & g_{13} \\
g_{12} & g_{22} & g_{23} \\
g_{13} & g_{23} & g_{33} 
\end{array}\right],
\eeq
where $\EE$, $\FF$, $\GG$, $\HH$, $\LL$ and $\AA$ are functions
linking the two coordinate systems through, 
\beqa \label{eqn:pderels}
\EE=\hxxdef, \;\;\;\; \FF=\hxydef, \;\;\;\; \GG = \hyydef,\nonumber \\
\HH=\hxzdef, \;\;\;\; \LL=\hyzdef, \;\;\;\; \AA = \hzzdef.  
\eeqa The associated (or inverse) metric
tensor, $g^{ij}$, is defined by $g_{ij}=\,|\mathbf{g}|\,g^{ij}$, where
$|\mathbf{g}|$ is metric tensor matrix discriminant.   The associated
metric tensor is given by,   
\beq \label{eqn:ametric}
  \left[g^{ij}\right] =\frac{1}{\left|\mathbf{g}\right|}
  \left[\begin{array}{ccc}
      \GG\AA -\LL^2 & \HH\LL-\FF\AA & \FF\LL-\HH\GG \\
      \HH\LL-\FF\AA & \EE\AA-\HH^2  & \FF\HH-\EE\LL \\
      \FF\LL-\HH\GG & \FF\HH-\EE\LL & \EE\GG-\FF^2 
    \end{array}\right],
\eeq
and has the following metric discriminant,
\beq \label{eqn:det}
  |\mathbf{g}| = \AA\,\JJ\, \lb 1 -
  \frac{\EE\LL^2+\GG\HH^2-2\FF\LL\HH}{\AA(\EE\GG-\FF^2)} \rb .
\eeq
A weighted metric tensor, $m^{ij}=\gdetf \, g^{ij}$, is also used
throughout the paper.   
%
\newpage
%
\APPENDIX{B}
%
The extrapolation wavenumber defined in equations~\ref{eqn:kzfinal}
and \ref{eqn:ametric2} cannot generally be implemented exactly in the
Fourier domain due to a simultaneous spatial dependence (i.e. a
function of both $x_1$ and $\kx$). This can be addressed using an
extended version of the split-step Fourier approximation
\cite{GEO55-04-04100421} that uses Taylor expansions to separate $\kz$
into two parts: $\kz \approx \kz^{PS} + \kz^{SSF}$.  Wavenumbers
$\kz^{PS}$ and $\kz^{SSF}$ represent a pure Fourier
($\omega-\mathbf{k_x}$) domain phase-shift and a mixed
($\omega-\mathbf{x}$) domain split-step correction, respectively.  
%
\par
The phase-shift term is given by,   
\beq
\kz^{PS} = 
- b_1 \kx
- b_2 \ky 
+ i b_3  
\pm
\lb
  b_4^2 \ww^2 
- b_5^2 \kx^2 
- b_6^2 \ky^2
- b_7 \kx\ky 
+ b_8 i\,\kx 
+ b_9 i\,\ky  
- b_{10}^2  
\rb^{\frac{1}{2}},
\eeq
where $b_i=b_i(\xi_3)$ are reference values of $a_i =a_i(\xi_1,\xi_2,\xi_3)$. 
The split-step approximation is developed by performing a Taylor
expansion about each coefficient $a_i$ and evaluating the results
at stationary reference values $b_i$.  Assuming that the stationary values of
$\kx$ and $\ky$ are zero, the split-step correction is  as follows,
\beq \label{eqn:kzssf}
\kz^{SSF}=
\left. \frac{\partial \kz}{\partial a_3} \right|_0 \lp a_3-b_3 \rp +
\left. \frac{\partial \kz}{\partial a_4} \right|_0 \lp a_4-b_4 \rp + 
\left. \frac{\partial \kz}{\partial a_{10}} \right|_0 \lp a_{10}-b_{10} \rp,
\eeq
where ``$0$'' denotes "with respect to a reference medium".  The
partial differential expressions in equation~\ref{eqn:kzssf} are,    
\beq
\left. \frac{\partial \kz}{\partial a_3} \right|_0 
= b_3,
\;\;\;\;\;
\left. \frac{\partial \kz}{\partial a_4} \right|_0
= \frac{b_4 \,\ww^2}{\sqrt{b_4^2 \, \ww^2 - b_{10}^2 }},
\;\;\;\;\;
\left. \frac{\partial \kz}{\partial a_{10}} \right|_0
= - \frac{b_{10}}{\sqrt{b_{10}^2 \, \ww^2 - b_{10}^2 }},
\;\;\;\;\;
\eeq
resulting in the following split-step Fourier correction wavenumber, 
\beq
\kz^{SSF} = i\,b_3 \, \lp a_3-b_3 \rp +
\frac{b_4 \, \ww^2 \, \lp a_4-b_4 \rp}{\sqrt{b_4^2 \, \ww^2 - b_{10}^2
}} - \frac{b_{10} \, \lp a_{10}-b_{10} \rp}{\sqrt{b_4\,\ww^2 -
    b_{10}^2 }}. 
\eeq
%
\newpage
%
\APPENDIX{C}
%
The extrapolation wavenumber developed in equation~\ref{eqn:kzfinal}
is appropriate for any non-orthogonal Riemannian geometry.  However,
there are a number of situations where symmetry or partial
orthogonality are present.  Moreover, one may wish to make a kinematic
approximation where all of the imaginary components of the wavenumber
are ignored.  These situations are discussed in this Appendix. 
%
\par
%
{\bf 3D Semi-orthogonal Coordinate Systems} - Semi-orthogonal
coordinate systems occur where one coordinate ($\xi_3$) is orthogonal
to the other two coordinates ($\xi_1$ and $\xi_2$)
\cite{Sava.geo.rwepub}.  In these cases the $m^{13}$ and $m^{23}$
components of the weighted metric tensor are identically zero, which
leads to the following extrapolation wavenumber,    
\beq 
\kz = 
i a_3
\pm   
\lb
  a_4^2 \ww^2 
- a_5^2 \kx^2 
- a_6^2 \ky^2
- a_7 \kx\ky 
+ i a_8 \kx 
+ i a_9 \ky  
- a^2_{10}  
\rb^{\frac{1}{2}},
\eeq
where, 
\beq
\mathbf{a} = 
\lb 
0 \;\;\;
0 \;\;\;
\frac{ n^3}{2 m^{33}} \;\;\;
\frac{ \ss }{ \sqrt{g^{33}} } \;\;\;
\sqrt{\frac{g^{11}}{g^{33}}} \;\;\;
\sqrt{\frac{g^{22}}{g^{33}}} \;\;\;
\frac{ 2 g^{12} }{ g^{33} }  \;\;\;
\frac{n^1}{m^{33}} \;\;\;
\frac{n^2}{m^{33}} \;\;\;
\frac{n^3}{m^{33}}
\rb^{\mathbf{T}}.
\eeq
which are identical to the coefficients recovered by \longcite{Sava.geo.rwepub}.
%
\par
%
{\bf 3D Kinematic Coordinate Systems} - Wave-equation migration
amplitudes are generally inexact in laterally variant media - even in
Cartesian-based systems.  Hence, one beneficial approximation that
reduces computational cost is to consider only the kinematic terms in
equation~\ref{eqn:kzfinal}.  This approximate generates the following
extrapolation wavenumber,   
\beq
\kz =
a_1 \kx +
a_2 \ky
\pm
\lb
  a_4^2 \ww^2 
- a_5^2 \kx^2 
- a_6^2 \ky^2
- a_7 \kx\ky 
- a_{10}^2  
\rb^{\frac{1}{2}},
\eeq
where, 
\beq
\mathbf{a} = 
\lb
-\frac{ g^{13} }{ g^{33} }\;\;\;
-\frac{ g^{23} }{ g^{33} } \;\;\;
0\;\;\;
\frac{\ss }{ \sqrt{g^{33}} }\;\;\;
\sqrt{\frac{g^{11} }{g^{33} } - \lp \frac{ g^{13} }{ g^{33} } \rp^2 }\;\;\;
\sqrt{\frac{g^{22}}{ g^{33} } - \lp \frac{ g^{23} }{ g^{33} } \rp^2}\;\;\;
\frac{ 2\,g^{12} }{ g^{33} } -\frac{ 2\, g^{13} g^{23}} {\lp
  g^{33}\rp^2 }\;\;\; 
0\;\;\;
0\;\;\;
\frac{n^3}{m^{33}}
\rb^{\mathbf{T}}.
\eeq
%
\par
%
{\bf 3D Kinematic Semi-orthogonal coordinate systems} - Combining the
two above restrictions yields the following extrapolation wavenumber,  
\beq
\kz =
\pm
\lb
  a_4^2 \ww^2 
- a_5^2 \kx^2 
- a_6^2 \ky^2
- a_7 \kx\ky 
- a_{10}^2  
\rb^{\frac{1}{2}},
\eeq
where,
\beq
\mathbf{a} = 
\lb
0 \;\;\;
0 \;\;\;
0 \;\;\;
\frac{ \ss }{ \sqrt{g^{33}} } \;\;\;
\sqrt{\frac{g^{11} }{g^{33} } } \;\;\;
\sqrt{\frac{g^{22}}{ g^{33} } }\;\;\;
\frac{ 2\, g^{12} }{ g^{33} }\ ;\;\;
0 \;\;\; 
0 \;\;\;
\frac{n^3}{m^{33}}
\rb^{\mathbf{T}}.
\eeq  
% 
\par
%
{\bf 2D Non-orthogonal coordinate systems} - Two-dimensional
situations are handled by identifying $\qy=0$.  All derivatives in the
associated metric tensor $g^{ij}$ with respect coordinate $\qy$ are
identically zero, and the resulting 2D non-orthogonal coordinate
system wavenumber is, 
%
\beq \label{eqn:2DNO}
\kz =
a_1 \kx+
i a_3 
\pm
\lb
  a_4^2 \ww^2 
- a_5^2 \kx^2 
+ a_8 i \kx 
- a_{10}^2  
\rb^{\frac{1}{2}},
\eeq
where,
\beq
\mathbf{a} = 
\lb
-\frac{ g^{13} }{ g^{33} }\;\;\;
0 \;\;\;
\frac{ n^3 }{ 2 m^{33} } \;\;\;
\frac{ \ss }{ \sqrt{g^{33}} } \;\;\;
\sqrt{ \frac{ g^{11} } {g^{33} } - \lp \frac{ g^{13} }{ g^{33} } \rp^2 } \;\;\;
0 \;\;\;
0 \;\;\;
\frac{n^1}{m^{33}} \;\;\;
0 \;\;\;
\frac{n^3}{m^{33}}
\rb^{\mathbf{T}}.
\eeq
%
\par
%
{\bf 2D Non-orthogonal Kinematic Coordinate Systems} -
Two-dimensional kinematic situations are handled through identity 
$\qy=0$.  Again, all derivatives in the associated metric tensor
$g^{ij}$ with respect coordinate $\qy$ are identically zero, and
the 2D non-orthogonal kinematic extrapolation wavenumber is 
\beq
\kz =
a_1 \kx
\pm
\lb
  a_4^2 \ww^2 
- a_5^2 \kx^2 
- a_{10}^2  
\rb^{\frac{1}{2}},
\eeq
where,
\beq
\mathbf{a} = 
\lb
-\frac{ g^{13} }{ g^{33} } \;\;\;
0 \;\;\;
0 \;\;\;
\frac{ \ss }{ \sqrt{g^{33}} } \;\;\;
\sqrt{ \frac{ g^{11} }{ g^{33} } - \lp \frac{ g^{13} }{ g^{33} } \rp^2 } \;\;\;
0 \;\;\; 
0 \;\;\;
0 \;\;\;
0 \;\;\;
\frac{n^3}{m^{33}}
\rb^{\mathbf{T}}.
\eeq
%
\par
%
{\bf 2D Orthogonal Coordinate Systems} - Two-dimensional situations
are handled with $\qy=\HH=0$. Accordingly, all derivatives in the
associated 
metric tensor $g^{ij}$ with respect coordinate $\qy$ are identically
zero, and the 2D non-orthogonal coordinate system is represented by 
\beq
\kz =
i a_3  
\pm
\lb
  a_4^2 \ww^2 
- a_5^2 \kx^2 
+ i a_8 \kx 
- a_{10}^2
\rb^{\frac{1}{2}},
\eeq
where,
\beq
\mathbf{a} = 
\lb
0 \;\;\;
0 \;\;\;
\frac{ n_3 }{ 2 m^{33} } \;\;\;
\frac{ \ss }{ \sqrt{g^{33}} } \;\;\;
\sqrt{\frac{g^{11} }{g^{33} } } \;\;\;
0 \;\;\; 
0 \;\;\;
\frac{n^1}{m^{33}} \;\;\;
0 \;\;\;
\frac{n^3}{m^{33}}
\rb^{\mathbf{T}}.
\eeq
%
\par
%
{\bf 2D Orthogonal Kinematic Coordinate Systems} - The above two
approximations can be combined to yield the following extrapolation
wavenumber for 2D orthogonal kinematic coordinate systems,  
\beq
\kz =
\pm
\lb
  a_4^2 \ww^2 
- a_5^2 \kx^2 
- a_{10}^2  
\rb^{\frac{1}{2}},
\eeq
where,
\beq
\mathbf{a} = 
\lb
0 \;\;\;
0 \;\;\;
0 \;\;\;
\frac{ \ss }{ \sqrt{g^{33}} } \;\;\;
\sqrt{\frac{ g^{11} }{g^{33} }} \;\;\;
0 \;\;\; 
0 \;\;\;
0 \;\;\;
0 \;\;\;
\frac{n^3}{m^{33}}
\rb^{\mathbf{T}}.
\eeq

\newpage
