Modeling and data-driven isolation of two-way wavefield constituents

Synthesizing individual wavefield constituents (such as primaries, first-order scattering, freesurface or internal multiples) is important in the development of seismic data processing algorithms, for instance for seismic multiple removal and imaging. A range of methods that allow for the computation of such wavefield constituents exist, but they are generally restricted to relatively simple, horizontally-layered media. For wave simulations on more complex models a straightforward and performant alternative are finite-difference methods. They are, however, generally not perceived as being capable of delivering isolated wavefield constituents. Based on recent advances, we show how this can be achieved for (non-horizontally) piecewise constant layered media. For example, we are able to accurately retrieve the isolated direct arrival of the transmission response (including tunneled waves), primary reflection data (without internal multiples), and all events related to a single (or multiple) interface(s) in a medium. The proposed methods require detailed knowledge of discretized medium parameters. If a medium is known only implicitly via recordings of reflection data, interface-related events can still be isolated through a combination of sub-domain related wavefields. We show how Marchenko redatuming can be used to derive these. Isolation of wavefield constituents Elison et al. 2 Isolation of wavefield constituents INTRODUCTION Full-waveform methods, such as finite-difference (FD) and finite-element (FE) methods, can be used for complex 2D and 3D wave propagation modeling. In contrast to other methods, where wavefield constituents such as upand downward parts can be naturally isolated, (e.g., like the reflectivity method by Kennett, 1983), FD or FE methods traditionally have been thought of as being capable of computing the full-wavefield only. However, recent work has shown that FD and FE methods can also be used to isolate certain parts of the two-way wavefield with numerical precision (Robertsson et al., 2015). Examples include the work by Amundsen and Robertsson (2014), which describes a method for up-down separation and direct wave removal, and the work by Vasmel et al. (2016) for surface-related multiple elimination. Simulation of individual wavefield constituents is of high interest. A fundamental processing step for seismic data is the removal of multiply scattered waves. This is because most seismic imaging algorithms rely on the Born approximation, such that the presence of multiply-scattered waves causes artifacts and ‘ghost’ events in the final image. The availability of a reference solution is therefore crucial to, for example, quality-control of internal demultiple algorithms. Moreover, the identification of target-related events is important for amplitude-versus-offset analysis in reservoir characterization and time-lapse seismic surveys. Event identification can also be used to recognize source-generated noise as well as prismatic or any other unknown events. In this work we introduce methods to isolate different wavefield constituents. First, we show how to isolate first-order wavefield constituents in a two-way wave propagation simulation (i.e., primary reflections and the direct transmitted arrival) by means of pseudo one-way wave propagation. Second, we show how to obtain only the reflections related to a single (or a stack of) interfaces, including both primaries and multiples. This is achieved by cloaking of these interfaces, followed by comparison with data obtained for the full medium. A data-driven as well as a model-driven approach for such cloaking is presented. We briefly outline the methods introduced in this work: Pseudo-one-way wave propagation. We achieve one-way wave propagation by isolating individual interfaces. Assuming sub-horizontally layered media, we split a wave simulation into a series of consecutive, smaller sub-simulations where every sub-simulation contains only the medium contrast of a single interface. Above and below every interface we use sound-transparent auxiliary surfaces to record and inject wavefields, using the method of multiple point sources (MPS, Amundsen and Robertsson, 2014). This approach assumes piecewise constant layered media. In a first step we downward propagate the transmission through the medium, computing an isolated direct arrival of the transmitted wavefield. This is particularly interesting for scenarios where the interfaces are separated by distances in the order of the signal’s wavelength or less. Since the direct arrival overlaps with its scattering coda in such cases, separations based on time-windowing are ruled out. We show that for such scenarios our method includes tunneled waves. In a second step we use the recorded primary reflections from the first step and return them to the surface. This yields isolated primary reflections at the acquisition surface. Comparison with conventional reflection data yields the multiply-scattered wavefield. All examples are given for FD solvers, but in principle the same approach could be implemented in FE solvers. Cloaking of interfaces. Scattering due to a single (or a stack of) interface(s) can be removed from a wavefield if we assume that a medium can be divided into submedia above and below the interface(s) of interest, for simplicity denoted here as target interfaces (but without making assumptions on the location of these interfaces). The key is to use only the direct arrival across the target to couple the domains above and below. In case of at least two closely spaced target interfaces, the Isolation of wavefield constituents Elison et al. 3 Isolation of wavefield constituents pseudo one-way wave propagation introduced in this work yields this direct arrival. The resulting data are free of any (primary and multiple) reflections related to the target interfaces. Subsequent comparison with data obtained in the full medium yields only these reflections. We introduce two methods that achieve such cloaking. First, we present convolutional expressions which allow us to express surface reflection data in terms of wavefields that are related to different sub-domains. In combination with Marchenko inversion (Wapenaar et al., 2014) these expressions can be used to convert measured surface reflection data of the full medium into data which do not contain (primary or multiple) reflections related to the cloaked interface(s). In addition to the measured data, this requires a background velocity model and information about transmission amplitudes. Second, we show how the same cloaking can be achieved in a FD simulation. Modified Immersive Boundary Conditions (IBCs, van Manen et al., 2007), placed in the vicinity of the target interface(s), couple FD simulations above and below and therefore cloak the target. The result is again a wavefield that contains the correct transmission across the target interface(s), but no reflections related to it. This model-driven approach requires detailed knowledge of the medium (including location of interfaces, seismic velocities and densities), but since it does not assume horizontal layering it significantly differs from many other existing methods for multiple removal (Weglein et al., 1997; Verschuur and Berkhout, 2005; Stork et al., 2006). THEORY The theory is divided according to the different wavefield constituents that are isolated. First, we outline how to achieve pseudo one-way wave propagation in layered media, yielding only primary reflections and an isolated direct transmission event. Thereafter, we show how one or several target interfaces can be cloaked, either by a series of convolutions of wavefields or in a FD simulation. To improve readability we define a short-hand notation to express multidimensional convolutions of Green’s functions according to


