The Marchenko algorithm can suppress the disturbing effects of internal multiples that are present in seismic reflection data. To achieve this, a set of coupled equations with four unknowns is solved. These coupled equations are separated into a set of two equations with two unknowns using a time window. The two unknown focusing functions can be resolved by an iterative or direct method. These focusing functions, when applied to reflection data, create virtual point-sources inside the medium. Combining individual virtual point-sources into a plane-wave leads to an efficient computation of images without internal multiples. In this study the internal multiples are eliminated in a redatuming step which is part of the imaging algorithm. To use the Marchenko algorithm with plane-wave focusing functions, the time window that separates the unknowns must be adapted. The design of the plane-wave Marchenko algorithm is explained and illustrated with numerically modeled and measured reflection data.

@en

Recently, the focus of reflection seismologists has shifted to applications where a high-resolution image of the subsurface is required. Least-Squares Reverse-Time Migration (LSRTM) is a common tool used to compute such images. Still, its high computational costs have led seismologists to use target-oriented LSRTM for imaging only a small target of interest within a larger subsurface block. Redatuming the data to the upper boundary of the target of interest is one approach to target-oriented LSRTM. Still, many redatuming methods cannot account for multiple scattering within the overburden. We apply a target-oriented least-squares reverse time migration algorithm that integrates Marchenko redatuming and double-focusing to a field dataset. This redatuming method accounts for all orders of multiple scattering in the overburden, thus improving the accuracy of target-oriented LSRTM. Moreover, we demonstrate the effectiveness of a double-focusing algorithm in reducing the data size by decreasing both spatial and temporal dimensions of the model and the data. The algorithm's performance is evaluated using field data acquired in the Norwegian Sea. The numerical results show that our target-oriented LSRTM algorithm can reduce the internal multiple effects and increase the resolution of the resulting image.@en

The Marchenko method is capable of estimating Green’s functions between the surface of the Earth and arbitrary locations in the subsurface. These Green’s functions are used to redatum wavefields to a deeper level in the subsurface. The Marchenko method enables the isolation of the response of a specific layer or package of layers, free from the influence of the overburden and underburden. In this study, we apply the Marchenko-based isolation technique to land S-wave seismic data acquired in the Groningen province, the Netherlands. We apply the technique for combined elimination of the overburden and underburden. Our results indicate that this approach enhances the resolution of reflection data. These enhanced reflections can be utilised for imaging and monitoring applications.@en

Waves in space-dependent and in time-dependent materials obey similar wave equations, with interchanged time- and space-coordinates. However, since the causality conditions are the same in both types of material (i.e., without interchangement of time- and space-coordinates), the solutions are dissimilar. We present a systematic treatment of wave propagation and scattering in 1D space-dependent and in 1D time-dependent materials. After formulating unified equations, we discuss Green's functions and simple wave field representations for both types of material. Next we discuss propagation invariants, i.e., quantities that are independent of the space coordinate in a space-dependent material (such as the net power-flux density) or of the time coordinate in a time-dependent material (such as the net field-momentum density). A discussion of general reciprocity theorems leads to the well-known source–receiver reciprocity relation for the Green's function of a space-dependent material and a new source–receiver reciprocity relation for the Green's function of a time-dependent material. A discussion of general wave field representations leads to the well-known expression for Green's function retrieval from the correlation of passive measurements in a space-dependent material and a new expression for Green's function retrieval in a time-dependent material. After an introduction of a matrix–vector wave equation, we discuss propagator matrices for both types of material. Since the initial condition for a propagator matrix in a time-dependent material follows from the boundary condition for a propagator matrix in a space-dependent material by interchanging the time- and space-coordinates, the propagator matrices for both types of material are interrelated in the same way. This also applies to representations and reciprocity theorems involving propagator matrices, and to Marchenko-type focusing functions.

@en

Coda-wave interferometry employs the sensitivity of multiply scattered waves to detect minute changes of the propagation velocity. In most applications the underlying assumption is that the velocity changes take place approximately uniformly in a large region. Time-lapse seismic reservoir monitoring, on the other hand, aims to infer local time-lapse changes from seismic reflection measurements at the surface. Hence, coda-wave interferometry, applied directly to time-lapse seismic data for reservoir monitoring, cannot infer such local changes in a reservoir. Recently, we developed a Marchenko-based methodology to isolate the response of a target zone from the full reflection response. We apply this methodology to Troll-field data and successfully predict the internal multiples of the target zone. This confirms that in principle the Marchenko-based isolation method can be used to prepare time-lapse seismic data for coda-wave interferometry, which in turn can be used to infer time-lapse changes in a reservoir caused by production or storage of resources in the subsurface.@en

