%% This file is automatically generated. Do not edit!
\documentclass[12pt,twoside,color]{geophysics}

\usepackage{amsmath,color}

\setfigdir{Fig}

\begin{document}
\chapter{RWE: Non-orthogonal coordinate systems} 
%\email{jeff@sep.stanford.edu}
%\author{Jeff Shragge}

% wavefield
\def\W{{U}}
% laplacian
\def\LAP{\Delta}
\def\dd {\nabla^2}
\def\bxi{{\bf \xi}}

% partial derivatives
\def\de#1{  \frac{\partial}{\partial #1} }
\def\dtwo#1#2{  \frac{\partial^2 #1}{\partial #2^2}   }
\def\done#1#2{  \frac{\partial   #1}{\partial #2  }   }
\def\eone#1#2{  \frac{\partial   }  {\partial #2  } #1}
\def\mtwo#1#2#3{\frac{\partial^2 #1}{\partial #2\partial #3}   }

\def\mmm#1#2{m^{#1#2}}
\def\ii{i}
\def\jj{j}
\def\z{z}
\def\x{x}
\def\y{y}

\def\qx{\xi_1}
\def\qy{\xi_2}
\def\qz{\xi_3}
\def\qi{\xi_i}
\def\qj{\xi_j}
\def\qk{\xi_k}
\def\kz{k_{\qz}}
\def\kx{k_{\qx}}
\def\ky{k_{\qy}}
\def\ki{k_{\xi_i}}
\def\kj{k_{\xi_j}}
\def\kk{k_{\xi_k}}
\def\xone{x_1}
\def\xtwo{x_2}
\def\xthr{x_3}

\def\DF#1#2{h_{#1#2}}
%%\def\hxx{\DF{\qx}{\qx}}
%%\def\hyy{\DF{\qy}{\qy}}
%%\def\hzz{\DF{\qz}{\qz}}
%%\def\hxy{\DF{\qx}{\qy}}

\def\EE{g_{11}}
\def\FF{g_{12}}
\def\GG{g_{22}}
\def\JJ{(\EE\GG-\FF^2)}
\def\HH{g_{13}}
\def\LL{g_{23}}
\def\PP{\Phi}
\def\AA{g_{33}}

\def\ww{\omega}
\def\ss{s}
\def\vv{v}

\def\xvec{\bf{x}}

\def\dxdex{ \frac{\partial x_k}{\partial \qx} }
\def\dxdey{ \frac{\partial x_k}{\partial \qy} }
\def\dxdez{ \frac{\partial x_k}{\partial \qz} }

% Partials
\def\ddex{ \frac{\partial}{\partial \qx} }
\def\ddey{ \frac{\partial}{\partial \qy} }
\def\ddez{ \frac{\partial}{\partial \qz} }
\def\ddexx{ \frac{\partial^2}{\partial \qx \partial \qx} }
\def\ddeyy{ \frac{\partial^2}{\partial \qy \partial \qy} }
\def\ddezz{ \frac{\partial^2}{\partial \qz \partiql \qz} }
\def\ddexy{ \frac{\partial^2}{\partial \qx \partial \qy} }
\def\ddexz{ \frac{\partial^2}{\partial \qx \partial \qz} }
\def\ddeyz{ \frac{\partial^2}{\partial \qy \partiql \qz} }

\def\gdet{ \left| \mathbf{g} \right| }
\def\gdetf{ \sqrt{\gdet} }


\def\hxxdef{\dxdex \dxdex}
\def\hyydef{\dxdey \dxdey}
\def\hzzdef{\dxdez \dxdez}
\def\hxydef{\dxdex \dxdey}
\def\hxzdef{\dxdex \dxdez}
\def\hyzdef{\dxdey \dxdez}

\def\c{c}
\def\cxx{\c_{\qx\qx}}
\def\cx {\c_{\qx   }}
\def\cyy{\c_{\qy\qy}}
\def\cy {\c_{\qy   }}
\def\czz{\c_{\qz\qz}}
\def\cz {\c_{\qz   }}
\def\cxy{\c_{\qx\qy}}

\def\ko{\; k_o}
\def\ofpo{ \lp \qz,\qx,\qy \rp}
\def\ofzo{ \lp \z ,\x ,\y  \rp}

% ray coordinates
\def\g{\gamma}
\def\t{\tau}
\def\hgg{\DF{\g}{\g}}
\def\htt{\DF{\t}{\t}}

\def\ofto{ \lp \t ,\g \rp}
\def\ctt{\c_{\t\t}}
\def\ct {\c_{\t  }}
\def\cgg{\c_{\g\g}}
\def\cg {\c_{\g  }}

\def\kt{k_\t}
\def\kg{k_\g}
\def\kn{k_{\nu}}
% discretization
\def\Wpo{\W^{\qx+1}_{\qz}}
\def\Woo{\W^{\qx}  _{\qz}}
\def\Wmo{\W^{\qx-1}_{\qz}}

\def\Wp {\W^{\qx+1}_{\qz+1}}
\def\Wo {\W^{\qx}  _{\qz+1}}
\def\Wm {\W^{\qx-1}_{\qz+1}}