INTRODUCTION
Full-waveform methods, such as finite-difference (FD) and finite-element (FE) methods, can be used for complex 2D and 3D wave propagation modeling.In contrast to other methods, where wavefield constituents such as up-and downward parts can be naturally isolated, (e.g., like the reflectivity method by Kennett, 1983), FD or FE methods traditionally have been thought of as being capable of computing the full-wavefield only.However, recent work has shown that FD and FE methods can also be used to isolate certain parts of the two-way wavefield with numerical precision (Robertsson et al., 2015).Examples include the work by Amundsen and Robertsson (2014), which describes a method for up-down separation and direct wave removal, and the work by Vasmel et al. (2016) for surface-related multiple elimination.
Simulation of individual wavefield constituents is of high interest.A fundamental processing step for seismic data is the removal of multiply scattered waves.This is because most seismic imaging algorithms rely on the Born approximation, such that the presence of multiply-scattered waves causes artifacts and 'ghost' events in the final image.The availability of a reference solution is therefore crucial to, for example, quality-control of internal demultiple algorithms.Moreover, the identification of target-related events is important for amplitude-versus-offset analysis in reservoir characterization and time-lapse seismic surveys.Event identification can also be used to recognize source-generated noise as well as prismatic or any other unknown events.
In this work we introduce methods to isolate different wavefield constituents.First, we show how to isolate first-order wavefield constituents in a two-way wave propagation simulation (i.e., primary reflections and the direct transmitted arrival) by means of pseudo one-way wave propagation.Second, we show how to obtain only the reflections related to a single (or a stack of) interfaces, including both primaries and multiples.This is achieved by cloaking of these interfaces, followed by comparison with data obtained for the full medium.A data-driven as well as a model-driven approach for such cloaking is presented.We briefly outline the methods introduced in this work: Pseudo-one-way wave propagation.We achieve one-way wave propagation by isolating individual interfaces.Assuming sub-horizontally layered media, we split a wave simulation into a series of consecutive, smaller sub-simulations where every sub-simulation contains only the medium contrast of a single interface.Above and below every interface we use sound-transparent auxiliary surfaces to record and inject wavefields, using the method of multiple point sources (MPS, Amundsen and Robertsson, 2014).This approach assumes piecewise constant layered media.In a first step we downward propagate the transmission through the medium, computing an isolated direct arrival of the transmitted wavefield.This is particularly interesting for scenarios where the interfaces are separated by distances in the order of the signal's wavelength or less.Since the direct arrival overlaps with its scattering coda in such cases, separations based on time-windowing are ruled out.We show that for such scenarios our method includes tunneled waves.In a second step we use the recorded primary reflections from the first step and return them to the surface.This yields isolated primary reflections at the acquisition surface.Comparison with conventional reflection data yields the multiply-scattered wavefield.All examples are given for FD solvers, but in principle the same approach could be implemented in FE solvers.
Cloaking of interfaces.Scattering due to a single (or a stack of) interface(s) can be removed from a wavefield if we assume that a medium can be divided into submedia above and below the interface(s) of interest, for simplicity denoted here as target interfaces (but without making assumptions on the location of these interfaces).The key is to use only the direct arrival across the target to couple the domains above and below.In case of at least two closely spaced target interfaces, the Isolation of wavefield constituents Isolation of wavefield constituents pseudo one-way wave propagation introduced in this work yields this direct arrival.The resulting data are free of any (primary and multiple) reflections related to the target interfaces.Subsequent comparison with data obtained in the full medium yields only these reflections.We introduce two methods that achieve such cloaking.First, we present convolutional expressions which allow us to express surface reflection data in terms of wavefields that are related to different sub-domains.In combination with Marchenko inversion (Wapenaar et al., 2014) these expressions can be used to convert measured surface reflection data of the full medium into data which do not contain (primary or multiple) reflections related to the cloaked interface(s).In addition to the measured data, this requires a background velocity model and information about transmission amplitudes.Second, we show how the same cloaking can be achieved in a FD simulation.Modified Immersive Boundary Conditions (IBCs, van Manen et al., 2007), placed in the vicinity of the target interface(s), couple FD simulations above and below and therefore cloak the target.The result is again a wavefield that contains the correct transmission across the target interface(s), but no reflections related to it.This model-driven approach requires detailed knowledge of the medium (including location of interfaces, seismic velocities and densities), but since it does not assume horizontal layering it significantly differs from many other existing methods for multiple removal (Weglein et al., 1997;Verschuur and Berkhout, 2005;Stork et al., 2006).

THEORY
The theory is divided according to the different wavefield constituents that are isolated.First, we outline how to achieve pseudo one-way wave propagation in layered media, yielding only primary reflections and an isolated direct transmission event.Thereafter, we show how one or several target interfaces can be cloaked, either by a series of convolutions of wavefields or in a FD simulation.
To improve readability we define a short-hand notation to express multidimensional convolutions of Green's functions according to where ω denotes the radial frequency and x i denotes spatial locations on a surface ∂D i .All wavefields in this work denote pressure-normalized monopole recordings due to dipole sources.
We introduce notations for Green's functions with different acquisition geometries, and examples can be found in Figure 3. Green's functions that are recorded with sources and receivers on a scattering-free surface ∂D i with an adjacent homogeneous half-space above or below are denoted as reflection responses R ∪ i,i and R ∩ i,i , respectively.If the medium is truncated with another homogeneous halfspace at ∂D j we denote this with an additional superscript in brackets R ∪(j) i,i [see an example for R ∪(1) 0,0 in Figure 3(k)].A reflection response recorded at ∂D 0 that contains only primary reflections is denoted it with P = R ∪ 0,0 (primaries), and if it contains only a primary reflection due to a single interface this is specified with a superscript, e.g., P (3) denotes a primary reflection due to the third interface.If sources and receivers are located on different surfaces ∂D i and ∂D j , respectively, with adjacent homogeneous halfspaces on both of these surfaces, we denote corresponding Green's functions as transmission responses T j,i , and if such a transmission contains only the direct (one-way) transmission it is denoted with D j,i .Ultimately, if at least one of the (source or receiver) surfaces is located inside the medium (i.e., without an adjacent homogeneous halfspace at that surface) we denote Green's functions with G j,i .

Isolation of wavefield constituents
Figure 1: Full-wavefield simulation and stepwise computation of the direct arrival through layered media.(a) A full-wavefield simulation of a source (red star) yields primary reflections at the top, and a direct transmitted wave at the bottom (thick arrows).These events may overlap with internal scattering (thin arrows), which make a separation of these wavefield constituents difficult; (b)-(d) stepwise computation of the direct arrival through the same layered medium: in three subsequent simulations the transmission is computed through each interface individually.

