Three‐dimensional Marchenko internal multiple attenuation on narrow azimuth streamer data of the Santos Basin, Brazil

ABSTRACT In recent years, a variety of Marchenko methods for the attenuation of internal multiples has been developed. These methods have been extensively tested on two‐dimensional synthetic data and applied to two‐dimensional field data, but only little is known about their behaviour on three‐dimensional synthetic data and three‐dimensional field data. Particularly, it is not known whether Marchenko methods are sufficiently robust for sparse acquisition geometries that are found in practice. Therefore, we start by performing a series of synthetic tests to identify the key acquisition parameters and limitations that affect the result of three‐dimensional Marchenko internal multiple prediction and subtraction using an adaptive double‐focusing method. Based on these tests, we define an interpolation strategy and use it for the field data application. Starting from a wide azimuth dense grid of sources and receivers, a series of decimation tests are performed until a narrow azimuth streamer geometry remains. We evaluate the effect of the removal of sail lines, near offsets, far offsets and outer cables on the result of the adaptive double‐focusing method. These tests show that our method is most sensitive to the limited aperture in the crossline direction and the sail line spacing when applying it to synthetic narrow azimuth streamer data. The sail line spacing can be interpolated, but the aperture in the crossline direction is a limitation of the acquisition. Next, we apply the adaptive Marchenko double‐focusing method to the narrow azimuth streamer field data from the Santos Basin, Brazil. Internal multiples are predicted and adaptively subtracted, thereby improving the geological interpretation of the target area. These results imply that our adaptive double‐focusing method is sufficiently robust for the application to three‐dimensional field data, although the key acquisition parameters and limitations will naturally differ in other geological settings and for other types of acquisition.


