Feasibility of reservoir monitoring in the Groningen gas field 1 using ghost reflections from seismic interferometry

28 Seismic interferometry (SI) retrieves new seismic responses, e.g., reflections, between either receivers 29 or sources. When SI is applied to a reflection survey with active sources and receivers at the surface, 30 non-physical (ghost) reflections are retrieved as well. Ghost reflections, retrieved from the correlation 31 of two primary reflections or multiples from two different depth levels, are only sensitive to the 32 properties in the layer that cause them to appear in the result of SI, such as velocity, density, and 33 thickness. We aim to use these ghost reflections for monitoring subsurface changes, to address 34 challenges associated with detecting and isolating changes within the target layer in monitoring. We 35 focus on the feasibility of monitoring pore-pressure changes in the Groningen gas field in the 36 Netherlands using ghost reflections. To achieve this, we utilise numerical modelling to simulate scalar 37 reflection data, deploying sources and receivers at the surface. To build up subsurface models for 38 monitoring purposes, we perform an ultrasonic transmission laboratory experiment to measure S-wave 39 velocities at different pore pressures. Applying SI by auto-correlation to the modelled datasets, we 40 retrieve zero-offset ghost reflections. Using a correlation operator, we determine time differences 41 between a baseline survey and monitoring surveys. To enhance the ability to detect small changes, we 42 propose subsampling the ghost reflections before the correlation operator and using only virtual sources 43 with a complete illumination of receivers. We demonstrate that the retrieved time differences between 44 the ghost reflections exhibit variations corresponding to velocity changes inside the reservoir. This 45 highlights the potential of ghost reflections as valuable indicators for monitoring even small changes. 46 We also investigate the effect of the sources and receivers’ geometry and spacing and the number of 47 virtual sources and receivers in retrieving ghost reflections with high interpretability resolution.

Subsurface activities, such as CO2 storage, oil and gas production, and geothermal energy production, involve substantial transportation of fluids, either into or out of geological formations.The changed pore-fluid pressure can result in potential risks such as induced seismicity (Bourne et al., 2014) or Time-lapse seismic studies compare an initial baseline study with subsequent monitoring studies (Lumley, 2001, Macbeth et al., 2020).For instance, Landrø (2013) explored the use of seismic amplitude analysis to assess changes in pressure and fluid saturation.Roach et al. (2015) employed two 3D time-lapse surveys at a CO2 storage site to monitor CO2 injection processes.Additionally, Hatchell and Bourne (2005) investigated the observation of compaction in reservoirs through seismic-attributes analysis.
Amplitude, travel time, and their combination can be used as seismic attributes for the purpose of timelapse studies (MacBeth et al., 2020;Van Ijsseldijk, et al., 2023;Trani et al., 2011).However, complex overburden structures can complicate the seismic response and obscure the desired temporal variations.
The presence of geological features such as faults or heterogeneities in the overlying layers can distort the seismic signals and make isolating the changes of the target layer challenging.Moreover, detecting small reservoir changes presents a significant challenge due to noise interference and limited sensitivity.
These factors can make detecting changes in the subsurface more difficult.In order to address those challenges, we aim to show the feasibility of using ghost reflections retrieved from seismic interferometry (SI) to monitor small temporal variations in subsurface characteristics.
SI refers to retrieving seismic responses through, e.g., cross-correlation, auto-correlation, or deconvolution of seismic observations at locations of seismic receivers or sources (Wapenaar andFokkema, 2006, Wapenaar et al., 2010).When using a dataset from active sources and receivers at the surface, ghost (non-physical) reflections are retrieved from SI because of insufficient destructive interferences (Snieder et al., 2006, Draganov et al., 2010, King and Curtis, 2012).Such ghost reflections propagate inside a specific layer or group of layers.This can be used advantageously for monitoring changes inside those specific layers, e.g., a gas reservoir or CO2 reservoir.Draganov et al. (2012) verified this concept by conducting numerical modelling and scaled laboratory experiments to monitor CO2 storage.Furthermore, Ma et al. (2022) showed that ghost reflections can be used to monitor the geotechnical behaviour of fluid mud using ultrasonic reflection measurements in a laboratory.In another study, Shirmohammadi et al. (2024) investigated the potential application of ghost reflections for characterising specific layers in the shallow subsurface.They demonstrated the technique's effectiveness using synthetic data for a subsurface model with a lateral change in velocity, a velocity gradient in depth, a thickness change, a velocity change of the target layer, and also a shallow field dataset.
To demonstrate the feasibility of using ghost reflections for monitoring reservoir pressure changes, we use the Groningen gas field as a well-known example of an onshore gas field in Europe, located in the Netherlands.The extraction of natural gas from this field since 1963 has led to a series of seismic events (Van Eijs et al., 2006, andMuntendam-Bos et al., 2021).The occurrence of such seismic activity has raised concerns regarding the need for effective monitoring methods.A number of studies suggest using SI for monitoring the Groningen subsurface.Brenguier et al. (2020) used a passive-seismic monitoring approach to detect velocity changes in the Groningen reservoir in ballistic waves recovered from seismic noise correlations.Their methodology requires dense arrays of seismic sensors.For the same gas reservoir, Zhou and Paulssen (2020) investigated the potential of passively recorded deep borehole noise data to detect temporal variations using SI by deconvolution.They showed the possibility of monitoring small temporal changes in the Groningen gas field if repeating noise sources are available.
Using the same approach, Zhou and Paulssen (2022) showed that the observed travel-time changes in P-wave and P-to-S converted waves could be related to fluctuations of the gas-water contact in the observation well.
We use a synthetic reflection dataset from the Groningen subsurface model to illustrate our approach of using ghost reflections for monitoring.To build up the subsurface model for the baseline survey and the monitoring surveys, we perform an ultrasonic laboratory experiment to measure the direct S-wave velocities for variations in the reservoir-rock pore pressures using in-situ conditions of pore pressure, stress, and confining pressure.We use the velocities measured in the laboratory experiment at different pore pressures to determine the effect of pressure depletion on the time-lapse changes of seismic velocity of the Groningen reservoir to see whether we could pick up those velocity changes using ghost reflections.Thus, after building up subsurface models, we apply SI by auto-correlation (SI by AC) to the synthetic reflection data to retrieve ghost reflections from inside the Groningen reservoir for the baseline survey and monitoring surveys.To determine the time difference in the reservoir, we use a correlation operator between the ghost reflection retrieved from the baseline survey and the monitoring surveys.To validate our approach, we also calculate the (relative) time difference of the ghost reflections from monitoring surveys with the baseline survey using the expected time difference.
Moreover, we discuss the source and receiver configuration for future practical applications of ghostreflection retrieval using the field dataset.
Below, we first present in section 2 the methodology of retrieving ghost reflections with SI and calculating the time differences between the baseline survey and the monitoring surveys, as well as the validation process.Then, in section 3.1, we show the results of SI by AC when applied to data from numerical modelling for the baseline survey.Subsequently, in section 3.2, we discuss time-lapse investigation using ghost reflections.In section 3.3, we also investigate the influence of the source and receiver configuration in retrieving high-resolution ghost reflections.This is followed by a discussion to address challenges, provide recommendations for the practical application of our approach, and draw conclusions.

