An analysis of acquisition-related subsampling effects on Marchenko focusing, redatuming, and primary estimation

Marchenko methods can retrieve Green’s functions and focusing functions from single-sided reflection data and a smooth velocity model, as essential components of a redatuming process. Recent studies also indicate that a modified Marchenko scheme can reconstruct primary-only reflection responses directly from reflection data without requiring a priori model information. To provide insight into the artifacts that arise when input data are not ideally sampled, we study the effects of subsampling in both types of Marchenko methods in 2D earth and data — by analyzing the behavior of Marchenko-based results on synthetic data subsampled in sources or receivers. With a layered model, we find that for Marchenko redatuming, subsampling effects jointly depend on the choice of integration variable and the subsampling dimension, originated from the integrand gather in the multidimensional convolution process. When reflection data are subsampled in a single dimension, integrating on the other yields spatial gaps together with artifacts, whereas integrating on the subsampled dimension produces aliasing artifacts but without spatial gaps. Our complex subsalt model indicates that the subsampling may lead to very strong artifacts, which can be further complicated by having limited apertures. For Marchenko-based primary estimation (MPE), subsampling below a certain fraction of the fully sampled data can cause MPE iterations to diverge, which can be mitigated to some extent by using more robust iterative solvers, such as least-squares QR. Our results, covering redatuming and primary estimation in a range of subsampling scenarios, provide insights that can inform acquisition sampling choices as well as processing parameterization and quality control, e.g., to set up appropriate data filters and scaling to accommodate the effects of dipole fields, or to help ensuring that the data interpolation achieves the desired levels of reconstruction quality that minimize subsampling artifacts in Marchenko-derived fields and images. INTRODUCTION By means of Marchenko methods, the full-waveform Green’s function to a virtual source in the subsurface can be retrieved from single-sided reflection data. Given a smooth velocity model (Wapenaar et al., 2014), this is accomplished by solving a set of coupled Marchenko equations by direct least-squares (LSQR) inversion (van der Neut et al., 2015a) or iterative substitution (van der Neut et al., 2015b). The method, in theory, retrieves internal multiples with correct amplitudes. One of the challenges of this method is that the multidimensional convolution operation of the reflection response with the focusing function requires a regularly sampled, dense data set with colocated sources and receivers. Similar multidimensional convolution operations are also used — and fairly well-understood — in surface-related multiple elimination (SRME) (Verschuur et al., 1992; Lopez and Verschuur, 2015). Here, we look into the effects of regular subsampling in the Marchenko method specifically — showing similarities and differences to the effects seen in SRME applications. Haindl et al. (2018) formulate the Marchenko equations with sparse inversion to alleviate the 2D irregular source subsampling effect. This study observed spatial gaps in the output redatumed Green’s function and the focusing function. Ravasi (2017) also studies the effect of irregular subsampling in sources Manuscript received by the Editor 28 December 2020; revised manuscript received 29 May 2021; published ahead of production 19 July 2021; published online 1 September 2021. Utrecht University, 3584 CB Utrecht, The Netherlands. E-mail: h.peng@uu.nl (corresponding author); i.vasconcelos@uu.nl; yanadet.sripanich@gmail.com. Delft University of Technology, 2628 CN Delft, The Netherlands. E-mail: l.zhang-1@tudelft.nl. © 2021 The Authors. Published by the Society of Exploration Geophysicists. All article content, except where otherwise noted (including republished material), is licensed under a Creative Commons Attribution 4.0 Unported License (CC BY-NC). See http://creativecommons.org/licenses/by/4.0/. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its digital object identifier (DOI). Commercial reuse is not permitted. The same license does not have to be used for derivative works. WC75 GEOPHYSICS, VOL. 86, NO. 5 (SEPTEMBER-OCTOBER 2021); P. WC75–WC88, 18 FIGS., 1 TABLE. 10.1190/GEO2020-0914.1 Downloaded from http://pubs.geoscienceworld.org/geophysics/article-pdf/86/5/WC75/5435069/geo-2020-0914.1.pdf by guest on 16 January 2022 in the Rayleigh-Marchenko method with LSQR inversion and observed artifacts, but no spatial gaps in the retrieved signals, thereby showing the Rayleigh approach’s flexibility in handling the source subsampling. Sripanich and Vasconcelos (2019) provide a comprehensive study on the effects of a limited aperture on Marchenko focusing functions and their radiation behavior of focusing. In synthetic and field data studies, respectively, Jia et al. (2019) and Staring andWapenaar (2020) present 3DMarchenko redatuming, in which they apply data interpolation to fill the gaps caused by subsampling. More recently, in a companion study to this paper, regular and irregular subsampling effects on 3D redatuming are discussed by Ravasi and Vasconcelos (2020), in the context of high-performance computing. To account for the imperfect source sampling, Wapenaar and van IJsseldijk (2020) develop a new discrete representation containing point-spread functions that can be multidimensionally deconvolved to deblur the focusing functions and Green’s functions. Van IJsseldijk and Wapenaar (2020) then alter the iterative Marchenko scheme to integrate this new representation and test it on synthetic data. All of the preceding studies alleviate the subsampling effects using different approaches, including sparse inversion, LSQR inversion, and/or data interpolation. Yet, there is still a need to build a basic understanding of these effects’ origin and patterns, which could not only inform acquisition sampling choices in the future but also crucially help in the design and quality control of preprocessing routines tailored toward Marchenko methods. Although the original form of the Marchenko system is designed to perform redatuming, i.e., retrieve subsurface wavefields, the framework can also be used for surface-to-surface internal multiple prediction — or more specifically — to predict primary-only data. To that end, an alternative, surface-projected form of the Marchenko equations presented by van der Neut and Wapenaar (2016) was recently used to eliminate multiples directly from reflection data without requiring model information (Zhang and Staring, 2018). In a comparative study, Zhang et al. (2019a) illustrate the benefits and shortcomings of Marchenko-based data-driven scheme for the internal multiple reflection elimination over the theory of inverse scattering series. The former method was further improved to compensate the transmission loss (Zhang et al., 2019b). In this context, i.e., Marchenko-based primary estimation (MPE), the effects of acquisition sampling remain largely unaddressed. Very recently, this projected scheme was also used in the augmented Marchenko method to account for short-period internal multiples in media with a horizontally layered overburden by Elison et al. (2020). In recent years, Marchenko-based methods have been developed in two related, yet distinct, contexts: depth-domain redatuming and data-domain primary estimation. In this paper, we choose to look into Marchenko approaches in both of these contexts because they typically pertain to different applications — namely, redatuming is tied to full-wavefield imaging at depth, whereas primary estimation is usually a preprocessing step prior to model building and imaging. To this end, building on the aforementioned studies that look into acquisition aspects of the Marchenko approaches in specific contexts, we intend to provide a basic understanding of acquisition sampling effects in focusing, redatuming, and primary estimation. In addition, two types of Marchenko methods are addressed in this paper because they require different inputs and temporal-windowing approaches. Albeit relying on the same general operations of multidimensional convolution and correlation, these two types of Marchenko approaches can behave differently in the presence of subsampling. There are two main reasons for this: (1) Marchenko redatuming does not perform time stepping and uses the full wavefield that has been retrieved by focusing for a given redatuming level, and primary estimation applies time stepping and just selects the samples at that time after applying the Marchenko method at a given time step and (2) the windowing operators themselves in Marchenko redatuming and primary estimation are different, leading to differences in the physical contributions of the convolution processes in either of the applications. We observe that subsampled input data cause artifacts for both methods and may further lead to divergence for the primary estimation, depending on the degree of subsampling and the algorithm used to solve the system, i.e., Neumann series expansion or LSQR inversion. Because we wish to analyze subsampling at varying levels of aliasing, we maintain a fixed bandwidth and only vary the sampling. In practice, sampling and aliasing are relative to the actual data bandwidth and, as a general rule, adequately sampling the integrals requires satisfying the general spatial Nyquist criteria for the integrands — similarly to well-known applications such as SRME or evaluating Kirchhoff-type integrals. Our study of the subsampling effects is structured in a simple, yet general, framework that could inform acquisition design or help in tuning and controlling data processing. With this goal in mind, we limit ourselves to 2D regular subsampling in just one dimension (source or receiver). We believe that the insight of this case shares similarities with the other cases, e.g., irregular subsampling in both dimensions. Because three dimensions and data interpolation are related, yet complex, topics, they are left for future research and therefore are not included in this study. First, we look at subsampling numerical examples in two dimensions to mimic subsampling caused by realistic acquisition design. Then, we analyze the results of the aforementioned two Marchenko methods in the context of subsampling without data interpolation, also by looking at specific yet simple measures to alleviate the effects of subsampling, such as invoking spatial reciprocity. In all of the current forms of Marchenko methods — whether in the context of redatuming or primary estimation — the theory clearly outlines that monopole (e.g., pressure-to-pressure) and dipole (e.g., particle velocity-to-pressure) fields are not only required but must be used in specific source and receiver configurations (Wapenaar et al., 2014; van der Neut et al., 2015b). When it comes to field data, e.g., in offshore acquisition settings, although monopole fields are easily available from streamer and ocean-bottom acquisitions, dipole fields are usually only available on the receiver side in the form of acceleration or velocity measurements from multicomponent sensors. Thus, in this paper, we show the effects of subsampling in the presence of dipole sources — which are generally not available in practice — as well as dipole receivers. We do this to provide insights on the importance and effects of dipole integration choices that may impact choices on how to use multicomponent sensor data for Marchenko applications or how to preprocess 1C data to be used in Marchenko schemes. The layout of this paper is as follows. We begin by presenting the Marchenko integral representations, while paying particular attention to the discrete form of the terms within the iterative solution of the coupled Marchenko system, followed by a brief review of the primary estimation method (Zhang and Staring, 2018). With a layered benchmark model, we then illustrate our observations on the subsampling effects on Marchenko focusing for multiple scenarios, including WC76 Peng et al. Downloaded from http://pubs.geoscienceworld.org/geophysics/article-pdf/86/5/WC75/5435069/geo-2020-0914.1.pdf by guest on 16 January 2022 the joint effects of subsampling versus the dipole source/receiver integration variable choice and the effects of focusing depth jointly with subsampling. Then, a complex subsalt model is used to illustrate subsampling effects together with limited aperture on Marchenko redatuming in a more realistic scenario. Subsequently, we show with numerical examples for both Marchenko systems that integrating on the monopole dimension instead of the (theoretically correct) dipole dimension actually yields negligible errors compared with the true results for media with mild lateral heterogeneity. Finally, for the primary estimation, we compare the joint effects of different subsampling schemes versus different approaches to solve the Marchenko system with a layered model and provide an operator norm analysis for the iterative substitution method. THEORY OF MARCHENKO REDATUMING AND PRIMARY ESTIMATION Discrete integration in the Marchenko framework According to Wapenaar et al. (2014), the coupled Marchenko representations for the acoustic medium can be defined in the frequency domain by integrating on the dipole-source dimension with the receiver being monopole: ĜðxF;xrÞ1⁄4 Z ΛR dxsR̂ðxr;xsÞf̂1 ðxs;xFÞ− f̂1 ðxr;xFÞ; (1) Ĝþ ðxF; xrÞ 1⁄4 − Z ΛR dxsR̂ ðxr; xsÞf̂1 ðxs; xFÞ þ f̂1 ðxr; xFÞ; (2) R̂ðxr; xsÞ 1⁄4 2 jωρðxsÞ ∂Ĝðxr; xsÞ ∂zs ; (3) where xF, xr, and xs are the spatial coordinates of the focal point, receiver, and source positions, respectively. The circumflex accent denotes the wavefields in the frequency domain. The terms f̂1 and f̂1 are the downand upgoing focusing functions, correspondingly, Ĝ and Ĝþ are the upgoing and the time-reversed downgoing Green’s functions, respectively, ΛR denotes an open boundary that corresponds to the surface acquisition level, and ∂zs represents the spatial derivative along the depth at the source position. As such, R̂ðxr; xsÞ is twice the pressure recorded at xr due to a vertical dipole source at xs, Ĝðxr; xsÞ is the pressure response from a monopole source, j is the imaginary unit, ω is the angular frequency, and ρðxsÞ is the mass density at the source position. Given only reflection surface data as a known input, the system in equations 1 and 2 can only be solved for the unknown focusing functions after using additional constraints. The more established version of the Marchenko redatuming scheme (e.g., van der Neut et al., 2015b) accomplishes this by introducing a known initial focusing function and a causality windowing operator leading to a discrete version of the Marchenko system that can be solved for the focusing functions: I 0 0 I − 0 WR WR⋆ 0 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} M f1 fþm 1⁄4 0 WR WR⋆ 0 0 f1d ; (4) where the focusing functions are discrete vectors corresponding to a fixed focal point and all points on the acquisition surface: fþm is the coda of the downgoing focusing function, f− is its upgoing response, and f1d is the input initial focusing function. Here, following van der Neut et al. (2015b), R and R⋆ are the discrete integral operators that apply multidimensional convolution and correlation in the time domain between the reflection response and the focusing functions, respectively. The term W is the operator enforcing causality by windowing in the time domain. Note that the Marchenko system in equation 4 has the same form either in the time or frequency domain — we hereafter refer to equations in the frequency domain for the sake of convenience. The system in equation 4 can be solved for fþm and f−, by using a Neumann series expansion of theM operator (under appropriate conditions, as described in Dukalski and de Vos, 2018), the so-called Marchenko iterative scheme is obtained. Once the focusing functions are solved for, the redatumed wavefields Ĝ and Ĝþ are then obtained by substituting the estimated focusing functions into the original system in equations 1 and 2. For our purposes, to facilitate the analysis of integration-related effects, we focus on the first operation Rf̂1d of the iterative substitution method in the form of matrix-vector multiplication (corresponding to the discrete spatial convolutions): 2 64 f̂1;K1⁄40ðx r ;xFÞ .. . f̂1;K1⁄40ðx r ;xFÞ 3 751⁄4 2 64 PNs i1⁄41 R̂ðxð0Þ r ;xðiÞ s Þf̂1dðx s ;xFÞdxs .. . PNs i1⁄41 R̂ðxðNrÞ r ;xðiÞ s Þf̂1dðx s ;xFÞdxs 3 75; (5) whereK 1⁄4 0 denotes the initial step of the iterative scheme, f̂1d is the initial focusing function, and Ns is the number of sources. An illustration of equation 5 is given in Figure 1a. Note that here we omit the action of the windowing operator because it does not affect the integrand within the spatial convolutions. This operation is interpreted as follows: The jth element on the left side corresponds to the trace at receiver position xðjÞ r , which is calculated by summing the convolution gather between R̂ðxðjÞ r ; xsÞ and f̂1dðxs; xFÞ for all the source positions and scaled by the (possibly subsampled) source interval. By using source-receiver reciprocity, the Marchenko equations can also be defined by integrating on the dipole-receiver dimension with the source being a monopole (van der Neut et al., 2015b): ĜðxF;xsÞ1⁄4 Z Λf dxrR̂ðxr;xsÞf̂1 ðxr;xFÞ− f̂1 ðxs;xFÞ; (6) Ĝþ ðxF; xsÞ 1⁄4 − Z Λf dxrR̂ ðxr; xsÞf̂1 ðxr; xFÞ þ f̂1 ðxs; xFÞ;