Pseudo one-way wave propagation
Conventional simulations of acoustic two-way wave propagation contain the full-wavefield.This is illustrated for a layered-medium example in Figure 1(a), which is assumed to be reflection free at its boundaries.Simulation of a source excited at the upper boundary ∂D 0 yields primary reflections P at the top and a directly transmitted wave D 3,0 at the bottom ∂D 3 , both illustrated with thick arrows and hereafter referred to as first-order wavefield constituents.In addition, a (theoretically infinitely long) series of internal scattering will be observed, illustrated with thin arrows for firstorder internal multiples.These latter events typically overlap with the primary reflections at the surface.They also overlap with the directly transmitted event at the bottom if at least one interface separation distance is on the order of the wavelength of the signal.Hence, separation of multiple scattering from first-order scattering is difficult in conventional two-way wave simulations.
To isolate first-order wavefield constituents we use repeated recordings and injections of twoway wavefields along arbitrarily curved (i.e., non-horizontal) auxiliary interfaces, assuming piecewiseconstant, layered media with continuous, non-crossing layer interfaces.Injection and recording of wavefields is achieved with the method of multiple point sources [MPS] (Amundsen and Robertsson, 2014) on surfaces that are aligned with the computational (FD) grid.This is illustrated in Figure 2 for an arbitrary gray surface that is discretized by the dashed red surface to align with pressure grid nodes.Other methods exist where this is different: Thomsen et al. (2018) have shown that an MPS surface can be expressed as the combination of two FD injection surfaces (Robertsson and Chapman, 2000), which are offset with respect to the grid nodes, and Mittet (2002) presented results for recording and injection surfaces which are not aligned to grid nodes.
Recordings of the pressure and the normal component of the particle velocity along a surface from one simulation can be injected in a subsequent simulation to reproduce the wavefield at that surface.The red background in Figure 2 highlights the area that contains associated MPS grid points.Along that surface the pressure p and the normal components of the particle velocity v are recorded, denoted with v x and v z for its horizontal and vertical components, respectively.The same nodes are used in a subsequent simulation for wavefield injection: the pressure recordings are distributed over grid nodes of the normal component of the particle velocity, and the recordings of the normal components of the particle velocity are averaged and injected in associated pressure Isolation of wavefield constituents Isolation of wavefield constituents nodes (Vasmel and Robertsson, 2016), summarized in the bottom of Figure 2. Grid nodes that are used by the MPS surface are illustrated in black in Figure 2, and if they are associated with a corner point they have a red center.In-line particle velocity nodes are not used for the MPS method, and hence they are illustrated in gray even if they lie within the red area.Red arrows in Figure 2 point in the normal direction of the MPS surface, and in the bottom we give examples of pressure and particle velocity injection on a straight section (bottom right) and on a corner point (bottom left).If normal vectors point in positive x or z direction, the normal components n x and n z in the corresponding equations take the values +1, and −1 otherwise.
Corner points of MPS surfaces are always located on pressure nodes, which have two normal components.Therefore their neighboring particle velocity nodes in both horizontal and vertical direction must be considered, which is summarized by the equations in the bottom left of Figure 2.For the corner point highlighted in that example the normal vector in x-direction points in negative direction, and hence n x takes the value −1.In the figure this is expressed by the minus signs next to the arrows that relate grid nodes in the horizontal direction.
Averaging of particle velocity measurements and distribution of pressure measurements are defined by the same FD coefficients that are used in the FD solver.The equations in Figure 2 are given for a second-order FD stencil, because this is what we we use in this work.General expressions were given by Vasmel and Robertsson (2016), who have shown that the approach for the second order stencil is accurate to numerical precision.The same authors have shown that for other FD stencils corner points introduce weak artifacts.However, these artifacts are orders of magnitude smaller than the actual signal, and for stencils of increasing order these artifacts become less.
MPS injection reproduces the previously recorded wavefield, given that the medium is identical along the MPS surface in subsequent recording and injection steps (i.e., the medium must be identical for all grid nodes within the red background in Figure 2).Everywhere else the medium may be different in two subsequent simulations.This is a key feature exploited by many applications of such wavefield injection, i.e., removing the complexity of the medium in the injection step enables the retrieval of individual wavefield constituents (e.g., compare with Vasmel et al., 2016).Here, it allows us to simulate reflections and transmission for individual interfaces in an isolated fashion.
In the following we show how first-order reflected and transmitted wavefield constituents can be expressed analytically, and how their isolation is achieved using a sequence of MPS injection and recording steps.

Direct transmission computation
We use auxiliary surfaces located between individual interfaces to express a one-way transmission in terms of a series of two-way transmissions.These auxiliary surfaces are scattering-free since we assume piecewise constant, layered media.With ∂D 0 defined as uppermost (acquisition) surface and auxiliary surface i located below the i th interface, we can express the direct arrival of the downward transmission response through a stack of n interfaces as where refers to a series of multidimensional convolutions according to Eq. ( 1), which we summarize with a one-way propagator W .For the three-interface example in Figure 1 two auxiliary Isolation of wavefield constituents Isolation of wavefield constituents Figure 2: Implementation of MPS along an arbitrarily curved surface in a second-order staggeredgrid FD scheme.The gray surface is discretized by the dashed red surface, and associated MPS nodes lie within the red band.Recordings of pressure (dots) and particle velocity (triangles) from the first simulation are injected in the associated other component's nodes, whos locations are defined by the normal direction (red arrows) of the discretized surface.The direction of the normal vectors define the normal components n x and n z which take the values ±1.Corner points (black and red) have two normal directions and need to take into account nodes in horizontal and vertical direction.
Equations for the injection in the second-order scheme are given in the bottom.
To compute the direct arrival of the transmission in a two-way wave solver, a series of MPS injections and recordings is used.We use the example in Figure 1 to explain how this is achieved.In the first injection [panel (b)] a source is excited on the (reflection-free) top of the model at ∂D 0 , similar to the simulation illustrated in panel (a).However, for this simulation only the medium contrast of the first interface is used, which is bounded by two homogeneous half-spaces (to ensure a reflection-free medium everywhere else).Consequently, only the downward traveling directly transmitted wave D 1,0 = T 1,0 is observed below the interface, where it is recorded along the scatteringfree auxiliary MPS surface ∂D 1 .Note that the source also causes a primary reflection P (1) above the interface.
A subsequent simulation uses a model that contains only the medium contrast of the second interface, see panel (c).The previously recorded (purely downgoing) wavefield is used as source wavefield and injected along the surface it was recorded on: the auxiliary recording surface ∂D 1 from the previous simulation is now the injection surface.Again, the transmission through the current interface is recorded on another auxiliary surface ∂D 2 just below (and now D 2,0 = T 2,1 T 1,0 is recorded).The upward (primary) reflection of this second interface does not cause an internal multiple since the overlying first interface is not present in this simulation [in contrast to the simulation of the full transmission response in Figure 1(a)].This procedure is repeated for every subsequent interface.For the example in Figure 1 this requires one additional simulation [panel (d)].The wavefield recorded below the last interface contains only the direct arrival D 3,0 across all interfaces, but does not contain any internal multiples.

