Activities underground, such as gas extraction or fluid injection, can disturb the natural stresses present and can cause human-induced earthquakes along pre-existing faults. Even though they are related to engineering, these earthquakes are currently unpredictable. Monitoring and understanding how these earthquakes occur are essential for a safe use of the subsurface and to progress with mitigation measures and earthquake forecasting.
Current monitoring relies on post-failure seismic recordings, emphasizing the need for advancements in monitoring and forecasting techniques. Detecting stress changes before seismicity (pre-failure) occurs allows for the timely implementation of mitigation measures. Active seismic monitoring methods have the potential to detect stress changes early and as such precursory information that can improve the forecasting methods and models. However, there is still much to discover regarding the relationship between precursors and the underlying physics. In general, the common fault mechanisms during the seismic cycle are well known. Initial stress build-up is followed by first slip instabilities where the local stress exceeds the fault strength, leading up to a seismic event, during which stress on the fault is released. However, robust and reliable predicting of fault failure and the resulting earthquake has proven to be challenging even for reactivating experimental faults in a controlled laboratory setting....@en

Seismic imaging and monitoring with reflected waves, originally used in the oil and gas industry to identify and assess potential hydrocarbon reservoirs and later monitor their exploitation, also have diverse applications in near-surface geophysics, mineral exploration, geothermal energy, and CO2 or H2 storage. Beyond revealing subsurface structures, these techniques enhance our understanding of how the subsurface responds to human activities, such as induced seismicity due to extraction processes. Seismic imaging and monitoring often focus on specific target layers within the subsurface, but challenges from interferences with surrounding layers and small changes within the specific layer(s) can distort the specific signals and lower the accuracy. Our research aims to address these challenges and provide practical solutions for more accurate and reliable seismic imaging and monitoring with reflected waves.
In this thesis, we aim to develop seismic data-driven methods for layer-specific imaging and monitoring, with a primary focus on advancing the technique of ghost-reflection retrieval using seismic interferometry (SI) and showing how the Marchenko method could be used with land seismic data.
SI often involves the cross-correlation of seismic observations at different receiver locations and the consecutive summation over the available sources, allowing the retrieval of new seismic responses from virtual sources located at the receiver positions. When using sources and receivers only at the surface, the virtual-source gathers retrieved by SI contain not only pseudo-physical reflections but also ghost (non-physical) reflections. These ghost reflections result mainly from the cross-correlation (CC) or auto-correlation (AC) of primary reflections from two different depths, representing reflections from inside specific subsurface layer(s), as measured with a virtual ghost source and a virtual ghost receiver positioned directly on top of the specific layer(s). Consequently, the ghost reflections can provide information about the specific layer(s) without the effects of the overburden and underburden layers.
We first explore the use of ghost reflections for layer-specific characterisation of the shallow subsurface using SI by AC, utilising numerically modelled data for a layered subsurface model down to 30 m depth, incorporating a lateral change in velocity, a velocity gradient with depth, a thickness change, and a velocity change in the target layer. Additionally, we present the first application of ghost reflections to shallow subsurface field data. Ghost reflections typically exhibit similar characteristics to other reflection events, appear close to or interfere with other events with only slight temporal differences. This makes their identification a significant challenge. To address this, we eliminate surface-related multiples and demonstrate how specific ghost reflections can be more efficiently retrieved by muting undesired reflections in the dataset before applying SI.
To extend the application of ghost reflections to deep structures, we focus on the feasibility of monitoring pore-pressure changes in the Groningen gas field in the Netherlands. We utilise numerical modelling to simulate scalar reflection data, deploying sources and receivers at the surface. We conduct an ultrasonic transmission laboratory experiment to measure S-wave velocities at different pore pressures. This data is used to create subsurface models, which are then utilized to simulate scalar reflection seismic data for monitoring purposes. We retrieve zero-offset ghost reflections by applying SI by AC to the modelled datasets. We then use a correlation operator to determine time differences between a baseline survey and monitoring surveys. Additionally, we investigate the effects of the sources and receivers' geometry and spacing, as well as the number of virtual sources and receivers, on retrieving ghost reflections with high interpretability and resolution. Besides observing time shifts in the ghost reflections, we also explore the feasibility of using the amplitude of ghost reflections for reservoir monitoring.
Having clear reflections from both the top and bottom of the specific layer(s) is crucial for retrieving ghost reflections, which can be challenging when using land seismic datasets due to the usual presence of strong surface waves. Conventionally, surface waves are suppressed during data processing using frequency-offset, frequency-wavenumber, or bandpass filters. However, these approaches can prove ineffective when the surface waves are scattered and/or overlap with the frequency regions of the reflected body waves that we intend to preserve. To overcome some of these challenges, we show the efficacy of the interferometric surface-wave suppression using a 2D seismic reflection dataset from Scheemda, Groningen province, the Netherlands. Interferometric surface-wave suppression can be used to effectively suppress surface waves by applying SI to first estimate the surface waves and second followed by their adaptive subtraction from the original data. We propose to apply these two steps recursively, i.e., several times, which yields better results than a single application in terms of clearer and more continuous reflections. This technique can function as a standalone technique or as part of a pre-processing flow.
When applying the seismic reflection method for monitoring purposes, specific reflections, e.g., from the top and bottom of the reservoir, are of interest. The reflections from both the top and bottom of the specific layer(s) can also be distorted by other events from the surrounding layers. To eliminate such distortions, the Marchenko-redatuming method was introduced. Several Marchenko-redatuming methods have been applied successfully to marine field data. We demonstrate, for the first time, the application of the Marchenko-based isolation technique to field land seismic data to isolate the target response by removing the overburden and underburden. Land data are intrinsically elastic, known for dominant surface waves and a low signal-to-noise ratio, posing a challenge for the Marchenko method, which requires high-quality reflection data. After we carefully apply several pre-processing steps, including recursive interferometric surface-wave suppression, we apply the Marchenko method twice: first, to remove the overburden effects by choosing a focal depth of 30 m, and then to remove the underburden effects by choosing a focal depth of 270 m. This process generates a new reflection response from the target area, providing clearer subsurface responses. The Marchenko method is particularly beneficial for data-driven techniques such as ghost-reflection retrieval, seismic imaging, and time-lapse studies using land seismic datasets.@en

