\title{Riemannian Wavefield Extrapolation: Non-orthogonal coordinate
  systems} 
\email{jeff@sep.stanford.edu}
\author{Jeff Shragge}

\input{mainmacros}
\input{macros} 

\begin{abstract}
Riemanninan wavefield extrapolation (RWE) is used to model one-way
wave propagation on generalized coordinate meshes.  Previous RWE
implementations assume that coordinate systems are defined by either
orthogonal or semi-orthogonal geometry.  This restriction leads to
situations where coordinate meshes suffer from problematic bunching
and singularities. In this paper, I develop a procedure to avoid many
of these problems by posing wavefield extrapolation on smooth, but
generally non-orthogonal and singularity-free, coordinate meshes. The
resulting extrapolation operators include additional terms that
describe non-orthogonal propagation.  These extra degrees of
complexity, however, are offset by smoother coefficients that are more
accurately implemented in one-way extrapolation operators.  I validate
my theory of non-orthogonal propagation with two analytic coordinate
system examples, and I present a method for eliminating any remaining
singularities from coordinate systems. I demonstrate non-orthogonal
RWE by numerical calculation of 2D and 3D Green's functions for
cylindrical and near-spherical geometry.  Results from 2-D benchmark
testing suggest that the computational overhead associated with the
RWE approach is roughly 35\% greater than Cartesian-based
extrapolation. 
\end{abstract}            
%
\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 large propagation angles and cannot
propagate turning waves by design (though extensions including
phase-shift table lookup \cite[]{Hale1992} or two-pass migration
\cite[]{Zhang} address 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.  Examples of this approach include
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 local propagation angle, reducing the need for expensive
global extrapolation operators and enabling imaging of turning waves.
\cite{Sava.geo.rwepub} followed this approach in developing
Riemannian wavefield extrapolation (RWE), a theory of one-way
wavefield propagation for 3D numerically generated,
  semi-orthogonal meshes.  This formulation specifies the
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 the wave-propagation
direction is aligned with 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 \cite[]{Sava.geo.rwepub}. 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 geometric restriction leads to wave-equation
dispersion relationships that contain a number of mixed spatial and
wavenumber domain  terms (i.e. a simultaneous dependence on ${\bf x}$ and
${\bf k}$) that encode coordinate system geometry.  In most examples,
modeled Green's function estimates interpolated into the Cartesian
domain are highly accurate at large propagation angles; however,
accuracy is compromised in certain situations exhibiting
unfavorable characteristics such as extensive mesh
compression/extension or in the presence of singularities.  
%
\par
%
Semi-orthogonal geometry, though, can be an overly restrictive
assertion.  One problematic example is illustrated by the coordinate
system singularities observed in Figure~\ref{fig:Figure1}a.  These
situations arise wherever a mesh is generated from a rayfield
that includes a crossing set of rays, which generates spatial
singularities and singular Jacobians that lead to zero-division during
wavefield extrapolation.  Although ray-coordinate singularities 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
\cite{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 (see
Figure~\ref{fig:Figure1}b). 
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.   
%
\par
%
In this paper, I argue that generalizing RWE to non-orthogonal
geometries resolves many of the issues associated with semi-orthogonal
meshes, affords the use of smoother computational meshes and leads to a
procedure for removing all coordinate singularities.  The non-orthogonal
theory results in an extrapolation wavenumber that consists of
two additional mixed-domain terms relative to the wavenumber
associated with semi-orthogonal RWE.  I then use this wavenumber to
implement a one-way extrapolation operator in the usual sense of
wave-equation imaging.
%
\par
%
A goal of this paper is to 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 singularity-free computational meshes.
My development follows that of \cite{Sava.geo.rwepub}; however, I
recast the theory in a more compact notation for 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 extrapolation operator used in the examples is derived in
Appendix~B.  Two analytic 2D non-orthogonal coordinate system examples
are then provided to validate the theory.  The final sections detail a
procedure for generating singularity-free coordinate systems, present
2D and 3D Green's functions estimates modeled in cylindrical and
near-spherical coordinates, respectively, and I discuss the relative
computational cost and memory overhead of the RWE method.  
%
\plot{Figure1}{width=7in}{Illustration of problems with semi-orthogonal
  RWE computational grids.  a) Singularities in a coordinate system
  generated by ray-tracing that lead to zero-division in RWE. b) Grid
  bunching for a topographically conformal coordinate system that
  cause significant lateral variation in each extrapolation step.}  
