Synthesizing individual wavefield constituents (such as primaries, first-order scattering, and free-surface 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 found how this can be achieved for (nonhorizontally) piecewise constant layered media. For example, we were 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. Our methods required detailed knowledge of discretized medium parameters. Alternatively, if a medium is known only implicitly via recordings of reflection data, interface-related events can still be isolated through a combination of subdomain-related wavefields. We found how Marchenko redatuming can be used to derive these, which enables data-driven identification (and removal) of interface-related events from surface data.


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., such as the reflectivity method by Kennett, 1984), 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 variation with 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 interface (or a stack of interfaces), including primaries and multiples. This is achieved by cloaking 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 are 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 subhorizontally layered media, we split a wave simulation into a series of consecutive, smaller subsimulations in which every subsimulation 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 (MPSs) (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 in which the interfaces are separated by distances in the order of the signal’s wavelength or less. Because 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 interface (or a stack of interfaces) 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 the case of at least two closely spaced target interfaces, the 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. They reveal the impact of an interface on seismic data, can be used to identify unknown events and strong multiple generators, or, as outlined above, can be used for reservoir characterization. We introduce two methods that achieve such cloaking. First, we present convolutional expressions that allow us to express surface reflection data in terms of wavefields that are related to different subdomains. 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 containing no (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 an 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 because 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).

The equivalency of the data-driven and the model-driven cloaking approaches will be demonstrated in the “Theory” section. We want to emphasize, however, that the two approaches solve very different problems. The model-driven approach is a purely numerical modeling process, whereas the data-driven approach can be understood as a data-processing method applied to (measured) reflection data.


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 an 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 xi denotes the spatial locations on a surface Di. All wavefields in this work denote pressure-normalized monopole recordings due to dipole sources, and we assume that proper scaling of these wavefields allows us to neglect scale factors.

We introduce notations for Green’s functions with different acquisition geometries. Green’s functions that are recorded with sources and receivers on a scattering-free surface Di with an adjacent homogeneous half-space above or below are denoted as reflection responses Ri,i and Ri,i, respectively. If the medium is truncated with another homogeneous half-space at Dj, we denote this with an additional superscript in brackets Ri,i(j). A reflection response recorded at D0 that contains only primary reflections is denoted with P=R0,0(primaries), and if it contains only a primary reflection due to a single interface this is specified with a superscript, for example, P(3) denotes a primary reflection due to the third interface. If sources and receivers are located on different surfaces Di and Dj, respectively, with adjacent homogeneous half-spaces on both of these surfaces, we denote the corresponding Green’s functions as transmission responses Tj,i, and if such a transmission contains only the direct (one-way) transmission, it is denoted with Dj,i. Ultimately, if at least one of the (source or receiver) surfaces is located inside the medium (i.e., without an adjacent homogeneous half-space at that surface), we denote the Green’s functions with Gj,i.

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 1a, which is assumed to be reflection-free at its upper and lower boundaries. Simulation of a source excited at the upper boundary D0 yields primary reflections P at the top and a directly transmitted wave D3,0 at the bottom D3, 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 first-order 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 numerically implement Kirchhoff-Helmholtz integrals, using repeated recordings and injections of two-way wavefields along arbitrarily curved (i.e., nonhorizontal) auxiliary interfaces, assuming piecewise-constant, layered media with continuous, noncrossing layer interfaces. Injection and recording of wavefields is achieved with the MPS method (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 in which this is different: Thomsen et al. (2018) show 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) presents 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 vx and vz 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 the associated pressure 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. Inline particle velocity nodes are not used for the MPS method; hence, they are illustrated in gray even if they lie within the red area. The 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 nx and nz 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 the horizontal and vertical directions 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 the x-direction points in negative direction; hence, nx 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 use in this work. General expressions were given by Vasmel and Robertsson (2016), who show that the approach for the second-order stencil is accurate to numerical precision. The same authors show 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 weaker.

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; that is, 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 because we assume piecewise constant layered media. With D0 defined as the uppermost (acquisition) surface and auxiliary surface i being located below the ith interface, we can express the direct arrival of the downward transmission response through a stack of n interfaces as (compare with Kennett, 1984)  
where refers to a series of multidimensional convolutions according to equation 1, which we summarize with a one-way propagator W. Note that subsequent terms are multiplied from the left because the convolutions are not commutative but the order matters. For the three-interface example in Figure 1, two auxiliary surfaces D1 and D2 are introduced; see Figure 1b1d. Consequently, according to equation 2, the direct downgoing transmission is given by D3,0=T3,2T2,1T1,0.

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 (Figure 1b), a source is excited on the (reflection-free) top of the model at D0, similar to the simulation illustrated in Figure 1a. 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 D1,0=T1,0 is observed below the interface, where it is recorded along the scattering-free auxiliary MPS surface D1. 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 Figure 1c). The previously recorded (purely downgoing) wavefield is used as the source wavefield and is injected along the surface it was recorded on: The auxiliary recording surface D1 from the previous simulation is now the injection surface. Again, the transmission through the current interface is recorded on another auxiliary surface D2 just below (and now D2,0=T2,1T1,0 is recorded). The upward (primary) reflection of this second interface does not cause an internal multiple because the overlying first interface is not present in this simulation (in contrast to the simulation of the full transmission response in Figure 1a). This procedure is repeated for every subsequent interface. For the example in Figure 1, this requires one additional simulation (Figure 1d). The wavefield recorded below the last interface contains only the direct arrival D3,0 across all interfaces but does not contain any internal multiples.