Reflection seismology aims to estimate the Earth's subsurface elastic parameters for further investigation by geologists and engineers. This involves generating elastic waves using seismic sources and recording the Earth's response with receivers. The subsurface model is typically considered a combination of a background model and a short-wavelength reflectivity model. There are two main paths to estimate these parameters: non-linear waveform inversion to directly compute the elastic parameters or depth migration to estimate a structural image or reflectivity of the subsurface.

Reverse-Time Migration (RTM) is a common depth migration technique that migrates recorded wavefields from the space-time domain to the space-depth domain. It utilizes the Born approximation and the adjoint of the Born operator to produce an RTM image. However, RTM can suffer from errors, such as noise, temporal and spatial limitations, and multiple reflections.

Least-Squares Reverse-Time Migration (LSRTM) is used to overcome some of these errors. LSRTM involves resolving the reflectivity model by least-squares inversion, which is computationally expensive. Gradient-based optimization algorithms are often employed to reduce the computational burden, but they still require solving the wave equation and its adjoint for a large model in multiple iterations. One way to reduce the computational cost is by limiting the computational domain to a target region of interest.

Target-oriented LSRTM, known as TOLSRTM, focuses on the wavefield just above the target by bypassing the overburden. This approach proves beneficial when the overburden generates strong internal multiple reflections that obscure the reflections from the target area. However, a redatuming method is required to predict all orders of multiples. Marchenko redatuming is a data-driven technique that predicts the Green's functions at the boundary of the target region, incorporating all orders of internal multiples. It allows for double-sided redatuming, considering both the source and receiver perspectives. By combining the LSRTM algorithm and Marchenko double-focusing, a target-oriented LSRTM algorithm is devised that can predict interactions between the target and overburden and remove the effects of the overburden in the image. Predicting these interactions results in an artifact-free image, a better convergence rate, and a high-resolution image of the target.