% LISTS
\def\bit{\begin{itemize}}
\def\eit{  \end{itemize}}

\def\ben{\begin{enumerate}}
\def\een{  \end{enumerate}}

\def\itt{\item}

% SEPARATORS
\def\lp{\left(}
\def\rp{\right)}

\def\lb{\left[}
\def\rb{\right]}

%\def\ll{\left|}
%\def\rl{\right|}

\def\lc{\left\{}
\def\rc{\right\}}

% EQUATIONS
%\def\begin{equation}{ \begin{equation} }

%\def\end{equation}{  \end{equation}}
\def\req#1{equation~(\ref{eqn:#1})}
\def\reqs#1{equations~(\ref{eqn:#1})} 
\def\reqo#1{(\ref{eqn:#1})}

% VECTORS
\def\iv{\vec{i}} % x unit vector
\def\jv{\vec{j}} % y unit vector

\def\abs#1{\left| #1 \right|}	% vector modulus
\def\vdp#1#2{\lp#1\cdot#2\rp}   % dot product

% SIDE-BY-SIDE objects
\def\sidetwo#1#2{     \matrix{ {#1} & {#2}               \cr}}
\def\sidethree#1#2#3{ \matrix{ {#1} & {#2} & {#3}        \cr}}
\def\sidefour#1#2#3#4{\matrix{ {#1} & {#2} & {#3} & {#4} \cr}}

% OVER-UNDER objects
\def\overtwo#1#2{     \matrix{ {#1} \cr {#2}                   \cr}}
\def\overthree#1#2#3{ \matrix{ {#1} \cr {#2} \cr {#3}          \cr}}
\def\overfour#1#2#3#4{\matrix{ {#1} \cr {#2} \cr {#3} \cr {#4} \cr}}

%\begin{abstract}
This chapter shows how Riemannian wavefield extrapolation (RWE) can be used to model one-way wave propagation on generalized coordinate meshes.  Previous RWE implementations assume that coordinate systems are defined by either orthogonal or partial geometry.  This restriction leads to situations where coordinate meshes suffer from problematic bunching and singularities. I develop a procedure to avoid 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 coordinate singularities. I demonstrate the accuracy of the non-orthogonal RWE approach by numerical calculation of 2D Green's functions.  Testing results in 3D analytic coordinates are performed using an elliptic cylindrical coordinate system.  Results from 2-D benchmark testing suggest that the computational overhead associated with this wavefield extrapolation implementation on numerically generated coordinates using mixed space and Fourier-domain extrapolation operators is roughly 35\% greater than the equivalent Cartesian-based implementation.   Wavefield propagation in 2D and 3D analytic coordinate systems, though, does not require storing the coordinate geometry in memory, resulting in significantly more efficient extrapolation operators.  The main results from this chapter have been published as \cite{ShraggeGEONO}.
%
%\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 multi-pathing in the presence of lateral velocity variation, they are of limited large-angle accuracy and cannot propagate turning waves by design [though extensions like two-pass migration \cite[]{Zhiming,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}, coherent states \cite[]{Albertin-SEG01}, beam-waves  \cite[]{BDEtgen-SEG03}  or tilted Cartesian meshes \cite[]{SEG2004sgj}.  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 with 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, partially-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
%
\cite{Sava.geo.rwepub} initially implemented RWE 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: partially 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_x}$) 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
%
Partial orthogonal geometry, though, can be an overly restrictive assertion.  One problematic example is illustrated by the coordinate system singularities observed in Figure~\ref{fig:Fig1}a. 
%
\plot{Fig1}{width=6in}{Illustration of problems with partially 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 causes significant lateral variation in each extrapolation step.}  
%
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 partially orthogonal geometry is illustrated in \cite{ShraggeTOPO}, 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:Fig1}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 partially orthogonal geometry.   
%
\par
%
In this chapter, I argue that generalizing RWE to non-orthogonal geometries resolves many of the issues associated with partially 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 consisting of two additional mixed-domain terms relative to the wavenumber associated with the partially orthogonal RWE wavenumber.  I then use this wavenumber to implement a one-way extrapolation operator in the usual sense of wave-equation imaging.
%
%\par
%
To these ends, I 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 that allows for a closer analytic connection of the generalized computation geometry with the underlying Cartesian grid. I begin 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 to compute this chapter's 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 elliptical-cylindrical coordinates, respectively, and discuss the relative computational cost and memory overhead of the RWE method.  
%
\section{Acoustic wave equation in 3D generalized Riemannian spaces}
%
To specify the acoustic wave-equation in a 3D Riemannian space, the physics of wave-propagation must be formulated 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).   I use a notation where a generalized coordinate system $\boldsymbol{\xi}=\{ \xi_1,\xi_2,\xi_3 \}$ is related to a Cartesian grid $\boldsymbol{x}=\{x_1,x_2,x_3 \}$ through transformation $x_j (\xi_k) = f_j, j,k=1,2,3$.  Provided these conditions are met, the monochromatic wave-equation for an acoustic wavefield, $\W$, in a generalized Riemannian space is,    
\begin{equation} \label{eqn:helm}
\dd_\bxi \W= - \ww^2\ss^2 \W,
\end{equation} 
where $\dd_\bxi$ is the Laplacian operator applied in coordinates $\bxi$, $\ww$ is frequency, and $\ss$ is the propagation slowness.
%
\par
%
A correct formulation of the wave-equation in the $\boldsymbol{\xi}$-coordinate system requires that Laplacian operator $\dd_\bxi$ be specified by differential geometry relationships.  (Appendix~A provides an overview of necessary differential geometry theory.)  The Laplacian operator in generalized coordinates is  \cite[]{tensor},  
\begin{equation} \label{eqn:laplac2}
  \dd_\bxi \W = 
  \frac{1}{\sqrt{|\mathbf{g}|}}\,
    \frac{\partial}{\partial \xi_j}\,
    \left(\,m^{jk}\,\frac{\partial \W}{\partial \xi_k}\right), \quad
    j,k=1,2,3, 
