{'fraction_non_alphanumeric': 0.02482758620689655, 'fraction_numerical': 0.007619047619047619, 'mean_word_length': 4.88084202394747, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Extracting small time-lapse traveltime changes in a reservoir using primaries and internal multiples after Marchenko-based target zone isolation

Geophysical monitoring of subsurface reservoirs relies on detecting small changes in the seismic response between a baseline and monitor study. However, internal multiples, related to the over- and underburden, can obstruct the view of the target response, hence complicating the time-lapse analysis. In order to retrieve a response that is free from over- and underburden effects, the data-driven Marchenko method is used. This method effectively isolates the target response, which can then be used to extract more precise time-lapse changes. Additionally, the method also reveals target-related multiples that probe the reservoir more than once, which further define the changes in the reservoir. To verify the effectiveness of the method, a numerical example is constructed. This test shows that when using the isolated target response, the observed time differences resemble the expected time differences in the reservoir. Moreover, the results obtained with target-related multiples also benefit from the Marchenko-based isolation of the reservoir. It is, therefore, concluded that this method has potential to observe dynamic changes in the subsurface with increased accuracy.

method has potential to observe dynamic changes in the subsurface with increased accuracy.

INTRODUCTION
Time-lapse seismic studies are concerned with detecting small changes in the seismic response between a baseline and a monitor study. These changes can either be a difference in amplitude (e.g. Landrø, 2001), a difference in traveltime (e.g. Landrø and Stammeijer, 2004) or a combination of both (e.g. Trani et al., 2011). These timelapse methods are essential for observing and monitoring subsurface reservoirs, with applications ranging from determining pressure and fluid saturation changes (Landrø, 2001) to monitoring CO 2 injection (Roach et al., 2015) or observing compaction in a reservoir (Hatchell and Bourne, 2005).
In order for these methods to work optimally, it is important that the reservoir response can be clearly identified in the seismic response. In practice this requirement is not always fulfilled, as multiple reflections from a (highly) scattering overburden can mask the response of the reservoir. It is, therefore, desirable to remove the overburden effects before applying any time-lapse analysis. The Marchenko method is able to redatum a wavefield from the surface of the earth to an arbitrary focal depth in the subsurface, while accounting for all orders of multiples Slob et al., 2014). This data-driven method can be used to remove all interactions from layers above the selected focal level, hence giving an unobstructed view of the reservoir response. From this new response the traveltime difference in the reservoir can more precisely be determined.
In addition to removing the overburden, the reservoir response can completely be isolated by also removing the underburden with the Marchenko method . Consequently, not only the primary response of the reservoir is uncovered, but internal multiples, which traversed through the reservoir more than once, will now also be clearly visible and unobstructed by primaries and multiples outside the target zone. Since these multiples have passed through the reservoir multiple times, the time-lapse traveltime change of the multiples will be larger, hence more sensitive. This is akin to coda-wave interferometry, which exploits the fact that time-lapse changes are exaggerated in the coda due to the longer paths traveled in the medium (Snieder et al., 2002;Grêt et al., 2005).
Inspired by this principle of coda-wave interferometry, Wapenaar and Van IJsseldijk (2020) show how correlation of multiples improve the ability to detect small changes in velocity compared to correlation of primaries. This method is then adapted to find changes in lateral varying media (van IJsseldijk and Wapenaar, 2021). In this work we further develop the method in order to account for time-lapse changes in the overburden. First, we revise the theory of isolating the reservoir response with the Marchenko method and we review how to extract traveltime changes from this isolated response. Subsequently, we present a numerical model that will be used to test the methodology. The reservoir response is then isolated from the modeled data, and the traveltime changes of the primary as well as the multiple reflections are calculated. Finally, we discuss the results and possible future improvements to the method. Figure 1 shows the principle of the method. Here the acoustic situation is considered, with a reservoir enclosed by two strong reflectors. Note how the primary from the first reflector, does not probe the reservoir, whereas the primary from the second reflector does. Moreover, the internal multiples generated by these reflectors will traverse the reservoir additional times, hence they travel the reservoir multiple times and experience a larger traveltime shift. In order to achieve this same situation from a regular reflection response measured at the surface, the medium is first divided into 3 parts: overburden "a", target zone "b" which contains the reservoir and two reflectors as in Figure 1, and underburden "c". The reflection response of the full medium is denoted by R abc (x R , x S , t), here x R , x S and t denote the receiver position, source position and time, respectively. Our first aim is to isolate the reflection response R b of the target zone with the help of the Marchenko method, which will briefly be discussed in the next section.