%
\section{Acoustic wave equation in 3D generalized Riemannian spaces}
%
To specify the acoustic wave-equation in a 3D Riemannian space, we
must formulate 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. singularity-free and one-to-one).   In this
paper, I use a notation where a generalized coordinate system $\bxi
=\{ \xi_1,\xi_2,\xi_3 \}$ is 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 monochromatic wave-equation
for an acoustic wavefield, $\W$, in a generalized Riemannian space is,    
\beq \label{eqn:helm}
\dd_\bxi \W= - \ww^2\ss^2_\mathbf{\bxi} \W,
\eeq 
where $\dd_\bxi$ is the Laplacian operator applied in coordinates
$\bxi$, $\ww$ is frequency, and $\ss_\bxi$ is the propagation slowness.
\par
%
A correct formulation of 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  \cite[]{tensor},  
\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 indicies (i.e. $i,j=1,2,3$) is assumed
throughout this paper.  Substituting equation~\ref{eqn:laplac2} into
\ref{eqn:helm} leads to the Helmholtz equation appropriate for
propagating waves through a 3D Riemannian space \cite[]{Sava.geo.rwepub},   
\beq \label{eqn:we1}
\frac{1}{\gdetf} \de\qi \lp \mmm\ii\jj \done\W\qj \rp = - \ww^2\ss^2 \W.
\eeq
%
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 derivation of the 3D RWE acoustic wave equation deviates here from
that found in \cite{Sava.geo.rwepub}, who represent the metric tensor
with only four independent coefficients.  The remaining two
independent coefficients are explicitly set to zero in order to
satisfy the semi-orthogonal geometry restriction.  In this
development, I follow a more general approach that represents the metric
tensor with the full six independent coefficients.  This is the key
extension developed in this paper that leads to the greater
flexibility in coordinate system design discussed below. 
%
\par
%
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}$.  Note that the use of these dual
operators is strictly accurate only for the case of constant
coefficients.  Situations where $\ss, m^{ij}, \gdet$ or $n^{j}$
spatially vary lead to a simultaneous spatial and Fourier wavenumber
dependence.  However, as discussed below, I handle this through
multi-coefficient extensions of standard approximations.   
%
\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 coefficients
$m^{ij}=0$ for $i \neq j$).  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}$:
\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 an evaluation of indicies and
rearranging the result 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}
Wavenumber $\kz$ is then isolated through a complete-the-square
  transform
\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 
+ {\rm i} a_8 \,\kx 
+ {\rm i} a_9 \,\ky  
- a_{10}^2  
\rb^{\frac{1}{2}},
\eeq
where the non-stationary coefficients, $a_i$ in
equation~\ref{eqn:kzfinal}, are presented in the following vector
$\mathbf{a}$,  
\begin{eqnarray}
 \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 }