I N T RO D U C T I O N
The Santos Basin in Brazil is known for its oil-bearing carbonate reservoirs below a highly reflective stratified salt layer (see Fig. 1). The salt layer generates strong internal multiples that [Corrections made on 18 June 2020, after first online publication: The horizontal axis labels have been added to  * E-mail: myrna.staring@gmail.com pose a problem for seismic imaging (Cypriano et al., 2015). Most imaging methods assume that the recorded wavefield was only reflected once and thus incorrectly interpret internal multiples as primaries from deeper reflectors. As a result, these methods create ghost reflectors that do not exist in reality. These ghost reflectors can interfere with the real reflectors in the target area and thereby corrupt the image. Therefore, we wish to attenuate internal multiples in order to obtain a reliable image of the target area. The attenuation of internal multiples is a challenge. Various methods have been proposed, ranging from filtering methods (e.g. Hampson, 1986;Foster and Mosher, 1992;Zhou and Greenhalgh, 1994) that transform the reflection response to an alternative domain in which the primaries and the internal multiples separate, to wave-equation-based methods that aim to predict the internal multiples by convolving and correlating the reflection response with itself (e.g. Jakubowicz, 1998;Weglein et al., 1997). The application of filtering methods is often challenging in settings with a complex overburden, since there is usually no distinct difference in properties between the primary reflections from the target area and the (strong) internal multiples generated in an overburden. Also, the application of wave-equation-based methods is not undisputed. Some wave-equation-based methods require the manual identification of internal multiple generators, thereby introducing bias and the risk of not correctly capturing all internal multiple generators in the process. In addition, some wave-equation-based methods predict internal multiples with incorrect amplitudes or use a layer stripping approach that results in error accumulation from shallow to deep reflectors. In order to attenuate internal multiples in a complex setting such as the Santos Basin, an alternative method is needed.
Marchenko methods (Ware and Aki, 1969;Broggini et al., 2012;Wapenaar et al., 2014) are data-driven and wave-equation-based methods that do not have these drawbacks. These methods have the ability to consider the entire overburden as a whole, instead of having to identify all individual internal multiple generators separately. In addition, Marchenko methods allow us to retrieve Green's functions including primaries as well as all orders of internal multiples at any desired depth level without having to resolve overlying layers first. When writing the retrieval of Green's functions using the coupled Marchenko equations as a Neumann series, Marchenko methods can be used for the prediction of internal multiples (van der Neut et al., 2015b). These predictions in principle have the correct amplitude and phase. However, minor amplitude and phase differences are usually present when applying the method to field data due to imperfect acquisition and preprocessing. A mild adaptive filter can be used to correct for these minor differences. We previously reported on the successful application of an adaptive Marchenko method (the adaptive double-focusing method) to two-dimensional (2D) synthetic data and a 2D line of streamer data of the Santos Basin in Brazil (Staring et al., 2018a). Internal multiples were predicted and adaptively subtracted from the target area, which improved the geological interpretation. In addition, we found that the adaptive double-focusing method was relatively robust for a sparse acquisition geometry in 2D and suitable for the application to large data volumes. In the hope that these properties also hold in three dimensions, we use this adaptive Marchenko method for the prediction and adaptive subtraction of internal multiples from three-dimensional (3D) narrow azimuth streamer data acquired in the Santos Basin.
The extension from 2D to 3D Marchenko methods may seem trivial in theory, but it is not the case in practice. Some aspects are similar, such as the data preparation requirements that include noise suppression, signature deconvolution, deghosting and the removal of surface-related multiples. However, aspects related to the sampling of the acquired data are different. In addition to the inline direction in 2D, there is in 3D also a crossline direction that typically has a limited aperture and less densely spaced sources and receivers. Also, streamers usually do not record responses at negative offsets, near offsets and far offsets in the inline direction. A thorough understanding of the effect of these acquisition limitations on the result of Marchenko internal multiple attenuation would allow us to estimate whether the application to any particular dataset is feasible. In addition, it would aid us in defining an interpolation strategy. Even though some researchers already applied a Marchenko method to 3D field data (Staring et al., 2018b;Pereira et al., 2018), they did not address the acquisition requirements and limitations of 3D Marchenko methods Cartoons showing (a) the starting acquisition geometry and (b) the final acquisition geometry for the synthetic decimation tests in this paper. The final acquisition geometry is based on our narrow azimuth streamer data. The stars represent sources, and the triangles represent receivers.
in detail. The objective of this paper is to gain a better understanding of the key acquisition parameters and limitations that affect the application of the adaptive double-focusing method to 3D data.
In this paper, we first revise the theory of the adaptive Marchenko double-focusing method. Second, we perform a series of 3D synthetic tests to study the effect of the acquisition parameters on the result of internal multiple prediction and adaptive subtraction using this method. Starting from a grid spacing of 50 m (inline direction) by 75 m (crossline direction) co-located sources and receivers with positive and negative offsets, near offsets, far offsets and a crossline aperture of 1.8 km (Fig. 2a), we step-by-step decimate the acquisition down to a narrow azimuth streamer geometry on which our 3D field data were acquired (Fig. 2b). Based on these tests, we identify the key limiting acquisition parameters and use these to design an interpolation strategy for the field data application. Next, we test the proposed interpolation strategy on 3D synthetic data. Finally, we apply the adaptive double-focusing method to 3D narrow azimuth streamer data. In the following discussion and conclusion section, we evaluate the performance of the adaptive double-focusing method.

M A R C H E N KO I N T E R NA L M U LT I P L E AT T E N UAT I O N B Y A DA P T I V E D O U B L E -F O C U S I N G
The adaptive Marchenko double-focusing method requires a preprocessed reflection response R(x R , x S , t ) acquired on a sufficiently dense grid of sources x S and receivers x R at the acquisition surface ∂D 0 . A smooth velocity model of the subsurface is needed to obtain the direct wave of the downgoing fo-cusing functionf + 0 . The direct wave is obtained by modelling and time-reversing the response from sources at the redatuming level ∂D i to receivers at the acquisition surface ∂D 0 using finite-difference modelling or an Eikonal solver (see Fig. 3). The· symbol indicates an user-specified wavelet that is convolved with the modelled wavefield. The direct downgoing focusing functionf + 0 initiates the iterative scheme that solves the coupled Marchenko equations. If the overburden were homogeneous, this initial wavefield would be sufficient to create a focus at the desired focal point at ∂D i . Otherwise, a coda for Figure 3 Cartoon illustrating the direct wavef + 0 , which we obtain using a smooth velocity model and an Eikonal solver. the downgoing focusing function has to be retrieved using the following series (van der Neut et al., 2015a): where i is the iteration number. Symbol x F denotes focal points at the redatuming level ∂D i that become virtual sources. Operators R and R perform a multi-dimensional convolution or correlation of the reflection response R with the wavefield that it acts upon. Window functions are tapered Heaviside step functions that separate the causal and the acausal wavefields (i.e. Green's functions and focusing functions) in time. See appendix A of Staring et al. (2018a) for details on the design of window function . The first update of the coda of the downgoing focusing functionf + 1 already contains many of the correct events to compensate for the inhomogeneous overburden, but with incorrect amplitude. Higher order estimates (i = 2, 3, 4, etc.) are needed to obtain the correct amplitude.
Using the downgoing focusing functionf + , we can also retrieve the receiver-redatumed upgoing Green's function: where mute = I − now selects the causal wavefield and symbol x F represents focal points at the redatuming level ∂D i that become virtual receivers. The iteration number is given by j. Initial estimateG − 0 is the standard receiver-redatumed upgoing Green's function at x F . The first updateG − 1 contains a first-order estimate of the receiver-side internal multiples generated in the overburden with incorrect amplitudes. Next updates (G − 2 ,G − 3 , etc.) contain higher order estimates that are necessary to obtain the correct amplitudes. An additional step is needed to also remove source-side and source-and-receiverside internal multiples generated by the overburden.
The retrieval of the upgoing Green's functionG − with a grid of sources at the acquisition surface ∂D 0 and a grid of virtual receivers at the redatuming level ∂D i is a single-focusing step. By creating double-focusing, we also remove other internal multiples generated by the overburden. To this end, we convolve the upgoing Green's functionG − at virtual receivers with the downgoing focusing functionf + at virtual sources van der Neut et al., 2018;Staring et al., 2018a) By applying this for many positions x F and x F at redatuming level ∂D i , a grid of downward radiating virtual sources and virtual receivers that measure the upgoing wavefield is created. The result is a redatumed Green's function̑G −+ in the physical medium. Internal multiples generated by the overburden (Fig. 4a) have been removed, but later arriving internal multiples generated by interactions between the target area and the overburden (Fig. 4b) and internal multiples generated by the target area ( Fig. 4c) remain. According to Cypriano et al. (2015), the main internal multiples that contaminate the image of the target area in the Santos Basin are generated between the water bottom and the top of salt (see Fig. 1). By using doublefocusing, we remove these internal multiples, while leaving some internal multiples below the target area behind. Note that one user-specified wavelet· has to be deconvolved from the redatumed Green's function̑G −+ . Also note that the integral over the acquisition surface ∂D 0 allows us to parallelize the implementation of the double-focusing method per pair of focal points, which makes this method particularly suitable for the application to large three-dimensional (3D) data volumes.
Using equations (1) and (2), we can write equation (3) as a series:G The first termG − 0 * f + 0 is the standard source-andreceiver redatumed Green's function including primaries and internal multiples. The second term −G − 1 * f + 0 contains firstorder predictions of receiver-side internal multiples generated Figure 4 Cartoons illustrating the redatumed sources at x F and the redatumed receivers at x F resulting from the adaptive double-focusing method. Internal multiples (a) generated by the overburden have been removed by double-focusing, while (b) later arriving internal multiples generated by interactions between the target area and the overburden and (c) later arriving internal multiples generated by the target area remain. by the overburden, while the third term −G − 0 * f + 1 contains first-order predictions of source-side internal multiples generated by the overburden and the fourth term −G − 1 * f + 1 contains first-order predictions of source-and-receiver-side internal multiples generated by the overburden. Subsequent terms contain higher order estimates of the predicted internal multiples that are needed to obtain the correct amplitude. Note that Marchenko methods in principle do not rely on an adaptive filter to accurately attenuate internal multiples. However, instead of retrieving all terms in the series in equation (4) by correlating and convolving the data with itself many times (see equation 2), we propose to retrieve only a few updates and use an adaptive filter as a substitute for higher order amplitude corrections. Also, an adaptive filter can correct for any minor amplitude and phase differences that are present in the internal multiple predictions due to imperfections in the data acquisition and preprocessing. Iteration numbers i and j that are needed to obtain predictions of all overburden internal multiples depend on the geological setting. In our case, new internal multiples were not predicted beyond the third term in equation (4), so we only use the terms −G − 1 * f + 0 and −G − 0 * f + 1 for the prediction of internal multiples in this particular setting. These predictions are treated as individual internal multiple predictions, which are orthogonalized to the data prior to simultaneous adaptive subtraction. We have chosen for an adaptive filter in the curvelet domain (Herrmann et al., 2008;Wu and Hung, 2015), since it can distinguish between primaries and internal multiples in space, time and dip. Naturally, care has to be taken not to subtract the primary reflections with the internal multiples.

S E N S I T I V I T Y T E S T S O N T H R E E -D I M E N S I O NA L S Y N T H E T I C DATA
We perform a series of three-dimensional (3D) synthetic tests to identify the key acquisition parameters that affect the result of the adaptive double-focusing method. In order to generate synthetic data that represent the geological contrasts in the area as realistically as possible, we use a velocity model (see a two-dimensional (2D) slice in Fig. 1) and a density model that are obtained from an acoustic inversion of field data based on the original seismic image and migration velocity. The grid size of these models is 18.75 m by 18.75 m by 10 m. Co-located sources and receivers are positioned with a spacing of 50 m in the inline direction and a spacing of 75 m in the crossline direction ( Fig. 2a), thereby simulating an inline spacing of 50 m and a sail line spacing of 75 m. The inline aperture is 20 km (offsets from −10 to 10 km), and the crossline aperture is 1.8 km. An acoustic finite-difference algorithm is used to model data up to 30 Hz, such that the dominant wavelength at the receivers is 50 m. The recording time is 8.5 s. Also, we generate an initial focusing functionf + 0 in the smooth velocity model using an Eikonal solver. Geometrical spreading is part of the simulation. In addition, we convolve the response with an Ormsby wavelet with tapers at the low and the high ends.
Starting from 24 lines of data modelled on this dense acquisition grid, we step-by-step decimate down to a realistic streamer acquisition geometry with a cable spacing of 150 m, a sail line spacing of 450 m and a cable length of 6 km. The inline source and receiver spacing remain 50 m. Inline offsets range from 250 to 6250 m, and the crossline aperture is 0.75 km (Fig. 2b). Throughout the decimation tests, we use the Marchenko double-focusing method to redatum to a grid of co-located virtual sources and virtual receivers below the overburden with a spacing of 25 m by 37.5 m. We have chosen the redatuming level to be just above the base of salt. The base of salt is the top of our reservoir and is therefore part of the target area. The main internal multiple generators in this geological setting, the water bottom and the top of salt, are part of the overburden. Internal multiple predictions are obtained by convolving the individual updates of the wavefields G − j andf + i with one another. We orthogonalize the predictions and the data before subtraction, but do not use a global leastsquares filter for pre-conditioning. Next, the internal multiple predictions are simultaneously subtracted from the data using a 3D curvelet filter (Wu and Hung, 2015). Parameters that need to be set are the number of scales in the transform, the number of angles in the transform, the window size and some sparsity parameters that control the inversion. We extensively test different filter settings and obtain the best results (the least damage of the primary reflections) using seven scales, eight angles and tapered windows of 768 ms by 256 traces. These settings are used for all synthetic examples shown here.

The complete dataset
First, we apply the adaptive double-focusing method to synthetic data generated on the dense grid in Fig. 2(a). Figure 5 shows redatumed common source gathers before and after internal multiple prediction and subtraction (G − 0 * f + 0 and̑G −+ from equation 4). The common source gathers are from a virtual source in the middle of the grid of focal points, as indicated by the red star in Fig. 2(a). A difference is visible, especially in the yellow ellipses and along the yellow lines. It seems that conflicting seismic events were resolved, resulting in a better continuity of the primary events.
Next, we deconvolve an user-specified wavelet· and migrate the result. Figure 6 shows reverse time migration (RTM) images of the reflection responseȒ at the acquisition surface (note that this image was truncated at the base of salt for comparison), the redatumed Green's function including primaries and internal multiples G − 0 * f + 0 and the redatumed Green's functionG −+ after internal multiple prediction and subtraction. Figure 6(a) and 6(b) are comparable, thereby demonstrating that source-receiver redatuming was correctly performed (according to the standard primary approach).  Fig. 2(a). The images show: a) the migrated reflection responseȒ, b) the redatumed and migrated Green's function G − 0 * f + 0 including primaries and internal multiples, and c) the redatumed and migrated Green's functionG −+ after prediction and adaptive subtraction of internal multiples. Below the RTM images are depth slices at 5900 m. The yellow ellipses indicate areas in which internal multiple attenuation is most visible.
A comparison between Fig. 6(b) and 6(c) shows a distinct difference. Internal multiples indicated by the yellow curved stripes in Fig. 6(a) and 6(b) are no longer visible in Fig. 6(c). Also, the yellow ellipses indicate areas where the removal of internal multiples is clearly visible. Overall, the continuity of the reflectors has improved. Below the vertical RTM images are depth slices of the 3D RTM volume at 5900 m depth, where the internal multiples are present in Fig. 6(a) and 6(b), but have been attenuated in Fig. 6(c). Based on these results, we conclude that adaptive double-focusing performs well in terms of redatuming and predicting and subtracting internal multiples when applying it to the initial dense acquisition geometry. Note that the image in Fig. 6(c) appears to have a lower frequency content compared to the images in Fig. 6(a) and 6(b). It seems that the internal multiple reflections tend to favor the high frequencies, possibly due to the generation mechanism in the stratified salt. The resulting images in other papers on internal multiple attenuation in the Santos Basin (e.g. Griffiths et al., 2011;Cypriano et al., 2015) confirm this observation.
In the following, we continue with synthetic data without negative offsets and use source-receiver reciprocity for reconstruction before applying the adaptive double-focusing method.