Extrapolated Marchenko representations
At the base of the Marchenko method are the focusing functions (f ± 1 ) that allow for retrieval of the Green's functions (G −,± ) between the acquisition surface S 0 and a focal level in the subsurface, the left superscript − denotes that the wavefield is upgoing at the receiver position and right superscript ± denotes a down-or up-going direction from the source position. Van der Neut and Wapenaar (2016) introduce modified functions that are extrapolated to the surface by convolution with the direct arrival of the transmission response (T d ). These extrapolated focusing functions (v ± ) are defined as follows: Here x F is the coordinate of the focusing point at focal depth S F , x S is a coordinate on the acquisition surface, and * denotes temporal convolution. Similarly, the extrapolated Green's functions (U −,± ) are defined as: Finally, the same convolutions are applied to the coupled Marchenko representations to find the extrapolated representations: The reflection response is denoted by R. In this paper this response will either be the response of the full medium R abc or the response after overburden removal R bc .
Moreover, these two equations have four unknowns. In order to solve this system a causality constraint is introduced, which exploits the fact that the focusing and Green's functions are separable in time. In order to apply this constraint, an estimate two-way travel time between the focal level and the surface is required. In our case, this is achieved by computing the direct arrival of the Green's function in a smooth velocity model with an Eikonal solver, and then convolving this response with itself to find the two-way travel time. A more elaborate derivation of the Marchenko method is beyond the scope of this paper, instead the reader is referred to Wapenaar et al. (2021), who give more background on both the regular and extrapolated expressions.

Isolation of the reservoir's response
Using the relations presented in the previous section, the focusing and Green's functions above and below the reservoir can now be retrieved. From these functions the reflection response of the target zone can be isolated. First, the overburden is removed, using the extrapolated Green's function at S 1 between the overburden and target zone : Here U −,± a|bc are the extrapolated Green's functions computed from Equation 3 and 4, where R abc is used as reflection response R. The vertical line in the subscript indicates the location of focal level, i.e. between the underburden "a" and the target-zone "b". Next, the newly acquired reflection response (R bc ) is used to compute the extrapolated focusing functions between the target zone and the underburden at S 2 . These focusing functions are then used to remove the underburden: The subscript b|c denotes that R bc was used to compute the focusing functions from Equation 3 and 4, with the focal level between target zone "b" and underburden "c".
Again, the reflection response of the target zone R b can be retrieved from Equation 6 by means of MDD.
Effectively, the target zone response has now been isolated, leaving a response analogous to the situation in Figure 1. During numerical tests it was found that the multiples in the response could be enhanced by applying a small, constant scaling factor to the reflection response R bc before applying the Marchenko scheme for the second time. However, this means that the Marchenko iterations are no longer unconditionally stable and hence some care is required when enhancing the amplitudes of the reflection response. Consequently, a new scheme is proposed to progressively enhance the scaling in the data. This is achieved by starting the operation with the regular scaling, which always produces a stable result, followed by a feedback loop with a gradually increasing scaling factor until the result becomes unstable or the multiples are no longer amplified. Figure 2 shows an overview of the process to isolate the target response. In this chart only a smooth velocity model of the baseline is used, it is assumed that this model can also be used for the monitor study, because the velocity changes are relatively small and only an approximation of two-way traveltime to the focal depth is needed. Next, the new response R b will be used to extract the traveltime shifts in the reservoir.