Geophysical monitoring of subsurface reservoirs relies on detecting small changes in the seismic response between a baseline and monitor study. However, internal multiples, related to the over- and underburden, can obstruct the view of the target response, hence complicating the time-lapse analysis. To retrieve a response that is free from the over- and underburden effects, the data-driven Marchenko method is used. This method effectively isolates the target response, which can then be used to extract more precise time-lapse changes. In addition, the method also reveals target-related multiples that probe the reservoir more than once, which further defines the changes in the reservoir. To verify the effectiveness of the method, a numerical example is constructed. This test finds that, when using the isolated target response, the observed time differences resemble the expected time differences in the reservoir. Moreover, the results obtained with target-related multiples also benefit from the Marchenko-based isolation of the reservoir. It is, therefore, concluded that this method has the potential to observe dynamic changes in the subsurface with increased accuracy.

@en

Minimum-phase properties are well-understood for scalar functions where they can be used as physical constraint for phase reconstruction. Existing scalar applications of the latter in geophysics include, for example the reconstruction of transmission from acoustic reflection data, or multiple elimination via the augmented acoustic Marchenko method. We review scalar minimum-phase reconstruction via the conventional Kolmogorov relation, as well as a less-known factorization method. Motivated to solve practice-relevant problems beyond the scalar case, we investigate (1) the properties and (2) the reconstruction of minimum-phase matrix functions. We consider a simple but non-trivial case of 2 × 2 matrix response functions associated with elastodynamic wavefields. Compared to the scalar acoustic case, matrix functions possess additional freedoms. Nonetheless, the minimum-phase property is still defined via a scalar function, that is a matrix possesses a minimum-phase property if its determinant does. We review and modify a matrix factorization method such that it can accurately reconstruct a 2 × 2 minimum-phase matrix function related to the elastodynamic Marchenko method. However, the reconstruction is limited to cases with sufficiently small differences between P- and S-wave traveltimes, which we illustrate with a synthetic example. Moreover, we show that the minimum-phase reconstruction method by factorization shares similarities with the Marchenko method in terms of the algorithm and its limitations. Our results reveal so-far unexplored matrix properties of geophysical responses that open the door towards novel data processing tools. Last but not least, it appears that minimum-phase matrix functions possess additional, still-hidden properties that remain to be exploited, for example for phase reconstruction.

@en

Green’s functions and propagator matrices are both solutions of the wave equation, but whereas Green’s functions obey a causality condition in time (G = 0 for t < 0), propagator matrices obey a boundary condition in space. Marchenko-type focusing functions focus a wave field in space at zero time. We discuss the mutual relations between Green’s functions, propagator matrices and focusing functions, avoiding up-down decomposition and accounting for propagating and evanescent waves. We conclude with discussing a Marchenko-type Green’s function representation, which forms a basis for extending the Marchenko method to improve the imaging of steeply dipping flanks and to account for refracted waves. @en

Marchenko-type integrals typically relate so-called focusing functions and Green's functions via the reflection response measured on the open surface of a volume of interest. Originating from one dimensional inverse scattering theory, the extension to two and three dimensions set in motion various new developments regarding imaging in complex materials. This extension, however, is based on wavefield decomposition inside the volume and a truncated medium state, i.e. a version of the medium that is reflection-free underneath the focusing location, suggesting that evanescent, refracted and diving waves cannot be included in the representation. We elaborate on a new derivation for Marchenko-like integrals that (i) extends the concept of wavefield focusing by using a generalised homogeneous Green's function, (ii) is based on partial differential equations and thereby allows for additional insights and a new physical intuition for Marchenko equations, (iii) unifies wavefield focusing for open and closed boundary systems, (iv) does not require wavefield decomposition or a truncated medium state, thus including the full wavefield Green's function, (v) enables using forward modelling to obtain, e.g., Marchenko-type, time-compact focusing functions. We place a particular focus on the latter point, illustrating and investigating how to solve the underlying partial differential equations for various types of focusing functions. This paves the way for a deeper understanding of focusing functions as well as advanced full wavefield Marchenko schemes. While the derivations are generally presented for the 3D case, we show numerical examples in 1D.

@en

Reservoir simulations for subsurface processes play an important role in successful deployment of geoscience applications such as geothermal energy extraction and geo-storage of fluids. These simulations provide time-lapse dynamics of the coupled poromechanical processes within the reservoir and its over-, under-, and side-burden environments. For more reliable operations, it is crucial to connect these reservoir simulation results with the seismic surveys (i.e., observation data). However, despite being crucial, such integration is challenging due to the fact that the reservoir dynamics alters the seismic parameters. In this work, a coupled reservoir simulation and time-lapse seismic methodology is developed for multiphase flow operations in subsurface reservoirs. To this end, a poromechanical simulator is designed for multiphase flow and connected to a forward seismic modeller. This simulator is then used to assess a novel methodology of seismic monitoring by isolating the reservoir signal from the entire reflection response. This methodology is shown to be able to track the development of the fluid front over time, even in the presence of a highly reflective overburden with strong time-lapse variations. These results suggest that the proposed methodology can contribute to a better understanding of fluid flow in the subsurface. Ultimately, this will lead to improved monitoring of reservoirs for underground energy storage or production.

