We implement the 3D Marchenko equations to retrieve responses to virtual sources inside the subsurface. For this, we require reflection data at the surface of the Earth that contain no free-surface multiples and are densely sampled in space. The required 3D reflection data volume is very large and solving the Marchenko equations requires a significant amount of computational cost. To limit the cost, we apply floating point compression to the reflection data to reduce their volume and the loading time from disk. We apply the Marchenko implementation to numerical reflection data to retrieve accurate Green's functions inside the medium and use these reflection data to apply imaging. This requires the simulation of many virtual source points, which we circumvent using virtual plane-wave sources instead of virtual point sources. Through this method, we retrieve the angle-dependent response of a source from a depth level rather than of a point. We use these responses to obtain angle-dependent structural images of the subsurface, free of contamination from wrongly imaged internal multiples. These images have less lateral resolution than those obtained using virtual point sources, but are more efficiently retrieved.

Seismic interferometry is a method used to calculate wavefields for sources and receivers that are located where only sources or only receivers are available. There are correlation- or deconvolution-based interferometric methods that can be used to reposition the seismic array from the earth's surface to an arbitrary datum at depth. Based on the one-way reciprocity theorems of convolution and correlation type, we have determined that interferometric redatuming can be achieved in a deconvolution-only procedure in three steps. The first two steps consist of separately retrieving, for sources at the earth's surface, the downward- and upward-propagating Green's functions at receivers at the datum, which are then used in the third step to reposition the sources to the datum. For the involved deconvolutions, transmitted and backscattered wavefields need to be modeled with a velocity model between the acquisition and datum levels. Our numerical experiments demonstrate that the method can help to reduce nonphysical events and other artifacts that commonly arise in purely correlation-based procedures. If a high-quality overburden-velocity model is available, it correctly accounts for inhomogeneities in the overburden medium. Because the method's sensitivity to the velocity model is mainly introduced by backscattering at overburden heterogeneities, a smooth model is sufficient when overburden scattering is weak.