Extraction of time differences
Before extracting the traveltime shifts in the reservoir, the different primaries and multiples are identified. Primary 1 and 2 are easily detected due to the isolation of the target zone (i.e. there are no interactions from the overburden to obscure the primaries). Subsequently, the arrival times of the internal multiples can be approximated based on the primaries, where the arrival time of the n-th multiple can be approximated by the arrival time of primary 2 plus n times the difference in time between the two primaries. The first step is now to eliminate any time shifts from a time-lapse change in the overburden. In order to do this the temporal correlation between primary 1 (i.e. the response that is not penetrating the reservoir) and primary 2 or an internal multiple (i.e. the responses that go through the reservoir) is computed. This gives the correlation between primary 1 and the target response below the reservoir: Here, C is the correlation of the two responses, and x 0 denote the position of the zero-offset traces in the data, where x S = x R . This correlation contains the timelag between the first reflector and P1, M1 or M2. Theta is a time window or mute function that isolates a specific primary or multiple as follows: In these equations the asterisk * can be replaced with P1, P2, M1 or M2, for primary 1, primary 2, multiple 1 and multiple 2, respectively. Hence t * is the travel time of one of the primaries or multiples. is a small constant that defines the window, and makes sure the whole waveform is included. Any traveltime differences in the overburden are removed by first calculating the time-lag with primary 1 in Equation 7, i.e. the time difference between primary 1 and either P2, M1 or M2 is free from overburden interactions. Next, from these time-lags for both the baseline and monitor study, the traveltime shifts in the reservoir can be determined with a second temporal correlation as follows: The bar denotes that the monitor correlation is used, thus the time lags, of primary 2 or one of the multiples with respect to primary 1, for the baseline and monitor are correlated. Next, the argument of the maximum of this correlation is taken to determine the traveltime differences in the reservoir. The process of extracting the time shifts is summarized in Figure 3. Note that there are a few additional practical considerations included in the chart, such as resampling and removing outliers. These will be discussed in more detail in the next section.

NUMERICAL EXAMPLE
A numerical experiment is designed, to show the viability of the method. The baseline velocity and density models are shown in Figure 4 A and B, respectively. between the surface and the focal points at 675m and at 1100m, are acquired using an eikonal solver in a smooth version of the baseline velocity model of Figure 4A. Lastly, a band-limited delta pulse is computed, which is used as initial focusing function (v + 0 ) in the iterative Marchenko scheme. Figure 5A shows a zero-offset section of the initial reflection response before applying the multiple internal removal (R abc ). Note that only the zero-offset data are

Results of target zone isolation
shown, but the data at all offsets are available and used to compute the isolated reservoir response R b from R abc . Due to the highly scattering overburden, the primaries (P1 in red and P2 in white) are nearly impossible to identify, and multiple 1 and 2 (in blue and orange) are completely obscured by the overburden and underburden interactions. After removing the overburden and acquiring R bc ( Figure 5B), the primaries are now clearly recognizable in the seismogram. However, the windows enclosing the multiples contain undesirable events from primary reflections from the underburden. The third panel (5C) reveals that these events are successfully removed after underburden removal. Also note how the multiple events are enhanced in the fourth display (5D), this is the result of the previously described scaling factor to retrieve the focus functions from R bc (which was scaled with a factor 2).
All the time-windows that select the primaries and multiples in Figure 5 are picked from this final response in panel C (R b ). First, the arrival time for both primaries are manually selected. As stated before these times are then used to provide an estimate for the multiple arrivals. Finally, the windows are manually adjusted to ensure they include the full response from each individual event. In the next section these windows will be used for the cross-correlations that compute the time differences.

