Marchenko Redatuming in Seismic Imaging

• Marchenko redatuming is a data-driven method that reconstructs multi-component, internal-multiple-free wavefields at arbitrary subsurface locations.
• It solves coupled integral equations for focusing functions, enabling the extraction of up- and downgoing Green's functions from reflection data.
• Its applications span high-resolution imaging, time-lapse reservoir monitoring, and efficient inversion in laterally heterogeneous and complex media.

Marchenko redatuming is a data-driven, single-sided method to reconstruct multi-component, multiple-free wavefields at arbitrary subsurface locations by utilizing surface reflection data and a smooth macro-model. It achieves this by solving coupled integral equations for so-called focusing functions, which, when used in Marchenko reciprocity integrals, yield up- and downgoing Green's functions corrected for all orders of internal multiples, including evanescent, refracted, and tunneling energy. The method is rooted in acoustic wave theory but extends to highly general settings, including laterally heterogeneous, multi-dimensional, and non-ideal acquisition scenarios. Its frameworks power advanced imaging workflows in seismology, reservoir characterization, and monitoring.

1. Theoretical Formulation and Integral Equations

At its core, Marchenko redatuming is governed by the acoustic wave equation with mixed Dirichlet-Neumann boundary conditions, enforcing a space-time focusing constraint. For a position $x = (x_H, z)$, let $f_\pm(x, t; x_f)$ denote the down- or up-going focusing function for a focal point $x_f = (x_{H, f}, z_f)$. In a lossless heterogeneous medium: $$\left[ \rho(x)\,\nabla \cdot \left(\frac{1}{\rho(x)}\nabla \right) - \frac{1}{c^2(x)}\frac{\partial^2}{\partial t^2} \right]f_\pm(x, t; x_f) = 0$$ These focusing functions satisfy:

• Homogeneous initial and final conditions (radiation):

$f_\pm(x, t \rightarrow \pm\infty; x_f) = 0$

• Focusing (Dirichlet) constraint at $z = z_f$:

$$f_\pm(x_H, z_f, t; x_f) = \delta(t)\delta(x_H - x_{H,f})$$

• One-way (Neumann/generalized impedance) constraint:

$$\left[\partial_z \mp \Lambda \right]f_\pm(x_H, z_f, t; x_f) = 0$$

where $\Lambda$ is the square-root spatial-temporal operator defined by the local dispersion relation.

These conditions guarantee the unique, physically admissible solution for the focusing function as a wavefield that focuses energy at $(x_f, t=0)$ from one boundary, without requiring sources/receivers in the subsurface. Marchenko redatuming then uses single-sided reflection data $f_\pm(x, t; x_f)$0 from the acquisition surface to construct the virtual wavefields and Green's functions between any focal point $f_\pm(x, t; x_f)$1 and receiver $f_\pm(x, t; x_f)$2 in the domain (Hajjaj et al., 2022).

2. Computational Algorithms and Inversion Strategies

The forward problem is discretized for practical implementation. The focusing function $f_\pm(x, t; x_f)$3 is approximated by a wavefield $f_\pm(x, t; x_f)$4 on a regular grid and is found by minimizing the cost functional: $f_\pm(x, t; x_f)$5 Here, $f_\pm(x, t; x_f)$6 is the discrete two-way wave operator, $f_\pm(x, t; x_f)$7 applies the spatial mask for the PDE, $f_\pm(x, t; x_f)$8 enforces the focusing and one-way conditions, $f_\pm(x, t; x_f)$9 encodes the delta focusing, and $x_f = (x_{H, f}, z_f)$0 is the Tikhonov parameter for regularization. This leads to the normal equations: $x_f = (x_{H, f}, z_f)$1 These are solved iteratively using methods such as conjugate gradients or LSQR (Hajjaj et al., 2022).

The resulting focusing functions, which capture evanescent, refracted, and tunneling waves as well as propagating modes, are used in Marchenko reciprocity integrals to reconstruct up- and downgoing Green’s functions at arbitrary subsurface points: $x_f = (x_{H, f}, z_f)$2 where $x_f = (x_{H, f}, z_f)$3 denotes time convolution over the acquisition surface (Hajjaj et al., 2022).

3. Multidimensional Effects, Evanescent Waves, and Stability

In higher spatial dimensions, handling evanescent waves requires explicit attention as they introduce severe ill-conditioning in the inverse problem, particularly manifesting at large horizontal wavenumbers and low frequencies. The spectrum of the discretized system spans many orders of magnitude, driven by strong evanescent modes that marginally satisfy the wave equation but do not respect the tight curvature requirements of the focusing constraints (Hajjaj et al., 2022, Wapenaar, 2020).

Mitigation strategies include:

• Careful selection of the Tikhonov parameter ($x_f = (x_{H, f}, z_f)$4 typical), balancing the residual norms of focusing and wave-equation constraints.
• Spectral or spatial-domain preconditioning, such as down-weighting large-$x_f = (x_{H, f}, z_f)$5 modes or frequency filters.
• Multiscale (low-frequency to high-frequency) incremental solvers.
• Variable-projection methods that partition and analytically eliminate portions of the model vector (Hajjaj et al., 2022, Wapenaar, 2020).

The downward-decaying evanescent solution ($x_f = (x_{H, f}, z_f)$6) can be retrieved stably over moderate depths, whereas the upward-decaying component ($x_f = (x_{H, f}, z_f)$7) is highly sensitive to small errors in the direct transmission estimate, with errors amplifying exponentially with depth (Wapenaar, 2020).

4. Extensions: Handling Imperfect Sampling, Plane Waves, and Wide-Angle Arrivals

Recent methodological advances extend classical Marchenko redatuming to handle realistic field data imperfection:

• Point-spread-function (PSF) based corrections in the Marchenko iteration enable accurate reconstruction with irregular or incomplete acquisition geometries by embedding multidimensional deconvolution after every Marchenko update (IJsseldijk et al., 2020, IJsseldijk et al., 2020).
• Plane-wave Marchenko methods leverage focusing on spatial planes, dramatically reducing computational cost in 3D and allowing efficient imaging of large datasets with internal multiple suppression, although with some limitations in dipping event illumination (Meles et al., 2017, Thorbecke et al., 2023).
• Dual-focusing and double-focusing schemes synthesize virtual sources and receivers at subsurface boundaries, providing the foundation for efficient, target-oriented least-squares RTM workflows and further reducing artifacts from overburden multiples (Shoja et al., 2023, Shoja et al., 2022).

Marchenko redatuming has also been extended to retrieve forward-scattered waveforms by using transmission data, thereby including essential non-reversed arrivals and their multiples in the reconstructed Green’s functions (Neut et al., 2022).

5. Applications: Imaging, Time-Lapse Monitoring, and Practical Implementation

Marchenko redatuming underpins a wide range of modern seismic imaging and monitoring procedures:

• Target-oriented LSRTM: Enables high-resolution imaging within selected subsurface zones by reducing both data and computational domain dimensions, while eliminating overburden-generated internal multiples (Shoja et al., 2022, Shoja et al., 2023).
• Time-lapse reservoir monitoring: Produces isolated, multiple-free reflection responses for arbitrary reservoir intervals, allowing accurate measurement of sub-millisecond differential traveltimes for fluid and pressure surveillance in CCS or hydrocarbon operations (IJsseldijk et al., 2023).
• Hard seafloor and seabed acquisition: Reciprocity-based formulations, such as Upside-Down Rayleigh–Marchenko, enable exact redatuming with coarsely sampled or irregular seabed receiver layouts by transferring all spatial integration to the denser source carpet (Wang et al., 2024).
• Large-scale HPC deployment: Distributed, matrix-free solvers and compressed representations (e.g., via Dask, Zarr, PyLops frameworks) allow scalable inversion and redatuming in 3D volumes, with demonstrated resilience to moderate aliasing and spatial decimation (Ravasi et al., 2020).

6. Mutual Relations with Propagator/Transfer Matrices and Full-Wave Extrapolation

Marchenko focusing functions and the wavefield propagator matrix are intimately linked—the focusing function is a linear combination of the propagator blocks (specifically, the pressure and pressure–velocity components), and the full propagator is reconstructible from knowledge of the focusing function via real and imaginary part separation and the full Helmholtz square-root operator. This linkage guarantees that Marchenko-based redatuming handles both propagating and evanescent waves, and that the two-way propagator inherits all internal-multiple elimination properties from the data-driven focusing function. Transfer matrices, used to decompose fields into up- and down-going waves, are also expressible in terms of the focusing functions, allowing for unified formulation of advanced redatuming and imaging operators without restrictive approximations (Wapenaar et al., 2023, Wapenaar et al., 2021).

7. Limitations, Assumptions, and Outlook

Marchenko redatuming requires:

• A sufficiently accurate, smooth overburden velocity model for estimating direct arrivals and defining focusing windows.
• Surface reflection response data with proper pre-processing (source deconvolution, surface-multiple removal, and amplitude control).
• Sufficient illumination from the acquisition geometry, although recent advances significantly relax regular sampling requirements.

Challenges persist in handling highly ill-conditioned scenarios (e.g., strong evanescent regimes, high-dip geology), extending to elastic and anisotropic media, and integrating robust data pre-processing in irregular field geometries. Nevertheless, the methodology continues to evolve, providing a universally consistent and physically complete approach to redatuming and multiple elimination for advanced seismic imaging and time-lapse studies (Hajjaj et al., 2022, Wapenaar et al., 2020, Wapenaar, 2020, Ravasi et al., 2020, Wapenaar et al., 2023).