Primary reflection computation
One-way propagators such as Eq. ( 2) can be used to compute primary reflections.Upward and downward one-way transmissions, combined with a reflection computed on single interfaces, yields an expression for the primary reflection due to interface n + 1 according to Implementation of Eq. ( 3) with MPS surfaces works similar as the computation of the direct transmission.As illustrated in Figure 1 (b)-(d) we use recording surfaces below the interfaces to compute the direct transmission across a medium.If the wavefield is additionally recorded above the interfaces we can also obtain the primary reflections.Subsequent upward propagation of these primary reflections in the same layer-by-layer fashion (i.e., using W 0,n ) allows us to retrieve primary reflection data at the acquisition surface.
To obtain all primary reflections we first need to compute the direct transmission across the entire medium downwards, storing the primary reflections obtained from all simulations individually.Subsequently, a single series of simulations in the upward direction suffices: we start by propagating the primary reflection from the bottom-most interface upwards through the second-last interface.In subsequent simulations the injected wavefield is always a combination of the primaries that have been propagated back upwards, and the primary reflection recorded in the initial transmission series (at the current surface).Every interface must therefore be simulated twice to achieve this.The first series of simulations yields the one-way transmitted wavefield, whereas the second series yields the Isolation of wavefield constituents Isolation of wavefield constituents transmission of the (total) first-order reflected fields.Alternatively every primary event can be propagated upwards individually, which yields isolated primary reflections at the acquisition surface, but requires considerably more computational effort.
Regarding the numerical implementation it is noteworthy that MPS injection and recording surfaces can not coincide, and hence we can not record the reflections on the same surface where the injection takes place.Therefore, in the downward transmission series, recording of the primaries must take place at surfaces that are offset with regard to the current injection surface, such that related grid nodes only correspond to one MPS surface (compare with the light red band in Figure 2 which may not overlap with a neighboring one).

Cloaking of interfaces
Next, we present two methods that allow us to derive a surface reflection response R∪ 0,0 that does not contain primary or multiple reflections related to a target interface (or a stack of target interfaces), denoted with a tilde.Subsequent comparison with 'conventional' reflection data R ∪ 0,0 yields all reflected events related to the target interfaces.The methods are based on dissecting a layered medium into three different sub-domains located on top of each other.The central sub-medium is defined as the target area, which contains one or more interfaces whose related (primary and multiple) reflections shall be isolated.

Convolutional expressions
The derivation of convolutional expressions is divided in three steps, illustrated in the three columns of Figure 3.The example that is used contains a single target interface, illustrated as solid red line for wavefields that contain the interface, and as dashed red line for wavefields where reflections related to it are absent.In this example the target is bounded by horizontal surfaces ∂D 1 and ∂D 2 at depth levels z 1 and z 2 , yet these surfaces need not be horizontal, but they should not cross any interfaces.The equations given here have similarities to the ones derived by Wapenaar and Staring (2018) and Elison et al. (2018).
When read from right to left, this equation is equivalent to the left column of Figure 3 read from top to bottom.A reflection response R∪ 1,1 is recorded with sources and receivers above the target area.It does not contain any reflected events related to it but only events related to the medium below.Note that for a single target interface D 2,1 = T 2,1 , however an expression such as Eq. ( 2) is needed for a target area containing more than a single interface.
In the subsequent step, illustrated in Figure 3 Eq.( 4) Eq.( 5)   4), (e) to (g) illustrate Eq. ( 5) and (h) to (k) illustrate Eq. ( 6).The solid red line in panels (a) and (c) denotes the interface of interest that should be removed from the data.It is drawn as dashed line for wavefields where its related events are absent.
Isolation of wavefield constituents Isolation of wavefield constituents to the target, according to (5) Powers of i indicate i-times repeated multidimensional convolutions of the respective terms, which corresponds to i-th order multiple interactions between the upper and lower sub-medium.Note that the first term on the right hand side always has downgoing source and upgoing receiver components, and vice verse for the second term.The third term consists of two similiar terms, one with upgoing and one with downgoing source and receiver components.
A third step ultimately enables synthesizing surface reflection data.Combining the previously obtained Green's function G1,1 with the overburden's transmission responses T 1,0 and T 0,1 [Figure 3(h) and (j)], it can be added to the reflection response of the overburden R ∪,(1) As a result a surface reflection response R∪ 0,0 where all reflections related to a single interface are absent is obtained.This is expressed by Here the additional superscripts −,+ denote only upgoing receiver components and downgoing source components of G1,1 , respectively, which can be obtained by evaluating only the first term in Eq. ( 5), i.e., G−,+ If data are available from a seismic reflection experiment, Marchenko inversion (Wapenaar et al., 2014) offers a way to derive the sub-domain-related wavefields used in Eqs.(4-6).Using the reflection response of the real medium R ∪ 0,0 and a smooth background velocity model as input, the Marchenko method enables the retrieval of Green's functions between virtual sources in the subsurface and receivers at the acquisition surface.Multidimensional deconvolution of the Marchenko solutions yields reflection responses R ∪,(1) 0,0 , R ∩ 1,1 and R ∪ 2,2 related to individual sub-domains.The overburden's transmission response T 0,1 and T 1,0 is obtained by inverting Marchenko focusing functions.However, additional information about the transmission amplitude is required to correctly scale these transmission responses as well as the direct arrivals across the interface of interest, D 1,2 and D 2,1 .Such information could be obtained for example from a well log.