\end{equation}
where $|\mathbf{g}|$ is the metric tensor discriminant with elements $g^{jk}$, and $m^{jk}=\sqrt{|\mathbf{g}|} g^{jk}$ is weighted metric tensor element that enables a more compact notation.  Unless otherwise stated, summation over all repeated indicies (i.e. $j,k=1,2,3$) is assumed throughout.  Note also that convention where superscript indicies (e.g. $m^{jk}$) on matricies indicates a contravariant tensor, while subscript indicies  (e.g. $g_{jk}$) on matricies denotes covariant tensors \cite[]{tensor}.  
%The relationship is further discussed in Appendix~A.
%
\par
%
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},   
\begin{equation} \label{eqn:we1}
\frac{1}{\gdetf} \de\qj \lp m^{jk} \done\W\qk \rp = - \ww^2\ss^2 \W.
\end{equation} 
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,    
\begin{equation} \label{eqn:helm1} 
\frac{\partial m^{jk}}{\partial_{\xi_j}} \done\W\qk + m^{jk} \mtwo\W\qj\qk
= - \gdetf \ww^2 \ss^2 \W.
\end{equation}
%
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 partially orthogonality 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 herein 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,  
\begin{equation}
n^j=\frac{\partial{m^{jk}}}{\partial \xi_k} = 
\done{m^{j1}}\qx + 
\done{m^{j2}}\qy +
\done{m^{j3}}\qz.
\end{equation}
Scalar 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, 
\begin{equation} \label{eqn:helm2}
n^j \done\W\qj + m^{jk} \mtwo\W\qj\qk
= -\gdetf \ww^2 \ss^2 \W.
\end{equation}
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}, 
\begin{equation} \label{eqn:disp}
\lp m^{jk} \kk - i n^j \rp \kj = \gdetf \ww^2\ss^2,
\end{equation}
where $-i \kj$ is the Fourier domain dual of differential operator $\frac{\partial}{\partial \xi_j}$.  Note that the use of these dual operators is strictly accurate only for the case of constant coefficients.  Situations where $\ss, m^{jk}, \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^{jk}$ in the first term, $m^{jk}\kj\kk$, is a measure of the dot product between wavenumber vectors in the $\kj$ and $\kk$ directions (i.e. orthogonal wavenumbers will have coefficients $m^{jk}=0$ for $j \neq k$).  Fields $n^j$ in the second term, $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^{jk}=\delta^{jk}$:
\begin{equation}
\kj \kj= \kx^2 + \ky^2 + \kz^2 = \ww^2 s^2.
\end{equation} 
%
\subsection{Extrapolation wavenumber isolation}
%
Specifying a one-way extrapolation operator requires isolating one of the wavenumbers in equation~\ref{eqn:disp}.  I associate the extrapolation direction with coordinate $\xi_3$.  Expanding equation~\ref{eqn:disp} and evaluating a complete-the-square transform yields an expression for the wavenumber $\kz$   
%\begin{small}
%\begin{equation} \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.
%\end{equation}
%\end{small} 
%Wavenumber $\kz$ is isolated through a complete-the-square transform
\begin{equation} \label{eqn:kzfinal}
\kz =
- a_1 \kx
- a_2 \ky 
+ 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_{10}^2  
\rb^{\frac{1}{2}},
\end{equation}
where the non-stationary coefficients, $a_j$ in equation~\ref{eqn:kzfinal}, are presented in vector $\mathbf{a}$,  
\begin{eqnarray}
 \label{eqn:ametric2}