INTRODUCTION
By means of Marchenko methods, the full-waveform Green's function to a virtual source in the subsurface can be retrieved from single-sided reflection data. Given a smooth velocity model (Wapenaar et al., 2014), this is accomplished by solving a set of coupled Marchenko equations by direct least-squares (LSQR) inversion (van der Neut et al., 2015a) or iterative substitution (van der Neut et al., 2015b). The method, in theory, retrieves internal multiples with correct amplitudes. One of the challenges of this method is that the multidimensional convolution operation of the reflection response with the focusing function requires a regularly sampled, dense data set with colocated sources and receivers. Similar multidimensional convolution operations are also usedand fairly well-understood in surface-related multiple elimination (SRME) (Verschuur et al., 1992;Lopez and Verschuur, 2015). Here, we look into the effects of regular subsampling in the Marchenko method specificallyshowing similarities and differences to the effects seen in SRME applications. Haindl et al. (2018) formulate the Marchenko equations with sparse inversion to alleviate the 2D irregular source subsampling effect. This study observed spatial gaps in the output redatumed Green's function and the focusing function. Ravasi (2017) also studies the effect of irregular subsampling in sources Manuscript received by the Editor 28 December 2020; revised manuscript received 29 May 2021; published ahead of production 19 July 2021; published online 1 September 2021. in the Rayleigh-Marchenko method with LSQR inversion and observed artifacts, but no spatial gaps in the retrieved signals, thereby showing the Rayleigh approach's flexibility in handling the source subsampling. Sripanich and Vasconcelos (2019) provide a comprehensive study on the effects of a limited aperture on Marchenko focusing functions and their radiation behavior of focusing. In synthetic and field data studies, respectively, Jia et al. (2019) and Staring and Wapenaar (2020) present 3D Marchenko redatuming, in which they apply data interpolation to fill the gaps caused by subsampling. More recently, in a companion study to this paper, regular and irregular subsampling effects on 3D redatuming are discussed by Ravasi and Vasconcelos (2020), in the context of high-performance computing. To account for the imperfect source sampling, Wapenaar and van IJsseldijk (2020) develop a new discrete representation containing point-spread functions that can be multidimensionally deconvolved to deblur the focusing functions and Green's functions. Van IJsseldijk and Wapenaar (2020) then alter the iterative Marchenko scheme to integrate this new representation and test it on synthetic data. All of the preceding studies alleviate the subsampling effects using different approaches, including sparse inversion, LSQR inversion, and/or data interpolation. Yet, there is still a need to build a basic understanding of these effects' origin and patterns, which could not only inform acquisition sampling choices in the future but also crucially help in the design and quality control of preprocessing routines tailored toward Marchenko methods.
Although the original form of the Marchenko system is designed to perform redatuming, i.e., retrieve subsurface wavefields, the framework can also be used for surface-to-surface internal multiple predictionor more specificallyto predict primary-only data. To that end, an alternative, surface-projected form of the Marchenko equations presented by van der Neut and Wapenaar (2016) was recently used to eliminate multiples directly from reflection data without requiring model information (Zhang and Staring, 2018). In a comparative study, Zhang et al. (2019a) illustrate the benefits and shortcomings of Marchenko-based data-driven scheme for the internal multiple reflection elimination over the theory of inverse scattering series. The former method was further improved to compensate the transmission loss (Zhang et al., 2019b). In this context, i.e., Marchenko-based primary estimation (MPE), the effects of acquisition sampling remain largely unaddressed. Very recently, this projected scheme was also used in the augmented Marchenko method to account for short-period internal multiples in media with a horizontally layered overburden by Elison et al. (2020).
In recent years, Marchenko-based methods have been developed in two related, yet distinct, contexts: depth-domain redatuming and data-domain primary estimation. In this paper, we choose to look into Marchenko approaches in both of these contexts because they typically pertain to different applicationsnamely, redatuming is tied to full-wavefield imaging at depth, whereas primary estimation is usually a preprocessing step prior to model building and imaging. To this end, building on the aforementioned studies that look into acquisition aspects of the Marchenko approaches in specific contexts, we intend to provide a basic understanding of acquisition sampling effects in focusing, redatuming, and primary estimation.
In addition, two types of Marchenko methods are addressed in this paper because they require different inputs and temporal-windowing approaches. Albeit relying on the same general operations of multidimensional convolution and correlation, these two types of Marchenko approaches can behave differently in the presence of subsampling. There are two main reasons for this: (1) Marchenko redatuming does not perform time stepping and uses the full wavefield that has been retrieved by focusing for a given redatuming level, and primary estimation applies time stepping and just selects the samples at that time after applying the Marchenko method at a given time step and (2) the windowing operators themselves in Marchenko redatuming and primary estimation are different, leading to differences in the physical contributions of the convolution processes in either of the applications. We observe that subsampled input data cause artifacts for both methods and may further lead to divergence for the primary estimation, depending on the degree of subsampling and the algorithm used to solve the system, i.e., Neumann series expansion or LSQR inversion. Because we wish to analyze subsampling at varying levels of aliasing, we maintain a fixed bandwidth and only vary the sampling. In practice, sampling and aliasing are relative to the actual data bandwidth and, as a general rule, adequately sampling the integrals requires satisfying the general spatial Nyquist criteria for the integrandssimilarly to well-known applications such as SRME or evaluating Kirchhoff-type integrals.
Our study of the subsampling effects is structured in a simple, yet general, framework that could inform acquisition design or help in tuning and controlling data processing. With this goal in mind, we limit ourselves to 2D regular subsampling in just one dimension (source or receiver). We believe that the insight of this case shares similarities with the other cases, e.g., irregular subsampling in both dimensions. Because three dimensions and data interpolation are related, yet complex, topics, they are left for future research and therefore are not included in this study. First, we look at subsampling numerical examples in two dimensions to mimic subsampling caused by realistic acquisition design. Then, we analyze the results of the aforementioned two Marchenko methods in the context of subsampling without data interpolation, also by looking at specific yet simple measures to alleviate the effects of subsampling, such as invoking spatial reciprocity.
In all of the current forms of Marchenko methodswhether in the context of redatuming or primary estimationthe theory clearly outlines that monopole (e.g., pressure-to-pressure) and dipole (e.g., particle velocity-to-pressure) fields are not only required but must be used in specific source and receiver configurations (Wapenaar et al., 2014;van der Neut et al., 2015b). When it comes to field data, e.g., in offshore acquisition settings, although monopole fields are easily available from streamer and ocean-bottom acquisitions, dipole fields are usually only available on the receiver side in the form of acceleration or velocity measurements from multicomponent sensors. Thus, in this paper, we show the effects of subsampling in the presence of dipole sourceswhich are generally not available in practiceas well as dipole receivers. We do this to provide insights on the importance and effects of dipole integration choices that may impact choices on how to use multicomponent sensor data for Marchenko applications or how to preprocess 1C data to be used in Marchenko schemes.
The layout of this paper is as follows. We begin by presenting the Marchenko integral representations, while paying particular attention to the discrete form of the terms within the iterative solution of the coupled Marchenko system, followed by a brief review of the primary estimation method (Zhang and Staring, 2018). With a layered benchmark model, we then illustrate our observations on the subsampling effects on Marchenko focusing for multiple scenarios, including the joint effects of subsampling versus the dipole source/receiver integration variable choice and the effects of focusing depth jointly with subsampling. Then, a complex subsalt model is used to illustrate subsampling effects together with limited aperture on Marchenko redatuming in a more realistic scenario. Subsequently, we show with numerical examples for both Marchenko systems that integrating on the monopole dimension instead of the (theoretically correct) dipole dimension actually yields negligible errors compared with the true results for media with mild lateral heterogeneity. Finally, for the primary estimation, we compare the joint effects of different subsampling schemes versus different approaches to solve the Marchenko system with a layered model and provide an operator norm analysis for the iterative substitution method.