A coarser sail line spacing
Since a sail line spacing of 75 m is not realistic, we study the effect of coarser sail line spacings on the result of our adaptive double-focusing method. Starting from the result in Fig. 6(c) with 75 m sail line spacing (here Fig. 7a), we compare RTM images showing the result of adaptive double-focusing when using a sail line spacing of 150 m (Fig. 7b), 300 m (Fig. 7c) and finally 450 m (Fig. 7d). The result obtained from data with a sail line spacing of 150 m looks very similar to the result obtained with 75 m sail line spacing; there are only some minor amplitude differences indicated by the arrows. A more significant difference becomes visible when decimating from 150 m sail line spacing to 300 m sail line spacing. Some internal multiples at the top of the image are no longer predicted and subtracted, probably because the traces that are necessary for the reconstruction of these multiples are missing. The realistic scenario of 450 m sail line spacing shows more internal multiples that could not be predicted and subtracted, now in the deeper part of the image as well. The depth slices confirm these observations: there is little difference between the results obtained with 75 m and 150 m sail line spacing, but internal multiple attenuation becomes less effective when moving to a sail line spacing of 300 and 450 m. Based on these tests, we conclude that the sail line spacing is a key acquisition parameter that affects our adaptive double-focusing method. Ideally, interpolation from 450 m sail line spacing to 150 m sail line spacing would be applied prior to the field data application in this geological setting. We remark that although we expect that the sail line spacing will also be a key acquisition parameter that affects the result of our adaptive double-focusing methods in other geological settings, the exact spacing at which the result is still acceptable will be different in every setting. In the following synthetic tests, we continue with a sail line spacing of 150 m, thereby assuming that the interpolation from 450 m sail line spacing to 150 m sail line spacing can be carried out correctly.

