\lefthead{Shragge} \righthead{RWE: NOC} \footer{SEP--124} %\usepackage{color} \title{Riemannian Wavefield Extrapolation: Non-orthogonal coordinate systems} \email{jeff@sep.stanford.edu} \keywords{} \author{Jeff Shragge} \sepinput{mainmacros} \sepinput{macros} \ABS{Riemanninan wavefield extrapolation (RWE) is a method for performing 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. \textcolor{red}{This paper develops a procedure that avoids many of these problems by posing wavefield extrapolation on smoother, but generally non-orthogonal and singularity-free, coordinate meshes}. \textcolor{red}{The resulting extrapolation operators include additional terms describing non-orthogonal propagation; however, these extra degrees of complexity are offset by smoother coefficients.} Two analytic coordinate system examples are used to validate the non-orthogonal propagation theory. \textcolor{red}{A method for eliminating any remaining singularities from coordinate systems is presented}, and is then used in calculating 2D and 3D Green's functions for cylindrical and near-spherical geometry. Results of 2-D benchmark testing suggest that the computational overhead associated with the RWE approach is roughly 35\% greater than Cartesian-based extrapolation.} % \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}. \textcolor{red}{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.} \longcite{Sava.geo.rwepub} followed this approach in developing Riemannian wavefield extrapolation (RWE), a theory of one-way wavefield propagation for \textcolor{red}{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 \textcolor{red}{the degree to which 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. 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 $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 \textcolor{red}{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 occurs in coordinate system singularities (see upper panel of figure~\ref{fig:Figure1}) that arise wherever a mesh is developed from a crossing set of rays. \textcolor{red}{This generates spatial singularities and singular Jacobians that lead to zero-divisions 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 \longcite{Shragge.seg.2005}, who formulate a wave-equation migration from topography strategy that poses wavefield extrapolation directly in locally orthogonal meshes conformal to the acquisition surface. This approach successfully generates subsurface images beneath areas exhibiting longer wavelength and lower amplitude relief; however, imaging results in situations involving more rugged acquisition topography degrade due to the grid compression/extension demanded by semi-orthogonality geometry (see lower panel of figure~\ref{fig:Figure1}). % \par % In this paper, I argue that the RWE approach can be made more broadly applicable by resolving the issues related to semi-orthogonal geometry using a combination of smoother computational meshes and the judicious removal of all coordinate singularities. In doing so, though, one must relax the semi-orthogonality geometry restriction and propagate wavefields on non-orthogonal grids. Fortunately, the existing RWE framework can be extended in a straightforward manner, leading to a wavenumber appropriate for extrapolating wavefields on non-orthogonal geometries. This wavenumber forms the basis of a one-way extrapolation operator used 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. The development presented herein naturally follows that of \longcite{Sava.geo.rwepub}; however, I recast the theory in a more compact notation that leads to a closer analytic connection of the generalized computation geometry with the underlying Cartesian grid. The paper begins with the formulation of the 3D Riemannian acoustic wave-equation and the corresponding non-orthogonal one-way wavefield extrapolation wavenumber. Appendix~A presents an overview of the required differential geometry theory, while the split-step Fourier 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 the procedure for generating triplication-free coordinate systems, present 2D and 3D Green's functions estimates modeled in cylindrical and near-spherical coordinates, respectively, and discuss the relative computational cost and memory overhead of the RWE method. % \plot{Figure1}{width=7in}{Illustration of problems with RWE computational grids. Top panel: Ray-coordinate system triplications. Bottom panel: Grid bunching for topographically conformal coordinate system.} % \section{Acoustic wave equation in 3D generalized Riemannian spaces} % Specifying the acoustic wave-equation in a 3D Riemannian space requires formulating the physics of wave-propagation in a generalized coordinate system framework. By definition, generalized Riemannian coordinates are related to the underlying Cartesian mesh by unique transformations (i.e. 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 acoustic wave-equation for a wavefield, $\W$, in a generalized Riemannian space is, \beq \label{eqn:helm} \dd_\bxi \W= - \ww^2\ss^2 \lp \mathbf{\bxi} \rp \W, \eeq where $\dd_\bxi$ is the Laplacian operator applied in coordinates $\bxi$, $\ww$ is frequency, and $\ss$ is the propagation slowness. % \par % Correctly formulating the wave-equation in coordinate system $\bxi$ requires that Laplacian operator $\dd_\bxi$ be specified by differential geometry relationships. (An overview of necessary differential geometry theory is provided in Appendix~A.) The Laplacian operator in generalized coordinates is, \beq \label{eqn:laplac2} \dd_\bxi \W = \frac{1}{\sqrt{|\mathbf{g}|}}\, \frac{\partial}{\partial \xi_i}\, \left(\,\sqrt{|\mathbf{g}|} g^{ij}\,\frac{\partial \W}{\partial \xi_j}\right) = \frac{1}{\sqrt{|\mathbf{g}|}}\, \frac{\partial}{\partial \xi_i}\, \left(\,m^{ij}\,\frac{\partial \W}{\partial \xi_j}\right), \quad i,j=1,2,3, \eeq where $g^{ij}$ is an element of the metric tensor $\mathbf{g}$, $|\mathbf{g}|$ is the metric tensor discriminant, and $m^{ij}=\sqrt{|\mathbf{g}|} g^{ij}$ is weighted metric tensor element that enables a more compact notation. Unless otherwise stated, summation over all repeated 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, \beq \label{eqn:we1} \frac{1}{\gdetf} \de\qi \lp \mmm\ii\jj \done\W\qj \rp = - \ww^2\ss^2 \W. \eeq % \par % The first step in developing a generalized RWE wave-equation dispersion relationship is to expand the derivative terms in equation~\ref{eqn:we1} and multiply through by $\gdetf$ to obtain, \beq \label{eqn:helm1} \done{\mmm\ii\jj}\qi \done\W\qj + \mmm\ii\jj \mtwo\W\qi\qj = - \gdetf \ww^2 \ss^2 \W. \eeq \textcolor{red}{The 3D RWE acoustic wave equation derivation here deviates from that in \longcite{Sava.geo.rwepub}, who introduce nine coefficients expressing the spatial metric tensor derivatives. In this development, I choose to express all metric derivatives in a notation that remains close to the underlying geometry.} 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 leads to a simultaneous spatial and Fourier wavenumber dependence; however, this is routinely handled through various approximations discussed below. % \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 evaluating indices 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} An expression for wavenumber $\kz$ is obtained through a complete-the-square transform, and isolating the contributions of wavenumber $\kz$ to yield %\begin{small} %\beqa %\,& m^{33} \lp \kz + %\frac{ \lp 2 m^{13}\kx + 2m^{23}\ky- {\rm i} \,n^3 \rp }{ 2m^{33} } %\rp^2 = \, %\gdetf\ww^2\ss^2 - %\kx^2 \lp m^{11}-\frac{ \lp m^{13} \rp^2 }{ m^{33}} \rp - \\ %\,& \ky^2 \lp m^{22}-\frac{ \lp m^{23} \rp^2 }{ m^{33}} \rp - %\kx\ky \lp 2 m^{12}-\frac{ 2 m^{13}m^{23}}{ m^{33}} \rp + %{\rm i} \, \kx \lp n^1 - \frac{m^{13}\,n^3}{ m^{33} } \rp + {\rm i} \, %\ky \lp n^2 %- \frac{m^{23}\,n^3}{ m^{33} }\rp - \frac{(n^3)^2 }{ m^{33}}. %\,\nonumber %\eeqa %\end{small} \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 of $a_i$ in equation~\ref{eqn:kzfinal} are presented in the following vector $\mathbf{a}$, \begin{tiny} \beq \label{eqn:ametric2} \mathbf{a} = \left[ \frac{g^{13}}{ g^{33}} \;\;\;\;\; \frac{g^{23}}{ g^{33} }\;\;\;\;\; \frac{ n^3 } {2 m^{33} } \;\;\;\;\; \frac{ \ss }{ \sqrt{g^{33}} } \;\;\;\;\;\; \sqrt{ \frac{g^{11}}{ g^{33}} - \lp \frac{g^{13}}{g^{33}} \rp^2 } \;\;\;\;\; \sqrt{\frac{g^{22}}{ g^{33} } - \lp \frac{ g^{23}}{g^{33}} \rp^2 } \;\;\;\;\; \lb \frac{2\,g^{12}}{g^{33}} - \frac{2\, g^{13}\,g^{23}}{ \lp g^{33} \rp^2} \rb \;\;\;\;\; \lb \frac{n^1}{m^{33}} - \frac{m^{13}\, n^3}{ \lp m^{33}\rp^2}\rb \;\;\;\;\; \lb \frac{n^2}{m^{33}} - \frac{m^{23}\, n^3}{ \lp m^{33}\rp^2}\rb \;\;\;\;\; \frac{ n^3 }{2 m^{33} } \right]^{\mathbf{T}}. \eeq \end{tiny} 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 % %\textcolor{red}{x} 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-\boldsymbol{k_\xi}$) phase-shift extrapolation scheme can be developed to recursively advance a wavefield from level $\xi_3$ to level $\xi_3+\Delta \xi_3$, % \beq \label{eqn:bulk} \W(\xi_3+\Delta \xi_3, \kx,\ky ; \omega) = \W(\xi_3, \kx, \ky ; \omega) {\rm e}^{ {\rm i} \kz \Delta \xi_3}. \eeq % If $\W$ represents a post-stack wavefield, an image $\mathcal{I}$ can be produced from the propagated wavefield through an imaging condition \cite{Claerbout85}, % \beq \label{eqn:ZOI} \mathcal{I}(\xi_3,\xi_1,\xi_2) = \sum_{\omega} \W(\xi_3,\xi_1,\xi_2; \omega). \eeq % \par % Situations where coefficients vary across an extrapolation step, though, require additional approximations. One straightforward approach is to use the multi-coefficient split-step Fourier (SSF) method \cite{GEO55-04-04100421} detailed in Appendix~B. This method uses Taylor expansions of the dispersion relation about a set of reference parameters to form a bulk phase-shift operator in the Fourier domain ($\omega-\boldsymbol{k_\bxi}$). Differences between the reference and true parameters then form a correction term applied in the mixed $\omega-\bxi$ domain. In addition, the phase-shift extrapolation scheme can be extended to incorporate multiple reference media followed by interpolation similar to the phase-shift plus interpolation (PSPI) technique of \longcite{Gazdag.1984}. All results in the following sections were generated with the combined PSPI plus SSF correction extrapolation operator. % \subsection{Discussion} % The accuracy of the multi-coefficient SSF approach is directly related to the degree to which coefficients in equation~\ref{eqn:ametric2} vary across a propagation step. At a first glance, one might expect that far too many expansions are required to make a PSPI approach practical. (For example, three reference expansions for each of the ten terms would seemingly require $3^{10}=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? In this paper, I calculate reference coefficients using a multi-dimensional Lloyd's algorithm developed by \longcite{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 \longcite{Clapp.Lloyds}. % \par % Second, numerous situations exists where some coefficients are zero or negligible. For example, the mixed-domain fields for a 3D weakly non-orthogonal coordinate system 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} % \textcolor{red}{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 left panel of figure~\ref{fig:Figure2}). This coordinate system is uniquely specified by one additional degree of freedom and is related to an underlying Cartesian mesh through the following transformation, % \beq \label{eqn:2Dred} \left[\begin{array}{c} \xone \\ \xthr \end{array}\right] = \left[\begin{array}{cc} 1 & \rm{cos}\, \theta \\ 0 & \rm{sin}\, \theta \end{array}\right] \left[\begin{array}{c} \qx \\ \qz \\ \end{array}\right], \eeq where $\theta$ is the shear angle of the coordinate system ($\theta=90^{\circ}$ is Cartesian). The metric tensor of this transformation is, \beq \label{eqn:2Dmet} \left[g_{ij}\right] = \left[\begin{array}{cc} \hxxdef & \hxzdef \\ \hxzdef & \hzzdef \end{array}\right] = \left[\begin{array}{cc} \EE & \HH \\ \HH & \AA \end{array}\right] = \left[\begin{array}{cc} 1 & \rm{cos}\, \theta \\ \rm{cos}\,\theta & 1 \end{array}\right], \eeq and has a discriminant $\gdet = \rm{sin}^2\,\theta$ and an associated metric tensor $g^{ij}$ given by, \beq \label{eqn:2Dmet6} \left[g^{ij}\right] = \frac{1}{\rm{sin}^2\,\theta} \left[\begin{array}{cc} 1 & - \rm{cos}\,\theta \\ -\rm{cos}\,\theta & 1 \end{array}\right]. \eeq % Note that because the tensor in equation~\ref{eqn:2Dmet6} is coordinate invariant, equation~\ref{eqn:helm2} simplifies to, \beq \label{eqn:laplac} \dd_\bxi \W = g^{ij} \,\frac{\partial^2 \W}{\partial \qi \partial \qj} = - \ww^2\ss^2\W, \eeq which generates the following dispersion relation, \beq g^{ij}\ki\kj = \ww^2\ss^2. \eeq Expanding out these terms leads to an expression for wavenumber $\kz$, \beq \label{eqn:kz2Dex} \kz = -\frac{ g^{13} }{ g^{33} } \kx \pm \sqrt{\frac{\ss^2\ww^2}{g^{33}}- \lp\frac{g^{11}}{g^{33}}-\lp \frac{g^{13}}{g^{33}} \rp^2\rp \kx^2}. \eeq Substituting the values of the associated metric tensor in equation~\ref{eqn:2Dmet6} into equation~\ref{eqn:kz2Dex} yields, \beq \kz = \rm{cos} \, \theta \;\kx \pm \rm{sin} \,\theta \sqrt{\ss^2\ww^2 - \kx^2}, \eeq which is appropriate for performing RWE on the 2D sheared Cartesian coordinate system shown in left panel in figure~\ref{fig:Figure2}. % \par % The right panel of figure~\ref{fig:Figure2} 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 4 flat plane-waves at times $t=0.2, 0.4, 0.6$ and $0.8$ s. Zero-offset migration results generated by equation~\ref{eqn:ZOI} show migrated reflectors at the expected depths of $z$=300, 600, 900, and 1200~m. 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. % \sideplot{Figure2}{width=7in}{Sheared Cartesian coordinate system test. Coordinate system (left panel) shear angle and velocity are $\theta=25^{\circ}$ and 1500 ms$^{-1}$, respectively. Zero-offset data consist of 4 flat plane-wave impulses at $t=0.2, 0.4, 0.6$ and $0.8$ s that are correctly imaged (right panel) at depths $z=300, 600, 900,$ and $1200~m$.} % \subsection{Polar Ellipsoidal Coordinates} % A second example is an ellipsoidal polar coordinate system (see upper left panel in figure~\ref{fig:Figure3}) 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 controling 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 - 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 \longcite{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. The upper left panel shows the polar ellipsoidal coordinate system mesh overlying a linear $v(z)=1500+0.35\,z$ m/s velocity function. The upper right panel presents the velocity model mapped into the RWE domain under the transformations defined in equation~\ref{eqn:2DPE1}. % \plot{Figure3}{width=6.5in}{Ellipsoidal polar coordinate system test example. Upper left: $v(z)=1500+0.35\,z$ velocity function overlain by a polar ellipsoidal coordinate system defined by parameter $a=1+0.2\,\xi_3-0.05\,\xi_3^2$. Upper right: Velocity model mapped in the RWE domain. Bottom right: Imaged reflectors in RWE domain. Bottom left: the RWE domain image mapped to a Cartesian mesh.} % \textcolor{red}{The test data consist of ten plane waves defined at the surface by a ray parameter $p_x=-0.5$ s/km. At greater depths the waves are no longer planar and propagate through a turning point before turning upward to the left.} 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. This problem is illustrated by the panels in figure~\ref{fig:Figure4}. The top panel 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}. Note that the 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 % The center panel 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 mesh generation procedure is to extract initial and final isochron surfaces from the Eikonal equation solution to form the inner and outer mesh boundaries. The mesh domain is then enclosed by interpolating between the edges of the inner and outer bounding surfaces. The interior mesh can then be formed through bi-linear interpolation methods, such as blending functions \cite{Liseikin,ShraggeDMG}. % \par % A corresponding singularity-free mesh is presented in the bottom panel of figure~\ref{fig:Figure4}. The grid is regularly-spaced on the outer isochron and has a couple of 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 mesh smoothing usually does not affect propagation accuracy because the coordinate system mesh forms only the skeleton on which wavefield extrapolation occurs. \textcolor{red}{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 triplication-free mesh generation. Top panel: Velocity model with three Gaussian velocity perturbations. Overlain is a coordinate mesh generated from ray-tracing. Note the triplication to either side of the shot-point, as well as the spreading beneath the shot point. Middle panel: Velocity model overlain by isochrons of an Eikonal equation solution for same shot-point. Bottom panel: Triplication-free computational mesh generated by through the approach specified in this section. } % \subsection{2D Green's Function Generation} % The third test uses RWE to model 2D Green's functions on coordinate systems constructed by the aforementioned approach. Figure~\ref{fig: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. The upper left panel shows the velocity model with an overlain coordinate system generated by meshing procedure discussed in the section above. The velocity model in the RWE domain is illustrated in the upper right panel. % \par % The lower right panel shows the impulse response tests in the RWE domain. The impulses conform fairly well to the travel-time steps, except where they enter the salt body in the lower left of the image. The migration results mapped back to Cartesian space are shown in the lower left panel. The complicated wavefield 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. Top left: Salt model in physical space with an overlain ray-based mesh. Top right: Velocity model in the RWE domain. Bottom right: Wavefield propagated in ray-coordinates through velocity model shown in the upper right. Bottom left: Wavefield in bottom right interpolated back to Cartesian space.} % Figure~\ref{fig: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 Cartesian and finite difference point-source versus the RWE plane-wave. 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, propagate 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 % \textcolor{red}{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. Finite differences better models amplitudes in the presence of velocity gradients in the propagation direction, and thus more accurately partitions energy at the top sediment-salt interface. The RWE wavefield underestimates the reflection contribution and allows significantly more energy to be transmitted into and through the salt body. This energy causes the additional multipathing noted in the RWE panel below the salt body, leading the more complicated wavefield behavior relative to the Cartesian wavefield example.} % \plot{Figure6}{width=6.5in}{Comparison between three different extrapolation methods. Top panel: two-way finite-difference; middle panel: Riemannian wavefield extrapolation; and bottom panel: Cartesian wavefield extrapolation. Note the fairly similar top-salt reflections in upper left of the top two panels. The lower two panels also have stronger amplitudes for the intra-salt phases.} \subsection{3D Point Source Coordinates} % The procedure discussed above can be extended to 3D point-source coordinate systems. \textcolor{red}{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 Cartesian domain. Herein, 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 % \textcolor{red}{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 (2 \xi_1 / \pi + 2 \, w\, \alpha {\rm cos} \, \xi_2 \, {\rm sin}(\xi_1 /2)) /2 \nonumber \\ \,&\,& x_2 = \xi_3 ( \xi_2 + 2\, 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 One then performs 3D wavefield extrapolation 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. % \plot{Figure8}{width=7in}{3D modeling test of Winkel-Tripel coordinate system using the Elf-IFP-CGG synthetic velocity. Eight wavefield snapshots are shown in the image. Note the lack of azimuthal anisotropy commonly found in 3D impulse response tests based on finite-difference splitting methods. Boundary reflection artifacts, though, are still present throughout the volume.} % 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 right creates a more complex geologic setting. Overlain are a number of impulse responses associated with different travel times. %Note that the wavefield does not exhibit azimuthally variable numerical dispersion commonly found in 3D impulse response tests based on finite-difference splitting methods \cite{3DSI}. Also, because the coordinate system does not cover the full velocity field, one can observe boundary reflection artifacts within the propagation domain (similar to the tilting Cartesian example). Overall, though, this extrapolation results demonstrates the stability of wavefield propagation on a 3D grid defined by the Winkel-Tripel spatial discritization. % \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. \textcolor{red}{A total of 82 frequencies were propagated a total of 511 extrapolation steps, thus requiring 41092 calls to the SSF operator. In addition, the tests involved 112996 calls to the phase-shift routine, or almost three per extrapolation step as this number varied with velocity model complexity.} \textcolor{red}{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. \textcolor{red}{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 Traynor 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{sep} \bibliography{SEP,MISC,GEOTLE,allrefs} \newpage \APPENDIX{A} % Geometry in a generalized 3D Riemannian space is described by a symmetric metric tensor, $g_{ij}=g_{ji}$, that relates the geometry in a general non-orthogonal coordinate system, $\{\xi_1,\xi_2,\xi_3\}$ , to an underlying Cartesian mesh, $\{x_1,x_2,x_3\}$ \cite{Guggenheimer}. In matrix form, the metric tensor is written, \beq \label{eqn:metric} \left[g_{ij}\right] = \left[\begin{array}{ccc} g_{11} & g_{12} & g_{13} \\ g_{21} & g_{22} & g_{23} \\ g_{31} & g_{32} & g_{33} \end{array}\right] = \left[\begin{array}{ccc} g_{11} & g_{12} & g_{13} \\ g_{12} & g_{22} & g_{23} \\ g_{13} & g_{23} & g_{33} \end{array}\right], \eeq where $\EE$, $\FF$, $\GG$, $\HH$, $\LL$ and $\AA$ are functions linking the two coordinate systems through, \beqa \label{eqn:pderels} \EE=\hxxdef, \;\;\;\; \FF=\hxydef, \;\;\;\; \GG = \hyydef,\nonumber \\ \HH=\hxzdef, \;\;\;\; \LL=\hyzdef, \;\;\;\; \AA = \hzzdef. \eeqa The associated (or inverse) metric tensor, $g^{ij}$, is 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} \eeq A weighted metric tensor, $m^{ij}=\gdetf \, g^{ij}$, is also used throughout the paper. % \newpage % \APPENDIX{B} % The extrapolation wavenumber defined in equations~\ref{eqn:kzfinal} and \ref{eqn:ametric2} cannot generally be implemented exactly in the Fourier domain due to a simultaneous spatial dependence (i.e. a function of both $x_1$ and $\kx$). This can be addressed using an 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_x}$) domain phase-shift and a mixed ($\omega-\mathbf{x}$) domain split-step correction, respectively. % \par The phase-shift term is given by, \beq \kz^{PS} = - b_1 \kx - b_2 \ky + i b_3 \pm \lb b_4^2 \ww^2 - b_5^2 \kx^2 - b_6^2 \ky^2 - b_7 \kx\ky + b_8 i\,\kx + b_9 i\,\ky - b_{10}^2 \rb^{\frac{1}{2}}, \eeq where $b_i=b_i(\xi_3)$ are reference values of $a_i =a_i(\xi_1,\xi_2,\xi_3)$. The split-step approximation is developed by performing a Taylor expansion about each coefficient $a_i$ and evaluating the results at stationary reference values $b_i$. Assuming that the stationary values of $\kx$ and $\ky$ are zero, the split-step correction is as follows, \beq \label{eqn:kzssf} \kz^{SSF}= \left. \frac{\partial \kz}{\partial a_3} \right|_0 \lp a_3-b_3 \rp + \left. \frac{\partial \kz}{\partial a_4} \right|_0 \lp a_4-b_4 \rp + \left. \frac{\partial \kz}{\partial a_{10}} \right|_0 \lp a_{10}-b_{10} \rp, \eeq where ``$0$'' denotes "with respect to a reference medium". The partial differential expressions in equation~\ref{eqn:kzssf} are, \beq \left. \frac{\partial \kz}{\partial a_3} \right|_0 = b_3, \;\;\;\;\; \left. \frac{\partial \kz}{\partial a_4} \right|_0 = \frac{b_4 \,\ww^2}{\sqrt{b_4^2 \, \ww^2 - b_{10}^2 }}, \;\;\;\;\; \left. \frac{\partial \kz}{\partial a_{10}} \right|_0 = - \frac{b_{10}}{\sqrt{b_{10}^2 \, \ww^2 - b_{10}^2 }}, \;\;\;\;\; \eeq resulting in the following split-step Fourier correction wavenumber, \beq \kz^{SSF} = i\,b_3 \, \lp a_3-b_3 \rp + \frac{b_4 \, \ww^2 \, \lp a_4-b_4 \rp}{\sqrt{b_4^2 \, \ww^2 - b_{10}^2 }} - \frac{b_{10} \, \lp a_{10}-b_{10} \rp}{\sqrt{b_4\,\ww^2 - b_{10}^2 }}. \eeq % \newpage % \APPENDIX{C} % The extrapolation wavenumber developed in equation~\ref{eqn:kzfinal} is appropriate for any non-orthogonal Riemannian geometry. However, there are a number of situations where symmetry or partial orthogonality are present. Moreover, 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 = i a_3 \pm \lb a_4^2 \ww^2 - a_5^2 \kx^2 - a_6^2 \ky^2 - a_7 \kx\ky + i a_8 \kx + i a_9 \ky - a^2_{10} \rb^{\frac{1}{2}}, \eeq where, \beq \mathbf{a} = \lb 0 \;\;\; 0 \;\;\; \frac{ n^3}{2 m^{33}} \;\;\; \frac{ \ss }{ \sqrt{g^{33}} } \;\;\; \sqrt{\frac{g^{11}}{g^{33}}} \;\;\; \sqrt{\frac{g^{22}}{g^{33}}} \;\;\; \frac{ 2 g^{12} }{ g^{33} } \;\;\; \frac{n^1}{m^{33}} \;\;\; \frac{n^2}{m^{33}} \;\;\; \frac{n^3}{2 m^{33}} \rb^{\mathbf{T}}. \eeq which are identical to the coefficients recovered by \longcite{Sava.geo.rwepub}. % \par % {\bf 3D Kinematic Coordinate Systems} - Wave-equation migration amplitudes are generally inexact in laterally variant media - even in Cartesian-based systems. Hence, one beneficial approximation that reduces computational cost is to consider only the kinematic terms in equation~\ref{eqn:kzfinal}. This approximate generates the following extrapolation wavenumber, \beq \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+ i a_3 \pm \lb a_4^2 \ww^2 - a_5^2 \kx^2 + i a_8 \kx - a_{10}^2 \rb^{\frac{1}{2}}, \eeq where, \beq \mathbf{a} = \lb -\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 + 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& 41092 & 4.2 & 0.051 \\ \hline RWE &Phase-shift & 112996 & 17.1 & 0.209 \\ \hline RWE &Interpolation & 112996 & 4.6 & 0.056 \\ \hline \hline Cart &Frequency Loop & 82 & 40.9 & 0.499 \\ \hline Cart &Split-step Fourier& 41092 & 2.1 & 0.026 \\ \hline Cart &Phase-shift & 112996 & 7.6 & 0.093 \\ \hline Cart &Interpolation & 112996 & 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}