Unified elimination of 1D acoustic multiple reflection

Migration, velocity and amplitude analysis are the employed processing steps to find the desired subsurface information from seismic reflection data. The presence of free‐surface and internal multiples can mask the primary reflections for which many processing methods are built. The ability to separate primary from multiple reflections is desirable. Connecting Marchenko theory with classical one‐dimensional inversion methods allows to understand the process of multiple reflection elimination as a data‐filtering process. The filter is a fundamental wave field, defined as a pressure and particle velocity that satisfy the wave equation. The fundamental wave field does not depend on the presence or absence of free‐surface multiples in the data. The backbone of the filtering process is that the fundamental wave field is computed from the measured pressure and particle velocity without additional information. Two different multiples‐free datasets are obtained: either directly from the fundamental wave field or by applying the fundamental wave field to the data. In addition, the known schemes for Marchenko multiple elimination follow from the main equation. Numerical examples show that source and receiver ghosts, free‐surface and internal multiples can be removed simultaneously using a conjugate gradient scheme. The advantage of the main equation is that the source wavelet does not need to be known and no pre‐processing is required. The fact that the reflection coefficients can be obtained is an interesting feature that could lead to improved amplitude analysis and inversion than would be possible with other processing methods.


I N T RO D U C T I O N
The ability to separate primary and multiple reflections in acoustic data is desired to find the subsurface velocity distribution and to create an image of the subsurface. Two of the older methods that perform one-dimensional (1D) multiple elimination are the forward recursion method (Kunetz, 1964), where a fundamental wave field was defined, and the backward recursion method (Robinson and Treitel, 1978) that uses the same fundamental wave field. They worked with the Goupillaud model and included free-surface multiples in their analysis with z-transforms. The goal of Kunetz was to remove * E-mail: e.c.slob@tudelft.nl overburden multiples from deeper reflection events. The goal of Robinson and Treitel was to recover the reflection coefficients as a function of travel time. Many velocity analysis and migration techniques were developed that assume all events in the data to be primary reflections (Weglein, 2016). The effect of the pressure-free surface on land and marine data is usually strong, and free-surface multiples can make it very difficult to identify primary reflections. For this reason, several methods have been developed to predict and remove free-surface multiples (Verschuur et al., 1992;van Borselen et al., 1996;Weglein et al., 1997;van Groenestijn and Verschuur, 2009;Amundsen, 2001). Internal multiples can mask primary reflections as well, and several methods have been proposed to attenuate them (Araújo et al., 1994;Weglein et al., 1997;Jakubowicz, 1998; 327 ten Kroode, 2002;Löer et al., 2016). Primary and multiple reflections can be used together in imaging (Brown and Guitton, 2005;Davydenko and Verschuur, 2017), or in a joint migration inversion scheme (Masaya and Verschuur, 2018).
The fact that the Marchenko equation and the Kunetz theory are connected was shown by Ware and Aki (1969) and recently illustrated in detail by Bardan and Robinson (2018). Slob et al. (2020) unified two bodies of thought. They showed that the fundamental wave field of the recursive methods of Kunetz (1964) and Robinson and Treitel (1978) composes the kernel of the Marchenko equation. This means that the Marchenko focusing functions, when projected back to the acquisition surface as proposed by van der Neut and Wapenaar (2016), and the fundamental wave fields are similar. The fundamental wave fields provide the ability to eliminate multiple reflections as part of a processing sequence. In this paper, we use the term attenuate when the theory is approximate and multiples cannot be totally removed even in a synthetic data test, while we use the term eliminate when the multiples are totally removed in a synthetic data test. Note that even in one dimension the finite bandwidth of the data results in incomplete removal of multiples in thin beds that are below the resolution limit posed by the bandwidth. In three dimensions, all known theories have principal limitations. Applied to field data, the distinction between attenuation and elimination may become irrelevant, because other effects could be larger than those caused by the differences between different methods.
The interest in Marchenko theory for seismic imaging was revived by understanding how subsurface reflection information is encoded in acoustic data, which has led to datadriven receiver redatuming requiring minimal knowledge of the subsurface (Broggini et al., 2012;Wapenaar et al., 2013). Several methods have been derived to create images without artefacts from internal multiples Slob et al., 2014;Wapenaar et al., 2014) with extensions to dissipative media (Slob, 2016;Cui et al., 2018) and elastic media (da Costa Filho et al., 2014;Wapenaar, 2014), and connections to reverse time migration . This has become known as the Marchenko method. Methods to create images while removing internal as well as free-surface multiples have been developed and are based on the same concept (Singh et al., 2015(Singh et al., , 2017Ravasi, 2017). A single Marchenko method does not exist, and we consider this understanding a concept rather than a method. Many other methods have been derived from this concept, such as fracture compliance characterization (Minato and Ghose, 2016), target-oriented velocity analysis (Mildner et al., 2017), extended imaging and inversion (van der Neut et al., 2017b;Cui, 2020) and monitoring , homogeneous Green's function retrieval , equations for inverse source problems (van der Neut et al., 2017a), virtual acoustics  and immersive wave simulation (Elison et al., 2018). Most field test results are with two-dimensional (2D) implementations on line data (Jia et al., 2018;Staring et al., 2018;Brackenhoff et al., 2019). The effect of line data on 2D methods and out-of-plane effects as well as sampling conditions for three-dimensional (3D) applications have been analysed (Jia et al., 2019;Lomas and Curtis, 2020). A 3D example of approximate double focusing can be found in Staring and Wapenaar (2020). Review, comparison and introductory publications on Marchenko redatuming and imaging are available in the literature Nowack and Kiraz, 2018;Lomas and Curtis, 2019).
The Marchenko redatuming concept has been used for internal multiple attenuation in combination with seismic interferometry (Meles et al., 2015;da Costa Filho et al., 2017;Liu et al., 2018). van der Neut and Wapenaar (2016) proposed to project the direct part of the down-going Marchenko focusing function back to a receiver at the acquisition surface. This idea eliminates the need to estimate the unknown initial part of the focusing function. In their version, limited model information is still necessary to find the offset-dependent truncation time instants. Zhang and Staring (2018) found that this information is not necessary, and the truncation time can be chosen independent of offset. Their work filters out primary reflections from acoustic data without using any model information and was successfully applied to field data (Zhang and Slob, 2020a). The scheme was adapted to remove free-surface and internal multiples simultaneously . Recently, a method for the retrieval of plane wave primary reflection responses with various angles from 3D data has been developed (Meles et al., 2020). These methods are relevant for future extensions beyond normal incidence treated in this paper. We unify, in 1D, the known Marchenko multiple elimination and include schemes for which redatuming equations exist when data are collected at the free surface (Singh et al., 2017) and in a marine setting (Ravasi, 2017). We start with the pressure and particle velocity field as data from which we want to retrieve only primary reflection events.
In this paper, we build our equations and the multiple elimination methods on the concept of the fundamental wave field. We show here that the fundamental wave field contains information of only the subsurface impulse reflection response. We derive new equations using two-way wave fields and reciprocity (Kang, 2018) and derive all known schemes from these equations. The paper is organized as follows: First, The model configuration for different settings with the receiver always at z = 0; (a) the marine setting, with the free surface at z f s < 0, where the acoustic pressure p(t ) and particle velocity v z (t ) can be measured, (b) as in (a) but now with up-and down-going pressure components, p − and p + , respectively, (c) as in (b) but now with the free surface with reflection coefficient r at z = 0, explicit expression of the initial downgoing impulse, δ(t ), and the reflected waves described by the impulse reflection response, R(t ), and (d) the setting without free surface and the up-going field is the subsurface impulse reflection response R 0 . we define the model underlying the acoustic reflection experiment and give expressions for the pressure and vertical component of the particle velocity. These two wave fields represent the data we work with and we show them after three different levels of data processing. Second, the fundamental wave field is defined and characterized as having a pressure and particle velocity field satisfying the wave equation. Third, four filter schemes are derived to eliminate multiple reflections from the data. Each scheme uses the data after a different level of processing, but uses always the same filters. These filters are the fundamental wave fields. Fourth, numerical aspects are discussed and illustrated.

DATA C O N F I G U R AT I O N S
In this section, we describe the acoustic experiment and the associated wave fields that are generated in a model configuration. The receiver level is at depth level z = 0 and the medium below the receiver level, for z > 0, is an arbitrary heterogeneous medium. In the marine setting, the pressure-free surface would be located at some depth level above the receiver level, at z f s < 0, and the source would be at the source level z s in between the free surface and the receiver level. The receivers could measure the pressure p(t ) and the particle velocity v z (t ), where t denotes the free time parameter that is used to measure time during an experiment. The pressure can be defined through its up-and down-going parts, given by p − and p + , respectively, as p = p + + p − . The particle velocity is then given by v z = Z −1 0 (p + − p − ), where Z 0 = ρ 0 c 0 is the acoustic impedance with mass density ρ 0 and velocity c 0 at the receiver level (see, e.g., Chapter 3 in Wapenaar and Berkhout, 1989).
Pressure and particle velocity form the measurable data of the first model. This configuration is depicted in Fig. 1(a).
The next step is to obtain the up-and down-going parts of the pressure from the data as This is the data of the second model for data where the measurement is decomposed into up-and down-going pressure fields. This configuration is depicted in Fig. 1(b). Note that this step can be carried out without information about the source time signature, but the medium properties at the receiver level must be known.
To continue with more processing steps, we need to know how the up-and down-going parts of the pressure are related to the action of the source and the reflection response of the medium. Let us assume that the initial pressure is generated by a monopole source with a source time function, W (t ). We can write the pressure at the receiver in its components as (see, e.g., Lindsey, 1960) p in which * denotes temporal convolution, r is the reflection coefficient of the free surface to an incoming pressure wave, R denotes the impulse reflection response including free-surface effects and t f s = |z f s |/c 0 is the one-way travel time from the receiver level to the free surface. The timing effects of the source being located above the receiver is incorporated in the effective source function S, which is given by where the arrival time of the direct wave from source to receiver is t − = |z s |/c 0 and of the direct wave via the free surface is t + = (z s − 2z f s )/c 0 . The factor Z 0 /2 in the effective source function is the scale factor for acoustic pressure generated by a monopole source. Note that in equations (3) and (4), R is the impulse reflection response including free-surface multiple reflections. It is defined as the up-going field at the receiver level and represents the response generated by a down-going impulse at the receiver level and at zero time. On the other hand, p − is a reflection response and, like R, defined at the receiver level. It represents a response that is generated by a general down-going wave field. We use R for an impulse reflection response and p − for a general reflection response throughout the paper. The third step is to remove the source time signature and ghost and to redatum the source and the receiver to the free surface, which we then choose to be at z = 0. In that case, we need to deconvolve the up-and down-going pressure by S(t ) and shift the impulse reflection response in p + by 2t f s . After these steps, equations (3) and (4) reduce to Equations (6) and (7) describe the down-and up-going pressure when the acquisition plane is at the free surface. The right-hand sides of equations (6) and (7) describe the data for the third model where the source time signature has been deconvolved, and, in the marine case the source and receiver ghosts have been removed as well. Notice that R(t ) in equation (7) is symbolically the same as in equation (4), because we have chosen the receiver level to be always at z = 0. This configuration is depicted in Fig. 1(c). It is well known that equation (7) can be written as p − (t ) = R 0 (t ) * p + (t ), (for a relation in 3D, see e.g., Amundsen, 2001). It is clear that we can deconvolve equations (6) and (7) by p + and obtain These equations represent the data for the fourth model and correspond to full data processing such that the resulting con- figuration is that of a reflection free acquisition surface and upper half space. This configuration is depicted in Fig. 1(d).
In all models, the wave speed and density of mass can be arbitrary piecewise continuous functions of depth below z = z 0 . We summarize the expressions for the down-and upgoing pressure wave fields for the three configurations in Table 1.

T H E F U N DA M E N TA L WAV E F I E L D
For convenience, we introduce a time parameter, ζ , that is related to propagation in the depth direction. The value of ζ indicates the two-way vertical travel time measured from the depth level z = 0 to any depth level below the receiver level. It is given by ζ = 2 z 0 c(ς ) −1 dς , which corresponds to the twoway travel time to the depth level z. In the derivations and examples that follow, we use t to denote time of an experiment, or recording time, and ζ to denote the two-way vertical travel time, which is a time that limits the existence of the fundamental wave field. Kunetz (1964) introduced the fundamental relations that remove the overburden multiples from the data. This relation makes recursive inversion possible on the acoustic reflection response of a discrete layered model. The fundamental relations of Kunetz are associated with the fundamental wave field. The importance of the fundamental wave field lies in the fact that when its down-going part is emitted into the subsurface, its up-going part will be the reflection response from which overburden effects are removed. More details on the fundamental wave field and the connection to the Marchenko focusing function can be found in the Appendix. Here we use the version that can be used in an arbitrary heterogeneous medium and denote the acoustic pressure of the fundamental wave field as δ(t ) + k(ζ , t ), with k(ζ , t ) = 0 for t ≤ 0 and t > ζ. It is a wave field with down-and up-going parts, δ(t ) + k + (ζ , t ) and k − (ζ , t ), respectively, and (10) The corresponding particle velocity of the fundamental wave field is then given by the down-going part, Illustration of the fundamental wave field in a truncated medium with three boundaries. The downward pointing arrows represent the down-going part, with the initial impulse in black and the coda in blue. The upward pointing arrows (red) represent the up-going part.
k + (ζ , t )], and up-going part, −k − (ζ , t )/Z 0 . We introduce the function h(ζ , t ) to represent the particle velocity, which is given by The down-going part of the fundamental wave field can be sent into the subsurface in which case its use can be seen as an acoustic reflection experiment with a complicated source time signature. The fundamental wave field is defined for a finite time duration for a fixed value of ζ , for example, ζ i . This corresponds to a medium that is reflection free for z < 0 and z > z i . We call this the truncated medium. In the Appendix, the acoustic reciprocity theorem of the time convolution type is used to derive expressions of the fundamental wave field in terms of impulse reflection response of the truncated medium, R i (t ). Here we give the equation as Equation (12) shows that k − (ζ , t ) is the reflection response of the truncated medium to the down-going part of the fundamental wave field. Here we give a simple illustration in a truncated medium with three reflectors in Fig. 2. The function k + (ζ i , t ) in equation (12) is the coda of down-going part, and k − (ζ i , t ) is the up-going part at the acquisition level. Figure 2 shows how the initial impulse (black line) and the coda (blue lines) propagate in the model and that when they bounce, they become part of the up-going field (red lines) which form k − (ζ i , t ) when they arrive at the acquisition plane. It can be seen that all downward travelling waves bounce only once and then become up-going waves that arrive at the acquisition plane. This means that all events in k − (ζ i , t ) are primary reflections. In this paper, the ones that can be traced back to the initial impulse are called physical primary reflections. The others are called non-physical primary reflections. In Fig. 2, four events can be seen in k − (ζ , t ) that arrive at the acquisition plane. Seen from the left, the first, third and fourth events are physical primary reflections. They have the amplitude of the local reflection coefficients r 0 , r 1 , r 2 , respectively. The second event is a non-physical primary reflection and has the amplitude r 0 r 1 r 2 , which is in terms of the order of multiplied reflection coefficients similar to the amplitude of a multiple. However, to call such event a multiple would be misleading, because the actual product of reflection coefficients does not correspond to a physical multiple and they have bounced only once. In this example, only one non-physical primary is present. In a discretely layered model with an arbitrary number of N + 1 reflectors, the number of physical primary reflections is N + 1, while the number of non-physical reflections is 2 N . All non-physical reflections have amplitudes that are a product of three of higher odd number of reflection coefficients. We refer to Slob et al. (2020) for a detailed discussion.