\mathbf{a}
= &
\left[
\frac{m^{13}}{ m^{33}} \;\;\;\;\;
\frac{m^{23}}{ m^{33} }\;\;\;\;\;
\frac{ n^3 } {2 m^{33} } \;\;\;\;\;
\sqrt{\frac{ \gdet \ss }{ m^{33} } } \;\;\;\;\;\;
\sqrt{ \frac{m^{11}}{ m^{33}} - \lp \frac{m^{13}}{m^{33}}  \rp^2 }
\;\;\;\;\; 
\sqrt{\frac{m^{22}}{ m^{33} } - \lp \frac{ m^{23}}{m^{33}}  \rp^2 }
\;\;\;\;\; 
... \right. \nonumber \\
\, & \left. ... 
\lb \frac{2\,m^{12}}{m^{33}} - \frac{2\, m^{13}\,m^{23}}{ \lp m^{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 globally positive terms $a_4, a_5, a_6$ and $a_{10}$ that are squared.  
%
\par
%
The special Cartesian case is again recovered from the two equations above by substituting $n^j=0$ and $m^{jk}=\delta^{jk}$ for the coefficients of equation~\ref{eqn:ametric2} 
\begin{equation}
\kz=\left[ s^2\ww^2-\kx^2-\ky^2  \right]^{\frac{1}{2}}.
\end{equation}  
%
\par
%
The dispersion relationship specified by equations~\ref{eqn:kzfinal} and \ref{eqn:ametric2} contains ten coefficients that represent mixed-domain fields.  For situations where all ten coefficients are constant, for example in Cartesian wavefield extrapolation through homogeneous media, a constant-coefficient Fourier-domain ($\omega-{\bf k_{\boldsymbol{\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}, 
%
\begin{equation} \label{eqn:bulk}
\W(\xi_3+\Delta \xi_3, \kx,\ky | \omega) = 
\W(\xi_3, \kx, \ky | \omega) {\rm e}^{ i \kz \Delta \xi_3}. 
\end{equation}
%
If $\W$ represents a post-stack wavefield, an image $I(\boldsymbol{\xi})$ can be produced from the propagated wavefield by evaluating an imaging condition \cite[]{Claerbout85}, 
%
\begin{equation} \label{eqn:ZOI}
I(\xi_3,\xi_1,\xi_2) = \sum_{\omega} \W(\xi_3,\xi_1,\xi_2|\omega). 
\end{equation}
%
\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_{\boldsymbol{\xi}}}$). Differences between the reference and true parameters then form a correction term applied in the mixed $\omega-\boldsymbol{\xi}$ 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}=59\,049$ separate wavefield extrapolations.)  However, three factors combine to greatly reduce the total number of required reference coefficient sets.  
%
\par
%
First, the $a_j$ coefficients in equation~\ref{eqn:kzfinal} are highly correlated because they are composed of similar metric tensor  elements $m^{jk}$.  Thus, the central issue is how accurately can we characterize these vector coefficient fields.  Concidentally, 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 \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.  One approximation is to set all terms containing imaginary numbers to zero, which largely affects only wavefield amplitudes.  This kinematic approximation can lead to a mixed-domain fields for a 3D weakly non-orthogonal coordinate system that contains only four coefficients.  A second approximation is to zero coefficients that are relatively small.  For example, in practice I use the following relationship to determine where non-orthogonal coefficients may be zeroed at any extrapolation step: 
\begin{eqnarray}
\,&\,&  \hat{m}^{jk}= 
\left\{ 
\begin{array}{l} 
0,  \quad \quad \quad m^{jk} < 0.01\left|{\rm min}
\{m^{11},m^{22},m^{33} \} \right| \\ 
m^{jk}  \quad \quad \quad {\rm otherwise} 
			\end{array} \right. 
\end{eqnarray}
where the circumflex accent $\hat{g}^{jk}$ 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 above 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{Sheared 2D Cartesian coordinates}
%
An instructive analytic coordinate system to examine is a sheared 2D
Cartesian grid formed by a uniform shearing action on a Cartesian
mesh (see Figure~\ref{fig:Fig2}a). 
%
\plot{Fig2}{width=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 four 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.}
%
This coordinate system is uniquely specified by one additional degree of freedom and is related to an underlying Cartesian mesh through  
%
\begin{equation} \label{eqn:2Dred}
  \left[\begin{array}{c}
  	\xone \\
  	\xthr
  	\end{array}\right] 
  	= 
  \left[\begin{array}{cc}
      1 & \rm{sin}\, \theta      \\
      0 & \rm{cos}\, \theta
  \end{array}\right]
  \left[\begin{array}{c}
    \qx \\
    \qz \\
  \end{array}\right],
\end{equation}
where $\theta$ is the shear angle of the coordinate system ($\theta=0^{\circ}$ is Cartesian).  The metric tensor of this transformation is, 
\begin{equation} \label{eqn:2Dmet}
\left[g_{jk}\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{sin}\, \theta \\
    \rm{sin}\,\theta & 1  
\end{array}\right],
\end{equation}
and has a discriminant $\gdet = \rm{cos}^2\,\theta$ and a weighted associated
metric tensor $m^{jk}$ given by, 
\begin{equation} \label{eqn:2Dmet6}
\left[m^{jk}\right]
=
\left[\begin{array}{cc}
    1  & - \rm{sin}\,\theta \\
    -\rm{sin}\,\theta & 1
\end{array}\right].
\end{equation}
%
Because the tensor in equation~\ref{eqn:2Dmet6} is coordinate invariant, equation~\ref{eqn:helm2} simplifies to,   
\begin{equation} \label{eqn:laplac}
   m^{jk} \,\frac{\partial^2 \W}{\partial \qj \partial \qk}
  = - \gdetf \ww^2\ss^2\W,
\end{equation}
which generates the following dispersion relation,
\begin{equation}
m^{jk}\kj\kk = \gdetf \ww^2 \ss^2.
\end{equation}
Expanding out these terms leads to an expression for wavenumber $\kz$, 
\begin{equation} \label{eqn:kz2Dex}
\kz =
-\frac{ m^{13} }{ m^{33} } \kx
\pm
\sqrt{\frac{ \gdet \ss^2\ww^2}{m^{33}}- \lp\frac{m^{11}}{m^{33}}-\lp
  \frac{m^{13}}{m^{33}} \rp^2\rp \kx^2}. 
\end{equation}
Substituting the values of the associated metric tensor in equation~\ref{eqn:2Dmet6} into equation~\ref{eqn:kz2Dex} yields, 
\begin{equation}
\kz =
\rm{sin} \, \theta \;\kx
\pm \rm{cos} \,\theta
\sqrt{\ss^2\ww^2 - \kx^2}, 
\end{equation} 
which is appropriate for performing RWE on the sheared 2D Cartesian coordinate system shown in Figure~\ref{fig:Fig2}a. 
%
\par
%
Figure~\ref{fig:Fig2}b shows the results of extrapolating plane waves in a Cartesian coordinate system sheared at $\theta=25^{\circ}$.  The background velocity model is 1500~ms$^{-1}$ 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 $x_3$=300, 600, 900, and 1200~m.  The propagation generates 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). 
%
\subsection{Stretched Polar coordinates}
%
A second example is a stretched polar coordinate system (see Figure~\ref{fig:Fig3}a) appropriate for migrating with turning waves. 
%
\plot{Fig3}{width=\textwidth}{Stretched polar coordinate system test
  example. a) Velocity function $v(x_3)=1500+0.35\,x_3$ overlain by a stretched polar 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.}
%
 A stretched polar coordinate system is specified by  
\begin{equation} \label{eqn:2DPE1a}
  \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],