Target-oriented migration algorithms typically consider only the upper horizontal boundary of the region of interest (ROI), neglecting wavefields entering the ROI from the medium beneath the lower boundary. To address this, a target-enclosed LSRTM algorithm is proposed, including both the ROI's upper and lower boundaries. Including the lower boundary provides transmission information and can improve inversion convergence. In addition, this algorithm is adopted for virtual receivers created by Marchenko redatuming. In the case of physical receivers at the boundaries of the target zone, the target-enclosed algorithm can incorporate the transmission information emanating from the lower boundary to the upper one. Consequently, when the initial model is far from the actual model, the resulting image partly recovers the long wavelength part of the model in agreement with the Born approximation criteria. Moreover, when an initial model closer to the actual model is used, the algorithm can partially recover the vertical interfaces of the perturbation. In the case of virtual receivers at the boundaries of the target zone, since the Marchenko redatuming is performed in the initial background model, the redatumed wavefields at the lower boundary suffer from kinematic errors. Therefore, the algorithm can not recover the long wavelength part of the model.

The thesis concludes with a discussion of the results obtained from applying the algorithms to marine datasets. The images resulting from the Marchenko double-focusing based target-oriented LSRTM algorithm show improvements in both resolution and artifact reduction by suppressing the overburden generated internal multiple effects. Moreover, the double-focusing enables the user to reduce the computational costs of the LSRTM algorithm and choose finer spatial sampling for the image.

An appendix proposes a formulation for integrating the target-oriented algorithms with non-linear inversion like Full Waveform Inversion (FWI). The results of this proposed algorithm show its effectiveness by reducing the internal multiple related artifacts and increasing resolution and faster convergence.@en

Seismic imaging is often used to interpret subsurface formations. However, images obtained by conventional methods are contaminated with internal multiples. The Marchenko method provides the means to obtain multiple-free subsurface images. Due to the high computational cost of the conventional pointsource Marchenko imaging method, the less expensive plane wave Marchenko imaging method can be used to produce subsurface images along planes. This method can be repeated for different incident angles to produce images that account for the variable dip of the subsurface structures. In this abstract, we present the results of applying the plane wave Marchenko imaging method to a 2D marine dataset over the Vøring basin, the North Sea. The results show that, in comparison to the conventional plane-wave image, the plane-wave Marchenko imaging method successfully suppressed internal multiples, resulting in improvements in both the amplitude and continuity of the seismic events.

Wave-Equation-Based Amplitude Vs Offset (WEB-AVO ) inversion solves the full elastic wave equation both for properties and for the total wavefield. It is a non-linear inversion technique that accounts for multiple scattering and mode conversions inside the target interval. When prior geological information interpreted from well logs is incorporated, stochastic inversion can be performed by honouring Bayes' theorem for probability density functions. The posterior function is proportional to the product of the likelihood function and the prior probability density function.

The prior probability function is built from well logs and is a complex mixture of Gaussians that account for thicknesses, property values and their corresponding standard deviations.

The likelihood function is built from the maximum likelihood estimator, the result of the deterministic inversion, and from the Hessian derived from the inversion kernel, scaled by the variance of the noise in the data.

The present work proposes that the best estimate of the noise in the data can be extracted from the residual of the seismic-to-well match.
The inaccuracy of the method can be quantified by taking the second derivative of the posterior function at the Maximum a Posteriori estimate. The present work also proposes that an additional source of inaccuracy is the intrinsic uncertainty, or non-uniqueness, of the method. It can be estimated with the help of random starting models on a perfect data set (synthetic data).