C O M P U T I N G T H E F U N DA M E N TA L WAV E F I E L D F RO M D I F F E R E N T DATA S E T S
In this section, we show that the fundamental wave field of a truncated medium does not depend on the data from which it can be computed. We either retrieve it as a two-way wave field as its pressure and particle velocity or in its up/down components of the pressure.

The marine setting
In the Appendix, the main equations are derived in detail that form the basis for all known schemes to compute the up-and down-going parts of the fundamental wave field. Here we repeat the general expressions of equations (A12) and (A13), which are only valid at the interval t ∈ (−t − , t − + ζ ). The reason for t − is that the source is above the receiver, which creates an extra delay of t − in the first arrival of the down-going wave. The equations are given by (Kang, 2018), Equations (13) and (14) can be solved for k(ζ , t ) and h(ζ , t ) using the pressure and particle velocity of the data. The righthand side of equation (13) is given in terms of the pressure and particle velocity, but their combination is equal to minus two times the up-going part of the particle velocity. When the only pre-processing step is up/down decomposition, these parts of the acoustic pressure as given in equa-tions (1) and (2) are known. We can substitute equations (10) and (11), which represent the up/down decomposition of the unknown part of the fundamental wave field, in equations (13) and (14), collect the terms for k ± and use equations (1) and (2) to arrive at These equations can be solved for k ± (ζ , t ) using the up-and down-going parts of the pressure of the data. Equations (15) and (16) are similar to the three-dimensional (3D) versions in Ravasi (2017), but he derived redatuming equations and used the up-and down-going parts of the particle velocity.
Equations (13) and (14) and (15) and (16) cannot be solved with a Neumann-type iterative scheme, because an initial estimate is not readily available for the unknown quantities. In the section on numerical aspects, we briefly discuss how a conjugate gradient method can be used.