Discrete integration in the Marchenko framework
According to Wapenaar et al. (2014), the coupled Marchenko representations for the acoustic medium can be defined in the frequency domain by integrating on the dipole-source dimension with the receiver being monopole: (2) where x F , x r , and x s are the spatial coordinates of the focal point, receiver, and source positions, respectively. The circumflex accent denotes the wavefields in the frequency domain. The termsf þ 1 and f − 1 are the down-and upgoing focusing functions, correspondingly, G − andĜ þÃ are the upgoing and the time-reversed downgoing Green's functions, respectively, Λ R denotes an open boundary that corresponds to the surface acquisition level, and ∂z s represents the spatial derivative along the depth at the source position. As such, Rðx r ; x s Þ is twice the pressure recorded at x r due to a vertical dipole source at x s ,Ĝðx r ; x s Þ is the pressure response from a monopole source, j is the imaginary unit, ω is the angular frequency, and ρðx s Þ is the mass density at the source position. Given only reflection surface data as a known input, the system in equations 1 and 2 can only be solved for the unknown focusing functions after using additional constraints. The more established version of the Marchenko redatuming scheme (e.g., van der Neut et al., 2015b) accomplishes this by introducing a known initial focusing function and a causality windowing operator leading to a discrete version of the Marchenko system that can be solved for the focusing functions: where the focusing functions are discrete vectors corresponding to a fixed focal point and all points on the acquisition surface: f þ m is the coda of the downgoing focusing function, f − is its upgoing response, and f þ 1d is the input initial focusing function. Here, following van der Neut et al. (2015b), R and R ⋆ are the discrete integral operators that apply multidimensional convolution and correlation in the time domain between the reflection response and the focusing functions, respectively. The term W is the operator enforcing causality by windowing in the time domain. Note that the Marchenko system in equation 4 has the same form either in the time or frequency domainwe hereafter refer to equations in the frequency domain for the sake of convenience. The system in equation 4 can be solved for f þ m and f − , by using a Neumann series expansion of the M operator (under appropriate conditions, as described in Dukalski and de Vos, 2018), the so-called Marchenko iterative scheme is obtained. Once the focusing functions are solved for, the redatumed wavefieldŝ G − andĜ þ are then obtained by substituting the estimated focusing functions into the original system in equations 1 and 2.
For our purposes, to facilitate the analysis of integration-related effects, we focus on the first operation Rf þ 1d of the iterative substitution method in the form of matrix-vector multiplication (corresponding to the discrete spatial convolutions): where K ¼ 0 denotes the initial step of the iterative scheme,f þ 1d is the initial focusing function, and N s is the number of sources. An illustration of equation 5 is given in Figure 1a. Note that here we omit the action of the windowing operator because it does not affect the integrand within the spatial convolutions. This operation is interpreted as follows: The jth element on the left side corresponds to the trace at receiver position x ðjÞ r , which is calculated by summing the convolution gather betweenRðx ðjÞ r ; x s Þ andf þ 1d ðx s ; x F Þ for all the source positions and scaled by the (possibly subsampled) source interval.
By using source-receiver reciprocity, the Marchenko equations can also be defined by integrating on the dipole-receiver dimension with the source being a monopole (van der Neut et al., 2015b): whereRðx r ; x s Þ is twice the vertical velocity recording at x r due to a pressure source at x s (compare to equation 3). The term ∂z r represents the spatial derivative along the depth, and ρðx r Þ is the density, now both at the receiver position instead. Analogous to equation 5, the initial-iteration discrete convolution now has the form where N r is the number of receivers. As a consequence, when integrating over the receiver locations, the jth-element output in the obtained upgoing focusing function now corresponds to a trace at source position x ðjÞ s , which is calculated by integrating the convolution gatherRðx r ; x ðjÞ s Þ andf þ 1d ðx r ; x F Þ over the receiver positions, scaled by the (possibly subsampled) receiver interval. Figure 1b illustrates equation 9. Here, we rely on equations 5 and 9 to interpret the observed effects of subsampling in the source and receiver dimensions.
One key point to note, here as well as throughout this paper, is that we refer to monopole and dipole sources and receivers so as to relate the numerical results to the fields in the preceding equations. However, in real life, practical acquisition, dipole fields are not necessarily always available. For example, in an offshore acquisition setting, sources are generally of the monopole type onlythis means that to obtain dipole sources for calculations one would need to apply, e.g., local obliquity factor corrections. On the receiver side, however, dipole observations correspond to either direct particle velocity or acceleration measurements, which are generally present in ocean-bottom systems and available in some multicomponent streamer systems.