The stochastic WEB-AVO inversion is a natural extension of the already existing deterministic WEB-AVO inversion workflow. The inversion result is constrained by the prior to honour the true geology observed in the wells.

Marine seismic acquisition is troubled with several factors of noise that can deteriorate seismic data. Marine seismic data are recorded with towed streamers that acquire the desired upgoing wavefields containing information of the geology beneath. The upgoing wavefield will travel past the streamer/hydrophone to be reflected at the sea surface and propagate back down to be collected as undesired downgoing (ghost) wavefields, before propagating further downward as a surface related multiple. Because up- and downgoing wavefields interfere, these ghost wavefields generate peaks and notches in the recorded amplitude spectrum that compromises the bandwidth, reducing the resolution and the interpretability of the seismic data. Two classes of ghosts exist - source and receiver - that can be removed (’deghosting’) through two main approaches; utilizing different acquisition strategies and/or computer based processing algorithms. Additional measurements may prove useful in acquiring broadband data but may be hampered by high costs and limited availability. Vast amounts of 2-D single streamer legacy data exist that can still benefit from enhanced deghosting techniques. Acquisition uncertainties, such as the unknown exact depth of sources and receivers, the unknown reflectivity of the free-surfaces, and the unknown propagation velocity of seismic waves in water, lead to an increased complexity in finding a solution to the deghosting problem. In this thesis, sensitivity analysis has shown that a variability in the propagation velocity of seismic waves in water has the greatest effect on the conventional 2D receiver deghosting result. The least sensitive parameter turned out to be the water surface reflectivity. A multitude of adaptive deghosting techniques - which optimize the parameter settings through data driven optimization - have been found in literature and are summarized to identify the shortcomings. A quantitative uncertainty analysis into these methods seems to be missing. A new adaptive deghosting method based on echo-deblending is introduced that incorporates an uncertainty analysis. Results look promising for determining receiver depths, however, the water surface reflectivities results are more challenging. Recommendations for future work to improve the method are given for future students pursuing the topic.

Marchenko Imaging is a new technology in geophysics which enables to retrieve Green's functions at any point in the subsurface having only reflection data. This method is based on the extension of the 1D Gelfand-Levitan-Marchenko equation to a 3D medium. One of the assumptions of the Marchenko method is that the medium is lossless. If the lossy reflection response is used in the Marchenko scheme, some artefacts in the Green's functions as well as in the seismic image are present. One way to circumvent this assumption is to find a compensation parameter for the lossy reflection series so that the lossless Marchenko scheme can be applied. The main tasks of this thesis are to: [1] use the Marchenko equation to estimate the attenuation in the subsurface, [2] find a compensation parameter for the lossy reflection series so that the lossless Marchenko scheme can be applied, and [3] to create an upscaling method for wave propagation. The Artefact Removal Method was created which makes it possible to calculate an effective temporal Q-factor of the medium between a virtual source in the subsurface and receivers at the surface. This method is based on the minimization of the artefacts produced by the lossless Marchenko scheme. The minimization was performed in three ways: [1] in the space-time domain, [2] in the frequency domain and [3] to the scales of the wavelet transform applied to the artefacts. This method can also be used to find the layers with high attenuation. The upscaling method which can be used to construct macro-scale homogenized viscoelastic properties of the medium from the micro-scale properties of a heterogeneous medium was developed. This is done through linking the macro- and micro- scale Lippmann-Schwinger equations which describe the wave field and the strain field scattering in an inhomogeneous medium, respectively. In this thesis, the macro-scale homogenized viscoelastic properties were calculated by using the T-matrix Approach and the Generalized Dvorkin-Mavko Attenuation Model. All theoretical results are supported by synthetic 1D modeling. The theoretical part of the thesis and the general work flow can be used for a very complex medium.