Using the free surface impulse reflection response
If we deconvolve the marine data for the source time signature and its ghost and remove the receiver ghost by moving the receiver level to the free surface, the up-and down-going parts of the pressure data are given by equations (6) and (7). Substituting these results in equations (15) and (16) yields where now the time interval where these equations are valid is 0 < t < ζ, because the source has been moved to the receiver level. These equations can be solved for the unknown up-and down-going parts of the pressure of the fundamental wave field from the impulse reflection response including free surface effects. For the pressure-free surface, we take r = −1 and add the equations to find for t ∈ (0, ζ ) (Ware and Aki, 1969), with k(t ) defined in equation (10). This equation can hardly ever be solved with the Neumann scheme . In the section on numerical aspects, we show a conjugate gradient scheme that can be used.

Using the subsurface impulse reflection response
When we remove all surface-related multiples from R(t ), the up-and down-going parts of the pressure data are given by equations (8) and (9). Substituting these results in equations (15) and (16) yields where these equations are valid for 0 < t < ζ. In the section on numerical aspects, we show that these equations can be solved with the Neumann scheme.

M U LT I P L E E L I M I NAT I O N
Equation (A8) is repeated here, because from this equation we explain how a multiple free dataset can be obtained once the fundamental wave field is known, Equation (22) has a right-hand side that is zero for t < ζ i and non-zero for t > ζ i . The non-zero part can be seen as an unmodified part of the up-going part of the data as explained in Slob et al. (2020). The time instant t = ζ i is a time instant where equation (22) can show a jump discontinuity. This happens when the time instant ζ i is associated with a reflector at z i . In that case, for t ↑ ζ i the equation gives the reflection event with its reflection coefficient in k − , but for t ↓ ζ i the equation gives the same event with the strength of what we call the physical primary. A physical primary is a primary reflection as it is recorded in the data. Apart from the local reflection coefficient at the reflector where it originates from, it contains all overburden effects, such as free surface and ghost effects, and two-way transmission effects. How to capture a primary in either k − or in the right-hand side of equation (22) is explained in this section.