Marchenko-based primary estimation
With the objective of retrieving primary-only data from internalmultiple interference, in addition to redatuming, van der Neut and Wapenaar (2016) propose that the coupled Marchenko equations could be implicitly projected to the acquisition surface by the direct wavefield, such that they can be solved without the need of a macrovelocity model. Based on this approach, they eliminate the primaries and multiples for the preceding medium a chosen depth level by adaptive subtraction. The projected focusing functions are defined in the frequency domain as follows:v where x P is the projection point at the acquisition surface, x F is the focusing point, and Λ f denotes the focusing level. The term G þ d is the direct wavefield, which does not need to be known a priori, v − is the projected upgoing focusing function, and v þ m is the coda of the projected downgoing focusing function. Likewise, the projected Green's functions are defined aŝ whereÛ − andÛ þ are the projected up-and downgoing Green's functions, respectively. The coupled Marchenko system after projection iŝ whereRðx r ; x s Þ is defined in equation 8. The band-limited 2D delta function δðx r − x P Þ in space and frequency appears due to that f þ 1d is defined as the inverse of the direct arrival G þ deven though in the practical redatuming scheme it is generally approximated by the time-reversed G þ d . Relying on the same projection in the time domain, Zhang and Staring (2018) propose to estimate multiple-free primaries directly from the input reflection data without prior model information or relying on adaptive subtraction. The discrete projected Marchenko equations in their theory are where W t 2 −ϵ ϵ is the windowing operator that removes all of the signals outside the interval ðϵ; t 2 − ϵÞ, where ϵ is, e.g., the half-length of the source wavelet and t 2 is the two-way traveltime between the acquisition surface and the The same for equation 9, in which the integration is conducted over the dipole receiver dimension, but with the source being monopole. For each frequency slice,R andR T are the 2D matrices and focusing functionsf þ 1d andf − 1;K¼0 are 1D vectors for a fixed focal point X F . fictitious focusing depth. The term δ contains band-limited 2D delta functions in space and time, R and R ⋆ are the same reflection-databased convolution and correlation integral operators as in the coupled Marchenko equations, and Rδ selects a given commonshot gather to be used for primary estimation. By inserting equation 16 into 17 and eliminating v − , we get This equation can be solved by iterative substitution or direct inversion to obtain v þ m and v − , which could then be used to calculate the projected upgoing Green's function by In an ideal case, when the focusing level coincides with an actual reflector in the subsurface, the first value in U − would be the primary reflection generated by that reflector with a two-way traveltime t 2 . If not, the value in U − at t 2 would be zero. By choosing t 2 as a variable and looping over all the time steps, the primaries can be retrieved by picking the values of U − at t 2 . The resulting algorithm is outlined in the following pseudocode: It is important to note here that although redatuming requires the numerical solutions only once per redatuming location, MPE requires equation 18 to be solved for every time sample of the preceding iterationfor primary prediction for a single chosen shot gather. Hereafter, we rely on this algorithm to analyze the subsampling effects in MPE. Because of the symmetry of the chosen model, the reflection responses are the same for colocated shot and receiver gathers, which allows us to freely integrate over either dimension. Figure 3 shows the focusing functions in the time domain obtained by integrating on different chosen variables, with source-or receiver-domain subsampled data. Here, we show results after a fixed number of iterations (13 in this case) of the Marchenko schemethe number of iterations is chosen arbitrarily to be more than sufficient for the scheme to achieve convergence, which in ideal sampling conditions occurs in a mere few (e.g., 3-5) iterations. Figure 3a shows the benchmark focusing function obtained without subsampling. Figure 3b displays the same focusing function obtained by subsampling in either the source or receiver and integrating on the same corresponding variable. Because of the subsampling, aliasing artifacts appear at far offsets but both focusing solutions in The sampling masks for regular subsampling in the source and receiver dimensions, respectively, where the white represents the missing data. Here, subsampling is achieved by zeroing out seven of every eight sources or receivers. The red stars and white triangles indicate the source and receiver locations, respectively. discussed by van der Neut et al. (2015b). Figure 3e-3g shows the magnified plots of Figure 3a-3c around the circles, respectively.