@en

For the elastodynamic wave equation discussed in Appendix A.4 in Ref. 1, the expressions for matrices ∼L 6 1 and ∼L 6 2 in Eqs. (A32) and (A33) must be multiplied by 61. In other words, the signs of ∼L _ ₁and ∼L _ ₂have to be changed, whereas the signs of ∼L 1 and ∼L 2 remain unchanged. With these corrections, matrix ∼L , defined in Eq. (A5), fulfills the symmetry property formulated by Eq. (A6). The corrections have no consequences for Eqs. (A36)-(A38).

@en

Previous studies indicate that scattering may pose a trade-off for the performance of seismic interferometry (SI) applications for retrieving body-wave reflections of a target reflector. While it has been demonstrated that a higher scattering strength of the overburden improves the Green's function estimated by cross-correlation SI, other theoretical and empirical studies showed that multiple scattering also gives rise to more artefacts. The implications of this trade-off are analysed in this numerical study for a lithospheric scenario with varying crustal scattering strength and passive illumination conditions. In this scenario, we apply SI by cross-correlation to elastodynamic responses to double-couple sources to reconstruct virtual Moho primary reflections. We include multidimensional deconvolution (MDD) methods in the analysis to investigate whether scattering-induced artefacts affect MDD methods in a similar way as was shown for the cross-correlation method. Our results show that there indeed exists a trade-off between the quality of the virtual primary reflection of the target that can be obtained by SI and the scattering strength of the overburden. Furthermore, we find that the full-field MDD method proves to be most resilient to the negative effects of multiple scattering for all illumination conditions and scattering strengths analysed.

@en

Many seismic imaging methods use wavefield extrapolation operators to redatum sources and receivers from the surface into the subsurface. We discuss wavefield extrapolation operators that account for internal multiple reflections, in particular propagator matrices, transfer matrices and Marchenko focusing functions. A propagator matrix is a square matrix that 'propagates' a wavefield vector from one depth level to another. It accounts for primaries and multiples and holds for propagating and evanescent waves. A Marchenko focusing function is a wavefield that focuses at a designated point in space at zero time. Marchenko focusing functions are useful for retrieving the wavefield inside a heterogeneous medium from the reflection response at its surface. By expressing these focusing functions in terms of the propagator matrix, the usual approximations (such as ignoring evanescent waves) are avoided. While a propagator matrix acts on the full wavefield vector, a transfer matrix (according to the definition used in this paper) 'transfers' a decomposed wavefield vector (containing downgoing and upgoing waves) from one depth level to another. It can be expressed in terms of decomposed Marchenko focusing functions. We present propagator matrices, transfer matrices and Marchenko focusing functions in a consistent way and discuss their mutual relations. In the main text we consider the acoustic situation and in the appendices we discuss other wave phenomena. Understanding these mutual connections may lead to new developments of Marchenko theory and its applications in wavefield focusing, Green's function retrieval and imaging.

@en

The data-driven Marchenko method is able to redatum wavefields to arbitrary locations in the subsurface, and can, therefore, be used to isolate zones of specific interest. This creates a new reflection response of the target zone without interference from over- or underburden reflectors. Consequently, the method is well suited to obtain a clear response of a subsurface reservoir, which can be advantageous in time-lapse studies. The isolated responses of a baseline and monitor survey can be more effectively compared; hence, the retrieval of time-lapse characteristics is improved. This research aims to apply Marchenko-based isolation to a time-lapse marine data set of the Troll field in Norway in order to acquire an unobstructed image of the primary reflections and retrieve small time-lapse traveltime difference in the reservoir. It is found that the method not only isolates the primary reflections but can also estimate internal multiples outside the recording time. Both the primaries and the multiples can then be utilized to find time-lapse traveltime differences. More accurate ways of time-lapse monitoring will allow for a better understanding of dynamic processes in the subsurface, such as observing saturation and pressure changes in a reservoir or monitoring underground storage of hydrogen and CO₂.

@en

Least-squares reverse-time migration (LSRTM) is a method that seismologists utilize to compute a high-resolution subsurface image. Nevertheless, LSRTM is a computationally demanding problem. One way to reduce the computational costs of the LSRTM is to choose a small region of interest (ROI) and compute the image of that region. However, finding representations that account for the wavefields entering the target region from the surrounding boundaries is necessary. This article confines the ROI between two boundaries above and below this region. The acoustic reciprocity theorem is employed to derive representations for the wavefields at the upper and lower boundaries of the target region. With the help of these representations, a target-enclosed LSRTM algorithm is developed to compute a high-resolution image of the ROI. Moreover, the possibility of using virtual receivers created by Marchenko redatuming is investigated.@en

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.