Transmission compensated Marchenko multiple elimination
It is well known that the up-going part of fundamental wave field contains all physical and non-physical primary reflections at the interval 0 < t < ζ (Robinson and Treitel, 1978).
In a piecewise homogeneous model, each reflector is characterized by a distinct reflection coefficient. All physical primary reflections in the data end up in the up-going part of the fundamental wave field with their reflection coefficient as strength. The up-going part of the fundamental wave field contains many more non-physical primaries. Those are reflections from the down-going coda of the fundamental wave field that is necessary to remove multiples generated in the time window of the fundamental wave field. For that reason, we cannot obtain all primary reflections from a single filtering step. Each time instant ζ must be used to compute the fundamental wave field for all t after which either a zero value or the value that corresponds to the reflection is found. In case the discrete truncation time, ζ i , corresponds to the two-way travel time to a reflector, k − (ζ i , ζ i ) = r i where r i denotes the local reflection coefficient of the reflector, otherwise a zero value is found. This forms the basis of the transmission-compensated Marchenko multiple elimination (T-MME) scheme . When k and h are computed, k − is obtained as . Independent of the data that is used to compute the result at each time instant ζ , we take the resulting value at this time instant and store it in a new primary reflections only dataset denoted R r (t ). Hence, This equation is independent of the type of data that is used to compute k − .

Marchenko multiple elimination
The Marchenko multiple elimination (MME) scheme filters the physical primary reflections from the data. Two different results can be obtained, but in all situations the source time signature is present in the result. The two situations in the marine setting, either no processing or only up/down decomposition, lead to the same result in equation (22). This is the primary reflections dataset in the setting of the actual acquisition. This means that the source ghost is present, and the time delay for the source being above the receiver together with the transmission effects. The situation where the acquisition surface is a free surface or a reflection free surface leads to the same result as well. This is the primary reflections dataset including transmission effects. The right-hand side of equation (22) is equal to an upgoing particle velocity scaled pressure given by Because we deal in most of the paper with the up-going part of the pressure, we capture the physical primary pressure reflection by rescaling the result of equation (24) and obtain The scale factor is Z 0 because this scaling occurs at the receiver level. This is similar to what was explained in Slob et al. (2014). We do not need to know the right-hand side of equation (24) to compute V − , we use it only to understand its meaning. Hence, we can write equation (22) for any value of ζ as Equation (27) is obtained from equation (26) by writing all fields under the integral in terms of their up-and down-going parts and collecting the results. We obtain the physical primary reflection by evaluating the right-hand side of equation (27). After source signature deconvolution, deghosting and moving source and free surface to the receiver level, the data can be described by the impulse reflection response R(t ) when the free surface is present and by R 0 (t ) after additional free surface multiple removal. The resulting up-going part of the data is an impulse reflection response similar to R(t ) or R 0 (t ) for t > ζ. Let us call these up-going parts of the pressure K − (ζ , t ) and K − 0 (ζ , t ), respectively. They can be written for ζ = ζ i as where G − and G − 0 are the Green's functions for waves arriving as up-going waves at z i generated by a monopole source in z = 0, which is the level of the free surface for G − and of the reflection free surface for G − 0 . We find the physical primary reflection from These two equations form the basis of  and Zhang and Staring (2018), respectively. Similar to the procedure for T-MME, also here we must evaluate the equations for each value of ζ separately and retrieve one sample from the result and store it in a new primary reflections only dataset. We call this dataset R t , and it includes all possible acquisition-related effects, depending on the data that were used to compute the results. Hence, R t is computed from equations (26) or (27) and the result is a primary reflections dataset with the acquisition geometry and medium the same as the data, or R t is computed from equations (30) or (31) and the result is a primary reflection dataset as if the source and receiver were placed at a depth level in a homogeneous upper half space.

N U M E R I CA L A S P E C T S
The integral equations of the kind we developed in the previous section can be solved in various ways. We look for a method that solves the problem with the least amount of numerical operations and memory use. For the kind of equations we solve here, these are iterative methods. Iterative schemes have the advantage that the system matrix is never computed but only its action on the unknown function. Equations (20) and (21) correspond to the situation when the subsurface impulse reflection response is available. In this situation, the equations can be solved with the well-known Neumann iterative scheme. A detailed exposition and examples of this scheme for the corresponding two-dimensional (2D) redatuming equations can be found in Thorbecke et al. (2017). For the situation when we use the pressure and particle velocity as measured, equations (13) and (14) must be solved. When these are decomposed in up-and down-going pressure fields, (15) and (16) must be solved. When deghosting, redatuming and source deconvolution have been carried out, (17) and (18) must be solved. In all these situations, the equations are solved with an unconditionally convergent iterative scheme. Here we discuss schemes that solve this symbolic equation in the least squares sense, in which u denotes the unknown fundamental wave field, f denotes the known field, which depends on the equations used, and L is the bounded linear operator that acts on u and contains only the known field. The associated error, r, is given by The operator has a bounded inverse that maps a Hilbert space H onto itself. The space is equipped with an inner product ·, · for real-valued functions, and a norm || · || given by u, v = u(t )v(t )dt, ||u|| = u, u 1/2 . The adjoint operator L † , associated with L, satisfies Lu, v = u, L † v . With these definitions, we give three iterative solutions of equation (32). One when the spectral radius σ = lim n→∞ ||(I − L) n || 1/n < 1, where I is the unit operator, and two when σ < 1. The schemes that follow minimize ||r||. Note that for a self-adjoint operator L, the well-known scheme of Hestenes and Stiefel (1952) minimizes u, 1 2 Lu − f rather than ||r||. In the situation where only internal multiples need to be removed, equations (20) and (21) need to be solved. In this case σ < 1, because all subsurface reflection amplitudes are less than one. This means the scheme can be taken as a Neumann scheme given by where the subscript n denotes the nth iterate and r 0 = f − Lu 0 , with u 0 as the initial estimate that we specify in the next section. This scheme is not only the simplest but has the best convergence rate of the schemes presented here. When the free surface reflections are present in the data, the spectral radius is easily larger than one. In those cases, the Neumann scheme does not converge (see, e.g., Dukalski and de Vos, 2018). In those situations, we distinguish between the case where we have L = L † and where L = L † . Either scheme is initialized as u 0 is arbitrary, r 0 = f − Lu 0 and w 0 = 0, where w 0 is the initial term of an intermediate result in the conjugate gradient (CG) scheme from which subsequent iterates w n are computed. Subsequent iterations are computed for n ≥ 1 as w n = w n−1 + T r n−1 T r n−1 , L † r n−1 , where T = I when L = L † and T = L † when L = L † . When T = L † , this scheme is equivalent to the scheme of Le Foll (1971). The error criterion to stop the iterations numerically is the normalized global root mean square error, given by err n = ||r n ||/|| f ||.
One Neumann iteration requires one convolution and one correlation. These can be computed efficiently using an fast Fourier transformation (FFT) routine. A correlation is the adjoint of convolution, and if we compare the Neumann scheme of equation (34) with the scheme for a self-adjoint operator of equation (35) we see that the latter scheme requires  two convolutions and two correlations to perform one iteration. The reason is that we force convergence through the normalizations. It will require more computation time and more iterations to achieve the same error value. When T = I, we still have fast convergence, but when L = L † we must use T = L † to have a norm in the update of w n , which ensures the error is decreasing at every iteration. This reduces the convergence rate compared with the scheme for a self-adjoint operator, while it costs the same amount of operations per iteration. A detailed derivation and discussion of the associated convergence properties of these two schemes can be found in van den Berg (1991). Because convolution and correlation are each other's adjoints, equation (19) can be solved with T = I in equation (35). The corresponding equation in 3D may not contain a self-adjoint operator because in general R(x 0 , x 0 , t ) = R(x 0 , x 0 , t ). Once k is known, we compute k − from equation (17) and K − from equation (30).
When the reflection response is not explicitly present in the data, such as in equations (13) and (14) and (15) and (16), we solve the equations with T = L † in equation (35). Once k ± or k and h are known, we also know k − and compute P − from equation (27) or (26).
We summarize the equations to solve for the type of input data and how and from which equation the primary reflection dataset is obtained in Table 2. The order in the table from top to bottom is with increasing complexity in the data, which is the order in which we show numerical examples in the next section.