Joint effects of subsampling versus integration variable choice
When integration is conducted over the subsampled dimensions (Figure 3b), the summation process in every element of the output upgoing focusing functionf − 1 in equations 5 and 9 becomes increasingly inaccurate as a result of the increased subsampling interval, but they remain well sampled in space because the row space of equations 5 and 9 is well sampled. This property carries over in the subsequent convolution operations needed in the iterative scheme or in any other iterative linear solver (e.g., LSQR), leading to error propagation at later iterations. As a consequence of the negligible scaling in the source-receiver reciprocity in terms of the horizontally independent obliquity factor present in the monopole-dipole configuration (in this case because the medium is laterally homogeneous), subsampling and integrating on the source or receiver dimension are equivalent. When the subsampling and integration are performed over different dimensions, the summation process in every nonzero element of the output upgoing focusing functionf − 1 in equations 5 and 9 will be accurate because the column space of equations 5 and 9 is well sampled. However, subsampling will result in spatial gaps in the output focusing function because equations 5 and 9 are now subsampled in the row space. At later iterations, the row-space gaps will remain and inaccuracies arise in the nonzero elements of the output focusing functions because the input focusing functions gaps carry over into the summation process, as shown in Figure 3d. A diagram of all four cases in integration versus subsampling is shown in Figure 4. The corresponding focusing functions in Figure 3a-3c are transformed to the frequency-wavenumber (f − k) domain, as shown in Figure 5a, 5b, and 5d, separately. Normalized by the maximum value of Figure 5a, 5c, and 5e shows the differences between the benchmark response (Figure 5a) and those in Figure 5b and 5d, respectively. By comparing all of the panels in Figure 5, we see that the aliasing effects caused by subsampling and integrating on the same dimension are relatively weakercompared with subsampling and integrating on different dimensions in which artifacts are generally stronger. Figure 6 illustrates how the spatially dependent artifacts at the far offset come into being in Figure 3b, as a result of the column-space subsampling in equation 5. Figure 6a, 6c, and 6e shows the convolution gathers (i.e., the integrands) of Rf þ 1d at different receiver positions, where the gaps represent the missing sources. Their summation (properly scaled by the subsampling interval), corresponding to different elements in the vector of the outputf − 1 ðx r ; x F Þ, is shown by the blue traces in Figure 6b, 6d, and 6f, respectively. For reference, we overlay the corresponding traces of the focusing function without source subsampling. As can be seen from Figure 6, the focusing function traces obtained by the first convolutional operation with subsampling are very close to the reference ones at the near offset because the traveltime varies slowly with respect to x s . Toward larger offsets, the rate of traveltime variations in the integrand increases, leading to steeper slopes in the integrand and thus aliasing-induced artifacts. Figure 7a-7d shows the errors caused by subsampling and integrating on the same dimension from the first to the fourth iterations. As can be seen, most of the error comes from the second iteration and then the error decreases in magnitude with increasing iteration numbersthis is expected assuming that the iterative scheme is applied within its convergence criteria of the spectral radius (the maximum magnitude eigenvalue) of the operator, and as such our observations should hold for any media that satisfy them (Dukalski and de Vos, 2018). This is the main reason that we focus our analysis here on the artifacts arising at the leading iteration of the Marchenko redatuming scheme.