\;\;\;\;\; 
... \right. \nonumber \\
\, & \left. ... 
\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 }{2 m^{33} }
\right]^{\mathbf{T}}.
\end{eqnarray}
Note that the coefficients contain mixed $m^{ij}$ and $g^{ij}$
terms and globally positive 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 substituting $n^j=0$ and $g^{ij}=\delta^{ij}$ for the
coefficients of 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-{\bf 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$  \cite[]{Gazdag}, 
%
\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 $I$ can
be produced from the propagated wavefield through an imaging condition
\cite[]{Claerbout85}, 
%
\beq \label{eqn:ZOI}
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, require further approximations.  One straightforward
approach is a multi-coefficient split-step Fourier (SSF)
method \cite[]{GEO55-04-04100421,Sava.geo.rwepub}.  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-{\bf k_\bxi}$). Differences between the reference and true
parameters then form a correction term applied in the mixed
$\omega-\bxi$ domain. For non-orthogonal coordinate systems described
by equations~\ref{eqn:kzfinal}~and~\ref{eqn:ametric2}, I modify the 
SSF approach of \cite{Sava.geo.rwepub} as detailed in Appendix~B.
%
\par
%
The accuracy of the multi-coefficient SSF approach is
directly related to the degree to which coefficients in
equation~\ref{eqn:ametric2} vary at each 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}=59049$ 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} are
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? To address this issue, I calculate reference
coefficients using a multi-dimensional Lloyd's algorithm developed by
\cite{Clapp.Lloyds}.  This iterative procedure represents the
multi-dimensional histogram of the coefficients with the sparsest
number of points within a specified error tolerance.   For further
information and examples the reader is directed to
\cite{Clapp.Lloyds}. 
% 
\par
%
Second, numerous situations exists where some coefficients are zero or
otherwise negligible.  For example, the mixed-domain fields for a 3D weakly
non-orthogonal coordinate system within a kinematic approximation
(i.e. zeroed imaginary terms) is accurately approximated with four
coefficients.  In practice, I use the following relationship to
determine where non-orthogonal coefficients may be zeroed at any
extrapolation step: 
\beqa
\,&\,&  \hat{g}^{ij}= 
\left\{ 
\begin{array}{l} 
0,  \quad \quad \quad g^{ij} < 0.01\left|{\rm min}
\{g^{11},g^{22},g^{33} \} \right| \\ 
g^{ij}  \quad \quad \quad {\rm otherwise} 
			\end{array} \right. 
\eeqa
where the circumflex accent $\hat{g}^{ij}$ denotes approximation.
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}
%
This section presents numerical modeling examples that help validate
the non-orthogonal RWE theory.  I begin with the two basic 2D analytic
examples of  sheared Cartesian and polar ellipsoidal coordinates.   I
then present a method for generating singularity-free coordinate
meshes and illustrate this approach with 2D and 3D Green's function
modeling. 

\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 Figure~\ref{fig:Figure2}a).  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 Figure~\ref{fig:Figure2}a. 
%
\par
%
Figure~\ref{fig:Figure2}b shows the results of
extrapolating plane waves in a Cartesian coordinate system sheared
25$^{\circ}$ from vertical.  The background velocity model is 1500~m/s 
and the zero-offset data consist of four 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.  The propagation
creates explainable boundary artifacts.  Those on the left are caused
by the common edge effect of waves reflecting off the boundary at
non-normal incidence.  Hyperbolic diffractions on the right arise from
propagating truncated plane waves and are independent of the
coordinate system.  Mitigating these types of artifacts is not
difficult, though, because existing techniques in Cartesian wavefield
extrapolation craft still apply (e.g. cosine tapers). 
%
\sideplot{Figure2}{height=0.8\textwidth}{Sheared Cartesian coordinate
  system test. a) Coordinate system shear angle and velocity are
  $\theta=25^{\circ}$ and 1500 ms$^{-1}$, respectively.  b) 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 at depths $z=300, 600, 900,$ and
  $1200~m$.}