N U M E R I CA L E X A M P L E S
We use a layered model without and with free surface to compute all responses that are used to filter out the primary reflections. The model parameters are given in Table 3. The first thickness is relative to the receiver level at z = 0. In the marine setting, we use z f s = −31.5 m and z s = −21 m. A 30-Hz Ricker wavelet, denoted W (t ), is used with time step of dt = 1 ms to model the data and to show the responses.
The data for this model for the various input datasets indicated in Table 1 are shown in Fig. 3. Figure 3(a) shows the subsurface reflection response R 0 (t ) * W (t ), and Fig. 3(b) the reflection response R(t ) * W (t ) at the free surface. Both plots have the same amplitude scale. All primary reflections occur within the 2.5-s time window shown. It can be understood that filtering out the primary reflections requires the elimination of more multiples from R(t ) than from R 0 (t ). Figures 3(c) and 3(d) show the down-and up-going parts of the pressure for the marine-type acquisition of equations (3) and (4), while 3(e) and 3(f) show the pressure and particle velocity, scaled with impedance. The amplitude scales are adapted to the fields but maximized at ±1. Comparing these figures, we see that from 3(a) to 3(f) more events are present in the data and less preprocessing has been done on the data. We use these fields and the corresponding expressions to filter the primary reflections from the data. The input data that is necessary as a separate term is always as displayed in Fig. 3, while the data that is necessary as an operator in equations (20) and (21) and (19) is deconvolved for the Ricker wavelet. This is done in the frequency domain by division after adding 10 −3 max(Ŵ ) to the Scaled velocity (f) (e) Figure 3 The wave fields at the receiver level in the model given in Table 3 as a function of time: (a) subsurface reflection response, (b) reflection response with free surface multiples, (c) up-and (d) down-going parts of the pressure, (e) pressure and (f) particle velocity scaled by impedance.
wavelet. We need to do this because these equations are derived using an impulse as the initial down-going part of the pressure, while in reality the pressure will have a finite bandwidth which is created by using the Ricker wavelet. Deconvolving the reflection response as operator ensures that the Ricker wavelet is automatically included in the fundamental wave field. In this way, we solve the bandlimited version of equations (20) and (21) and (19). Equations (15)and (16) and (13) and (14) are not derived with any assumption about the source, and we can use the data as they appear in the equations, both as input and as operator. In those cases, we convolve the fundamental wave field with the Ricker wavelet for display purposes.