\end{equation}
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 stretch that has curvature parameters
$b=\frac{{\rm d} a}{ {\rm d} \xi_3}$ and $c=\frac{ {\rm d}^2
  a}{ {\rm d} \xi_3^2}$.  The metric tensor $g_{jk}$ for the stretched polar
 coordinate system defined in equation~\ref{eqn:2DPE1a} is, 
\begin{equation} \label{eqn:2DPE2}
\left[g_{jk}\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],
\end{equation}
and has a metric discriminant given by $\gdet = a^4 \xi_1^2$. The weighted associated metric tensor is given by,  
\begin{equation} 
\quad \label{eqn:2DPE4}
\left[m^{jk}\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].
\end{equation}
%	
Tensor $m^{jk}$ is 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^j$ 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),     
\begin{equation} \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 }.
\end{equation}
The kinematic approximation of equation~\ref{eqn:2Dex5} (see equations
C.9 and C.10) is 
\begin{equation} \label{eqn:2Dex6}
\hat{\kz} = \xi_1 \lb \frac{b}{a} \kx \pm \sqrt{ a^2 s^2 \omega^2 -  \kx^2 }\rb, 
\end{equation}
and further restricting to the orthogonal polar case that is a circular geometry, where $a=1$ and
$b=0$ (see equations C.13 and C.14), yields,  
\begin{equation} \label{eqn:2Dex7} 
\hat{\kz} = \pm \xi_1 \sqrt{ s^2 \omega^2 -\kx^2 },
\end{equation}
which is examined in \cite{Nichols}.
%
\par
%
Figure~\ref{fig:Fig3} shows a wavefield extrapolation example for an polar-ellipsoidal coordinate system in equation~\ref{eqn:2DPE1a} defined by stretch parameter $a(\xi_3)=1+0.2\,\xi_3-0.05\,\xi_3^2$.  The upper and lower panels of Figure~\ref{fig:Fig3} 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~\ref{fig:Fig3}a shows the stretched polar coordinate system mesh overlying a linear $v(x_3)=1500+0.35x_3$ ms$^{-1}$ velocity function.  Figure~\ref{fig:Fig3}b presents the velocity model mapped into the RWE domain under the transformations defined in equation~\ref{eqn:2DPE1a}.     
%

%
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$ skm$^{-1}$.  The waves, propagated to greater depths, are no longer planar and pass through a turning point before moving upward to the left (panels \ref{fig:Fig3}c-d).  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 stretch parameter $a=a(\xi_3)$, then a stretched polar mesh could form an effective coordinate basis for plane-wave migration.   
%
\section{Generating singularity-free coordinate meshes}
%
A computational mesh design challenge is generating a RWE coordinate system fairly conformal to the wavefield propagation direction yet unconditionally singularity-free.  Panel~\ref{fig:Fig4}a shows a $v(x_3)$ 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. 
%
  \plot{Fig4}{width=\textwidth}{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, but weakly non-orthogonal, computational mesh generated by Eikonal mesh smoothing. }     