Effects of focusing depth in conjunction with subsampling
Because most data sets have fixed source and receiver apertures, the behavior of focusing functions and resulting redatumed fields changes with the focal-point depth. Hence, we complement our observations by studying the effects of focusing depth together with subsampling. Figure 8a and 8b shows the focusing functions calculated by integrating on source without and with subsampling at a shallower focusing depth ðX F ¼ 0 m; Z F ¼ 1125 mÞ. Comparing  Figures 3b and 8b, we see that the artifacts at the far offset brought by the same subsampling appear to be stronger for the shallow focusing depth than those for the deep focusing depth. This effect can be explained by comparing the traveltime slopes of the two input focusing functions as shown in Figure 8c and 8d because the shallower depth focusing function in Figure 8c has steeper slopes at far offsets, the integrand events at far offsets will be more aliased than those from a deeper focal point (Figure 6), given a fixed subsampling rate.

Source-receiver reciprocity for Marchenko methods in laterally inhomogeneous media
The Marchenko methods have a requirement of monopoles and dipoles, either on the source or receiver side, as we show in equation 3. Because of the symmetrical geometry in previous examples, the source-receiver transpose yields the same reflection record; therefore, integrating on the dipole or the monopole dimension yields the same result in the Marchenko method for laterally homogeneous media, as shown in Figure 3. However, this may no longer be true for laterally inhomogeneous media because, e.g., the angles of rays coming in and out at a fixed source-receiver pair may not be equal anymore. In practice, real dipoles are typically only available on the receiver side, for example, through particle velocity or acceleration measurements. On the source side, they are generally unavailable on offshore experiments because the source is monopole in its nature. The common practice in industry when one needs a dipole on the source side is to apply model-based corrections that would require knowledge of the local wavenumber and impedance at the source locations. Here, we show with numerical examples that, for smoothly laterally varying 2D media, the error caused by integrating over the monopole dimension is generally negligible so long as the near-surface properties in the vicinity of sources and receivers do not have pronounced variations. In the case in which the near-surface properties vary rapidly laterally, in addition   Figure 6. For the layered model in Figure 2a with the focal position ðX F ¼ 0 m; Z F ¼ 2670 mÞ corresponding to Figure 3b, i.e., subsampling and integrating over the same dimension: (a) the convolution gather from Rf þ 1d for the output near-offset receiver position X r ¼ −40 m. (b) Integration of (a), with scaling by the subsampled source interval. The blue and red lines represent the results calculated with and without subsampling, respectively. (c and d) as well as (e and f) are counterpart to (a and b), but for X r ¼ −400 m (mid-offset) and X r ¼ −960 m (far-offset), respectively. to obliquity factor, differences in local impedance at source and receiver locations must be taken into account. The 2D model used in our test is shown in Figure 9. The full reflection data have 201 colocated sources (monopole) and receivers (dipole) with spacing intervals of 10 m lying within the range of x ¼ ð3000 m; 5000 mÞ at the acquisition surface with the time sampling rate of dt ¼ 2.4ms. The focal point is located at ðX F ¼ 4000 m; Z F ¼ 1000 mÞ. The coupled Marchenko system is solved with 13 iterations. Figure 10a-10c shows the Green's function obtained from SEM modeling, Marchenko redatuming by integrating on the dipole dimension and the monopole dimension, respectively. These Green's functions correspond to waves generated from a monopole source sitting at the focal point and recorded at receivers at the acquisition surface, or vice versa, due to the source-receiver reciprocity. The focusing functions calculated with dipole and monopole integration as well as their difference are shown in Figure 11a-11c, respectively. For primary estimation, Figure 12a shows the reference reflection data with the source position X s ¼ 4000 m. Figure 12b-12d shows the result from dipole and monopole integration as well as their difference, respectively. Note that our observations are entirely consistent with wellestablished experience in SRME: amplitude effectsthose related to the incidence-angle-dependent obliquity factorsare second order with respect to kinematic errors that result from sampling deficiencies.