Removal of internal multiples
In the situation where the acquisition surface is in the homogeneous upper half space, we solve equations (20) and (21) with the iterative scheme of equation (34). The resulting primary reflections only datasets, R t and R r , are shown in Fig. 4. The top plot shows the primary reflection events convolved with the wavelet as the expected outcome in black and the corresponding Marchenko multiple elimination (MME) result in red. The middle plot shows the reflectivity convolved with the wavelet as the expected outcome in black and the corresponding transmission-compensated Marchenko multiple elimination (T-MME) result in red. The bottom plot shows the difference between the black and red lines in the top and middle plots.
As discussed at the beginning of the multiple elimination section, the time instant t = ζ is an instant of jumpdiscontinuity in the functions when the data have infinite bandwidth. Using finite bandwidth data requires making a choice. Each choice leads to a different dataset with only primary reflections. Instead of taking a single limit, we must now avoid overlapping time windows. This can lead to resolution problems when layers become too thin as discussed in Slob et al. (2014), and a way to reduce this problem can be found The modelled (black) and retrieved (red) primary reflections as a function of time, with reflection amplitudes including two-way transmission effects, R t (top plot), and reflection coefficients as amplitudes, R r (middle plot); the difference between modelled and retrieved results with T-MME in black and MME in blue (bottom plot).
in Elison et al. (2020). To demonstrate how the result for R t is obtained, we show the results for K − 0 and k − as a function of one-way vertical travel time, ζ /2, and time, t, in Figs. 5(a) and 5(b). In each plot, the horizontal black lines indicate the depth of the reflectors in one-way vertical travel time. They are presented in the figure for illustration purposes only. The slanted black line indicates the truncation line t = ζ for the integrals, and we take the values along this line and put them in R t (t ). To solve equations (20) and (21) for k ± at the interval t ∈ ( , ζ ± ), is a value determined by the size of the wavelet. The 30-Hz Ricker wavelet in our examples has a time duration of 60 ms down to 1% in amplitude, and we take = 0.030 s. The minus-sign is used in the time interval for MME, and the plus-sign for T-MME. In this way, we obtain the results shown in Fig. 5. Figure 5(a) illustrates MME, and we can see that the function K − 0 is still non-zero just below the diagonal black line that marks the truncation time. From this solution, we obtain R t as indicated in the top row and last column of Table 2. The result is shown in the top plot of Fig. 4. We show the result for k − (ζ , t ) as a function of ζ /2 and t in Fig. 5(b) to illustrate T-MME. We can see that the function k − is already non-zero just above the black line that marks the truncation time. From this result, we obtain R r as indicated in the top row and third column of Table 2 and shown in the middle plot of Fig. 4.
Finally, it is good to remark that the functions k ± (ζ , t ) do not change when we increase the value of ζ by ζ when the depth associated with the vertical travel time ζ + ζ is inside the same layer as ζ . In that case, no new primary reflections occur in the data and, hence, no new events can occur in the fundamental wave field (Zhang and Slob, 2020b). The existing events are unchanged because we always start with an initial impulse at zero time in k + . This can be seen in Fig. 5(b) for k − . When ζ is small enough, the evaluation of equation (20) using k + (ζ , t ) with padded zeros, for ζ < t < ζ + ζ , as the estimate for k + (ζ + ζ , t ) will result in k − (ζ , t ) with the same padded zeros and the scheme will terminate. When ζ is large enough to move across a reflector, all previous results are  , t ) with the aim to obtain R t from K − 0 (ζ , t = ζ ); (b) the fundamental wave field k − (ζ , t ) with the aim to obtain R r from k − (ζ , t = ζ ); in both plots the wave fields are shown as a function of one-way vertical travel time, ζ /2, and time t, the line t = ζ is indicated by the diagonal solid black line, and the horizontal black lines indicate one-way vertical travel time to the reflectors. unchanged and new events are found in k ± (ζ + ζ , t ). Once an event is created, it remains unchanged for all larger values of ζ .
We illustrate the effect of using previous results as initial estimate for a new value of ζ with the same model example with 11 reflectors as used for Figs. 4 and 5. Figure 6 shows the number of iterations required to obtain the same quality in the result. The dashed line corresponds to the result shown for obtaining R r in Fig. 4, which is obtained with a zero initial estimate for the up-and down-going parts of the fundamental wave field for every new value of ζ and a stop criterion of err = 10 −3 . The number of iterations increases rapidly with increasing values of ζ . Using the previous result for k ± as initial estimate for a new value of ζ does not require any additional iterations unless the next value of ζ is connected to a new reflector within the time duration of the source time signature. For most reflectors, the scheme has converged in less than three iterations as can be seen in the solid line in Fig. 6. The fundamental wave field is computed for 1251 different  values of ζ with dζ = 2 ms. When a zero value is used as initial estimate, the total number of iterations is 41,768, whereas the total number of iterations is 283 when the previous result is used as initial estimate. A stop criterion of err = 10 −2 is sufficient to achieve the same overall quality in reflection amplitude retrieval and multiple suppression when the previous result is used as initial estimate in this example. Dukalski and de Vos (2018) demonstrated that the scheme of equation (34) applied to similar redatuming equations as equations (17) and (18) does not converge in most situations. We solve equation (19) with equation (35) with T = I. It is noted that the discrete evaluation of the correlation and convolution in equation (19) can be done simultaneously using their discrete Fourier transformed (DFT) counterparts. The discrete convolution and correlation after DFT together amount to a single product of the fundamental wave field k and the real part of the impulse reflection response. The time domain result is then obtained after inverse DFT and time truncations as indicated by the integral limits. The DFTs can be computed efficiently using a fast Fourier transform algorithm. In all configurations and all other schemes, the time integrals are computed in this way, but then separately for convolution and correlation. The iterative scheme of equation (34) works down to the fourth reflector and starts diverging right after that. For the conjugate gradient scheme, we use the same criteria as before. The results for K − and k − up to 5 s are shown in detail in Figs. 7(a) and 7(b). The function K − looks like the one shown in Fig. 5 but has the free surface multiples, whereas the fundamental wave field k − is the same because it does not depend on the presence or absence of the free surface in the data. Figure 7(a) shows many more events than Fig. 5(a). For vertical travel times below the bottom reflector, the function K − should be zero and we can see in the figure that it is virtually zero. The colour scale gives saturated red for positive and blue for negative values that are 10% of the maximum value in the data. This scale maximum is five times larger than the amplitude of the reflection from the last reflector. Visible white occurs for values that are 0.65% of the colour maximum. The zone where the Green's function should be zero is visible white. We conclude that the method performs well and is robust to recording times well beyond the arrival time of the last reflection. The reason to show this result for k − to 5 s is because Dukalski and de Vos (2018) argue that removing free surface and internal multiples simultaneously is not necessarily a good idea in Marchenko redatuming. The fact that we do not redatum does not change the numerical challenge.

Joint removal of free surface and internal multiples
Here we show that the conjugate gradient scheme performs very well and is quite stable to times that relate to large depths. Several reflection coefficients we use are quite high as can be seen from R r in Fig. 4. These results correspond to the second row and the third column for R r and the last column for R t in Table 2. They are both exactly the same as those shown in Fig. 4, but now obtained from the data in Fig. 3(b). Figure 8 shows the number of iterations required to solve the problem with err = 10 −3 with u 0 = 0 and with err = 10 −2 when using the previous results for every new value of ζ /2. Considering the whole 5 s window (2501 ζ -values), the total number of iterations with zero values as initial estimate is 148,912 and it drops to 1165 when we use the previous results as initial estimate. For the 2.5-s time window of the previous example, the numbers are 46,134 and 959, respectively. Compared with the required number of iterations in the situation without free surface shown in Fig. 6, the increase is less than 10% using zero initial estimate and it is a factor 3.4 with previous result as initial estimate. This shows that the method is stable, robust and accurate. Using the previous result as initial estimate keeps the number of iterations well behaved. It is clear that using previous results for every new value of ζ is beneficial from a computational point of view. This benefit cannot be exploited easily in the known redatuming schemes, because there the focus is always at zero time. This implies that previously obtained results should shift in time when a new focus point is chosen, and the quality of the estimate of the shift depends on the quality of the background model.

Joint removal of ghost, free-surface and internal multiples
In a marine setting, the presence of the source and receiver ghosts need to be taken into account. When source and receiver are not at the same depth, the schemes that assume an initial down-going impulse at the receiver level are not feasible. We use the data either as recorded or after up/down decomposition. This means that we do not make assumptions about the reflection strength of the free surface, and the source time function can remain unknown because we keep it in the data. We solve equations (15) and (16) when the upand down-going parts of the pressure, shown in Figs. 3(c) and 3(d), are available. This leads to the desired dataset with only primary reflections in R r according to the third row and column in Table 2. We find the physical primaries with the source wavelet and the ghost present in P − according to the third row and last column in Table 2. When up/down decomposition is not performed but pressure and particle velocity are available, we use equations (13) and (14) with the data shown in Figs. 3(e) and 3(f). From this result, we compute R r and we compute P − as indicated in the fourth row, third and fourth columns, respectively. Using an error criterion of err = 10 −3 and the same time steps as before, the CG scheme requires a total of 22,900 iterations to solve equations (15) and (16) and 20,950 to solve equations (13) and (14) for 1251 ζ -values in the 2.5 s time window when the previous filter result is used for the next truncation time value. The distribution of the iterations for solving the latter two equations is shown in Fig. 9. It can be seen that up to almost 300 iterations are required at the stronger reflectors. The reflectivity is quite well retrieved, although the amplitude of the last reflector is 16% too small.
The global error criterion is too relaxed at 10 −3 , and the smallest reflection amplitude is recovered to less than 1% error at its maximum by using 10 −4 as stopping criterion, at the cost of a total of 109,000 iterations. In that case, the largest number of iterations for a single ζ -value is 1160, which occurs The iteration counts as a function of ζ /2 to reach the error of 10 −3 with pressure and particle velocity data from the marine setting to retrieve the primary reflections.
at the eighth reflector. When p ± are used instead of p and v z the total number of iterations is 120,000 and the largest number of iterations for a single ζ -value is 1300. Figure 10 shows the primary reflections in the up-going pressure, R t , (top plot) and in the up-going part of the fundamental wave field, R r (middle plot). These results are obtained with 10 −4 as stopping criterion. The most striking difference in these result is the ghost effect that is present in R t and absent in R r . This is because in R t we retrieve the primary reflection as they are in the data, hence in the physical domain with the free surface, and the source and receiver at the levels used in acquisition. The primary reflections in R r are obtained from the fundamental wave field which has source and receiver at the reflection free boundary. For the same reason, a slight time shift can be seen in the events in the top and bottom traces. The bottom plot of Fig. 10 shows the differences in the expected and retrieved primary reflections with the MME and T-MME methods in blue and black lines, respectively. It can be seen that the MME result has its most visible errors at the tail of the physical primary reflections, whereas the T-MME result has its largest errors just before the arrival times of the primary reflections. The latter are visible in the middle plot of Fig. 10, especially just before the third to sixth events.