@en

In a surface-seismic setting, Distributed Acoustic Sensing (DAS) is still not a widely adopted method for near-surface characterisation, especially for reflection seismic imaging, despite the dense spatial sampling it provides over long distances. This is mainly due to the decreased broadside sensitivity that DAS suffers from when buried horizontally in the ground (that is when the upgoing wavefield (e.g. reflected wavefield) is perpendicular to the optical fibre). This is unlike borehole settings (e.g. zero-offset Vertical Seismic Profiling), where DAS has been widely adopted for many monitoring applications. Advancements in the field, like shaping the fibre to a helix, commonly known as helically wound fibre, allow better sensitivity for the reflections. The promise of spatially dense seismic data over long distances is an attractive prospect for retrieving the local variations of near-surface properties. This is particularly valuable for areas impacted by induced seismicity, as is the case in the Groningen Province in the north of The Netherlands, where near-surface properties, mostly composed of clays and peats, play an essential role on the amount of damage on the very near-surface and the structures built on it. Installing fibre-optic cables for passive and active measurements is valuable in this situation. We installed multiple cables containing different fibre configurations of straight and helically wound fibres, buried in a 2-m deep trench. The combination of the different fibre configurations allows us to obtain multi-component information. We observe differences in the amplitude and phase information, suggesting that these differences can be used for separating the different components of the wave motion. We also see that using enhanced backscatter fibres, reflection images can be obtained for the helically wound fibre as well as the straight fibre, despite the decreased broadside sensitivity for the latter. @en

Standard Marchenko redatuming and imaging schemes neglect evanescent waves and are based on the assumption that decomposition into downgoing and upgoing waves is possible in the subsurface. Recently we have shown that propagator matrices, which circumvent these assumptions, can be expressed in terms of Marchenko focusing functions. In this paper we generalize the relation between the propagator matrix and the Marchenko focusing functions for a 3D inhomogeneous dissipative medium. Moreover, for the same type of medium we discuss a relation between the transfer matrix and the Marchenko focusing functions.@en

Green’s functions in an unknown elastic layered medium can be retrieved from single-sided reflection data by solving a Marchenko equation. This methodology requires a priori knowledge of all forward-scattered (non-converted and converted) waveforms. Moreover, the medium should satisfy stringent monotonicity conditions, which are often not met in realistic scenarios. In this contribution, we show that the situation is significantly less cumbersome if two-sided reflection and transmission data are recorded (for instance in laboratory settings). A novel methodology is presented to retrieve elastodynamic Green’s functions from such data. Apart from the two-sided reflection and transmission responses, our methodology requires knowledge of the direct non-converted PP- and SS-transmissions (a priori knowledge of forward-scattered converted waveforms is not needed). We demonstrate the success of our methodology by conducting a numerical experiment in an elastic layered medium that violates the monotonicity conditions of the Marchenko equation for single-sided reflection data. The limitations of the methodology and the sensitivity to errors in our initial estimates require further investigation.@en

Seismic interferometry (SI) is a method that retrieves new seismic traces from the cross-correlation of existing traces, where one of the receivers acts as a virtual seismic source whose response is retrieved at other receivers. When using sources only at the surface, and so-called one-sided illumination of the receivers occurs, not only desired physical reflections are retrieved, but also non-physical (ghost) reflections. These non-physical reflections are caused by internal reflections inside subsurface layers. They are thus particularly interesting to use for monitoring changes in the specific subsurface layer that causes them to appear in the SI result because they could provide valuable information about the physical properties of the subsurface. We illustrate the potential of SI with active-source data from numerical acoustic modelling using the Groningen subsurface model. This model describes the natural gas field located in the Groningen province in the northeastern part of the Netherlands. The reservoir of the Groningen gas field is located at depths between 2600 m and 3200 m, the total thickness ranges from approximately 100 m to 300 m. The Groningen field is cut by several fault systems, subdividing the field into a large number of fault blocks, and it is a clear example of induced seismicity by gas production. We investigate the utilization of non-physical reflections retrieved from surface active-source data using SI by cross-correlation and auto-correlation. With multi-offset gathers, besides physical reflections, we retrieve non-physical reflections as well; by muting undesired reflections, we can retrieve better target-related non-physical reflections. To illustrate the potential of the non-physical reflections for monitoring purposes, we apply velocity changes in the Groningen reservoir. With zero-offset gathers, which are retrieved from SI by auto-correlation, we show that in case of velocity changes, the non-physical reflections show a clear change; furthermore, they show a good agreement with the geometry of specific subsurface layers, specifically with the faulted structure. Thus, we can utilize non-physical reflections for imaging and monitoring in the Groningen reservoir.@en