Cloaking with immersive boundaries
In contrast to the data-driven method introduced in the previous section, cloaking of interfaces can also be achieved in a wave propagation simulation.Wave simulations carried out in an upper and a lower sub-domain are coupled via the direct transmission response D of the intermittent target medium.This results in an accurate reconstruction of interactions between the upper and lower submedia (as observed in the two-way modeling for the full medium), but does not include any (first or higher order) reflections related to the target medium.The coupling between the wave simulations is based on IBCs as introduced by van Manen et al. (2007); however the IBC configuration used here resembles the one introduced by Elison et al. (2018).The direct transmission D across the target medium can be obtained with the MPS method presented earlier in this work.We refer to this approach as model-driven cloaking, because the wave simulation requires detailed knowledge of the medium.In the following we first give a brief description of conventional IBCs, followed by Isolation of wavefield constituents Isolation of wavefield constituents IBC extrapolation surface IBC injection boundary introducing the alternative implementation that we use, which enables cloaking of a target area in layered media.
IBCs are time-dependent boundary conditions which enable the immersion of a wavefield propagating on a truncated sub-domain into an arbitrary larger domain (van Manen et al., 2007;Vasmel et al., 2013).They are based on an implementation of non-reflecting boundaries by Ting and Miksis (1986), and although the concept of IBCs has been extended to physical experiments (Becker et al., 2018;Börsing et al., 2019), we restrict our explanations to numerical wave simulations.IBCs require an inner extrapolation surface ∂D e and an outer injection boundary ∂D i , where the latter constitutes the actual boundary of the sub-domain D. An arbitrary IBC sub-domain of a layered model is illustrated in Figure 4(a).The wave simulation takes place only in D and is coupled with the outside medium D such that the resulting wavefield is equivalent to a simulation performed on the entire medium (D and D ).Using a Kirchhoff-Helmholtz-like integral extrapolation the coupling is achieved with pre-computed extrapolation Green's functions Γ: they are used to predict the boundary condition on the injection boundary ∂D i at every time step of the simulation by means of measurements of the pressure p and particle velocity v on the (sound-transparent) extrapolation surface ∂D e .Note that the extrapolation Green's functions Γ, computed in the full medium (D and D ), must be available for all combinations of points between the extrapolation surface and injection boundary.A key feature of IBCs is that the medium can be arbitrarily perturbed in D without the need to recompute the extrapolation Green's functions, which allows for efficient wavefield modeling of a large set of scenarios.
We use an alternative IBC configuration that enables isolation of the effects of a single or multiple interfaces.To explain this configuration we use an auxiliary setup illustrated in Figure 4(b) which uses two truncated domains at the top and bottom, denoted with D 1 and D 2 (and with corresponding injection and extrapolation surfaces D [e,i], [1,2] ).They are chosen such that the remaining medium D consists of three interfaces 'sandwiched' between the two truncated domains.Following the theory given by Vasmel (2016) wave simulations in multiple sub-domains can be linked to each other by evaluating surface integrals for all combinations of extrapolation surfaces and injection Isolation of wavefield constituents Isolation of wavefield constituents boundaries.As sketched in Figure 4(b) for two subdomains this results in two sets of Green's functions for each individual sub-domain (Γ 1,1 and Γ 2,2 ) as well as two sets of cross-relating Green's functions (Γ 2,1 and Γ 1,2 ).Where the IBC boundaries coincide with the model boundaries their contribution to the extrapolation integral is zero, and they can hence be neglected.To link the two truncated domains it is therefore sufficient to use only the two IBC pairs that dissect the model in the center.This is the configuration used in this work, shown in Figure 4(c).
The IBC configuration we propose allows for removing wavefield constituents related to D from the full-wavefield if only the direct arrivals are used in the extrapolation Green's functions.We exploit that any arbitrary medium can be used in D when computing the extrapolation Green's functions (Thomson, 2012;Vasmel, 2016;Elison et al., 2018).If we hence replace D 1 and D 2 by homogeneous half-spaces in this preliminary step, the extrapolation Green's function Γ 2,1 contains the direct arrival from sources at ∂D e,1 to receivers at ∂D i,2 , and a series of internal multiples caused by the three interfaces in D .Instead, we use only the direct arrival D across D from ∂D e,1 to ∂D i,2 , that can be obtained with the MPS method as outlined in this work.With this modification interactions between the wave simulations in D 1 and D 2 contain the correct transmission across D , but do not contain any reflections related to the interfaces located therein.The same reasoning applies to Γ 1,2 .Moreover, our configuration allows us to also reduce Γ 1,1 and Γ 2,2 to the direct arrivals.This ensures that the injection boundaries are non-reflective for outgoing waves [compare with Ting and Miksis (1986)].
Another way of interpreting our IBC configuration is the cloaking application presented by van Manen et al. (2015).If the two extrapolation surfaces and injection boundaries in Figure 4(c) are combined and closed, the configuration is reversed with respect to the original IBC formulation in Figure 4(a): the extrapolation surface becomes the outer surface and the injection boundary the inner surface, which entirely surrounds the domain D over which the extrapolation takes place.Instead of such a closed configuration as used by van Manen et al. (2015), however, here we deploy an open configuration which splits extrapolation surfaces and injection boundaries apart, because it is better suited for geologic scenarios with layered models.Such a modification requires special tapering at the lateral IBC ends to avoid artifacts due to edge effects (Elison et al., 2018).

Equivalence of the two methods
We demonstrate the equivalence of Eqs.(4-6) and the IBC-based cloaking.Consider simulation of a reflection experiment with sources and receivers located at the top of the medium ∂D 0 using the IBC setup outlined in the previous subsection [Figure 4(c)].A schematic illustration of such a simulation is shown in Figure 5, where blue and green arrows denote FD propagation and IBC extrapolation, respectively.Moreover, for simplicity the extrapolation surfaces ∂D e and injection boundaries ∂D i for the upper and lower subdomains have been summarized to ∂D 1 and ∂D 2 .For this example, we only consider two interactions of the upper sub-domain with the lower sub-domain.
The IBC simulation can be expressed analytically using Eqs.(4-6).If we combine these three equations, using i = 2 in Eq. ( 5) and evaluating only the first term of that equation to ensure purely down-and upgoing source and receiver components, respectively, we obtain Isolation of wavefield constituents Isolation of wavefield constituents  which again should be read from right to left.Comparison of this equation with Figure 5 shows the equivalence of both cloaking approaches, provided that the direct arrival D is used for the extrapolation Green's functions Γ in Eq. ( 7).Every term on the right hand-side of that equation corresponds to one arrow reaching the receiver surface ∂D 0 in Figure 5.

NUMERICAL EXAMPLES
We now present numerical examples for the methods we have introduced.In the first subsection pseudo one-way wave propagation in a finite-difference scheme is demonstrated.We compute the direct transmitted wave and the primary reflection events.The model is designed to cause shortperiod reflections for a 20 Hz-Ricker source wavelet.Therefore, a conventional time-windowing approach could not be used to distinguish between the different primary wavefield constituents.
In the second subsection we demonstrate how convolutions of sub-domain-related wavefields can be combined to cloak a certain target.In a first step we use finite-difference reference data to verify our equations, and in a second step we demonstrate the same approach using Marchenko-derived wavefields.The model used for this example contains only well-separated interfaces, since the conventional Marchenko method breaks down for short-period reflections.Although the Marchenko method has been extended to handle such short-period internal multiples (Dukalski et al., 2019; Isolation of wavefield constituents Isolation of wavefield constituents   Elison et al., 2019), this can currently only be solved for a horizontally layered overburden.
The third subsection demonstrates cloaking of an interface pair in an IBC simulation, using the same model that is introduced for the pseudo one-way propagation example.