%
\subsection{Polar Ellipsoidal Coordinates}
%
A second example is an ellipsoidal polar coordinate system (see
Figure~\ref{fig:Figure3}a) appropriate for turning wave migration.  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 coordinate $\xi_1$ is the radius from the center focus, $\xi_3$
is polar angle, and $a=a(\xi_3)$ is a smooth function controlling
coordinate system ellipticity that 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
equations C-7 and C-8),     
\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 - {\rm i} \kx \xi_1 \lp \frac{a^2+2b^2-ac}{a^2} \rp }.
\eeq
The kinematic approximation of equation~\ref{eqn:2Dex5} (see equations
C-9 and C-10) is given by, 
\beq \label{eqn:2Dex6}
\hat{\kz} = \xi_1 \lb \frac{b}{a} \kx \pm \sqrt{ a^2 s^2 \omega^2 -  \kx^2 }\rb, 
\eeq
and further restricting to the orthogonal polar case where $a=1$ and
$b=0$ (see equations C-13 and C-14) yields,  
\beq \label{eqn:2Dex7} 
\hat{\kz} = \pm \xi_1 \sqrt{ s^2 \omega^2 -\kx^2 },
\eeq
which is examined in \cite{Nichols}.
%
\par
%
Figure~\ref{fig:Figure3} 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:Figure3} 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.  Figure~3a shows the
polar ellipsoidal coordinate system mesh overlying a linear
$v(z)=1500+0.35\,z$ m/s velocity function.  Figure~3b
presents the velocity model mapped into the RWE domain under the
transformations defined in equation~\ref{eqn:2DPE1}.     
%
\plot{Figure3}{width=\textwidth}{Ellipsoidal polar coordinate system test
  example. a) Vsselocity function $v(z)=1500+0.35\,z$ 
  overlain by a polar ellipsoidal coordinate system defined  by
  parameter $a=1+0.2\,\xi_3-0.05\,\xi_3^2$.  b) Velocity
  model mapped in the RWE domain. c) Imaged reflectors in
  RWE domain. d) RWE domain image mapped to a Cartesian 
  mesh.}  
%
The test data consist of ten plane waves defined on the surface
between 1000~m and 9000~m by ray parameter $p_x=-0.5$ s/km.  After
propagation to greater depths, the waves are no longer planar and
pass through a turning point before moving upward to the left (Figure
3c and 3d).  The wave tops, though, travel through slower material and
have not yet overturned.  One observation is that if propagating
wavepaths can be well represented by a single ellipticity parameter
$a=a(\xi_3)$, then a polar elliptical mesh could form the coordinate
basis for a plane-wave migration strategy that is potentially more
dynamic than tilted Cartesian coordinates.   
%
\section{Generating singularity-free coordinate meshes}
%
A computational mesh design challenge is finding a RWE coordinate
system that is fairly conformal to the wavefield propagation direction
yet unconditionally singularity-free.  Figure~\ref{fig:Figure4}a shows
a v(z) velocity model with three Gaussian anomaly inclusions overlain
by a ray-coordinate system calculated by Huygens' ray-front tracing 
\cite[]{Sava.hwt}.  These anomalies cause both mesh singularities to the
left and right of the model as well as a grid rarefaction directly
beneath the shot-point.    
%
\par
%
Figure~4b shows the single-valued isochrons of the
first-arrival Eikonal equation solution for the same shot-point
presented in the top panel. Note that isochrons generally conform to
the propagation direction and can be used to construct the
extrapolation steps of a RWE computational mesh.  The first step in
the mesh generation procedure is to extract the 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,ShraggeDMG}.  
%
\par
%
Figure~4c presents the corresponding singularity-free mesh.  The
grid is regularly-spaced on the outer isochron and has dimples at the
locations of removed singularities.   These discontinuities have been
reduced by applying a 
smoothing operator to the Eikonal equation solution before calculating
the mesh.  Importantly, coordinate smoothing usually does not
affect propagation accuracy because the coordinate system mesh forms
only the skeleton on which wavefield extrapolation occurs.
However, for meshes exhibiting rough and/or
discontinuous boundaries, even excessive local smoothing cannot
generate coefficients that are smooth enough to be accurately
represented with standard extrapolation techniques. 
%
\plot{Figure4}{width=6.5in}{Example of singularity-free mesh
  generation. a) 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. b)
  Velocity model overlain by isochrons of an Eikonal equation 
  solution for same shot-point. c) Singularity-free
  computational mesh generated by Eikonal mesh smoothing. }   
%
\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:Figure5} presents a slice through the SEG-EAGE salt
velocity model used for the test.  Importantly, the contrast between
the salt body and sediment velocities leads to complicated wavefield
propagation including triplication and multipathing.  Figure~5a 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 Figure~5b.     
%
\par
%
Figure~5c 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 Figure~5d.  The complicated wavefield to the left of the
shot point advances through the salt body and subsequently refracts
upward.  Note also the presence of wide-angle reflections from the 
top-salt/sediment interface.     
%
\plot{Figure5}{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. a) Salt model in physical
  space with an overlain ray-based mesh. b) Velocity model in
  the RWE domain.  c) Wavefield propagated in
  ray-coordinates through velocity model shown in b).
  d) Wavefield in c) interpolated back to
  Cartesian space.}   