%
\par
%
Panel~\ref{fig:Fig4}b 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
%
Panel~\ref{fig:Fig4}c presents the corresponding singularity-free, but weakly non-orthogonal mesh.  The grid is regularly spaced on the outer isochron and has dimples at the locations of the 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. 
%
\subsection{2D Green's function generation}
%
The third test uses RWE to model 2D Green's functions on coordinate systems constructed by the smoothed Eikonal meshing approach. Figure~\ref{fig:Fig5} presents a slice through the SEG-EAGE salt velocity model used for the test.  
%
 \plot{Fig5}{width=\textwidth}{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-coordinate mesh. b) Velocity model in the transform domain.  c) Wavefield propagated in ray coordinates through velocity model shown in b).  d) Wavefield in c) interpolated back to Cartesian space.} 
%
Importantly, the contrast between the salt body and sediment velocities leads to complex wavefield propagation including triplication and multi-pathing.  Panel~\ref{fig:Fig5}a shows the velocity model with an overlain coordinate system generated by the smoothed Eikonal meshing procedure.  The velocity model in the RWE domain is illustrated in panel~\ref{fig:Fig5}b.     
%
\par
%
Panel~\ref{fig:Fig5}c 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 panel~\ref{fig:Fig5}d.  The complex 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.     
%
  
%
Figure~\ref{fig:Fig6} presents a comparison test between two-way finite-difference modeling, RWE and Cartesian extrapolation. 
%
 \plot{Fig6}{width=\textwidth}{Comparison between three different
  extrapolation methods.  a) Two-way finite-difference modeling. b) Riemannian 
  wavefield extrapolation. c) Cartesian wavefield 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 in panels~\ref{fig:Fig6}a and \ref{fig:Fig6}c versus the RWE plane-wave in panel~\ref{fig:Fig6}b.  [See \cite{Hudson} for a complete explanation of the phase differences associated with line- versus point-source modeling].  These phase-changes were also observed in the polar coordinate examples of \cite{Nichols}.  However, significant differences are noted to the left of the shot-point.  Panels~\ref{fig:Fig6}a-b 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 the finite-difference modeling and one-way wavefield extrapolation implementations. Finite-difference propagation 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 sub-salt multi-pathed arrivals.  The RWE approach 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 multi-pathing 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 modeling and RWE results.
%
\subsection{3D analytic coordinates}
%
There are a number of analytic coordinate systems that can be used in various 3D RWE applications (i.e. cylindrical polar, elliptic cylindrical, ellipsoidal coordinates).   The optimal choice of coordinate system generally depends on the acquisition geometry of the wavefield to be propagated.  For example, line-sources are more conformal in either cylindrical polar or elliptical cylindrical coordinates, while point sources are better matched to ellipsoidal coordinates.
%
\par
%
To provide a 3D RWE example, I examine wavefield extrapolation in elliptical-cylindrical coordinates (ECC) (see figure~\ref{fig:TECgeom}).  This computational mesh forms the basis of the inline delayed-shot migration algorithm detailed in Chapter~5.    
%One advantage of ECC-based wavefield extrapolation is that the breadth of the coordinate system at the surface allows multiple streamers of a full sail line to be positioned directly on a single mesh.  Hence, this geometry is applicable to both narrow- and wide-azimuth acquisition.  A second advantage is that one direction of high-angle propagation can be handled by coordinate system tilting (discussed in Chapter~5), while the other is naturally handled by the ellipticity of the mesh.  
%
%\par
%
I set up the extrapolation geometry of the elliptical-cylindrical mesh as follows:
\begin{itemize} 
\item $\xi_3 \in [0,\infty]$ is the extrapolation direction, where surfaces of constant $\xi_3$ form concentric elliptic cylinders, shown in Figure~\ref{fig:TECgeom}a;
\item $\xi_2 \in [0,2\pi)$ is the cross-line direction, where surfaces of constant $\xi_2$ are folded hyperbolic planes, shown in Figure~\ref{fig:TECgeom}b; and 
\item $\xi_1 \in [-\infty,\infty]$ is the inline direction, where surfaces of constant $\xi_1$ are 2D elliptic coordinate meshes, shown in Figure~\ref{fig:TECgeom}c.
\end{itemize}
%
\plot{TECgeom}{height=0.8\textheight}{Constant surfaces of the elliptic-cylindrical coordinate system   Cartesian coordinate axes are given to  a) Five constant $\xi_3$ surfaces forming confocal elliptic cylindrical shells.  This represents the direction of extrapolation direction.  b) Five constant $\xi_2$ surfaces representing folded hyperbolic planes.  c) Five constant $\xi_1$ surfaces representing 2D elliptic meshes.}
%
\par%
The mapping relationship between the two coordinate systems is \cite{Arfken}
\begin{equation}
\label{eqn:2DPE1}
  \left[\begin{array}{c}
  	x_1 \\
  	x_2 \\
  	x_3 
  	\end{array}\right] 
  	= 
  \left[\begin{array}{c}
\xi_1 \\
a \,\rm{cosh} \, \xi_3 \, \rm{cos}\,\xi_2  \\ 
a \,\rm{sinh} \, \xi_3 \, \rm{sin}\,\xi_2 \\  
  \end{array}\right],  