Effects of acquisition aperture and subsampling in redatuming: Complex model
To illustrate how acquisition sampling can impact results in a more realistic scenario, we study the effect of a limited aperture, together with source integration and subsampling in the redatumed Green's function with a salt model, which has a constant density and varying wavespeed ( Figure 13). This model has also been used by Vargas and Vasconcelos (2020) in the context of scattering-based focusing. We choose this model such that the effects of density-related impedance changes do not play a role, which allows us to focus on the kinematics and reconstruction of the eventsin any case. Because sampling-related inaccuracies are controlled by the event slope (kinematics) and bandwidth, we expect that tests on a variable-density version of this model would yield the same conclusions as those we present later. The full reflection data have 201 colocated sources and receivers with spacing intervals of 40 m in the subsalt model. The acquisition spread lies within the range of x ¼ ð4000 m; 12;000 mÞ at the acquisition surface with the time sampling rate of dt ¼ 4ms. The focal point is located ðX F ¼ 8130 m; Z F ¼ 4400 mÞ, as shown in Figure 13. The same subsampling scheme is used as shown in Figure 2b. Figure 14a shows the Green's functions obtained with the full aperture. As observed in Figure 14b, subsampling by a ratio of eight brings strong artifacts for the complex salt a) b) Figure 9. (a and b) Velocity and density profiles of the 2D model used in the source-receiver transpose study for the Marchenko focusing and primary estimation, respectively. The magenta dot denotes the location of focal point (x ¼ 4000 m, z ¼ 1000 m). The red stars and white triangles indicate the source and receiver locations, respectively.  Figure 10. For the 2D model shown in Figure 9: (a-c) the Green's function from forward modeling, from Marchenko redatuming by integrating on the dipole dimension, and from Marchenko redatuming by integrating on the monopole dimension, respectively.  Figure 14a and 14c illustrates that the acquisition with smaller aperture also increases apertureinduced artifacts in the Green's function. Figure 14d shows that a smaller aperture and subsampling compound into more complicated artifacts. Here, all of the Green's functions are calculated by solving the coupled Marchenko equations by inversion with the LSQR solver (Paige and Saunders, 1982;Ravasi and Vasconcelos, 2020).

Effects of subsampling in MPE
Completing our analysis of subsampling effects on redatuming, here we study the receiver subsampling effects on MPE in a fourlayer model with colocated 201 sources and receivers, as shown in The redatumed Green's function calculated by subsampling and integrating over the source and integrating on the source. The same subsampling scheme is used as shown in Figure 2b. (c and d) Difference plots between the true and estimated responses use two-third of the original data aperture: (c) is from fully sampled data, whereas (d) is from subsampled data. In (c and d), the missing traces indicate the part of the data that was removed for these tests.
Initialize Rδ with a chosen common-shot gather; for t 2 ¼ 0 to t max do compute the windowing operator W t 2 −ϵ ϵ ; solve equations 18 and 6 for v þ m and v − ; use v þ m and v − to compute U − (equation 19); store the value of U − at t 2 ; Figure 15. Figure 16a shows the reference reflection data with the source location X s ¼ 0 m, where four strong primaries and the ensuing relatively weak multiples can be observed. Figure 16b shows the primary prediction resultusing the algorithm presented in Algorithm 1without subsampling, where only four primaries are visible and the multiples are suppressed significantly. Figure 16c-16e shows the results of primary estimation with three different subsampling schemes (evenly zeroing out 66%, 75%, and 80%), calculated with 20 iterative substitutions, respectively. Here, as previously, the choice of 20 iterations is arbitrary, chosen to be sufficiently higher than our observed number of iterations to achieve convergence on fully sampled data. Figure 16f shows the result using the same subsampling scheme as shown in Figure 16e but calculated with LSQR instead of explicit Neumann iterations. As expected, the strength of the artifacts increases with increasingly subsampled data. The errors in the first three iterations for the 75% subsampling scheme are given in Figure 17a-17c, which show that the major error comes from the first iteration and can decrease in subsequent iterations. Figure 16e shows that the iterative substitution method fails to converge for this example when the subsampling degree increases to 80%. Based on the theory of Neumann series expansion, in Figure 18, we examine the operator norm of A ¼ W t 2 −ϵ ϵ R ⋆ W t 2 −ϵ ϵ R in equation 18 for different subsampling schemes as a function of time. In general, the operator norm first increases and then plateaus with time t 2 . This occurs as more finite-energy arrivals are included by the windowing operator W t 2 −ϵ ϵ as the time increases, but the contribution of an increasing number of arrivals is offset by the decreasing energy of high-order multiples that arrive late (Zhang and Staring, 2018). As shown in Figure 18, the operator norm with full data for primary estimation with 20 iterations lies below that of the Marchenko redatuming for the chosen focal depth (0.8094), which satisfy the convergence condition jjAjj < 1. It is important to note that jjAjj < 1 guarantees convergence, whereas for jjAjj > 1, it may diverge not guaranteed, but true for our examplefor details on Neumann series convergence, we refer readers to Dukalski and de Vos (2018). Subsampling evenly by 66% and 75% brings the operator norm to approximately 1.2 and 2.1, respectively, where the iterations diverge slowly but lie within an acceptable range (see Figure 16c and 16d). Subsampling evenly by 80% further pushes the operator norm to approximately 3.1, resulting in a much faster divergence as shown in Figure 16e. Here, we point out again that this operator norm analysis does not carry the same implications when the problem is solved without relying on explicitly evaluating the Neumann-series-based solution.