The removal of the near offsets
The responses at near offsets are typically not recorded by streamers, so we study the effect of removing the first 250 m of inline offsets. Figure 8 shows a comparison of RTM images with and without near offsets, both after internal multiple prediction and subtraction. The removal of the near offsets deteriorates the result somewhat in terms of a few remnant internal multiples (at the ellipse and at the arrows), but not as much as expected. The depth of the first reflector influences how much the near offset responses contribute to the image. Since the water bottom in this setting is very deep (see Fig. 1), most reflections originating from this depth would have simply not been recorded by the first 250 m of receivers. Even though the near offset responses do not have a large effect on the result of internal multiple prediction and removal in this very deep marine setting, we will interpolate the field data for the miss-ing offsets in order to predict as many internal multiples as accurately as possible.

The removal of the far offsets
Next, we assume that we could correctly reconstruct the responses at near offsets and we study the effect of removing the far offsets (the inline offsets 6250-10,000 m). Figure 9 shows the RTM images of the result of internal multiple prediction and subtraction using adaptive double-focusing, where Fig. 9(a) shows the result when including far offsets in the reflection response, and Fig. 9(b) shows the result when excluding far offsets from the reflection response. Only minor differences are visible in the result, mostly in terms of amplitudes. Surprisingly, the far offsets seem to have little impact on the result of adaptive double-focusing, similar to the near offsets. Verschuur (2013) reports that missing offsets have a particularly large effect on multiple prediction methods in a shallow water setting. Since we are in a very deep marine setting, missing offsets seem to only have a minor effect.