Computation of the direct transmission and primary reflections
We compute the direct transmission response through a model with thin layers.The velocity model is shown in Figure 6 and the corresponding densities are computed using Gardner's relation (Gardner et al., 1974).Figure 7(a) shows a snapshot of a regular second-order staggered-grid finitedifference simulation (Virieux, 1984) for a dipole source excited at the top of the model.Besides the downward propagating direct wave one can observe a complex pattern of primary reflections and internal multiples propagating between the different interfaces.This is manifested in the wavefield recorded with receivers below the model, shown in subfigure (d), where the direct arrival overlaps with a series of internal multiples.In contrast to that, Figure 7(b) illustrates a composition of multiple subsequent MPS simulations in which the wavefield is propagated from one dashed red surface to the next.Note that this subfigure actually shows a composition of eight subsequent finitedifference runs.No interactions occur between different interfaces, e.g., the primary reflection from the third last interface at about 300 m depth [highlighted with a green arrow in Figure 7(b)] does not 'reach' the overlying interface and therefore does not cause another downward propagating internal multiple.Consequently, the recordings at the model bottom illustrated in Figure 7(e) exhibit a single event only.To support these observations the right panels of Figure 7 [panels (c) and (f)] show the difference between the left and middle panels.Note that absorbing boundaries on either side of the medium (Roden and Gedney, 2000;Komatitsch and Martin, 2007)   recording and injection surfaces in order to reduce edge effects (Elison et al., 2018).In the given example 200 m wide absorbing boundaries are applied on either side of the simulation, denoted with dashed black lines in Figure 7.
We observe that our composite modeling scheme indeed propagates only the first event through the series of interfaces.Figure 7(f) highlights the series of internal multiples that overlap with the direct arrival in the original simulation (panel (d)), and which is not present in panel (e).The traveltime of the direct event, computed with an Eikonal solver, is highlighted with a dashed green line in panels (d) to (f).Comparing this traveltime with the events in panel (f) emphasizes that a time-windowing approach would not be possible in such a scenario with thin layers.At x = 400 m a reduction in energy can be observed in the isolated direct arrival (panel (e)), which is caused by tunneling through thin layers: critically refracted waves cause evanescent waves which propagate along interfaces.Due to sub-wavelength interface separation distances evanescent waves translate again into propagating waves at subsequent interfaces.Ultimately, some slight artifacts visible in panel (e) are caused by scattering off the staircase-like discretization of curved interfaces, although this effect has been significantly suppressed using the approach described by Moczo et al. (2004).
In addition to computing the direct transmission, we can use the recordings from the previous simulations to obtain primary reflection data of the model in Figure 6.This is shown in Figure 8 using two different source signatures: the top row uses a Ricker wavelet with 100 Hz central frequency and the bottom row uses a 20 Hz central frequency.We show results for wavelets of different Isolation of wavefield constituents Isolation of wavefield constituents fequencies to demonstrate the method with and without the presence of short-period internal scattering.Panel (a) shows surface reflection data using a conventional (full-wavefield) FD scheme.Consequently, the data contain primary and multiple scattering.Panel (b) shows only the primary events obtained by stepwise simulation of one-way wave propagation in downward and upward direction.Since we have used a high-frequency source wavelet one can nicely observe eight primary events, related to the eight interfaces in the medium.The difference between the full reflection data and the primary events yields all (first and higher order) multiple scattering (panel (c)).
A more realistic input wavelet for a seismic reflection experiment is the Ricker wavelet with 20 Hz central frequency.The corresponding results are shown in Figure 8(d)-(f).The combination of this source wavelet and the layer thicknesses in the current model clearly leads to overlapping primary reflections.Whilst it might be possible to pick individual events in the full data shown in panel (a), this is certainly not possible for the data in panel (d).Our method is able to correctly simulate the primary reflections only, even when using a wavelet that is in the order of the separation distance of the interfaces.Note that we can even obtain each primary reflection individually [Eq.( 3)].