Time differences inside the reservoir
From the isolated response the traveltime changes can be estimated. First, the data are subsampled from 4ms to 1ms to achieve a better time resolution. Second, the primary enclosed in the red window is selected from the data, similarly either the second primary, first multiple or second multiple is also extracted using their respective window (as shown in Figure 5). For both the baseline and monitor studies, these responses are then correlated to effectively remove the time differences developed in the overburden. Finally, the baseline and monitor correlations are mutually correlated to find the time differences in the target zone.
The results of these final correlations are shown in Figure 6, note that these results were acquired with the true Marchenko scaling without the previously described multiple enhancement. Here, the response of the full medium, the response after overburden removal and the response of the target zone are displayed with blue, orange and green lines, respectively. The light-blue area marks the 1D-based zero-offset traveltime difference, which was computed by multiplying two times the reservoir thickness with the slowness change in the reservoir. However, this reference solution does not take into account lateral variations. The red line gives a second reference solution made by cross-correlating "ideal" data. This data was acquired in a model with a smooth overburden (A), the same target zone (B) as the actual model and a transparent underburden (C). The zero-offset response is then modeled with finitedifferences, providing an accurate isolated response of the target zone R b for both the monitor and baseline response. Subsequently, the primaries and multiples are identi-fied in the modeled zero-offset response, and then correlated as described in flowchart in Figure 3. Figure 6A shows the results for primary 2. Although, all three responses capture some differences in the reservoir, the response of the full medium still contains overburden events as shown by the time-shifts greater than 0ms. These positive time-shifts are almost fully removed after overburden removal. Note that, on the one hand, the correlations do not match the 1D reference very well, because of the lateral variations in the model. On the other hand, the match with correlation of the ideal data is a lot better, which implies that the Marchenko based isolation was successful.
Based on these results it is concluded that the expected time differences are smaller than 15ms for P1, M1 and M2. This observation is used to achieve more accurate results, by removing outliers that give a time difference greater than 15ms at any lateral distance (i.e. they are removed before applying the Gaussian smoothing along the lateral direction in Figure 3).
Next, the procedure is applied to multiple 1 and multiple 2, the results of which are shown in Figure 6B and 6C, respectively. This time none of the results match either reference, and seemingly no meaningful information can be acquired from the multiples. However, these results are achieved without any multiple enhancement (i.e. the feedback loop in Figure 2 is ignored). In the next section we will explore how the results can be improved by using the multiple enhancement.

Results after multiple enhancement
To improve the correlations of the multiples the feedback loop from Figure  Especially, the results for primary 2 and multiple 1 accurately match the correlation of the "ideal" data. The results of multiple 2 do not have the same match, this is due to the fact that events from multiple 1 are interfering within the correlation window of multiple 2, and leaving an imprint on the reference solution. From this match with the "ideal" data, it can be concluded that the Marchenko method succeeded in correctly isolating the target response, as the results for multiple 1 coincide with the reference solution. This also highlights the importance of isolating the response, since the correlations of R abc and R bc do not come close to the reference solution at all. Finally, the results outside the reservoir compartments should show a time differences equal to 0ms, and time difference greater than 0ms indicate that the result is contaminated by effects from the overburden.The correlation results of the isolated response R b display less of these positive time differences compared to the results of the full response R abc . Consequently, this suggests that the isolation process has successfully eliminated the effects of the overburden.