The removal of the outer cables
Lastly, we study the effect of removing the outer cables. Instead of a crossline aperture of 1800 m, as used in the previous tests, we now use a crossline aperture of 750 m. The RTM images in Fig. 10 show that removing the outer cables has a significant effect on the adaptive double-focusing result. The quality of the image in Fig. 10(b) has deteriorated, and some internal multiples were not predicted and subtracted. Although the missing outer cables have a large effect on the result of adaptive double-focusing, the image in Fig. 10(b) is still of acceptable quality. This becomes especially clear when comparing it to the standard redatumed Green's function in Fig. 6(b), which is constructed from a dense and wide azimuth grid of sources and receivers at the acquisition surface. Compared to this image, Fig. 10(b) still shows a significant reduction in internal multiple energy. This is promising for the field data application, since we cannot compensate for missing outer cables during preprocessing. We remark that the effect of removing the outer cables is expected to become more severe in geological settings with strongly dipping reflectors in the crossline direction. In those cases, the missing outer cables can be a limiting factor that hinders the application of the adaptive double-focusing method. This observation is supported by reports on the performance of similar multiple prediction  and removal methods (Moore and Dragoset, 2008;Wang and Hung, 2014).

The combination of all effects
Although the results of the synthetic tests in the previous sections are encouraging, they do not provide an indication on the feasibility of the application of the adaptive double-focusing method to our field dataset, where all acquisition restrictions are imposed simultaneously. The negative offsets, the near offsets, the far offsets, some sail lines and the outer cables are all missing. Therefore, we model 32 lines of 3D synthetic data based on the acquisition geometry of our narrow azimuth streamer data (Fig. 2b). Next, we reconstruct the negative offsets (by applying source-receiver reciprocity) and the near offsets (by interpolation) and perform interpolation for the sail line spacing (from 450 to 150 m). Figure 11 shows the RTM images of the reflection responseȒ, the standard redatumed Green's function G − 0 * f + 0 with primaries and internal multiples and the redatumed Green's functionG −+ after internal multiple prediction and subtraction, zoomed in at the target area. We observe an unexpected difference in illumination when comparing Fig. 11(a) and 11(b), especially on the right side of the images. Figure 11(a) is constructed by applying an RTM method to the reflection responseȒ, which uses a finite-difference method to back-propagate the wavefield from the acquisition surface to the redatuming level. In contrast, Fig. 11(b) is constructed by first back-propagating the reflection response R using convolutions with the modelled direct downgoing focusing function, according tof + 0 * R * f + 0 (see equation 3), to obtainG − 0 * f + 0 , which is subsequently back-propagated from the redatuming level into the target using the RTM method. In principle, back-propagation using a multi-dimensional convolution is equivalent to Figure 11 RTM images of the result after applying the adaptive double-focusing method to 3D synthetics with a field data acquisition geometry. The images show (a) the migrated reflection responseȒ, (b) the redatumed and migrated Green's function G − 0 * f + 0 including primaries and internal multiples and (c) the redatumed and migrated Green's functionG −+ after prediction and adaptive subtraction of internal multiples. Below the RTM images are depth slices at 5900 m. back-propagation using an RTM method (Esmersoy and Oristaglio, 1988). However, in practice, these are only equivalent when the same numerical method is used. We use an Eikonal solver to model the direct wavef + 0 , which is different from the finite-difference method used in the RTM method. As a result, there are slight differences in illumination between the two images.
Next, we evaluate the effectiveness of internal multiple attenuation. A comparison of Fig. 11(b) and 11(c) shows that the adaptive double-focusing method succeeded in predicting and subtracting internal multiples from the standard redatumed Green's function. Especially inside the yellow ellipses, the internal multiple energy is significantly reduced, resulting in a better continuity of the reflectors. Again, we observe that the internal multiples seem to mainly have a highfrequency content. We conclude that the adaptive doublefocusing method appears to be sufficiently robust for the prediction and adaptive subtraction of internal multiples from narrow azimuth streamer data in this geological setting.