S U M M A RY A N D D I S C U S S I O N
The fundamental wave fields are computed in a simple filter procedure that uses the input data and a time truncation that is a free parameter. For each new truncation time instant, we can use the fundamental wave fields obtained from the previous truncation time instant as a first estimate. This reduces the required number of iterations by at least an order magnitude. The truncation time step does not have to be the same as the sampling time, but can be adapted such that it satisfies the sampling criterion, which reduces the computation time further. If for display purposes a smaller time step is desirable to present the final results, proper zero padding in the frequency domain will produce exact sinc-interpolation results. Four schemes are presented for which the input data has received various levels of pre-processing before being used to filter out the primary reflections. All schemes filter the subsurface reflection response, but in forms that depend on the pre-processing. They have been tested in a model with 11 reflectors, with local reflection coefficients given in Table 3. All primary reflection events are retrieved with errors less than 1% in their maximum values. These results show that the required number of iterations to filter out the primary reflections with reflectivity as amplitudes depends strongly on the amount of events in the data that is used as input. This should be no surprise because the more events are present in addition to the primary reflections the more multiples must be eliminated to compute the fundamental wave field.
The first example used the subsurface reflection response, in which only internal multiples are present in the data to retrieve the primary reflections. All primary reflection events are obtained with less than four iterations per time instant per reflector. In the second example, free-surface multiples are present as well in the reflection response. Here we assumed a known free-surface reflection coefficient r = −1. The presence of the free surface has two consequences. First, the number of events in the input data increases, which makes it harder to filter the primary reflections from the data. Second, the resulting equation cannot be solved with the Neumann iterative scheme. The proposed conjugate gradient scheme seems a better alternative. It solves the same equation and minimizes the number of iterations while ensuring convergence. In this example, we split the down-going wave field in an initial impulse and the free-surface effect. This leads to an operator equation with a self-adjoint operator.
The third and fourth examples come from a computed marine acquisition dataset. Considering the number of iterations necessary to find the primary reflections, it does not really matter whether up/down decomposition is carried out. In either case, the wavelet is not deconvolved for, the ghosts and all multiples are present in the data. This leads to an equation with an operator that is not self-adjoint. The number of numerical operations per iteration compared with the previous The modelled (black) and retrieved (red) primary reflections as a function of time in the marine setting, with reflection amplitudes including twoway transmission effects, R t (top plot), and reflection coefficients as amplitudes, R r (middle plot); the difference between modelled and retrieved results with T-MME in black and MME in blue (bottom plot) with an error criterion of err = 10 −4 . example does not increase, but the convergence rate is less. In addition, the number of events and their amplitudes in the input data have dramatically increased due to the source receiver distance and the fact that both source and receiver are below the free surface. This alone makes the problem already much harder to solve, because the global root mean square error criterion should be set one to two orders of magnitude smaller to accurately retrieve also the weakest reflection event.
In combination with the decreased convergence rate, we see that solving this problem requires almost two orders of magnitude more iterations. In this example, no assumption is made about the properties of the free surface, the primary reflections in R t have retained their ghost and the arrival time is the same as in the data, while in R r the ghost is removed and the source is redatumed to the receiver. Not having to assume anything about the free-surface reflection properties is an advantage, as is the absence of needing to know the source time signature. Equations (13)-(15) are possibly the best equations to remove free-surface and internal multiple reflections from the data together with the ghost. These equations assume the pressure and particle velocity as measured to be mutually well balanced in amplitude, the acoustic impedance to be known at the receiver level and the medium to be lossless.
All the results presented here in 1D can be extended to 3D with various levels of assumptions that are similar to those known from the concept of Marchenko redatuming. The general 3D situation is beyond the scope of this paper. Recently, an interesting combination of 1D multiple elimination and 3D data has been developed by Meles et al. (2020). The 2D version of this scheme takes line data and sums over all sources to create a plane wave input reflection dataset. The fundamental wave field is also integrated over all source locations, and a plane wave equation results where the reflection response as operator is still the original dataset with all sources and receivers. When the plane wave is a normal incidence plane wave and the earth would be horizontally layered, this scheme would be identical to the scheme presented in this paper. However, the scheme of Meles et al. (2020) can work with plane waves with non-zero angles of incidence and be applied to general data. In the latter case, the computational cost is, however, very similar to that of the 1D scheme. The method was shown to work well for a variety of angles of incidence of the plane waves. Numerical tests on 2D models showed that the method is fully automated, able to remove multiples and that a combination of 11 different angles could be used to make a 2D image with very good results. This has been tested only on reflection data with the acquisition surface in a homogeneous upper half space.
Whether the results presented in this paper should lead to a change in the workflow remains to be seen. Existing free-surface multiple removal techniques are well established. Solving the free-surface multiple problem before the internal multiple problem can be more advantageous from a computational point of view. Problems with estimated source amplitudes and remnant free-surface multiples can be mitigated using a small adaptive filter in the method that filters the subsurface reflection response R 0 . An important advantage of removing all multiples and ghost simultaneously is possibly the fact that it can be done with minimum assumptions on data and model. Because the multiple elimination methods presented here are, in principle, non-recursive, multiples can be eliminated in a target-oriented manner, starting at any desired truncation time, which is another advantage. This will come at the cost of many iterations that will be necessary for the initial time instant, after which the previous results can be used again as initial estimates for consecutive time instants. Target-oriented demultiple on time slices could contribute to a target-oriented workflow including pre-processing and imaging (e.g. Karrenbach, 1990). How these methods perform with noisy data is an important topic for the overall performance assessment, but is outside the scope of the present paper. Advantages of removing all multiples simultaneously have to be weighed with the disadvantage of having to run many iterations, which may be problematic on noisy data. This remains to be investigated.