DISCUSSION
Our study with the four-layer model shows that for Marchenko redatuming the subsampling effects jointly depend on (1) the choice of integration variable and (2)  receiver) is subsampled. For the same reason, when the sampling of reflection data is symmetric with respect to the source and receiver dimensions, integrating on either dimension yields the same result, regardless of subsampling. When the data are regularly subsampled on one dimension, integrating on the other dimension yields inaccurate focusing results, with focusing function gaps at surface locations corresponding to the missing data. However, in this same sampling scenario, integrating on the subsampled dimension yields focusing functions with artifacts but no spatial gaps, with their strength affected by the chosen focusing depth. The analysis of the integrand gather shows that the mechanism behind lies at the integration within the multidimensional convolution between the reflection response and the focusing function formulated in the form of a matrix-vector multiplication in the frequency domain. Naturally, we expect that subsampling in sources and receivers will compound and cause stronger artifacts together with gaps, compared with subsampling in a single dimension. In addition, limiting acquisition aperture will also bring artifacts (Sripanich and Vasconcelos, 2019), which can be compounded by the subsampling effects, as we show with our salt model example.
Furthermore, we show with a model of mild lateral heterogeneity that conducting integration over the monopole dimension instead of the dipole brings negligible errors to both Marchenko systems investigated here, which may facilitate our choice of integration variable in combination with the subsampling scheme. This also reveals that the Marchenko methods are not highly sensitive to relatively small errors contained within the reflection data. Therefore, it implies that in practice interpolating the gaps of the input reflection data should be preferred over postprocessing the output focusing functions or Green's functions, as shown by Jia et al. (2019) and Staring and Wapenaar (2020) for 3D redatuming.
One important take away from our results with different combinations of dipole and monopole field quantities is that, although theory strictly requires the use of dipole sourceswhich are not practically availableit is entirely feasible to rely on dipole data instead, e.g., those acquired with multicomponent sensors. Although replacing dipole sources by receivers may violate spatial reciprocity in heterogeneous media, we believe that, in most offshore conditions, the associated errors are likely negligible. In the absence of any dipole data, e.g., the case of conventional streamer data, we recommend that a monopole-to-dipole filter be applied in the best-sampled domain, likely to be the receiver domain in streamer-based data.
In the case of ocean-bottom nodes, which are multicomponent but often highly sparse, it may be better to apply dipole corrections on the more densely sampled source domain and perform convolutions over sources.
In the context of focusing and redatuming, it is important to point out that the behavior of subsampling-related errors in Marchenko operations is analogous to those occurring in SRME. In SRME literature and practice (Dragoset et al., 2006(Dragoset et al., , 2010Verschuur, 2013), acquisition-related subsampling has been studied by many authors in many contexts -2D, 3D, synthetic, and field data. Throughout the literature in multiple elimination, subsampling effects have been extensively investigated, together with several approaches suggested to condition data for multiple suppressionincluding interpolation/regularization techniques. For example, in closed-loop SRME (Lopez and Verschuur, 2015), one may use a focal transformation and a sparse norm regularization to remove the undersampling noise. Although our study begins to highlight the connections between these methods, we expect that more about Marchenko can be learned from what has been studied in the context of SRMEe.g., by taking inspiration from data conditioning approaches used in that context. However, appropriately translating SRME lessons into Marchenko practice will require further testing and research, particularly in the context of field data examples.
Although here we base our discussions on 2D examples, the arguments that we build regarding the integration over monopole fields are general, and as such in principle apply to 3D media, i.e., having source and receiver integrals each act on a plane of source/receiver coordinates. We point out that the work of Ravasi and Vasconcelos (2020) contains examples of subsampling 3D Marchenko operatorsindeed in agreement with this study. Their example does further highlight that although aliasing effects in three Operator norm Marchenko redatuming with 100% data Primary estimation with 100% data Primary estimation with zeroing out 66% data Primary estimation with zeroing out 75% data Primary estimation with zeroing out 80% data Figure 18. (a) Comparison of the operator norm of iterative substitution in primary estimation with different subsampling schemes (100% data, regularly zeroing out 66%, 75%, and 80%, respectively), as well as in the Marchenko redatuming with full data for the same model. Here, we note that the Marchenko redatuming operator norm corresponds to a fixed focal point at depththus not being one-to-one comparable with primary estimation on the acquisition surfacebut we choose to show it here merely to indicate that the operator norm for redatuming is not time-dependent.
dimensions follows the same patterns as we present here, they can be lower in magnitude when either source or receiver integration is subsampled in only one of the two spatial coordinates. This observation is likely to depend on the degree of subsampling and on the degree of complexitydue to the subsurface heterogeneitypresent in the 3D wavefields that make up the reflection operator.
A study of subsampling effects on the primary estimation based on the projected Marchenko system shows a similar pattern of artifacts that appear in the original coupled Marchenko system, which can be explained by the common multidimensional convolution and correlation between reflection data and focusing functions. As expected, the artifacts increase with the degree of subsampling and may cause divergence issue with iterative substitution when subsampling to a certain extent, which is illustrated by the analysis of the time-dependent operator norm lying in the core of the iteration. This problem can be mitigated using more robust iterative linear solvers, such as LSQRbut this would likely come at a greater computational cost because such solvers tend to convergence more slowly than the Neumann series, thus requiring a larger number of iterations (i.e., of operator evaluations; Dukalski and de Vos, 2018). In our case, LSQR yields consistent results, as opposed to Neumann-series-based iterative substitution, when presented with data with high degrees of subsampling. When the data are well sampled, the Neumann-series iterative scheme and the LSQR solver perform equally well. As in using the coupled Marchenko system for redatuming, the errors caused by subsampling decrease rapidly with the number of iterationsprovided that the Neumann series remains stable.

CONCLUSION
By relying on representative numerical cases to mimic acquisition scenarios over varying subsurface complexity, we provide a study of the acquisition-related subsampling effects on the 2D Marchenko focusing, redatuming (original coupled Marchenko system), and primary estimation (projected Marchenko system) to shed light on the underlying physical and numerical mechanisms. For Marchenko focusing and redatuming, we show that the artifacts brought by subsampling are rooted in the integration within the multidimensional convolution between the reflection response and the focusing function, which also applies to primary estimation because of how similar the underlying operations are. These subsampling effects depend jointly on the integration variable choice and the subsampled dimension (the dipole source or receiver). Artifacts appear without spatial gaps when these two are identical, and spatial gaps are present with artifacts for different dimensions. Our study also shows that integrating on the monopole dimension also gives very close approximation to the results of both Marchenko systems for smoothly laterally varying media. In this regard, future studies are expected regarding the sampling requirement for convergence, the interpolation accuracy, etc. The artifacts in redatumed Greens' functions depend on the reflection data sampling and focusing depths. In imaging, these artifacts will very likely cause small errors, but for full-waveform inversion, time lapse, and reservoir characterization, these could lead to larger errors.
Our study focuses only on regularly subsampling in one dimension of the input reflection data, i.e., source or receiver dimension. Meanwhile, the conclusions drawn from our study are limited to the 2D case. For three dimensions, the irregular subsampling effects are likely more complicated, due to the 2D integration over the source and receiver dimension, particularly in the case of highly complex media. Our results can inform acquisition design, but more importantly they tell us how to handle single and multicomponent data for Marchenko applications, while providing insight that can help in ensuring the quality of data conditioning approaches for Marchenko applications.