T H E T H R E E -D I M E N S I O NA L F I E L D DATA A P P L I CAT I O N
Based on the results of the synthetic tests, we continue with the field data application. We have 24 lines of narrow azimuth streamer data, acquired using six flat streamers with a cable spacing of 150 m and a sail line spacing of 450 m. The length of the cables is 6000 m, covering offsets from 250 to 6250 m. The source and receiver spacing is 50 m in the inline direction. The crossline aperture is 750 m. Deghosting is performed in the frequency-p (f-p) domain (Wang et al., 2013), and a designature filter is obtained from the water bottom reflection. Shot and near offset reconstruction are performed using a par-tial normal moveout (NMO) correction of traces per common depth point (CDP) (e.g. Dragoset et al., 2010). In addition, the data are projected on a regular grid using a τ − p transform (Wang and Nimsaila, 2014). We also remove surface-related multiples, the evanescent wavefield and noise. After preprocessing, we obtain a dataset with a sail line spacing of 150 m. Similar to the synthetic tests, we have chosen the redatuming level just above the base of salt. A smoothed version of the velocity model in Fig. 1 and an Eikonal solver are used to model our direct downgoing focusing functionf + 0 (including geometrical spreading), which we subsequently convolve with a 30-Hz Ormsby wavelet. Convergence is tracked by computing the L 2 norm of the updates of the downgoing focusing functionf + i . After convolving the individual updates of the wavefieldsf + i and G − j and migrating them, we obtain the internal multiple predictions in the image domain. Extensive testing shows us that primary reflections are better preserved when Figure 13 Depth slices corresponding to the RTM images in Fig. 12 at 5900 m. The slices show (a) the migrated reflection responseȒ, (b) the redatumed and migrated Green's function G − 0 * f + 0 including primaries and internal multiples and (c) the redatumed and migrated upgoing Green's functionG −+ after prediction and adaptive subtraction of internal multiples.
subtracting the internal multiple predictions in the image domain instead of in the redatumed domain. Prior to subtraction, the predictions are orthogonalized. We use the full curvelet transform (for all scales) and tapered windows of 768 ms by 256 traces. Figures 12 and 13 show the result of applying the adaptive double-focusing method to predict and adaptively subtract internal multiples. First, we observe that there is still a slight illumination difference between the RTM migrated image of the reflection response in Fig. 12(a) and the image of the RTM migrated redatumed reflection response in Fig. 12(b) (see the yellow ellipse in the top right), which has been explained above. However, the difference is not as pronounced as in Fig. 11. Second, a clear difference is visible between Fig. 12(b) and 12(c), indicating that the adaptive double-focusing method succeeded in predicting and adaptively subtracting events which are likely internal multiples. The yellow arrow, stripes and ellipses indicate areas in which events were attenuated, while the yellow boxes indicate what we interpret to be an improvement in fault definition. The white circle shows conflicting events that are being resolved, but it also shows events that we believe to be remnant internal multiples. These were most likely generated by the base of salt, which is not part of our overburden. In order to also remove these events, we can obtain additional internal multiple predictions by repeating our procedure while placing the redatuming level just below the base of salt. The RTM depth slices in Figure 13 demonstrate the attenuation of events in the form of ellipses. This example shows that the adaptive doublefocusing method is sufficiently robust for the application to our three-dimensional (3D) narrow azimuth streamer dataset. Naturally, this result can be improved by the use of more lines of data.