Method
For an active-source reflection seismic survey, where the sources (the red stars in Fig. 1) and receivers (the blue triangles in Fig. 1) are restricted to the surface, the frequency-domain response  ( ,  , ) and its complex conjugate at  from a virtual source at  can be obtained from the equation (Halliday et al., 2007)  * ( ,  , ) +  ( ,  , ) ∝  * ( ,  , )  ( ,  , ) .
The right-hand side of this equation corresponds in the time domain to a cross-correlation between two observations at positions  and  , both originating from active sources located at  at the surface.
The symbol (*) shows the complex conjugate in the frequency domain, while  representing the total number of active sources at the surface.Given the source-receiver reciprocity theorem,  and  can also represent the positions of the active sources, and the obtained response on the left-hand side of Eq. (1) 1 corresponds to the frequency domain response at virtual receivers.In this case, we turn the active sources into virtual receivers, while  represents the number of receivers at the surface.
The theory behind SI requires having sources which effectively surround the receivers completely.
When the positions  are at the surface, like for Eq. 1 (the red stars in Fig. 1) as in a typical activesource exploration survey, one-sided illumination of the receivers occurs (Wapenaar, 2006).A consequence of the one-sided illumination is that the application of Eq. 1 will retrieve not only the desired (pseudo-) physical reflections but also ghost reflections.(Note that we labelled the retrieved reflection arrivals as pseudo-physical reflections because they exhibit kinematics coinciding with those of reflections in active-source reflection data, but the amplitudes and phases are not directly comparable (Loer et al., 2013)).Ghost reflections are retrieved primarily from the correlation of two primary reflections from two different depth levels.For example, the ghost reflection inside the reservoir (red lines in Fig. 1) can be retrieved by correlation of the primary reflection from the top and bottom of the reservoir (the blue and the purple dotted lines, respectively, in Fig. 1).These ghost reflections propagate only inside the reservoir without any kinematic effect of the overburden and underburden; they are equivalent to reflections by a virtual ghost receiver (the green triangle in Fig. 1) due to a virtual ghost source (the black star in Fig. 1).
If we substitute the response at  instead of the response at  in the right-hand side of Eq. 1,  acts as a collocated virtual source and receiver, which means we perform SI by AC and the result represents a zero-offset reflection trace at  .
In order to improve the accuracy of retrieving a specific ghost reflection, e.g., the ghost reflection inside the reservoir, we implement a process of muting observations before the reflection from the top of the reservoir and after the reflection from the bottom of the reservoir, which are used in the right-hand-side of Eq. 2 ( / ( ,  , )).As a result, the Green`s functions on the left-hand of Eq. 2 contains only the ghost reflection from inside the reservoir, which we refer to as  / ( , ) for simplification in further equations.We apply Eq. 2 for a baseline survey and a monitoring survey.So, a subscript "b" would denote the baseline survey, while a subscript "m" would denote the monitoring survey.If we apply Eq. 2 to all the active sources, a ghost zero-offset section is retrieved directly (the red lines of Fig. 1b) corresponding to both the baseline survey and the monitoring survey.Subsequently, by obtaining these zero-offset sections for both surveys, we determine the time difference ( ) using a correlation operator between the ghost reflection of the baseline survey and monitoring survey: The right-hand side of Eq. 3 represents multiplication between two terms in the frequency domain: * ( , ), which presents the ghost reflection retrieved from the baseline survey at the source (or receiver) location  , and  ( , ), which represents the ghost reflection retrieved from the monitoring survey at the source (or receiver) location  .The symbol (*) denotes the complex conjugate.Thus, we use an inverse Fourier transform of this correlation as described by ( −1 ), then determine the time difference by identifying the argument of the absolute maximum value within this correlation. (2) (3) Since our technique is being applied to a numerically modelled dataset, we can evaluate the retrieved time differences.First, we determine the expected time difference for a monitoring survey: Here,  is the calculated arrival time of the ghost reflection for the baseline survey and  is the calculated time of the ghost reflections for the monitoring survey.Assuming a constant thickness for the reservoir and homogeneous lateral velocity changes within the reservoir in our subsurface model for numerical modelling, the ghost reflections exhibit the same arrival time at all virtual receiver positions.Consequently, we can apply averaging of the retrieved time differences over the positions to simplify the validation process.To accomplish this, we compare the retrieved time difference from the ghost reflection ( in Eq. 3) with the expected time difference (∆ in Eq. 4).Given that the time difference inside the reservoir may exhibit varying scales across different monitoring survey, we compute the average relative time difference ( ) between the retrieved time differences and expected time difference: where  is the number of the virtual receivers.
Up to this point, we have been considering the virtual receivers individually.Given our assumption about the subsurface model, we can also calculate the time difference for the stacked ghost reflection.
This involves stacking all the ghost reflections retrieved from all virtual receivers using Eq.2: Here,  / represents the ghost reflection for each virtual receiver, and  represents the total number of virtual receivers.This calculation results in a single trace representing the ghost reflection inside the reservoir for both the baseline survey and the monitoring survey ( / ).
We apply the same procedure to calculate the time difference ( ) and the relative time difference ( ) using the stacked ghost reflection: (5) (4)

Results
In this section, we first show the retrieval of ghost reflections from SI for the baseline survey, utilising the Groningen subsurface model.We then proceed to demonstrate the retrieved ghost reflections for the monitoring surveys, with a specific emphasis on analysing the time differences between the baseline survey and the monitoring surveys.Furthermore, we delve deeper into understanding the influence of source and receiver configurations on successfully retrieving high-resolution ghost reflections.

Numerically modelled data using the Groningen subsurface model
Our objective is to retrieve and analyse ghost reflections from inside the Groningen reservoir to monitor changes occurring within it.One of the key factors contributing to these changes is the pressure depletion resulting from gas extraction.To investigate this phenomenon, we conducted a laboratory experiment to measure velocity variations caused by pressure depletion.
To simulate the conditions of the Groningen gas reservoir, we utilise a Red Felser sandstone cylindrical sample, known as an analogue to the Rotliegend sandstones found in the reservoir.The sample has a porosity of 19% and dimensions of 30 mm in diameter and 60 mm in height.Our experimental setup involves active-source ultrasonic transmission measurements, employing two S-wave transducers integrated into the pistons of the loading system.The transducers are positioned such that one serves as the source at the top while the other acts as the receiver at the bottom.By utilising a centre frequency of 1 MHz, we determine the S-wave velocities (a schematic illustration of the experimental setup can be found in Veltmeijer et al., 2022 andNaderloo et al., 2023).During the experiment, we systematically reduced the pore pressure from 30 MPa to 1 MPa, employing decrement steps of 5-10 bar.Throughout this process, the axial stress (65 MPa) and confining pressure (33 MPa) are kept constant to investigate the specific impact of pore-pressure depletion.These specific values were adopted from Spiers et al.Using the subsurface model, we generate a seismic reflection dataset employing a finite-difference modelling code (Thorbecke and Draganov, 2011) in a scalar mode.We use S-wave velocities in our numerical modelling because, in a 2D field survey, it is possible to use S-wave sources and horizontalcomponent receivers oriented in the direction perpendicular to the line.As a result, the recorded horizontal S-waves are completely decoupled from the P-and vertical S-waves.
Fig. 3 shows the location of the Groningen region and the S-wave velocity model for the baseline survey.The Groningen reservoir is located at a depth of around 3000 m, and we consider a constant thickness of 267 m in our modelling to be able to interpret the result as monitor velocity changes inside the reservoir and to avoid the velocity/thickness ambiguity.Additionally, in order to avoid using an excessively small grid size in our finite-difference modelling, we have taken into account a higher velocity for the North Sea formation (the top layer in the velocity model).This adjustment aims to reduce computational costs associated with the modelling process without having any conceptual influence on the results.The fixed receiver positions for our numerical modelling range from 3000 m to 7000 m (the blue triangles in Fig. 3b); the sources are placed from 2001.25 m to 8001.25 m at the surface (the red stars in Fig. 3b).The receivers and sources are regularly sampled with 1.25 m (the grid length) and 40 m spacing, respectively.We use a Ricker wavelet with a peak frequency of 20 Hz as a source wavelet.We also use an absorbing boundary at the surface to remove free-surface multiples in the numerically modelled data to better retrieve ghost reflections (Shirmohammadi et al., 2024).Fig. 4a shows the modelled common-source gather for a source at 5001.25 m for the subsurface model for 8 MPa pore pressure as a baseline survey.In this shot gather, we can see primary reflections from subsurface layers, including the reflections from the top and bottom of the reservoir (the blue and the purple arrows, respectively).We apply SI by AC by turning active sources into virtual receivers, which means we correlate each trace with itself for each common-receiver gather, and then we stack all correlated gathers along the receivers to retrieve a zero-offset section.Fig. 4b shows the results of SI by AC while all events in the common-receiver gathers are used.The result contains several retrieved ghost reflections from inside different layers of the subsurface model.They result from the correlation of all primaries and internal multiples in the source gathers.Note that all the events in the result of SI in our numerical modelling are ghost events because we used an absorbing boundary at the surface, and there are no surface-related multiples in the source gathers to create pseudo-physical reflections.
As indicated in Fig. 4b by the magenta arrows, it becomes challenging to distinguish specific ghost reflections from target layers, such as the ghost reflection from inside the reservoir.To improve the retrieval of ghost reflections propagating specifically inside the reservoir, which is our target layer, it is essential to correlate only the reflection from the top and the bottom of the target layer (Eq.2).
Therefore, we mute by manual picking all undesired events arriving earlier than the reflection from the top of the reservoir and later than the reflection from the bottom of the reservoir before applying SI by AC.In this condition, we only correlate the primary reflection from the top and the bottom of the reservoir (the blue and the purple arrows, respectively, in Fig. 4a).As a result, we only retrieve the ghost reflection from inside the reservoir.
Fig. 4c shows the retrieved result for the muted common-receiver gathers.The event indicated by the red arrow is the ghost reflection from inside the reservoir, which propagates only inside the reservoir.
The arrival time of the retrieved ghost reflection (0.27 s in Fig. 4c) corresponds to the travel time of a reflection inside the reservoir, which is equivalent to a reflection that would be recorded at a virtual ghost receiver placed directly at the top of the reservoir from a virtual ghost source at the same position.
(Because such direct recordings are impossible with our acquisition geometry, we call them virtual ghost receiver and virtual ghost source, respectively.)We aim to use the retrieved ghost reflections from inside the Groningen reservoir for monitoring.So, we continue with retrieving ghost reflections for monitoring surveys.

Time-lapse investigation using ghost reflections
For retrieving ghost reflections from the subsurface models for the monitoring surveys for pore pressures of 30, 20, 10, and 5 MPa, we apply a procedure similar to the one for the base survey.To better compare the results of SI, we extract one trace for a virtual receiver at 5001.25 m from the retrieved zero-offset sections for all different pore pressures (Fig. 5a) and a virtual receiver at 3801.25 (Fig. 5b) and the stacked ghost reflection using all virtual receivers using Eq.6 (Fig. 5c).
As we can see in Fig. 5, the retrieved ghost reflections show changes in time and amplitude for the different pore pressures.Our goal is to assess the time difference between the ghost reflections retrieved in the baseline survey and those in the monitoring surveys.It is challenging to extract the exact time differences in the ghost reflections by comparing them, as depicted in Fig. 5.To address this challenge, we employ a correlation technique between the ghost reflections obtained from the baseline survey and the corresponding ghost reflections from the monitoring surveys.By applying this correlation operator across all virtual receivers, we can determine the time at which the correlation yields its maximum amplitude (Eq.3).This time value directly corresponds to the time difference between the ghost  In Fig. 6(a-d), we present as the blue areas the retrieved time differences between the ghost reflections from the baseline survey and the monitoring survey for different pore pressure conditions: 30 MPa (Fig. 6a), 20 MPa (Fig. 6b), 10 MPa (Fig. 6c), and 5 MPa (Fig. 6d), considering all virtual receivers denoted by the blue highlighted region.The thick black line in this figure represents the expected time difference, calculated using Eq. 4.
As observed in Figs 6a and 6b, for monitoring surveys for pore pressures of 20 and 30 MPa, the time differences for all virtual receivers are more or less similar to the expected one, especially for those virtual receivers in the middle of the line of the sources.These virtual receivers are characterised by a complete illumination of the actual receivers from both sides.As also illustrated in Figs 5a and 5b, it becomes evident that the ghost reflections at the virtual receivers positioned in the middle present a better signal-to-noise ratio.So, we further zoom in on the virtual receivers positioned in the middle of the line, as shown in Fig. 6(e-h), depicting the time differences specifically for this subset of virtual receivers.For the monitoring surveys for pore pressure of 10 and 5 MPa (Figs 6c and 6d, respectively), the time difference is consistently zero, indicating that the ghost reflections cannot accurately predict time differences in these scenarios.respectively.The blue highlighted area shows the time difference using the original dataset, and the orange one shows the time difference using the subsampled dataset.The horizontal black lines show the calculated time differences.
To increase our ability to estimate time differences between ghost reflections, we employ time series subsampling, which, in this context, involves increasing the number of samples through the use of interpolation (Mikesell et al., 2015).Consequently, we subsample the retrieved ghost reflections before applying the correlation operator for calculating time differences.The time differences derived from the subsampled dataset are depicted in the highlighted orange regions in Fig. 6.We observe a significant improvement in the time-difference estimation, particularly for pore pressures of 10 MPa and 5 MPa (Figs 6c, 6d, 6g, and 6h), which were previously challenging to predict using the original dataset.
Furthermore, there is still a noticeable enhancement in the time-difference estimation for pore pressures of 30 MPa and 20 MPa (Figs 6a, 6b, 6e, and 6f).
As mentioned in the methodology section, we validate our retrieved time differences by computing the average relative time difference for the monitoring surveys using Eq. 5. Fig. 7 shows the average relative time differences using all virtual receivers (Fig. 7a) and specifically using selected virtual receivers positioned in the middle of the line (Fig. 7b) derived from both the original dataset (the blue circles) and the subsampled dataset (the purple triangles).Note that we have chosen to show the average relative time difference due to our assumption of no thickness change or lateral velocity change in our subsurface model.This procedure simplifies the validation process of our technique.
As we can see in Fig. 7, the relative time differences for pore pressures of 5 MPa and 10 MPa are both equal to 1 (the blue circles for the pore pressure of 5 MPa and 10 MPa).This indicates we are unable to estimate time differences for these monitoring surveys, given the very small changes inside the reservoir using the original dataset without subsampling.Contrary to this, after the application of subsampling, we can efficiently estimate time differences for these two pore pressures, approaching an average relative time difference close to zero.Conversely, for pore pressures of 20 MPa and 30 MPa, the relative time differences for both the original dataset and the subsampled dataset remain relatively consistent and close to zero.
Fig. 7b, illustrating the selected virtual receivers, exhibits a similar trend to Fig. 7a for both the original dataset and the subsampled dataset (the blue circles and purple triangles, respectively).However, the relative time difference for the selected virtual receivers approaches zero, signifying that they provide better estimates compared to using all virtual receivers.This improvement in results for the virtual receivers positioned in the middle of the source line is because of the complete illumination from both sides.
To mitigate the influence of virtual receivers positioned at the sides, we apply stacking as introduced in Eq. 6.The stacked ghost reflections are shown in Fig. 5c.Subsequently, we calculate the relative time difference for both the original dataset (the black plus symbols in Fig. 7) and the subsampled (the grey plus symbols) using Eq. 7. Comparing these relative time differences, it becomes evident that the stacked dataset yields superior results when all virtual receivers are utilised, while there is no significant distinction between the relative time difference of stacked ghost reflections and the average of the time differences of the individual virtual receivers when we use the selected virtual receivers.This implies that stacking effectively mitigates the influence of the virtual receivers positioned at the sides.
In conclusion, the subsampled dataset yields the most favourable outcome for all monitoring surveys (the grey plus symbols in Fig. 7).Consequently, we opt to present the retrieved time differences from the stacked subsampled dataset as the final time difference for different pore pressure conditions of our monitoring surveys.
Fig. 8 presents the time difference between the baseline and monitoring surveys, employing the stacked subsampled ghost reflections.In this context, the baseline survey corresponds to a pore pressure of 8 MPa, while the monitoring survey at 5 MPa can be interpreted as representing pore pressure depletion within the reservoir.On the other hand, the monitoring surveys at pore pressures of 10, 20, and 30 MPa can be viewed as examples of an injection mode within the reservoir.In both scenarios, we observe that the time difference remains detectable, even in the presence of small changes within the reservoir.It is worth noting that, as depicted in Fig. 2, we observe velocity variations in the reservoir ranging from 0.3% to 6%, and the retrieved time differences depicted in Fig. 8 exclusively reflect changes inside the reservoir that correspond to these velocity variations inside the reservoir.and the baseline survey (pore pressure of 8 MPa) using the stacked ghost reflections.
To recap the procedure for time-lapse monitoring using ghost reflections, we initially apply SI by AC using Eq. 2.Then, we apply subsampling to the retrieved ghost reflections, and we extract the time differences between the monitoring survey and the baseline survey using the correlation operator in Eq.
3 for each individual virtual source, same as shown in Fig. 6.Note that stacking and averaging are only used to simplify the validation of our technique (in Figs.7 and 8), which is allowed due to the constant velocity and uniform thickness within the target layer.However, it can be applied to models with homogeneous changes within the target layers or when utilizing very localized receivers, assuming that the changes inside the target layer are homogeneous.

Sources and receivers configuration for retrieval of ghost reflections
In the previous subsections, we looked at the result of SI by AC, where we turned actual sources into virtual receivers using SI.We opted for this approach because the spacing between the receivers was 1.25 m, which is shorter than the spacing between the sources.In this subsection, we will show the result of SI by turning receivers into virtual sources.Additionally, we demonstrate the effect of the source and receiver configuration, such as the number of sources and receivers, the spacing between them, and their geometry, in the retrieval of ghost reflections.This will offer valuable insights into the future practical applications of ghost reflections using field datasets.Note that in this subsection, we only look at the virtual zero-offset section for the baseline survey.
Fig. 9a shows the zoomed-in virtual zero-offset section resulting from SI by AC from Fig. 4c when turning the active sources into virtual receivers, which means we stack the correlated responses over all receivers with 1.25 m spacing.Fig. 9b shows the obtained virtual zero-offset section by turning receivers into virtual sources.In this situation, we stack the correlated responses for each receiver over all active sources with 40 m spacing.Table 1 shows all information regarding the spacing between virtual sources and receivers and spacing between stacked traces for the virtual zero-offset sections in Fig. 9.
Comparing Fig. 9a with Fig. 9b, we can see we have similar ghost reflections despite different numbers of traces in the zero-offset sections and different spacing in responses used for stacking in SI.This is an important finding, especially for the applications of ghost reflections using seismic field data when the number of receivers and sources is different.To better understand the effect of optimal spacing and the number of traces, we reduce the number of active sources and then apply SI by AC.
Fig. 9c shows the virtual zero-offset section when we use half of the sources as virtual receivers but the same number of receivers for stacking as in Fig. 9a.By comparing Fig. 9a with Fig. 9c, we can see that retrieved ghost reflections using optimal spacing between responses for stacking allow high-resolution interpretation of the target, even with the limited number of virtual receivers.Note that the optimal spacing is determined by half of the dominant wavelength.Fig. 9d shows a virtual zero-offset section using the same dataset as in Fig. 9c but now after turning the receivers into virtual sources.Comparing this result with the one in Fig. 9b, we can see that the retrieved ghost reflections allow interpretability with lower resolution compared to Fig. 9b, but we can still observe the desired ghost reflection.This means that if we have sub-optimal spacing between traces for stacking, it is necessary to have sufficient virtual sources.
Figure. 9.The result of SI by AC (a, c, and e) after turning sources into virtual receivers and (b, d and f) after turning receivers into virtual sources for different spacing between the virtual sources and receivers and different spacing between traces used for summation in SI, see Table 1 for exact numbers.Fig. 9e shows the virtual zero-offset section for virtual receivers retrieved using receivers ranging from 3000 m to 7000 m with spacing of 40 m and sources from 2001.25 m to 8001.25 m with spacing of 40 m.It is evident that the interpretation resolution of the ghost reflections is lowered compared to the one observed in Fig. 9a, particularly for virtual receivers at lateral distances before 4000 m and after 6000 m.However, for virtual receivers located in the central interval around 5000 m, we achieve the required interpretation resolution.This difference can be attributed to the fact that we utilise for stacking receivers ranging from 3000 to 7000 m, ensuring complete illumination from both sides of the virtual receivers (i.e., the sources to be turned into virtual receivers).It is crucial to note that, for the retrieval of pseudo-physical or ghost reflections, only a limited number of points fall into the stationary-phase region for a laterally homogenous medium and thus significantly contribute to the Green's function estimate.Consequently, for virtual receivers outside the receiver arrays, there is a higher chance for a limited number of receivers within the stationary-phase region and incomplete illumination (Draganov et al, 2012 andBalestrini et al., 2020).
Fig. 9f shows the virtual zero-offset section for virtual sources.We used sources from 2001.25 m to 8001.25 m for while the receivers are placed from 3000 m to 7000 with the same 40 m spacing.
Comparing this result to the one in Fig. 9b, we can see that we have achieved a comparable resolution, but due to the courser distance between the virtual sources, the section now appears as if spatially lowpass-filtered.
In summary, achieving higher interpretation resolution for ghost reflections requires a sufficient number of sources (or receivers) for stacking within the stationary-phase zone, a sufficient number of virtual sources (or virtual receivers), and a geometry of sources (or receives) that provides a complete illumination from both sides of the virtual receivers (or sources).

Discussion
We showed that ghost reflections are retrieved from SI by AC by using two primaries from two different depth levels from active sources and receivers at the surface.Notably, ghost reflections exhibit sensitivity exclusively to changes occurring within the target layer and remain kinematically unaffected by the overburden and underburden layers.This study provides further support for the existing studies regarding the use of ghost reflections in monitoring CO2 reservoirs through numerical modelling and laboratory experiments (Draganov et al., 2012), as well as in monitoring and characterising the shallow subsurface using numerical modelling and a field dataset (Shirmohammadi et al., 2024).
Although the results indicate that the application of ghost reflections for monitoring subsurface changes could be very practical, specifically for monitoring pore-pressure changes within the Groningen reservoir, several issues require a more in-depth discussion, especially for the future application of this technique across various subsurface models or with the utilization of field datasets.
First, our findings indicate that the ghost reflections exhibit sensitivity exclusively to velocity changes inside the reservoir due to pore-pressure changes.We showed that even small changes occurring within the reservoir can be detected, as ghost reflections are similar to isolated reflections, propagating only within the reservoir.To quantify the time differences between the baseline survey and the monitoring surveys, we employed the correlation operator between the ghost reflections retrieved from SI by AC.
We used subsampling before the correlation operator to detect small differences.Furthermore, as depicted in Fig. 5, it is evident that the maximum amplitudes of retrieved ghost reflections vary between the baseline survey and the monitoring surveys.The amplitudes can be relatively correct in SI; however, when specific conditions are not met, amplitude errors can be significant (Wapennar and Fokkema, 2006).Therefore, a more in-depth investigation is required to assess these amplitude variations.
Second, in our numerical-modelling approach, we considered a uniform reservoir thickness to highlight the sensitivity of the ghost reflections to the velocity changes.However, in reality, the ghost reflections will also be sensitive to other types of changes within the reservoir, such as thickness change, density change or a combination of these factors.It is important to note that these changes are not limited to being constant; they can also involve gradient changes, as demonstrated by Shirmohammadi et al. (2024).In these scenarios, we should use ghost reflection for each virtual receiver (or source) to investigate the changes, as illustrated in Fig 6.
Third, we used ghost reflections propagating exclusively inside the reservoir, which was our target layer.These specific ghost reflections can be retrieved more efficiently by muting undesired reflection arrivals earlier and later than the target arrivals, which requires those events to be clearly interpretable, especially in field data.However, we propose this technique for monitoring applications in which the target layers are already clearly identified from the baseline survey.It is worth mentioning that under certain cases, such as when dealing with a thin target layer, it becomes challenging to accurately separate the reflection from the top and bottom of the target layer before SI.In such cases, the retrieved ghost reflection would represent propagation inside multiple layers, including the thin target layer.
Thus, the vertical resolution of our methodology for retrieving the ghost reflections relies on the resolution of the reflection dataset as input.Using S-wave surveys increases the vertical resolution compared to P-wave surveys due to the shorter wavelength of the S-waves for the same frequencies.
Additionally, thin layers above or below the target layer, such as an anhydrite layer above the reservoir in our subsurface model, make it challenging to separate reflections from the top of the reservoir for far offsets.Nonetheless, this issue is mitigated to some extent through the stacking using the dataset from different offsets in SI (Eqs. 1 and 2).
To mitigate the presence of undesirable reflections that could overlap with the target ghost reflections in the result of SI, such as pseudo-physical reflections (i.e., reflections with physical kinematics), we used an absorbing boundary at the surface in our numerical modelling.The pseudo-physical reflections are the result of the correlation between primary reflections and their surface-related multiples.By using an absorbing boundary, we effectively prevented the occurrence of surface-related multiples and surface waves, thereby enhancing the quality of the SI results in numerical modelling.It is thus important to suppress surface waves as in, e.g., Balestrini et al. (2020) and surface-related multiples as in, e.g., Ghose and Goudswaard (2004) prior to applying SI when working with field datasets.These two preprocessing steps would ensure the accuracy and reliability of the SI outcomes by eliminating any potential confusion arising from retrieved pseudo-physical reflections but also interference from surface waves.
Finally, the interpretability resolution of the ghost reflections relies not only on the resolution of the input data but also on the source and receiver configuration.We explored two approaches: turning sources into virtual receivers and turning receivers into virtual sources.Determining which approach is more reliable depends on the spacing between the traces utilised for stacking.The optimal spacing is determined by half of the dominant wavelength.We demonstrated that it is possible to retrieve highresolution ghost reflections even with a limited number of virtual sources or receivers when employing the optimal spacing.However, if achieving the optimal spacing is not feasible, more traces are required to aid the interpretation, ideally at least two per wavelength.The geometry of the sources and receivers plays a crucial role in obtaining ghost-reflection results with high interpretability resolution.Note that for a laterally homogeneous medium, only a few sources positioned in line with the receivers remain stationary, which means rays from such source positions are nearly parallel, and interfere constructively in the summation (Draganov et al., 2012 andBalestrini et al., 2020).In the case of other media, determining the exact stationary-phase zone region is challenging.Therefore, it is recommended to maintain optimal spacing between receivers (or sources) and utilize an inline array with a two-sided uniform distribution of the sources (or receivers) with an extended length for sources (or receivers) on both sides to ensure the retrieval of high-resolution ghost reflections, at least for the virtual receivers (or sources) located in the middle of the array as illustrated in Figure 9b.

Conclusion
We focused on exploring the practical applicability of ghost reflections in monitoring subsurface activities, specifically in relation to pore-pressure changes, within the gas reservoir in Groningen, the Netherlands.Applying seismic interferometry (SI) to surface reflection data results in the appearance of ghost reflections in the retrieved virtual-source or virtual-receiver gathers.The ghost reflections result primarily from the correlation of primary reflections from different depth levels.
We determined the time differences between the ghost reflection retrieved from a baseline survey and several monitoring surveys using a correlation operator.The monitor surveys represented data after To detect minor changes effectively, we highlighted the importance of employing subsampled ghost reflections and considering virtual receivers with sources at the surface illuminating from both sides of the line.We demonstrated that the retrieved time differences between the ghost reflections exhibit variations corresponding to velocity changes within the reservoir.Notably, the retrieved ghost reflections are only sensitive to changes occurring in the target layer, effectively eliminating the kinematic influence of the overburden and underburden.
We also investigated the factors that contribute to obtaining ghost-reflection images with high interpretability resolution.We found that the geometry of the sources and receivers, the number of virtual sources and receivers, as well as the spacing between traces used for stacking for SI all play significant roles in ensuring clear ghost-reflection images.Additionally, by muting undesired reflection arrivals earlier and later than the target arrivals to correlate, we were able to enhance the clarity and robustness of the retrieved ghost reflections, which is crucial when working with field datasets.

O
Downloaded from https://academic.oup.com/gji/advance-article/doi/10.1093/gji/ggae099/7625094 by guest on 16 March 2024 Therefore, detecting temporal variations in subsurface characteristics is crucial for risk mitigation.Numerous time-lapse seismic studies have demonstrated the feasibility of detecting temporal variations in subsurface characteristics.

Figure. 1 .
Figure. 1. Schematic representation of seismic interferometry by (a) cross-correlation (CC) and (b) auto-correction (AC).The ghost reflections (red lines) are retrieved from the correlation of the primary reflections from the top and the bottom of the reservoir (the blue and the purple dotted lines, respectively) from the active sources (the red stars) and recorded at the receivers (the blue triangles) at the surface.The black stars and the green triangles indicate virtual ghost sources and receivers, respectively.

O
Downloaded from https://academic.oup.com/gji/advance-article/doi/10.1093/gji/ggae099/7625094 by guest on 16 March 2024 and represent the stress regime of the Groningen reservoir.Fig.2shows the measured velocity changes due to the pore pressure changes.We use the measured S-wave velocities for the pore pressures of 30, 20, 10, 8, and 5 MPa (the blue circles in Fig.2) for the Rotliegend reservoir in the Groningen subsurface model (derived from the Groningen Velocity Model 2017 by Nederlandse Aardolie Maatschappij (NAM, 2017), which is one of the best-known realistic velocity models for the Groningen subsurface).We consider a subsurface model with a pore pressure of 8 MPa as the baseline survey and the others as monitoring surveys.Our choice of 8 MPa for the baseline survey is dictated by the current pore pressure in the reservoir (NAM, 2021).

Figure. 2 .
Figure. 2. Measured velocity in the laboratory experiment as a function of changes in pore pressure.The blue circles indicate the chosen pore pressures for the baseline and monitoring surveys.

Figure
Figure. 3. (a) The location of the velocity model, (b) The Groningen velocity model and the geometry of the active sources (red stars) and the receivers (blue triangles) at the surface used in our numerical modelling.

OFigure
Figure. 4. (a) Common-source gather for an active shot at 5001.25 m.(b) The result of SI by AC when all events are used.(c) Same as (b) but for muted shot gathers.The blue and purple arrows indicate the reflections from the top and the bottom of the reservoir, the magenta arrows indicate the possible ghost reflections from inside the Rotliegend reservoir or other ghost reflections, and the red arrow shows the specific ghost reflection from inside the Rotliegend reservoir.

O
Downloaded from https://academic.oup.com/gji/advance-article/doi/10.1093/gji/ggae099/7625094 by guest on 16 March 2024 base and monitoring surveys.Note that these ghost reflections exclusively propagate within the reservoir, and, thus, the extracted time difference can be interpreted as a direct measure of the time difference within the reservoir itself.

Figure. 5
Figure. 5. (a) Retrieved zero-offset ghost reflections for a virtual receiver at 5001.25 for the baseline survey and the monitoring surveys; (b) same as (a) but for a virtual receiver at 3801.25 m.(c) Results of stacking the retrieved ghost reflections over virtual receivers.

Figure. 6 .
Figure. 6.The retrieved time differences in the reservoir using ghost reflections for the monitoring surveys with pore pressure of (a) 30 MPa, (b) 20 MPa, (c) 10 MPa, and (d) 5 MPa for all virtual receivers.(e), (f), (g), and (h) Zoomed in results from inside the black rectangles in (a), (b), (c), and (d),

Figure. 7 .
Figure. 7. The average relative time differences for the monitoring surveys (a) using all virtual receivers and (b) specifically for selected virtual receivers positioned in the middle of the line derived from the original dataset (the blue circles) and the subsampled dataset (the purple triangles).The plus symbols show the relative time differences using the stacked ghost reflections from the original dataset (black) and the subsampled dataset (grey).
The implications of our findings are notable, as ghost reflections can serve as valuable indicators for monitoring both near-surface structures and deeper formations such as fluid reservoirs or storage sites for H2 and CO2.

Table 1 .
Information about the source and receiver configuration used for retrieving the ghost reflection section in Fig.9.