A seismic survey should be designed such that imaging of the acquired data leads to a sufficiently accurate subsurface image. For that purpose, methods for acquisition geometry analysis and design are available. These methods are used to judge whether an acquisition geometry is suited for the specified objectives. Conventional 2D and 3D acquisition geometry analysis methods are largely based on common-midpoint processing that assumes a horizontally layered earth model. Consequently, the influence of the inhomogeneity of the subsurface on the data quality and image quality is not taken into account. However, in practice it has been shown that the data and image quality can suffer significantly from complex geology. Therefore, the inhomogeneity of the subsurface must be taken into account in the acquisition analysis methods. This can be done by the use of a macro model of the subsurface. Additionally, recent developments in seismic imaging and reservoir characterization use the multiple-reflections in the data to extend the illumination in areas that cannot be reached by the primary-reflections. The use of multiples yields a better vertical resolution as well as to suppress migration artefacts caused by crosstalk of multiple-reflections. The seismic value chain suggests that the new developments in seismic imaging and reservoir characterization should lead to new opportunities in the seismic acquisition. Therefore, the goal of this research is the development of a method that meets all the above mentioned requirements. The method presented in this thesis is based on the previously developed $\it{focal}$ $\it{beam}$ $\it{analysis}$ concept. This concept emphasizes the separate analysis of the source geometry and the detector geometry, leading to two outputs: the focal source beam and the focal detector beam. This gives the opportunity to separately judge and adjust the configuration of the sources and the configuration of the detectors. These beams are inspected in the space-frequency domain and in the Radon-frequency domain. In the spatial domain, it is visible whether the wavefield has properly been focused into a point. In the Radon domain, it is visible how the angle-dependent amplitudes are affected by the acquisition geometry and overburden structures. This approach provides thorough understanding of the cause of image deficiencies. The source and detector beams can be multiplied to compute a migrated image. The multiplication can be carried out in two domains to assess the different quality parameters of a seismic image: • Multiplication of the source and detector beam in the space-frequency domain yields the resolution function, which represents the image of a unit point diffractor. • Multiplication in the Radon-frequency domain yields the AVP-function, which is the angle-dependent image of one reflection point on an angle-independent reflector. However, so far this method assumed primary reflections only as signal, leaving out multiple reflections as noise. In this thesis, I discussed that multiples can be considered as signal in the seismic imaging if they are handled correctly. Therefore, I extended this focal beam method to the multiples as well. In the extended focal beam method, the full wavefield propagation between a subsurface point and the acquisition surface is simulated using a recursive full wavefield modeling engine. It uses a macro velocity model for the wavefield extrapolation from one depth level to another depth level and a reflectivity model to include all the reflection and transmission properties related to the same depth level. Subsequently, the data are focused to the target point using the source and the detector geometry. This full wavefield data is a complex wavefield which includes the effects of all multiples (i.e., propagation as well as reflection and transmission effects). Therefore, the focusing is achieved by a minimization scheme. I have demonstrated with the help of some examples, that the gap in the Radon transformed focal source beam due to sparse sampling issues can be filled with the use of all multiples. It means more angle-dependent information can be obtained. Imaging of subsalt sediments is challenging in practice, because of the high velocity contrasts and irregular shapes. My analysis shows that in such situation, image quality varies strongly with the position of the target point with respect to salt. In that case, even if one has a perfect source distribution, primary illumination may be limited due to geology. The image quality cannot be improved any further by adjusting the acquisition geometry in this case. I show that utilizing all multiples may be the part of the solution. The illumination from below is important in such cases. The method presented in this thesis offers opportunities for the investigation of the added value of the surface-related and the internal multiples. To summarize, this new method meets all the following criteria: • Dealing with the complexity of the subsurface model. • Considering all the multiple reflections as useful information. • Illumination from below in the complex sub-salt scenarios. • Development of geophysical-based infill specifications to assess the impact of coverage holes on data quality@en