Coordinate systems created by ray tracing in a background medium often well represent wavefield propagation. In this context, we effectively split wave propagation effects into two parts: one part accounting for the general trend of wave propagation, which is incorporated in the coordinate system, and the other part accounting for the details of wavefield scattering due to rapid velocity variations. If the background medium is close to the real one, the wave-propagation can be properly described with low-order operators. However, if the background medium is far from the true one, the wavefield departs from the general direction of the coordinate system and the low-order extrapolators are not enough for accurate description of wave propagation.
For coordinate system describing a geometrical property of the medium (e.g. migration from topography), there is no guarantee that waves propagate in the direction of extrapolation. This situation is similar to that of Cartesian coordinates when waves propagate away from the vertical direction, except that conformal mapping gives us the flexibility to define any coordinates, as required by acquisition. In this case, too, low-order extrapolators are not enough for accurate description of wave propagation.
Therefore, we need to develop higher-order Riemannian wavefield extrapolators in order to handle correctly waves propagating obliquely with the coordinate system. Usually, the high-order extrapolators are implemented as mixed operators, part in the Fourier domain using a reference medium, part in the space domain as a correction from the reference medium. Many methods have been developed for high-order extrapolation in Cartesian coordinates. In this paper, I explore some of those extrapolators in Riemannian coordinates. In particular, I concentrate on high-order finite-differences solutions, and methods from the pseudo-screen family Huang et al. (1999) and Fourier finite-differences family Biondi (2002); Ristow and Ruhl (1994). In theory, any other high-order extrapolator developed in Cartesian coordinates can have a correspondent in Riemannian coordinates.
In this paper, I implement the finite-differences portion of the high-order extrapolators with implicit methods. Such solutions are accurate and robust, but they face difficulties for 3D implementations because the finite-differences part cannot be solved by fast tridiagonal solvers anymore and require more complex and costlier approaches Fomel and Claerbout (1997); Rickett et al. (1998). The problem of 3D wavefield extrapolation is addressed in Cartesian coordinates either by splitting the one-way wave-equation along orthogonal directions Ristow and Ruhl (1997), or by explicit numerical solutions Hale (1991). Similar approaches can be envisioned for 3D Riemannian extrapolation. The explicit solution seems more appropriate, since splitting is difficult due to the mixed terms of the Riemannian equations. I do not address this subject in this paper, and concentrate on developing higher-order kernels with implicit methods.