D I S C U S S I O N A N D C O N C L U S I O N
In this paper, we identified the key acquisition parameters that affect the application of our adaptive Marchenko internal multiple attenuation method to narrow azimuth streamer data. Tests on three-dimensional (3D) synthetic data evaluated the effect of removing sail lines, near offsets, far offsets and the outer cables. The results of these tests show that the aperture in the crossline direction and the sail line spacing have the strongest effect on the quality of the result. Typically, the sail line spacing can be interpolated, but the aperture in the crossline direction can possibly be a limiting factor for our method. Surprisingly, the missing near offsets and the far offsets only had a modest effect on the result of our method, possibly due to the very deep target area. In addition, we found that the responses at the negative offsets and the near offsets could be accurately reconstructed. We remark that these tests are only valid for this particular dataset, but they give an impression of the possibilities and limitations of the adaptive Marchenko double-focusing method. For an ocean bottom node acquisition geometry, these tests imply that our method will be most sensitive to the node separation (especially in the direction of the strongest geological variation).
Based on the decimation tests, we defined an interpolation strategy that was first tested on a realistic synthetic dataset. We reconstructed the negative offsets and the near offsets and interpolated the sail line spacing from 450 to 150 m. When applying the 3D adaptive double-focusing method, an important aspect that was not visible in earlier two-dimensional (2D) applications became visible, thereby showing that the extension of a method from 2D to 3D is not always trivial. In 3D, when using an Eikonal solver for the modelling of the direct wave and a finite-difference-based RTM method, a slight difference in illumination between the RTM image of the reflection response and the RTM image of the redatumed response occurs. Nevertheless, the doublefocusing method predicted and adaptively subtracted internal multiples, thereby improving the image of the target area.
Next, we applied the adaptive double-focusing method to 24 lines of narrow azimuth streamer data. We reconstructed the negative and the near offsets and interpolated the sail line spacing. We interpret that internal multiples were predicted and adaptively subtracted, which resulted in an improved geological interpretation of the target area. Therefore, we conclude that 3D Marchenko internal multiple attenuation using an adaptive double-focusing method is sufficiently robust for the application to narrow azimuth streamer data in a deep marine setting, provided that there is sufficient aperture in the crossline direction and that the sail lines are interpolated.
Note that redatuming is optional for Marchenko methods. The adaptive double-focusing method used in this paper includes source-receiver redatuming, which is particularly useful when the aim is to attenuate internal multiples in the target area and at the same time reducing the data volume for a next processing step (e.g. a target-oriented full waveform inversion). In contrast, when the aim is only to attenuate internal multiples as part of a larger workflow, the adaptive double-focusing method might not be the Marchenko method of choice. A direct quality check of the input common source gathers and the redatumed common source gathers is not possible, which is a disadvantage in a general processing workflow. In addition, a quality check on the resulting images is only possible when the same numerical method is used to obtain the direct wave for the Marchenko method as for the migration of the original data. Therefore, for the purpose of internal multiple elimination only, we propose the use of other Marchenko methods that do not include redatuming and thus allow for an easier quality check. An example is the adaptive overburden elimination method (van der Neut and , as shown in papers by Pereira et al. (2018), Krueger et al. (2018) and Pereira et al. (2019). A modified version of the adaptive double-focusing method is presented by Staring et al. (2019). Other alternatives are the Marchenko multiple elimination scheme  and the primary-only method proposed by Meles et al. (2016).
The Marchenko method used in this paper is acoustic. The synthetic data are acoustic, but naturally the field data are elastic. A suggestion for further research is to evaluate the effect of the presence of mode conversions on the acoustic Marchenko method in this geological setting, as was done by Reinicke et al. (2019) for the offshore Middle East. By applying an acoustic Marchenko method to elastic synthetic data, and comparing it to the result of applying an acoustic Marchenko method to acoustic synthetic data, Reinicke et al. (2019) evaluated whether the acoustic approximation is valid for structural imaging in the region. They concluded that the acoustic approximation may be sufficient when used for structural imaging in 1.5D geological settings.

AC K N OW L E D G E M E N T S
This research was performed in the framework of the project 'Marchenko imaging and monitoring of geophysical reflection data', which is part of the Dutch Open Technology Programme with project number 13939 and is financially supported by NWO Domain Applied and Engineering Sciences. The research of K. Wapenaar has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 742703). We thank Roberto Pereira, Peter Mesdag and Adel Khalil from CGG for the close collaboration and for providing the field data. Also, we would like to thank two anonymous reviewers, Lele Zhang, Giovanni Meles and Jan Thorbecke for valuable discussions.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data sharing is not applicable to this article.