Primary reflection computation

One-way propagators such as equation 2 can be used to compute primary reflections. Upward and downward one-way transmissions, combined with a reflection computed on single interfaces, yield an expression for the primary reflection due to interface n+1 according to (e.g., compare with Berkhout, 1980)  

Implementation of equation 3 with MPS surfaces works similar to the computation of the direct transmission. As illustrated in Figure 1b1d, 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 W0,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 upward through the second-to-last interface. In subsequent simulations, the injected wavefield is always a combination of the primaries that have been propagated back upward 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 transmission of the (total) first-order reflected fields. Alternatively, every primary event can be propagated upward 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 cannot coincide; hence, we cannot 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 R0,0 yields all reflected events related to the target interfaces. The methods are based on dissecting a layered medium into three different subdomains located on top of each other. The central submedium 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 a 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 D1 and D2 at depth levels z1 and z2; 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).

The first step consists of deriving a downward reflection response across the target area (Figure 3d). This requires a downward reflection response just below the target (Figure 3b) and, similar to the previous subsection, the direct transmission responses between D1 and D2 (Figure 3a and 3c), summarized by  
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 D2,1=T2,1; however, an expression such as equation 2 is needed for a target area containing more than a single interface.
In the subsequent step, illustrated in Figure 3e3g, two single-sided reflection responses are combined to yield a Green’s function G˜1,1 (Figure 3g) that does not contain reflections related to the target, according to  
Powers of i indicate i-times repeated multidimensional convolutions of the respective terms, which corresponds to ith order multiple interactions between the upper and lower submedium. Note that the first term on the right side always has downgoing source and upgoing receiver components, and vice versa for the second term. The third term has upgoing source and receiver components, and the fourth term has downgoing source and receiver components.
A third step ultimately enables synthesizing surface reflection data. Combining the previously obtained Green’s function G˜1,1 with the overburden’s transmission responses T1,0 and T0,1 (Figure 3h and 3j), it can be added to the reflection response of the overburden R0,0,(1) (Figure 3k). As a result, a surface reflection response R˜0,0 in which 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 G˜1,1, respectively, which can be obtained by evaluating only the first term in equation 5.

If data are available from a seismic reflection experiment, Marchenko inversion (Wapenaar et al., 2014) offers a way to derive the subdomain-related wavefields used in equations 46. Using the reflection response of the real medium R0,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 R0,0,(1), R1,1, and R2,2 related to individual subdomains. The overburden’s transmission response T0,1 and T1,0 is obtained by inverting the 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 D1,2 and D2,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 subdomain 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 introducing the alternative implementation that we use, which enables cloaking of a target area in layered media.