C O N C L U S I O N S
We have shown that the fundamental wave field is related only to the subsurface impulse reflection response. This relation can be written in various forms depending on the processing performed on the measured data before multiple elimination. The main equations are derived from the reciprocity theorems using the pressure and particle velocity fields in a marine setting. From these equations, all known equations follow directly by applying different levels of pre-processing. These are up/down decomposition; up/down decomposition, deghosting, source wavelet deconvolution, and source and free-surface redatuming; and, all the before with additional free-surface multiple reflection removal.
The new scheme is derived for a source above the acquisition surface but below the free surface. The reflection properties of the free surface and the source wavelet remain unknown. The fundamental wave fields are retrieved as impulse responses within the bandwidth of the data. Up-down decomposition of the pressure and the vertical component of particle velocity data leads to a similar and known scheme. Other known scheme that follow from this scheme are when the free surface coincides with the acquisition surface and the scheme where the acquisition surface is in a homogeneous half space. In the latter two schemes, it is assumed that the source wavelet is known and removed from the data.
For each of the four different datasets, a method of solution is presented for the retrieval of primary reflections from input data. Less pre-processing means more work to obtain the fundamental wave fields and the desired primary reflections. Retrieving the reflectivity is beneficial, because the transmission effects are removed and true reflection amplitude events are obtained that in 3D can be used for amplitude versus offset (AVO) or other types of inversion. It remains to be seen whether retrieving primary reflections directly from the measured data is the most advantageous way to obtain a dataset that can be used for velocity analysis, imaging and inversion.

DATA S H A R I N G
The data that is used to make any of the plots in this paper can be obtained from the first author on request.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data to reproduce the plots in this paper can be obtained from the first author on request.

A P P E N D I X : S I N G L E -S I D E D F U N DA M E N TA L WAV E F I E L D R E P R E S E N TAT I O N S R e c i p r o c i t y t h e o r e m s
The starting point for the single-sided fundamental representations lies in the reciprocity theorems of the time-convolution and time-correlation types (de Hoop, 1995). Here we use the theorems involving two states of wave fields in one and the same medium and with sources outside the region where reciprocity is applied. Then in a 1D model they are given by in which * denotes temporal convolution and the expression should be understood to be evaluated at both depth levels in the limits of z ↑ 0 and z ↓ z i . The book can be found online for personal use: http://www.atdehoop.com/Books/ DeHoopAT_book_2008/DeHoopAT1995Handbook.index. html. Equation (A1) is the 1D version of equation (7.2-7) in the book and equation (A2) of equation (7.3-7). To derive the desired representations, we need to define the medium in states A and B only in the region 0 < z < z i . If the medium outside this region is reflection free, we call this the truncated medium. The truncated medium and the medium that represents an experiment are the same medium in the region 0 < z < z i but can be different outside this region. For the choices we make a concise derivation of the reciprocity relations, we use can be found in Eisner (1981).

D e f i n i t i o n o f t h e f u n d a m e n t a l wave f i e l d
To define the fundamental wave field, we use the truncated medium in both states. The truncated medium is characterized by its impulse reflection R i (t ) and transmission T i (t ) responses, where the subscript i indicates the truncated medium. At z = 0, we define the pressure in state A as p A (0, t ) = δ(t ) + R i (t ) where the impulse is the initiating down-going part p + and the reflection response is the up-going part p − . The particle velocity is then given by v z;A (0, t ) = Z −1 0 [δ(t ) − R i (t )], where Z 0 is the acoustic impedance of the medium at z = 0. In state B, we define the pressure to be the fundamental wave field given by p B (0, t ) = δ(t ) + k + (ζ i , t ) + k − (ζ i , t ), where the impulse and k + together form the initiating down-going part p + and k − is the up-going part p − . The corresponding velocity field is then given by v z;B (0, t ) = Z −1 0 [δ(t ) + k + (ζ i , t ) − k − (ζ i , t )]. The depth level z i is indicated by the two-way vertical travel time ζ i in the argument of the fundamental wave field. It is given by ζ i = 2 z i 0 1/c(z)dz. At depth level z i , the fields in state A are the impulse transmission responses in the truncated medium given by p A (z i , t ) = T i (t ) and v z;A (z i , t ) = T i (t )/Z i . The fundamental wave field in state B is defined such that at z i only the direct part of the impulse transmission response arrives, T d;i , such that p B (z i , t ) = T d;i (t ) and v z;B (z i , t ) = T d;i (t )/Z i . Substitution of these fields in equation (A1) results in When the down-going part of the fundamental wave field is emitted into the truncated medium, the transmission response is the direct part of the impulse transmission response. This can be written as T d;i (t ) = T i (t ) + T i (t ) * k + (ζ i , t ). Hence the down-going part of the fundamental wave field is the multiple eliminator for the truncated medium. The inverse of a multiple eliminator is a multiple generator. The impulse transmission response can be written as T i (t ) = T d,i (t ) * M i (t ), where M :i (t ) is the multiple generator . This leads to the following expression: which states that the down-going part of the fundamental wave field is the multiple reflection eliminator for the truncated medium. It consists of an impulse and a following coda k + (ζ i , t ) such that k + (ζ i , t ) = 0 for t ≤ 0 and t > ζ i . We convolve both sides of equation (A3) with M i (t ) and use the result of equation (A4) to find which shows that the impulse reflection response of the truncated medium is equal to the up-going part of the fundamental wave field convolved with the multiple generator. Therefore, the up-going part of the fundamental wave field is zero for t ≤ 0 and for t > ζ i . To summarize, we have found that the fundamental wave field satisfies the wave equation. At the acquisition surface, it has down and up-going parts that exist in the same time window. This window is bounded by zero time and the two-way vertical travel time to a depth level z i , below which the truncated medium is reflection free. The down-going part is always an impulse at zero time followed by a coda that together form the multiple reflection eliminator of the truncated medium. The up-going part is the reflection response to the down-going To find the Marchenko multiple elimination (MME) method, let us look again at equation (A8) and assume we computed the fundamental wave fields from equations (A12) and (A13). In that case, the entire left-hand side of equation (A8) is known and we rewrite the equation as where we have used the fact that p = p + + p − and v z = (p + − p − )/Z. The function p − on the left-hand side of equation (A15) can be interpreted as the up-going pressure just below the depth level z i that is generated by the source used in the experiment. The function T d;i /Z i can be interpreted using reciprocity as the direct wave transmission from the depth level z i to the receiver level z = 0. Their convolution represents an unmodified part of the up-going field as measured in the data. It contains all the effects of heterogeneities above the receiver level that are in the down-going field, but all multiples of the overburden are eliminated such that the first arrival is the physical primary reflection of a possible reflector whose two-way travel time to the receiver is equal to ζ + t − that is free from multiple reflections. Without further specification of the medium, we can only write it as This expression is the general 1D form of the MME method.