%
Figure~\ref{fig:Figure6} presents a comparison test between two-way
finite-difference, RWE and Cartesian extrapolation.  The three
wavefields are fairly similar beneath and to the right of the
shot-point except for a $90^{\circ}$ phase-change associated with
differences between modeling the finite difference and Cartesian 
point-source (Figures~6a and 6c) versus the RWE plane-wave (Figure~6b).
However, significant differences are noted to the left of the
shot-point.  The upper two panels contain strong reflections from the
salt-sediment that are fairly well matched in location.
Cartesian-based extrapolation, though, propagates wavefields
laterally neither with the same accuracy nor upward at all. Hence,
this energy is absent from the propagating wavefield in the lower
panel.  
%
\par
%
Differences in the modeled amplitudes at and above the salt interface
in the upper two panels are attributed to differences between
finite difference and one-way wavefield extrapolation implementations.
Finite differences better models amplitudes in the presence of
velocity gradients in the propagation direction.  Thus, 
incident energy is more accurately partitioned at the top sediment-salt 
interface leading to the lower (and more correct) amplitudes of
subsalt multi-pathed arrivals.  The RWE wavefield underestimates the
reflection contribution and allows significantly more energy to be
transmitted into and through the salt body.  This modeling inaccuracy
leads to the more pronounced multipathing below the salt body and the 
more complicated wavefield behavior relative to the Cartesian
wavefield example.  Incorporating higher-order terms
\cite[]{ZhangFactor} into the RWE formalism, though, would likely
diminish the differences between the finite difference and RWE results.
%
\plot{Figure6}{height=0.75\textheight}{Comparison between three different
  extrapolation methods.  a) Two-way finite-difference. b) Riemannian 
  wavefield extrapolation. c) Cartesian wavefield extrapolation.}    

\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 (i.e. vertically downward).  This requires appropriately
discritizing the spherical shell making up an extrapolation step in
the presence of an ever diminishing solid angle toward the polar
regions.  One approach is 
to rotate the poles such that the small arc-length regions are at the
surface rather than vertically downward \cite[]{Schneider}.  A second
approach is to simply 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 spherical shell discritization  that
has a fairly uniform distribution in the Cartesian domain.  I use a
3D point source gridding technique based on the Winkel-Tripel
cartographic projection \cite[]{Winkel}.  Figure~\ref{fig:Figure7}
illustrates this projection for a constant spherical surface.  Note
that the mapping does not eliminate distortions in area, direction or
distance; rather, it minimizes the sum of all three.   
%
\par
%
A 3D point-source coordinate system can be
  constructed by using a set of concentric shells with spherical
  Winkel-Tripel gridding of increasing radii (or equally by the
Eikonal equation solution method described above).  This results in
the following analytic 3D point-source coordinate system: 
\beqa
\,&\,& x_1 = \xi_3 (\xi_1 / \pi + \, w\, \alpha {\rm cos} \, \xi_2
\, {\rm sin}(\xi_1 /2)) \nonumber \\ 
\,&\,& x_2 = \xi_3 ( \xi_2 + 2\, w \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
latitude parallels and longitude meridians, 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
Wavefield extrapolation is performed by stepping a point-source
wavefield outward in the radial direction $\xi_3$.   
%
\plot{Figure7}{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.  A 3D coordinate system is generated by
  constructing many concentric  Winkel-Tripel spherical shells of
  increasing radii.}  
%
\par 
%
Figure~\ref{fig:Figure8} presents an example of a wavefield
extrapolated on 3D point source coordinate mesh (Figure~8a) and, for
comparison, on a Cartesian grid (Figure~8b).    
%
\plot{Figure8}{height=0.9\textheight}{3D modeling example using the
  Elf-IFP-CGG synthetic velocity. a) RWE on Winkel-Tripel
  coordinate system. b) Cartesian wavefield extrapolation.} 