IBCs are time-dependent boundary conditions that enable the immersion of a wavefield propagating on a truncated subdomain into an arbitrary larger domain (van Manen et al., 2007; Vasmel et al., 2013). They are based on an implementation of nonreflecting 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 De and an outer injection boundary Di, where the latter constitutes the actual boundary of the subdomain D. An arbitrary IBC subdomain of a layered model is illustrated in Figure 4a. 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 precomputed extrapolation Green’s functions Γ: They are used to predict the boundary condition on the injection boundary Di 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 De. Note that the extrapolation Green’s functions Γ, computed in the full medium (DandD), must be available for all combinations of points between the extrapolation surface and the 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 4b, which uses two truncated domains at the top and bottom, denoted with D1 and D2 (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 subdomains can be linked to each other by evaluating surface integrals for all combinations of extrapolation surfaces and injection boundaries. As sketched in Figure 4b for two subdomains, this results in two sets of Green’s functions for each individual subdomain (Γ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; hence, they can 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 4c.

The IBC configuration that 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 the fact 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 D1 and D2 by homogeneous half-spaces in this preliminary step, the extrapolation Green’s function Γ2,1 contains the direct arrival from sources at De,1 to receivers at Di,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 De,1 to Di,2, that can be obtained with the MPS method as outlined in this work. With this modification, interactions between the wave simulations in D1 and D2 contain the correct transmission across D, but they 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 nonreflective for outgoing waves (compare to 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 4c are combined and closed, the configuration is reversed with respect to the original IBC formulation in Figure 4a: The extrapolation surface becomes the outer surface, and the injection boundary becomes the inner surface, which entirely surrounds 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 that splits the 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 equations 46 and the IBC-based cloaking and consider simulation of a reflection experiment with sources and receivers located at the top of the medium D0 using the IBC setup outlined in the previous subsection (Figure 4c). A schematic illustration of such a simulation is shown in Figure 5, where the blue and green arrows denote FD propagation and IBC extrapolation, respectively. Moreover, for simplicity, the extrapolation surfaces De and injection boundaries Di for the upper and lower subdomains have been summarized to D1 and D2. For this example, we only consider two interactions of the upper subdomain with the lower subdomain.

The IBC simulation can be expressed analytically using equations 46. If we combine these three equations, using i=2 in equation 5 and evaluating only the first term of that equation to ensure purely down- and upgoing source and receiver components, respectively, we obtain  
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 equation 7. Every term on the right side of that equation corresponds to one arrow reaching the receiver surface D0 in Figure 5.


We now present numerical examples for the methods we have introduced. In the first subsection, pseudo one-way wave propagation in an FD scheme is demonstrated. We compute the direct transmitted wave and the primary reflection events. The model is designed to cause short-period 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 subdomain-related wavefields can be combined to cloak a certain target. In a first step, we use FD 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 because 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; 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 7a shows a snapshot of a regular second-order staggered-grid FD 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 Figure 7d, where the direct arrival overlaps with a series of internal multiples. In contrast to that, Figure 7b illustrates a composition of multiple subsequent MPS simulations in which the wavefield is propagated from one dashed gray surface to the next. Note that this subfigure actually shows a composition of eight subsequent FD runs. No interactions occur between different interfaces, for example, the primary reflection from the third-to-last interface at an approximately 300  m depth (highlighted with a green arrow in Figure 7b) 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 7e exhibit a single event only. To support these observations, Figure 7c and 7f shows 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) overlap with the open MPS recording and injection surfaces 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 the 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 7f highlights the series of internal multiples that overlap with the direct arrival in the original simulation (Figure 7d), which is not present in Figure 7e. The traveltime of the direct event, computed with an eikonal solver, is highlighted with a dashed green line in Figure 7d7f. Comparing this traveltime with the events in Figure 7f 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 (Figure 7e), which is caused by tunneling through thin layers: Critically refracted waves cause evanescent waves that propagate along the interfaces. Due to subwavelength interface separation distances, evanescent waves translate again into propagating waves at subsequent interfaces. Ultimately, some slight artifacts visible in Figure 7e are caused by scattering off the staircase-like discretization of curved interfaces, although this effect has been significantly suppressed using the smoothing 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: In Figure 8a8c a Ricker wavelet with a 100 Hz-central frequency is used, and in Figure 8d8f a 20 Hz-central frequency is used. We show results for wavelets of different frequencies to demonstrate the method with and without the presence of short-period internal scattering. Figure 8a shows surface reflection data using a conventional (full-wavefield) FD scheme. Consequently, the data contain primary and multiple scattering. Figure 8b shows only the primary events obtained by stepwise simulation of one-way wave propagation in the downward and upward directions. Because we have used a high-frequency source wavelet one can 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 (Figure 8c).

A more realistic input wavelet for a seismic reflection experiment is the Ricker wavelet with a 20 Hz-central frequency. The corresponding results are shown in Figure 8d8f. The combination of this source wavelet and the layer thicknesses in the current model clearly leads to overlapping primary reflections. Although it might be possible to pick individual events in the full data shown in Figure 8a, this is certainly not possible for the data in Figure 8d. 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 (equation 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, in which we wish to isolate events due to the target interface at x=800  m. To this end, two auxiliary depth levels z1 and z2 at z=600  m and z=1000  m (denoted with the 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 subdomains are computed with FDs and used in equation 46 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 recombined 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 equation 5, which is sufficient for the recording length of interest.

We obtain reference data by synthesizing FD wavefields. The required inputs are the reflection responses R2,2, R1,1, and R0,0,(1) as well as the transmission responses D12, D21,T0,1, and T1,0 (compare with Figure 3). Using these in equations 46 yields the surface reflection data R˜0,0 shown in Figure 10d, which does not include the target interface. Additionally, we compute R1,1 and use it to replace R˜1,1 (from equation 4), which allows us to synthesize surface reflection data of the full medium R0,0 (including the target interface) as shown in Figure 10c. The difference between these two reflection responses R˜0,0 and R0,0, shown in Figure 10e, contains all reflections related to the target interface. To validate our equations, we compare the synthesized reflection data of the full medium (Figure 10c) with conventional FD data computed in the original model (Figure 10a). Inspecting the difference between the two in Figure 10b gives us confidence in our method.

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 z0 to two horizontal depth levels at z1 and z2 above and below a target interface, respectively. Subsequent multidimensional deconvolution of the Marchenko-derived wavefields yields redatumed reflection responses R1,1 and R2,2 (see Figure 3e and 3b) 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 R0,0,(1) (see Figure 3k). Apart from the surface reflection data R0,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 resynthesize surface data are transmission responses. More precisely, these are the transmissions across the overburden T1,0 and T0,1 (see equation 6 and Figure 3h and 3j) as well as the direct arrivals of the transmissions across the target interface D2,1 and D1,2 (see equation 4 and Figure 3a and 3c). The former transmission responses are derived by inverting the downgoing Marchenko focusing function, commonly denoted as f1+. For 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 correct 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 Figure 10f10h. The results are similar to the reference wavefields shown in Figure 10c10e 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 a significantly reduced offset. This is due to the different wavenumber content of different events, 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 because these early times do not contain events related to the target area.

A pronounced difference between Figure 10e and 10h is highlighted with the black arrows in Figure 10h, which are remnants of the primary reflections of the two bottom-most interfaces 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 such as 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 10f (where the transmission is implicitly contained in the Marchenko-derived R2,2), but we are using the correct transmission response for the derivation of Figure 10g. Note that this problem cannot be resolved by comparing the surface data without the target interface Figure 10g directly with the input reflection data Figure 10a because the erroneous offset-dependent amplitudes in the Marchenko redatumed data are also affecting all other wavefields used to synthesize the data in Figure 10g (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.

For quality control of the Marchenko-based results, we are comparing individual traces in Figure 10i10k. Comparing the reflection data of the full medium R0,0 in Figure 10i confirms the previous observations: Synthesized FD reference data and conventional FD data plot on top of each other (the black trace is not visible in this panel), whereas Marchenko-derived data cannot recover offsets at early times. Therefore, the first event has a too-small amplitude, whereas all later events are recovered correctly. The red traces in Figure 10j and 10k follow the blue traces closely. However, there are discrepancies that increase with time (and offset). They are related to relative amplitude errors as discussed above, to deconvolution artifacts, as well as to refracted waves that 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 the 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 that encompass a target reflector pair that dips from left to right for 180  m<z<250  m. Precomputed Green’s functions at these IBC surfaces enable cloaking of the target interfaces and yield the reflection data recorded at z=0  m shown in Figure 11b. 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 Figure 11a. The difference between the two reflection data, shown in Figure 11c, contains all events related to the two target interfaces, including primaries and multiples. Note the reduced offset in these data with respect to the model in Figure 6, due to the absorbing boundaries at the lateral model sides.


Stepwise computation of first-order wavefield constituents can be interpreted as pseudo one-way wave propagation in an FD 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, the auxiliary MPS surfaces must also 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 because the injection and recording surface may not coincide. The number of grid nodes associated with an MPS surface depends on the FD stencil; 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. Because this is achieved by a single series of transmission simulations (but now in the 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 because 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 that 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 recombined 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.

The use of an open IBC configuration with two separated pairs of IBC surfaces enables cloaking applications in layered media in numerical wave simulations. 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 a layer 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 because the computational domain can be limited to the medium between z1 and z2. Nevertheless, considerable memory resources are required for IBC computations (Broggini et al., 2017).

Our choice of model dissection makes the IBC-based cloaking application especially suited for geologic applications. Scattering of a circular object can be isolated by “removing” the object from a simulation. Alternatively, surrounding such an object by closed IBC surfaces can also achieve cloaking (van Manen et al., 2015), even beyond numerical simulations (Börsing et al., 2019). However, in contrast to the cloaking method presented here, both approaches break down for cloaking of interfaces in layered media because an interface can neither be removed, nor can it be entirely surrounded by IBC surfaces.

Using open MPS and IBC surfaces comes at the cost of reduced precision because both methods are only numerically accurate for closed surfaces (van Manen et al., 2007; Amundsen and Robertsson, 2014). Nevertheless, 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 cannot recover accurate wavefields.


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 that allow simulations of wave propagation across individually interfaces. This results in pseudo one-way wave propagation in an FD 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 subdomains. 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 because subdomain-related wavefields can be obtained using the Marchenko method. However, a background velocity model and the 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. The key to 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.


This work has received funding from the European Union’s Horizon 2020 research and innovation program under Marie Skłodowska-Curie grant agreement no. 641943, as well as from the European Research Council under grant agreement no. 694407. We thank the three anonymous reviewers for their constructive comments that helped improve this paper.


Data associated with this research are available and can be obtained by contacting the corresponding author.

Freely available online through the SEG open-access option.