DISCUSSION
Although the results show that the reservoir response can successfully be isolated and accurate time difference can be retrieved from both the primary reflection and the internal multiples, there are a number of issues that require a more in-depth discussion. In addition to this discussion, future improvements and practical considerations for the implementation on real data will be considered as well.
First, the scaling factor introduced to amplify multiple events was found experimentally. Also, it is noted that the application of this scaling factor only works for Hence the MDD of the focusing functions produces constructed multiples that can benefit from a scaling factor, but the MDD of the Green's functions uses the recorded multiples, therefore a scaling factor does not provide any benefit. The authors conducted numerous 1D experiments to get a better understanding of this behavior, and found that multiple enhancement is always possible for relatively simple models with two (strong) reflectors in the target zone. Hence, this method lends itself perfectly for our application, where only two reflectors remain after overburden removal. These findings were confirmed on a handful of 2D models (e.g. the model used in van IJsseldijk and Wapenaar, 2021). However, the experiments also demonstrated that if the scaling factor becomes too large the solution is no longer stable. To deal with these ambiguities, it is important to iteratively apply the scaling, in order to ensure a stable solution. In the first iteration the correct scaling is applied, therefore the solution is always stable. Subsequently, the scaling can be increased with small steps until the solution is no longer stable or no longer enhances the multiples.
A second point of discussion is the order of operations used to isolate the target response. In theory it does not matter whether the underburden is removed before or after overburden removal. Numerical experiments indeed showed that removing the underburden before the overburden is also a viable approach to isolate the target response, by first retrieving R ab from a MDD of v ± ab|c and then R b from U −,± a|b . However, the previously discussed multiple enhancements are no longer available when the method is applied in this order, because the MDD of the focusing functions would be applied before the MDD of the Green's functions that require proper scaling.
Therefore, it was decided to start with overburden elimination followed by removal of the underburden.
Third, the method is designed to use as little a priori information as possible, needing solely three prerequisites: the baseline reflection response (R abc ), the monitor reflection response (R abc ) and a smooth version of the baseline velocity model. The smooth velocity model is used to approximate the two-way travel time between the surface and the focal depth. The same model can be used for both the baseline and monitor reflection response, because it is assumed that the velocity changes in the medium are relatively small. Additionally, the extrapolated Marchenko functions are also less sensitive to velocity variations compared to regular functions. If the velocity changes are larger, a separate model is required to isolate the target zone from the monitor reflection response.
Next, the primary results before the target zone isolation already match the actual results quite closely, and only minor deviations due to overburden effects are present (i.e. the parts where ∆t > 0ms). However, these results indirectly benefit from the isolated result, because the windows, which are used to identify the primaries, are selected from the isolated results ( Figure 5C). Moreover, previous results by van IJsseldijk and Wapenaar (2021) show that the correlations from R abc are insufficient to get accurate time differences in less complex models with overburden events interfering with the primaries.
Application of the Marchenko algorithm to field data can be quite cumbersome.
Especially the MDD that is used to remove the overburden effects tends to be very sensitive to errors in the amplitude of the data. To overcome this limitation, Staring et al. (2018) introduce a double-focusing method, which is more stable but leaves some remaining interactions of the overburden. A similar approach is envisioned to acquire R bc , when applying this method to real data.
Finally, we would ideally find the velocity change of the reservoir rather than the time differences. For very simple situations a similar approach as coda-wave interferometry can be considered, which finds the velocity perturbation from the change in traveltime and initial velocity (Snieder, 2006). However, this only holds when the relative velocity perturbation is constant at every location. In our case, the perturbation is different outside the reservoir, hence the relation does not hold. Alternatively, the velocity perturbation can be found by inversion of cross-correlations at all offsets (instead of just the zero-offset data used here). Compared to the traveltime differences, the velocity changes can more directly be related to physical processes due to flow in the reservoir. Currently, this is subject to ongoing research.