%
The 3D test example uses the Elf-IFP-CGG synthetic velocity model
generated 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 left creates a more complex geologic
setting.  Overlain are a number of impulse responses associated with
different travel times.  Note that the RWE approach models turning
waves to the left of the shot point, whereas Cartesian propagation
cannot by design.  However, because the coordinate system does not
cover the full velocity field, boundary artifacts are present within
the RWE propagation domain and need to be suppressed.  Overall,
though, the extrapolation result demonstrates the stability of
wavefield propagation on 3D Winkel-Tripel meshes and an improved
accuracy over one-way propagators on 3D Cartesian grids. 
%
\section{Implementation Costs}
%
The introduction of additional mixed-domain coefficients into the
dispersion relationship leads to both increased computation costs and
memory requirements.  To give an example of the cost overhead of the
RWE approach, relative to Cartesian, the algorithm was benchmarked on
the 2D computational grid (512x512 samples) used to generate
Figure~\ref{fig:Figure6}.  Tests were conducted on two codes that
differed only in the phase-shift and split-step Fourier subroutines.
The RWE code implemented the 2D non-orthogonal extrapolation operator
in equation~\ref{eqn:testtest}, while the Cartesian implementation
used the regular expression (i.e. $g_{33}=g_{11}=1$ and $g_{13}=n_3=0$
in equation~\ref{eqn:testtest}). 
%
\par
%
Table~1 presents the results of the benchmark test.
A total of 82 frequencies were propagated a total of 511 extrapolation
steps requiring 41 092 calls to the SSF operator.  The tests involved
112 996 calls to the phase-shift routine, or almost three per
extrapolation step as this number varied according to velocity model
complexity.  The most significant observation is that the RWE
algorithm is roughly 1.35 times slower than the equivalent Cartesian 
code.  Most of the overhead occurs in the phase-shift and SSF
subroutines that are roughly 2.25 and 2.0 times slower,
respectively. Whether these costs may be reduced by implementing
look-up tables remains an unresolved question.  An additional
computational overhead is the time required to calculate the
geometrical factors $a_i$ in equation~\ref{eqn:kzfinal}. This cost,
though, can usually be spread over the total number of shots for
stationary geometries.  Furthermore, the extra cost
  of non-orthogonal propagation, relative to that on semi-orthogonal
  mesh, is $<$5\% since this affects only phase-shift operation and
  occurs outside of the more costly square-root calculation. 
%
\par
%  
A second major implementation issue is the memory required to store
the non-stationary $a_i$ coefficients.  Holding each additional
coefficient in core requires allocating memory equivalent to a
velocity model, which can become the limiting issue for large 3D
models.  (For example, a 3D non-orthogonal grid requires an additional
20\% memory to store coefficients relative to a semi-orthogonal mesh.)
Unfortunately, the alternatives to allocating memory, recalculating
the $a_i$ coefficients locally each time or reading them from disk,
are inefficient.  An alternate approach is to consider analytic
coordinate systems \cite[]{ShraggeSEG2007} similar to those found in
Figures~\ref{fig:Figure2}~and~\ref{fig:Figure3}.  Their main advantage
is that they result in analytically defined extrapolation  operators
that avoid most problems associated with additional computational and
memory overhead costs.  Further discussion on this point, however, is
beyond the scope of this paper.  
% 
\section{Concluding Remarks}
%
This paper addresses existing issues with Riemannian wavefield
extrapolation theory by extending RWE to smoother, but non-orthogonal,
coordinate systems.   The paper demonstrates that one can generate an
acoustic wave equation in general 3D Riemannian spaces, and that the
corresponding extrapolation wavenumber decouples from the other
wavenumbers.  This wavenumber can then  be incorporated into a one-way
wavefield extrapolation operator appropriate for propagating
wavefields.  An approximate PSPI plus multi-coefficient SSF solution
of the one-way wavefield extrapolation operator is then derived.  A
method for generating computational meshes that are unconditionally
singularity-free is detailed, and used to generate examples that
illustrate wavefield propagation on non-orthogonal coordinate meshes
using RWE operators.  Accordingly, this opens up a range of imaging
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 Traynin for enlightening discussions.
I also acknowledge the reviewers of this manuscript reviewer and thank
them for their helpful comments.  The sponsors of the SEP consortium
provided funding for this work.  
%
\bibliographystyle{seg}
\bibliography{SEP,MISC,GEOTLE,allrefs}