Cloaking using Marchenko-derived wavefields
We demonstrate the convolutional expressions to cloak interfaces using the velocity model shown in Figure 9 and a density distribution derived from Gardner's relation (Gardner et al., 1974).The model consists of five well-separated (i.e., beyond the dominant wavelength), laterally varying interfaces, where we wish to isolate events due to the target interface at x = 800 m.To this end, two auxiliary depth levels z 1 and z 2 at z = 600 m and z = 1000 m (denoted with dashed lines in Figure 9) are used to derive the individual wavefields required for the convolutional expressions given in the theory section.In a first step, all wavefields related to sub-domains are computed with finite-differences and used in Eqs.(4-6) to obtain synthesized reference data.In a second step, the same wavefields are derived with the Marchenko method from surface data, after which they are re-combined again to isolate events related to an individual interface.The results are shown in Figure 10.For all results presented here we have used i = 4 in Eq. ( 5), which is sufficient for the recording length of interest.
We obtain reference data by synthesizing finite-difference wavefields.The required input are reflection responses R ∪ 2,2 , R ∩ 1,1 and R ∪,(1) 0,0 as well as transmission responses D 12 , D 21 ,T 0,1 and T 1,0 (compare with Figure 3).Using these in Eqs.(4-6) yields the surface reflection data R∪ 0,0 shown in Figure 10(d  Modified surface data can be obtained from conventional surface data using the Marchenko method (Wapenaar et al., 2014).In a first step, Marchenko Green's functions are derived from the acquisition surface at z 0 to two horizontal depth levels at z 1 and z 2 above and below a target interface, respectively.Subsequent multidimensional deconvolution of the Marchenko-derived wavefields yields redatumed reflection responses R ∩ 1,1 and R ∪ 2,2 [see Figure 3(e) and (b)] that illuminate the medium above and below the target interface, respectively.In a similar way, multidimensional deconvolution of Marchenko focusing functions yields the overburden reflection response R ∪,(1) 0,0 [see Figure 3(k)].Apart from the surface reflection data R ∪ 0,0 these steps require a background model of the seismic velocities to obtain traveltimes from the acquisition surface to the dashed surfaces at depth.
The additional wavefields required to re-synthesize surface data are transmission responses.More precisely, these are the transmission across the overburden T 1,0 and T 0,1 [see Eq. ( 6) and Figures 3(h .For an accurate overall scaling additional knowledge of the overburden's transmission amplitude is assumed to be available.The latter transmissions (across the target interface) are obtained using the same background velocity model required for the Marchenko method.However, these (arbitrarily scaled) transmission wavefields require proper overall and offset-dependent scaling.If available, this can be achieved using transmission data or well-log information.Alternatively, the correct scaling may be obtained by comparison with the input reflection data and a matching filter.For the example in this section we have used the location and impedance contrast of the actual target interface.
Synthesized results using Marchenko-derived wavefields are shown in the third row of Figure 10 [i.e., panels (f)-(h)].The results are similar to the reference wavefields shown in the same figure (i.e., the second row) and contain all relevant events.Differences can be observed in the early part of the wavefields (t < 1 s), where the Marchenko-derived reflection data have significantly reduced offset.This is due to tapering that is required in the Marchenko scheme and, more importantly, the multidimensional deconvolution, which is not able to retrieve large offsets at early times.For the application presented here, however, this does not pose a significant limitation since these early times do not contain events related to the target area.
A pronounced difference between panels (e) and (h) is highlighted with black arrows in panel (h), which are remnants of primary reflections of the two bottommost interfaces and that do not cancel out properly for large offsets.The reason for this difference is an inherent limitation of the conventional Marchenko method: as a result of an erroneous estimate of the direct arriving focusing function (that neglects transmission losses) Marchenko Green's functions suffer from erroneous relative amplitudes as a function of offset.For laterally invariant media this error drops out upon multidimensional deconvolution, but for a laterally heterogeneous medium like the one we are investigating, this causes offset-dependent amplitude errors in the Marchenko-derived redatumed data.Consequently, we are using a transmission across the target interface that contains erroneous relative amplitudes when computing the data shown in Figure 10(f) (where the transmission is implicitly contained in the Marchenko-derived R ∪ 2,2 ), but we are using the correct transmission response for the derivation of panel (g).Note that this problem cannot be resolved by comparing the surface data without the target interface [panel (g)] directly with the input reflection data [panel (a)], since the erroneous offset-dependent amplitudes in the Marchenko redatumed data are also affecting all other wavefields used to synthesize the data in panel (g) (and hence other events would not cancel correctly when taking the difference).Moreover, the fact that a smooth velocity model is used to constrain the Marchenko solutions (by using it to estimate the early focusing function) causes additional phase errors in the Marchenko solutions.

Isolation of wavefield constituents Isolation of wavefield constituents
For quality-control of the Marchenko-based results we are comparing individual traces in .Comparing the reflection data of the full medium R ∪ 0,0 in panel (i) confirms the previous observations: synthesized finite-difference reference data and conventional finite-difference data plot on top of each other (the black trace is not visible in this panel), whereas the Marchenkoderived data cannot recover offsets at early times.Therefore the first event has a too small amplitude, while all later events are recovered correctly.The red traces in panels (j) and (k) follow the blue traces closely.However, there are discrepancies which increase with time (and offset).They are related to relative amplitude errors as discussed above, to deconvolution artifacts as well as to refracted waves which are not properly handled in the Marchenko scheme.Nevertheless, a satisfactory match of the data synthesized from Marchenko wavefields with the reference is observed.

IBC-based cloaking
The same model that was used to demonstrate the computation of first-order wavefield constituents (Figure 6) is used to isolate the effects of two interfaces.A reflection experiment with a dipole source and pressure receivers located at the upper surface is simulated.In Figure 6 we have used solid and dashed lines to denote the locations of two injection and extrapolation surfaces which encompass a target reflector pair which dips from left to right for 180 m < z < 250 m.Pre-computed Green's functions at these IBC surfaces enable cloaking of the target interfaces and yield the reflec-Isolation of wavefield constituents Isolation of wavefield constituents tion data recorded at z = 0 shown in Figure 11(b).To obtain the extrapolation Green's functions across the target the MPS-based method introduced in this work has been used.For comparison we show conventional reflection data computed in the full model (containing all interfaces) in panel (a).The difference between the two reflection data, shown in panel (c), contains all events related to the two target interfaces, including both primaries and multiples.Note the reduced offset in these data with respect to the model in Figure 6 which owes to the absorbing boundaries at the lateral model sides.

DISCUSSION
Stepwise computation of first-order wavefield constituents can be interpreted as pseudo one-way wave propagation in a finite-difference scheme.Although the wave simulation propagates the full (two-way) wavefield, the fact that we treat each discontinuity individually in adjacent homogeneous half-spaces allows us to record the one-way transmission and primary reflections.Our approach is limited to piecewise constant layers that are continuous and not crossing.Therefore, also the auxiliary MPS surfaces must not cross any inhomogeneities (although there is no such restriction for general MPS applications).This has implications on the minimal separation distance of interfaces: for the computation of the direct transmission interfaces must be offset far enough to accommodate one MPS surface (compare with the red band in Figure 2).To compute primary reflections, two MPS surfaces must fit between two adjacent interfaces, since injection and recording surface may not coincide.The number of grid nodes associated with an MPS surface depends on the FD stencil, and hence for higher order stencils a larger separation distance is required.
The total number of simulations required to compute the direct arrival of the transmission response is equal to the amount of interfaces in the medium.Note, however, that the computational domain for each simulation is restricted to the area between two subsequent auxiliary MPS surfaces, reducing the computational requirement for each simulation significantly.Yet, to each simulation a 'frame' of absorbing boundaries must be added, which leads to an increased total computational cost as opposed to one 'normal' simulation of the entire medium.The results of the transmission computation can subsequently be used to retrieve primary reflections.Since this is achieved by a single series of transmission simulations (but now in upward direction), this only doubles the computational cost.
Modified reflection data synthesized from Marchenko wavefields suffer from erroneous amplitudes, which limits the offset range for which we can accurately determine events related to a target interface.This is because conventional Marchenko wavefields never account for transmission losses at interfaces, since these would have to be known a-priori.This is different in the augmented Marchenko scheme (Dukalski et al., 2019;Mildner et al., 2019), in which, in addition to the well-known coupled Marchenko equations, energy conservation of the focusing state is evaluated.Deviations from the Marchenko focusing condition enable correction of the Marchenko wavefields such that they properly account for transmission losses.Currently, however, this correction is only possible for media which are horizontally layered between the acquisition surface and the focusing depth level.
Besides seismic reflection data, additional transmission amplitude information is required when Marchenko wavefields are re-combined to surface data.For the numerical example given in this work, we assumed knowledge of such transmission information, which in a real application may be available from well logs.Alternatively, such amplitude information could be inferred from comparison with the initial surface reflection data, using (adaptive) subtraction.

Isolation of wavefield constituents Isolation of wavefield constituents
The use of an open IBC configuration with two separated pairs of IBC surfaces enables cloaking applications in layered media.As outlined in the theory section this is not possible using the original IBC-cloaking theory introduced by van Manen et al. ( 2015), because it is neither possible to entirely surround an interface that spans the model horizontally, nor to remove a single interface.In this work we therefore split the extrapolation and injection surfaces into two separated IBC pairs.For cloaking of the intermediate medium in such a configuration the use of extrapolation Green's functions containing only the direct arrival becomes key.This leads to a significant speed-up in computation time for the computation of these Green's functions, since the computational domain can be limited to the medium between z 1 and z 2 .Nevertheless, considerable memory ressources are required for IBC computations (Broggini et al., 2017).
Our choice of model dissection makes the IBC-based cloaking application especially suited for geological applications.Scattering of a circular object can be isolated by simply 'removing' the object from a simulation.Alternatively, surrounding such an object by closed IBCs surfaces can also achieve cloaking (van Manen et al., 2015), even beyond numerical simulations (Börsing et al., 2019).However, both approaches break down for cloaking of interfaces in layered media, since an interface can neither simply be removed, nor can it be entirely surrounded by IBC surfaces.
Using open MPS and IBC surfaces comes at the cost of reduced precision, since both methods are only numerically accurate for closed surfaces (van Manen et al., 2007;Amundsen and Robertsson, 2014).Neverthless, applying absorbing boundary conditions (Komatitsch and Martin, 2007) at the lateral ends allows us to successfully suppress artifacts caused by edge effects.However, in that boundary region we therefore can not recover accurate wavefields.

CONCLUSIONS
By combining conventional two-way FD solvers for acoustic wave propagation with advanced methods for wavefield injection and extrapolation we have shown how to obtain isolated wavefield constituents.Given a layered medium we can introduce auxiliary MPS surfaces which allow simulations of wave propagation across individually interfaces.This results in pseudo one-way wave propagation in a finite-difference scheme.We have shown how this can be used to obtain the direct arrival of the transmission response and primary reflections only.The transmission includes tunneled waves that occur if interfaces are separated by less than the dominant wavelength of the signal.Surface reflection data, for which a single (or several) target interface(s) are cloaked, can be expressed as combinations of measured reflection and transmission responses, related to different sub-domains.Comparison with data obtained in the full medium allows isolation of events related to selected interfaces.Implicit model knowledge via measured reflection data is sufficient, since subdomain related wavefields can be obtained using the Marchenko method.However, a background velocity model and transmission amplitudes must be known.
Cloaking of a target area can also be achieved directly in a two-way wave propagation simulation, using an IBC setup with a configuration that is reversed with respect to the original IBC definition.They key of this method lies in the use of the direct arrivals only for the extrapolation Green's functions.If the target contains closely-spaced interfaces the direct transmission can be obtained with the MPS-based pseudo one-way wave simulation introduced in this work.

Isolation of wavefield constituents
The first step consists of deriving a downward reflection response across the target area].This requires a downward reflection response just below the target [Figure3(b)] and, similar as in the previous subsection, the direct transmission responses between ∂D 1 and ∂D 2 [Figure 3(a) and (c)], summarized by (e)-(g), two single-sided reflection responses are combined to yield a Green's function G1,1 [Figure 3(g)] which does not contain reflections related Isolation of wavefield constituents Isolation of wavefield constituents

Figure 3 :
Figure 3: Illustration of Eqs.(4-6) to derive reflection data without wavefield constituents related to a single interface or a stack of interfaces.(a) to (d) illustrate Eq. (4), (e) to (g) illustrate Eq. (5) and (h) to (k) illustrate Eq. (6).The solid red line in panels (a) and (c) denotes the interface of interest that should be removed from the data.It is drawn as dashed line for wavefields where its related events are absent.

Figure 4 :
Figure 4: Different IBC configurations.(a) conventional configuration for a single truncated domain as introduced by van Manen et al. (2007); (b) two truncated IBC domains as described by Vasmel (2016), arranged such that the outside domain D contains a package of three reflectors; (c) open IBC configuration deployed in this work, where the outer surfaces with zero contribution to the extrapolation integral are removed.This setup allows for a straightforward cloaking of the intermittent medium D .Green arrows illustrate combinations of extrapolation surfaces and injection boundaries.

Figure 5 :
Figure 5: Schematic illustration of IBC-based cloaking in an FD simulation of surface reflection data.The schematic contains two interactions of the overburden with the underburden.Blue and green arrows denote FD propagation and IBC extrapolation, respectively.

Figure 6 :
Figure 6: Velocity model used for the model-driven finite-difference examples.The dashed white and solid pink lines denote the IBC extrapolation surfaces and injection boundaries, respectively.

Figure 7 :
Figure 7: Computing the direct transmission.(a) Snapshot at t = 0.18 s for a dipole source located at the model top.The pressure recordings at the model bottom is shown in panel (d); (c) composition of eight subsequent MPS simulations across the individual interfaces which yield only the direct arrival at the model bottom shown in panel (e).The dashed red lines indicate the auxiliary MPS recording and injection surfaces (compare with Figure 1); (c) difference between (a) and (b); (f) difference between (d) and (e).The vertical dashed black lines in (b), (c), (e) and (f) mark the onset of 200 m wide absorbing boundaries which reduce edge effects.The dashed green line in panels (d)-(f) marks the travel time of the direct transmission, obtained from an Eikonal solver.All results are displayed with a t 2 -gain, panels (a)-(c) are clipped to 20 %, panels (d)-(f) to 5 % of the maximum amplitude.
overlap with the open MPS Isolation of wavefield constituents Isolation of wavefield constituents

Figure 8 :
Figure 8: Computation of primary reflections for the model shown in Figure 6.(a) Full-wavefield surface reflection data for a source with a 100 Hz Ricker wavelet; (b) Primary reflections only; (c) difference between (a) and (b); (d)-(f) are like (a)-(c), but for a source with a 20 Hz Ricker wavelet.Vertical dashed black lines denote the onset of the absorbing boundaries, which are five times broader for the low-frequency example due to a coarser discretization.All results are clipped to 3 % of the maximum amplitude and displayed with a t 2 -gain.

Figure 9 :
Figure 9: Velocity model used for the Marchenko-based interface removal.The dashed black lines indicate the two Marchenko focusing depth levels at z 1 and z 2 that encompass the target interface.
), which does not include the target interface.Additionally, we compute R ∪ 1,1 and Isolation of wavefield constituents Isolation of wavefield constituents
) and (j)] as well as the direct arrivals of the transmissions across the target interface D 2,1 and D 1,2 [see Eq. (4) and Figure3(a) and (c)].The former transmission responses are derived by inverting the downgoing Marchenko focusing function, commonly denoted as f + 1

Figure 11 :
Figure 11: Synthetic reflection data for the model shown in Figure 6 for a dipole source and pressure receivers at the surface.(a) Data computed with a conventional finite-difference scheme; (b) data obtained using IBC-based cloaking of two central interfaces removed; (c) difference between (a) and (b).All results are displayed with a t 2 -gain and are clipped to 5 %.
Isolation of wavefield constituents