CONCLUSION
A good understanding of fluid flow, temperature variations and mechanical changes in subsurface reservoirs is essential for a large variety of geoscientific methods. These dynamic changes can be observed with seismics by identifying amplitude changes, time shifts or both, between a baseline and a monitor study. However, the response of a subsurface target can be obscured by interferences from reflectors in the overburden and/or underburden, making it harder to detect the time-lapse effects.
The Marchenko method can be used to remove primaries as well as internal multiples above or below an arbitrary focal level in the subsurface from the reflection response.
Hence, this method can be used to isolate the reservoir response in both the baseline and monitor response, providing an unobstructed examination of changes in the target zone.
A twofold methodology has been proposed to extract time differences. With this methodology, first the target response is isolated, by overburden removal, followed by underburden removal. This new response is then used to identify the primary and multiple reflections in the target zone. Second, time differences are retrieved by cross-correlating the different reflections of the baseline and monitor studies. By first correlating the response with primary 1 above the reservoir, all possible time shifts of the overburden are removed.
A numerical example with a strongly scattering overburden was created to test the methodology. The isolation of the target zone revealed the primary responses of the reservoir, allowing their extraction from the data. Next, the time differences of the reservoir could be approximated from correlations with a primary reflection below the reservoir. Furthermore, additional traveltime changes were retrieved from the first and second order multiples, created by the two reflectors enclosing the reservoir.
These multiples confirmed the earlier observations, but also improved specific blind spots in the original approximation of the time changes.
The proposed methodology provides a new means to extract traveltime differences, especially for situations where complex overburden and underburden interactions mask the target response. Future developments should also make it possible to invert for velocity changes in the reservoir, rather than time differences. The method will then enable us to more accurately observe dynamic changes in the subsurface.

DATA AVAILABILITY
Data associated with this research are available and can be accessed via the following URL: https://gitlab.com/geophysicsdelft/OpenSource. The figures in this paper can be reproduced using the demo in the "vmar" folder of the repository.

LIST OF FIGURES
1 Graphic displaying the principle of the method. Note how the reservoir layer is located in between two reflectors. Primary 1 from reflector 1 does not propagate through the reservoir, whereas primary 2 from reflector 2 does. The multiples (1 and 2) are propagating through the reservoir twice or thrice, hence experiencing double or triple the traveltime changes compared to primary 2. Target zone "b" is located in between focal depths S 1 and S 2 , overburden "a" and underburden "c" are above and below the target zone, respectively.
2 Flowchart depicting how the reservoir response is isolated with the Marchenko method.
3 Flowchart to get the time differences in the reservoir from the isolated response (continuation of Figure 2). 4 Velocity (A) and density (B) model of the baseline study for the numerical example. The black dashed lines define the focal levels above (S 1 ) and below the reservoir (S 2 ), used for the Marchenko method. The solid white contour depicts the three different reservoir pockets. C shows the difference in velocity between the baseline and monitor study. The density inside the reservoir pockets is also increasing with 100 kg/m 3 for the monitor study, with no density changes outside the reservoir.
Primary 1 and Primary 2 originate from the green to blue contrast at 700m to 900m and the blue to green contrast at 1000m, respectively.

5
Zero-offset gathers for the baseline reflection response of the entire medium (C), after isolation and multiple enhancement of R b (D). The windows used for crosscorrelation for primary 1, primary 2, multiple 1 and multiple 2 are highlighted in red, white, blue and orange, respectively.

6
Results showing the estimated time difference in the reservoir from primary 2 (A), multiple 1 (B) and multiple 2 (C), no scaling factor has been applied to enhance the multiples. First, the time-lag with primary 1 is computed for the baseline and monitor study. 2) are propagating through the reservoir twice or thrice, hence experiencing double or triple the traveltime changes compared to primary 2. Target zone "b" is located in between focal depths S 1 and S 2 , overburden "a" and underburden "c" are above and below the target zone, respectively.
Extracting small time-lapse differences from a reservoir

Marchenko-based isolation of the reservoir response
Smooth baseline velocity model  Figure 2: Flowchart depicting how the reservoir response is isolated with the Marchenko method.    the three different reservoir pockets. C shows the difference in velocity between the baseline and monitor study. The density inside the reservoir pockets is also increasing with 100 kg/m 3 for the monitor study, with no density changes outside the reservoir.

Marchenko-based isolation of the reservoir respons
Primary 1 and Primary 2 originate from the green to blue contrast at 700m to 900m and the blue to green contrast at 1000m, respectively.       Figure 7: Same results as Figure 6, but with multiple enhancement by applying a scaling factor of 2 to R ab before retrieving R b with the Marchenko method. Note that the multiple correlations are a lot closer to the ideal result compared to the results without multiple enhancement.