\newpage

\append{}
%
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 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, $\left| \mathbf{g}\right|$,
\beq \label{eqn:det}
  |\mathbf{g}| =g_{11}g_{22}g_{33} - g^2_{12}g_{33} - g^2_{23}g_{11} -
  g^2_{13}g_{22} + 2g_{12}g_{13}g_{23}.ss 
\eeq
A weighted metric tensor, $m^{ij}=\gdetf \, g^{ij}$, is also used
throughout the paper.   
%
\newpage
%
\append{}
%
The extrapolation wavenumber defined in equations~\ref{eqn:kzfinal}
and \ref{eqn:ametric2} generally cannot be implemented exactly in the
Fourier domain due to a simultaneous spatial dependence (i.e. a
function of both $\xi_1$ and $\kx$). This can be addressed using an
multi-coefficient 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_\xi}$) domain phase-shift and a mixed
($\omega-\mathbf{\xi}$) domain split-step correction, respectively.  
%
\par
The phase-shift term is given by,   
\beq
\kz^{PS} = 
- b_1 \kx
- b_2 \ky 
+ {\rm i} b_3  
\pm
\lb
  b_4^2 \ww^2 
- b_5^2 \kx^2 
- b_6^2 \ky^2
- b_7 \kx\ky 
+ {\rm i} b_8 \,\kx 
+ {\rm i} b_9 \,\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 i  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} = {\rm 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
%
\append{}
%
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, kinematic approximations can be
made where one ignores all imaginary wavenumber components.  This
Appendix discusses these situations.  
%
\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 = 
{\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 
+ {\rm i} a_8 \kx 
+ {\rm 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}{2 m^{33}}
\rb^{\mathbf{T}}.
\eeq
which are identical to the coefficients recovered by
\cite{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
\hat{\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}{2 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
\hat{\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}{2 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+
{\rm i} a_3 
\pm
\lb
  a_4^2 \ww^2 
- a_5^2 \kx^2 
+ {\rm i} a_8 \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}} - \frac{m^{13}n^{3}}{(m^{33})^2} \;\;\;
0 \;\;\;
\frac{n^3}{2 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 \label{eqn:testtest}
\hat{\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}{2 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 orthogonal coordinate system is represented by 
\beq
\kz =
i a_3  
\pm
\lb
  a_4^2 \ww^2 
- a_5^2 \kx^2 
+ {\rm 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}{2 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
\hat{\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}{2 m^{33}}
\rb^{\mathbf{T}}.
\eeq

\newpage

\begin{table}
\begin{tabular}{|c|c|c|c|c|}
\hline
Extrapolation Type & Operation & Number of Calls & Total Time (s) &
Time per Frequency (s) \\ 
\hline
\hline  RWE &Frequency Loop  & 82     & 55.4 & 0.676 \\
\hline  RWE &Split-step Fourier& 41 092  &  4.2 & 0.051 \\
\hline  RWE &Phase-shift     & 112 996 & 17.1 & 0.209 \\
\hline  RWE &Interpolation   & 112 996 &  4.6 & 0.056 \\
\hline
\hline  Cart &Frequency Loop & 82     & 40.9 & 0.499 \\
\hline  Cart &Split-step Fourier& 41 092  &  2.1 & 0.026 \\
\hline  Cart &Phase-shift    & 112 996 &  7.6 & 0.093 \\
\hline  Cart &Interpolation  & 112 996 &  4.4 & 0.054 \\
\hline
\end{tabular}
\caption{Comparison of computational costs of the split-step Fourier
  and phase-shift subroutines for the RWE and equivalent Cartesian
  implementations.  Test results cited here were for calculating the
 2D example shown in Figure~\ref{fig:Figure6}. }\label{Efficiency}  
\end{table}
\newpage