\end{equation}
where parameter $a$ controls the coordinate system breadth.  
%Panels~\ref{fig:TEC}a and %\ref{fig:TEC}b show the TEC coordinate system at $0^\circ$ and  $25^\circ$ tilt angles, respectively. 
%
%\plot{TEC}{width=\textwidth}{Four extrapolation steps in $\xi_3$ of an TEC coordinate system, where the $\xi_1$  and $\xi_2$ coordinate axes are oriented in the inline and crossline directions, respectively. a) 0$^\circ$ tilt angle.  b) 25$^\circ$ tilt angle.}   
%
%\par
%
%A metric tensor $g_{jk}$ can be specified from the mapping relationship given in equations~\ref{eqn:2DPE1}:   
%\begin{equation}
%\left[g_{jk}\right] =
%\left[\begin{array}{ccc}
%    1 & 0 & 0\\
%    0 & A^2 & 0 \\
%    0 & 0 & A^2 \\
%\end{array}\right],
%\end{equation}
%where $A=a \sqrt{ {\rm sinh}^2 \,\xi_3 + {\rm sin}^2 \,\xi_2 }$.  The determinant of the metric %tensor is: $|\boldsymbol{g}| = A^4$, and the corresponding inverse weighted metric tensor, %$m^{jk}$, is given by:  
%\begin{equation}
%\left[m^{jk}\right]
%=
%\left[\begin{array}{ccc}
%    A^2 & 0 & 0 \\
%    0 & 1 & 0 \\
%    0 & 0 & 1 \\
%\end{array}\right].
%\end{equation}
%Note that even though the elliptic cylindrical coordinate metric tensor varies spatially, the local curvature parameters ($n^j=\frac{\partial   m^{jk}}{\partial\xi_k}$) remain constant: $n^1=n^2=n^3=0$.   The corresponding extrapolation wavenumber, $k_{\xi_3}$, can be generated by inputting tensor $m^{jk}$ and fields $n^j$ into the general wavenumber expression for 3D non-orthogonal coordinate systems in equation \ref{eqn:kzfinal} 
%
%\begin{equation} \label{eqn:EWN}
%k_{\xi_3} = \pm \sqrt{ A^2 s^2 \omega^2 - A^2 k_{\xi_1}^2 -  k_{\xi_2}^2 }.
%\end{equation}
%
\par
%
I conducted elliptic coordinate imaging test results using a 3D Gulf of Mexico velocity model consisting of typical $v(x_3)$ velocity gradients and a salt body.  The source wavefield consisted of five line sources of zero inline dip at times 1.0, 1.5, 2.0, 2.5 and 3.0s. 
%(I use line sources because their impulse responses are well matched to the geometry of the elliptical-cylindrical coordinate system.)  
Figure~\ref{fig:NewPoint} presents the ECC and CC imaging test results superposed over the a velocity model with the inline (crossline) axis in the front (side) panel.  The ECC image (panel~\ref{fig:NewPoint}a) has energy that propagates at steep angles through the salt body in the crossline section, and even overturns to the right-hand side.  The CC image (panel~\ref{fig:NewPoint}b) contains similar propagation in most locations; however, the high-angle and turning-wave features present in the ECC wavefield are absent because of one-way propagation limitations.
%
\plot{NewPoint}{width=\textwidth,height=0.8\textheight}{Green's function computation test in the 3D elliptical-cylindrical analytic coordinate system for a set of five flat plane waves at time 1.0 1.5, 2.0, 2.5, and 3.0 s.  Top: Elliptical-cylindrical coordinate imaging result. Bottom: Cartesian coordinate imaging result.}
%
%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 discretizing 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 discretization leads to highly variable sampling in Cartesian space.  
%
\par
%
%Another solution is to choose a spherical shell discretization 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:Fig7} illustrates this projection for a constant spherical surface.  This mapping does not eliminate distortions in area, direction or distance; rather, it minimizes the sum of all three.   
%
%  \sideplot{Fig7}{width=\textwidth}{The Winkel-Tripel projection of unit hemisphere.  Cells on the unit hemisphere surface have nearly consistent area indicating fairly uniform sampling.  Not all parts of the unit hemisphere are sampled.  A 3D coordinate system is generated by concatenating concentric Winkel-Tripel spherical shells of  increasing radii.}  
%
%\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 smoothed Eikonal meshing method described above).  This results in the following analytic 3D point-source coordinate system: 
%\begin{eqnarray}
%\,&\,& 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
%\end{eqnarray}
%where $\xi_3$ is in the radial direction, $\xi_1$ and $\xi_2$ are latitude parallels and longitude %meridians, and
%\begin{eqnarray}
%\,&\,& \alpha = {\rm cos}^{-1} \, 
%\left[
%{\rm cos} \, \xi_2 \, {\rm cos} (\xi_1/2)
%\right] ,\nonumber \\ 
%\,&\,& w = \left\{ \begin{array}{l}  
%			0, \quad \quad \quad {\rm sin} \,\alpha = 0, \\
%			({\rm sin}\, \alpha)^{-1}, \quad {\rm otherwise}. 
%			\end{array} \right. 
%\end{eqnarray}
%Extrapolation is performed by stepping a point-source wavefield outward in the radial direction %$\xi_3$.   
%
%\par 
%
%Figure~\ref{fig:Fig8} presents an example of a wavefield extrapolated on 3D point source coordinate %mesh (panel~\ref{fig:Fig8}a) and on a Cartesian grid (panel~\ref{fig:Fig8}b).    
%
%\plot{Fig8}{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(x_3)$ 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 extrapolates 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 that need to be suppressed.  Overall, 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 generally non-orthogonal Riemannian 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, I benchmarked the algorithm on the 2D computational grid (512x512 samples) used to generate Figure~\ref{fig:Fig6}.  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. $m^{33}=m^{11}=1$ and $m^{13}=n^3=0$ in equation~\ref{eqn:testtest}). 
%
\par
%
Table~\ref{tab:Efficiency} presents the results of the benchmark testing. A total of 82 frequencies were propagated a total of 511 extrapolation steps requiring 41 092 calls to the SSF operator.  
%
\begin{table}
\begin{center}
\begin{tabular}{|c|c|c|c|c|}
\hline
Extrapolation & Operation & Number of &    Total  & Time per      \\  
Type          & \,        & Calls     &  Time (s) & Call (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 RWE and equivalent Cartesian implementations.  Results computed for the 2D example shown in Figure~\ref{fig:Fig6}. }\label{tab:Efficiency}  
\end{center}
\end{table}
%
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.35x 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_j$ 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 partially 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_j$ coefficients.  Holding each additional coefficient in core requires allocating memory equivalent to that of the 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 partially orthogonal mesh.) Unfortunately, the alternatives to allocating memory, recalculating the $a_j$ coefficients locally each time or reading them from disk, are inefficient.  
%
\subsection{Analytical Coordinates}
%
An alternate approach, discussed in greater detail in Chapters~3-5, is to consider analytic coordinate systems similar to those illustrated in Figures~\ref{fig:Fig2}, \ref{fig:Fig3} and \ref{fig:NewPoint}.  The main advantage of these coordinate systems that they result in analytically defined extrapolation operators that avoid most problems associated with additional computational and memory overhead costs.  The key reason is that  all geometric factors can be calculated locally knowing just the local coordinate position using just a few floating point operations.  To illustrate this point, I computed four different 3D analytical coordinate migration volumes: one Cartesian and one elliptic cylindrical coordinate image where rays are pre-computed and stored in memory, and one Cartesian and elliptical cylindrical coordinate image where the geometry is computed at run time.  The migration geometries were all of the mesh size dimensions: 400x300x500 grid points. 
%
\par
%
Table~\ref{tab:RayDifference} summarizes the results of the analytic coordinate test migration runs.   Overall, the Cartesian coordinate migrations using geometry stored in memory were the fastest.  The Cartesian test run using a migration code where the geometry was explicitly recalculated each time ran about 2\% slower, but used only 78\% of the memory.  The elliptical-cylindrical geometry migrations overall ran a little bit slower (2\% and 6\% for stored and recalculated geometries, respectively).  However, the migrations used little extra memory relative to the Cartesian geometry runs.  Overall, the test results illustrate the significant implementation advantages afforded by analytic coordinate systems with respect to run time and memory requirements.  In the following chapters, I will build upon these advantages to help implement 2D and 3D prestack migration algorithms.
%
\begin{table}
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|}
\hline
Mig.  & Geometry & Run & Relative Run &Max Memory  & Relative \\
Type & \,   & Time (s) & Time & Usage (Gb) & Memory Usage \\
\hline
\hline  CC &Stored  &  1304  &  1.00 & 3.3  & 1.000 \\
\hline  CC & Calculated &  1333 & 1.022 & 2.6 & 0.788 \\
\hline
\hline  EC & Stored  & 1330  & 1.020 &  3.4 & 1.030  \\
\hline  EC & Calculated &  1380 & 1.058 &   2.7 & 0.812\\
\hline 
\end{tabular}
\caption{Comparison between performing migration with stored or locally computed coordinate system geometry information for Cartesian (CC) and Elliptic Cylindrical (EC) coordinates.}
\label{tab:RayDifference}  
\end{center}
\end{table}
% 
\section{Conclusions}
%
This chapter addresses existing issues with Riemannian wavefield extrapolation theory by extending RWE to smoother, but non-orthogonal, coordinate systems.   I demonstrate that acoustic wave equations can be generated in general 3D Riemannian spaces, and that the corresponding extrapolation wavenumber decouples from the other wavenumbers.  Using extensions to the PSPI and SSF approximation techniques, I incorporate this wavenumber into a one-way extrapolation operator appropriate for propagating wavefields.   I present a method for generating computational meshes that are unconditionally singularity-free, which is used to generate examples illustrating wavefield propagation on non-orthogonal coordinate meshes using RWE operators.  Finally, I argue that analytic coordinate systems afford the possibility of more accurate extrapolation operators implementations, which opens up a range of imaging possibilities including shot-profile migration in elliptic coordinates (Chapter~3), the calculation of angle-domain common-image gathers in generalized coordinates (Chapter~4) and inline delayed-shot migration in tilted elliptical-cylindrical coordinates (Chapter~5).   
%
\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 the Geophysics manuscript and thank them for their helpful comments.

\cleardoublepage
\end{document}
