Induced seismicity, referring to earthquakes triggered by human activities, has become a pressing concern, particularly in regions with extensive subsurface resource extraction, such as the Groningen gas field in the Netherlands. Decades of natural gas production have altered regional subsurface stress and pressure, reactivating faults and leading to frequent earthquakes, which pose significant structural and societal concerns. To better understand the nucleation and mechanisms behind these earthquakes, source characterization has become a fundamental scientific approach. Accurate characterization not only provides insights into earthquake dynamics but also plays a crucial role in developing subsurface models, optimizing extraction techniques, and helping to mitigate future seismic risks. This research focuses on designing an efficient characterization method by developing a probabilistic workflow for estimating induced earthquake source parameters. Specifically, it integrates advanced probabilistic inversion techniques with geological knowledge to enhance accuracy while reducing computational costs. The proposed workflow leverages a variant of the Hamiltonian Monte Carlo (HMC) algorithms, which offers advantages over a traditional Markov Chain Monte Carlo (MCMC) algorithm, particularly in handling high-dimensional parameter spaces. The workflow estimates key earthquake source parameters, including the centroid, moment tensor, and origin time, even when prior information is limited or uncertain. Initial validation through a synthetic earthquake scenario demonstrates the workflow’s ability to recover source characteristics with high confidence. To further evaluate its validity, the workflow is first tested on a geomechanically simulated synthetic event before being applied to real earthquake data from Groningen. A geomechanical simulation of the 2018 M$_L$ 3.4 Zeerijp event evaluates whether a point source assumption sufficiently represents induced earthquakes compared to finite fault models. The results indicate that, despite simplifications, estimated earthquake parameters remain robust, suggesting that computationally expensive finite fault models may not always be necessary. Building on this, the workflow is then applied to real earthquake data from the 2019 M$_L$ 3.4 event below Westerwijtwerd village, where geological prior knowledge is incorporated to enhance computational efficiency. Using fault-based prior distributions instead of uniform priors significantly reduces computational time while maintaining accuracy, demonstrating the benefits of integrating geological constraints into probabilistic inversion. Expanding its applicability, the workflow is further employed to characterize ten additional earthquakes (M$_L > 2$) in Groningen, where it is combined with rupture directivity analysis to provide insights into rupture dynamics. The results show that normal faulting is the dominant mechanism, aligning with regional geological structures. While rupture directivity effects are minor, they exhibit small correlations with known fault orientations. The study concludes with recommendations for future improvements, such as incorporating S-wave analysis, which could enhance the accuracy of the estimated parameters. Additionally, a more rigorous treatment of uncertainty in both observed and synthetic seismograms is suggested to improve probabilistic estimates. All in all, this research contributes to the advancement of seismic inversion techniques, offering valuable tools for assessing and mitigating the risks associated with induced seismicity. The proposed workflow provides a framework applicable beyond Groningen, aiding in the characterization of induced earthquakes in other regions affected by subsurface activities.

Seismic survey design deals with determining the acquisition parameters that lead to the best possible imaging and characterization of the subsurface. The design of the survey is constrained by health, safety and environmental considerations and the available budget, seeking for a balance between quality and cost. Because seismic exploration is a widely used geophysical method for revealing underground resources, information about the subsurface is available in many areas. Therefore, it can potentially be used for purposes supplementary to exploration such as the monitoring of producing fields and fluids injection. However, the available budget for these purposes is usually lower than for exploration, and it becomes a priority to maximize the benefits derived from a potentially cheaper acquisition. In this thesis, we propose new methods for the analysis and design of seismic surveys that are based on previous knowledge from existing subsurface models and that aimto maximize image quality with the lowest acquisition efforts.

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....

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.

Monitoring time-lapse changes inside the subsurface is of great significance to many geotechnical applications, such as storage of gasses in underground geological formations. Minute differences in the seismic wavefield between an initial baseline and a subsequent monitor survey have to be detected in order to observe fluid flow inside subsurface reservoirs. This problem becomes even more challenging when the reservoir is situated underneath a series of complex, highly reflective layers. Such an overburden will generate strong multiple reflections that will interfere with the reflections of the target zone. Ideally, a methodology is designed in order to remove these internal multiples to allow a clear view of the reservoir response for time-lapse analysis. The Marchenko method can redatumthe seismic wavefield to arbitrary depth levels or points in the subsurface, while accounting for all orders of internal multiple reflections. This method, therefore, has great potential to solve some of the time-lapse issues, as it is able to closely examine specific zones of interest in the subsurface without distortions from surrounding layers. Time-lapse studies are often hampered by irregular or imperfect sampling, whereas the Marchenko method relies on densely sampled, co-located sources and receivers. It is, therefore, important that the Marchenko method is able to handle more complex acquisition geometries. This can either be achieved by interpolating the reflection data as a pre-processing step or by correcting for errors inside the Marchenko scheme. Here, point-spread functions are introduced that describe the imperfections in the reflection data. These imperfections distort the focusing and Green’s functions retrieved from the Marchenko method. Next, each iteration of theMarchenko scheme is extended to deblur the imperfect focusing and Green’s functions by multidimensional deconvolution with these point-spread functions. Additionally, a slight modification is required to ensure stability of the new scheme. This new iterative Marchenko scheme is computationally more expensive, but removes all sampling artifacts. Finally, the migrated images of the target zone show significant improvements, when using either the new scheme or interpolation as pre-processing step...

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.

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.

Monitoring seismic wavefields caused by induced seismicity in the subsurface is a difficult process. Ideally, it requires physical receivers in the subsurface, which is unpractical. Frequently, only measurements at the surface of the Earth are available, which give a limited amount of information about the subsurface. One way to improve the monitoring of the subsurface is through the use of virtual sources and receivers, which are not physically present but are created from the measured reflection data at the surface. This can be achieved through the use of the classical homogeneous Green's representation, however, this method requires two Green's functions measured on an enclosing boundary, which is an unrealistic requirement. Instead, a single-sided representation of the homogeneous Green's function can be used, where a focusing function, which is a wavefield that focuses from a single-sided boundary to a focal position in the subsurface without artifacts related to the internal multiples, is employed together with a Green's function. To obtain the Green's function and focusing function that are needed for this representation, the Marchenko method is used. This method employs reflection data, without free-surface multiples, at the surface of the Earth and an estimation of the first arrival, which can be modeled in a macro velocity model.

To test whether induced seismicity in the real subsurface can be monitored using the single-sided representation, synthetic data are first considered, which include a synthetic reflection response and macro velocity model. The Marchenko method is used in combination with these data to obtain the focusing functions and Green's functions that are required for the homogeneous Green's function representations. The classical representation and the single-sided representation of the homogeneous Green's function employ the Green's functions and focusing functions to obtain the homogeneous Green's function of the medium. The homogeneous Green's function is visualized by creating snapshots of the homogeneous Green's function and these snapshots are compared to a directly modeled reference wavefield. This demonstrates that the classical representation, when applied to data at an open acquisition boundary, yields significant artifacts in the results, while the single-sided representation obtains accurate results. It is also shown that the radiation pattern of a double-couple source can be included in the retrieval of the homogeneous Green's function. The synthetic reflection data are truncated by limiting the offsets and sampling distance and applying attenuation to simulate field conditions. These truncations show that the single-sided homogeneous Green's function contains artifacts and lacks physical events if the reflection data are not ideal. 2D field reflection data and a macro velocity model from the V\o ring basin are considered and pre-processed to account for these truncations. The classical and the single-sided homogeneous Green's function representation are both applied to the field data and the results show that the retrieval of the homogeneous Green's function is possible for 2D field data using point sources while employing the single-sided representation. The results of the classical representation contain a large amount of errors. It is also shown that a homogeneous Green's function can be retrieved that has a virtual source with a double-couple radiation pattern.

Next, the application of the single-sided representation is considered in greater detail. The representation is used to forecast a wavefield in the subsurface as well as to monitor a wavefield in the subsurface. For the monitoring of the wavefield, it is assumed that a physical source in the subsurface causes a wavefield which is measured at the surface of the Earth. The Marchenko method is used to create virtual receivers inside the subsurface, which are used in combination with the physical measurement in the single-sided representation. This is a one-step process, because the Marchenko method is only used to create the virtual receivers. The single-sided representation of the homogeneous Green's function requires the source wavelet to be symmetric in time, which is unlikely for physical sources. Hence, a different single-sided representation can be used, which retrieves the causal Green's function and does not require a symmetric source wavelet. The single-sided representation of the causal Green's function can retrieve a majority of the correct events, however, the results contain anti-symmetric artifacts when the physical source is located above the virtual receiver. To forecast a wavefield in the subsurface, given a specific source configuration, the single-sided representation of the homogeneous Green's function can be used. In this case, a two-step process is applied, where both the source and the receiver in the subsurface are created by the Marchenko method and are therefore both virtual. After the homogeneous Green's function is obtained, it can be convolved with a non-symmetric wavelet. To demonstrate the difference between the one-step monitoring process and the two-step forecasting process, 2D synthetic reflection data are utilized. For the source configuration, a rupture plane is considered, which is modeled by superposing and time-shifting point sources, which contain a double-couple radiation pattern and are all scaled differently to simulate the heterogeneity of the rupture plane. The total wavefield created by this rupture plane is monitored using the single-sided representation of the causal Green's function. There are anti-symmetric artifacts present in the result, related to each point source, however, the correct wavefield is retrieved above the shallowest source location and below this source location after the first arrivals of all sources. The single-sided representation of the homogeneous Green's function is applied to forecast a virtual rupture plane, by retrieving the homogeneous Green's function for each source separately. The retrieved homogeneous Green's functions are transformed to causal Green's functions, shifted in time and superposed to forecast the total wavefield, which is free of the anti-symmetric artifacts at any depth. Both the monitoring approach and the forecasting approach are tested on 2D field data and the retrieved wavefields show similar results as were seen when the synthetic data were used. When the total wavefield is forecasted, there are no anti-symmetric artifacts present and when the wavefield is monitored, there are artifacts, however, they are only present in part of the result, below the sources before and during the first arrival of each source.

To test the application of the single-sided representation in 3D, a 3D implementation of the Marchenko method is required. The implementation is straightforward from a theoretical standpoint, as the surface integrals are performed over two dimensions instead of just one. The practical implementation is more difficult, however. The Marchenko method requires that the reflection data are well sampled in both space and time for sources and receivers, hence, the 3D reflection data are of a large size. As a result, not only a large amount of storage space is required, but the loading time of the reflection data is high, both of which are unpractical for efficient computation. We limit these problems by pre-transforming the reflection data to the frequency domain and compressing the data using floating point arrays, which reduces the storage space and loading time. Two datasets are considered, one modeled in a simple four layer model and the other in a subsection of the complex 3D Overthrust model. For both models, a Green's function inside the medium is retrieved, using a first arrival in the Marchenko method that was modeled in the exact medium, and compared to a reference Green's function that was directly modeled. The results for both models are accurate for the single Green's function. Next, imaging is performed for the models, however, instead of modeling the first arrivals, they are estimated using an Eikonal solver, because the modeling time of all the first arrivals is too high. The results of the imaging using the Marchenko method are compared to the results of conventional imaging, which demonstrates that artifacts, related to the internal multiples, are attenuated.

The 3D implementation of the Marchenko method is used to retrieve the Green's functions and focusing functions in 3D using 3D synthetic reflection data modeled in the Overhtrust model. The classical homogeneous Green's function representation and the single-sided representation of the causal Green's function and the homogeneous Green's function are all applied using these data, for three different combinations of a virtual source and a virtual receiver. The results are compared to a directly modeled wavefield, which shows that the result obtained by using the classical representation is contaminated by artifacts and lacks physical events. The result of the single-sided representation of the causal Green's function contains anti-symmetric artifacts related to the focusing function when the virtual receiver is located below the virtual source. The result of the single-sided representation of the homogeneous Green's function shows a good match to the reference result. The single-sided representation of the homogeneous Green's function is also applied using an Eikonal solver to obtain the first arrival that is required for the Marchenko method. The homogeneous Green's function that is obtained in this way shows a small decrease in quality for the result, however, this approach is more computationally feasible. The single-sided representation is used in combination with the Eikonal solver to retrieve a large amount of virtual receivers, so that the propagation of the wavefield in the subsurface can be visualized in time through the use of snapshots. This reveals that the part of the wavefield that is traveling at angles that are close to the normal of the surface is retrieved properly, while the part of the wavefield that is traveling at greater angles to the normal is reconstructed with less accuracy. This lack of proper retrieval is caused by the limited aperture of the reflection data. A rupture plane in 3D is considered and constructed in a similar way as is done for the 2D synthetic data. Point sources are used to model wavefields, which are time-shifted and superposed, however, to further represent the heterogeneity of the rupture plane, each wavefield is modeled using an unique causal wavelet. Both monitoring, using the single-sided causal Green's function representation, and forecasting, using the single-sided homogeneous Green's function representation, are performed on the rupture plane configuration. The two-step forecasting approach yields accurate results, for a given distribution of sources. The one-step monitoring approach retrieves accurate results above the shallowest source location, however, the result contains artifacts at the locations below the shallowest source, before and during the first arrival of each source.

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.

Curiosity regarding what we cannot see has always driven research. Science has helped us to uncover many of those hidden secrets. In particular, geophysics has helped us to image the inside of the Earth. By sending a seismic signal into the Earth and recording the signal that comes back, geophysicists can characterize the layers of the subsurface. Nowadays, geophysics is used for many purposes, for example, the localization of fossil fuels, the characterization of the subsurface for the construction of wind farms and the evaluation of reservoirs for geothermal energy. In order to decrease the risk and cost involved in these activities, we need images of the subsurface that are as accurate as possible.

These images can only be obtained if we fully understand the propagation of the seismic signal in the subsurface. A long-standing problem in geophysical imaging is the presence of internal multiple reflections. When imaging the subsurface, we assume that the signal only reflects once when there is a contrast in velocity and/or density (for example, when changing from sand to rock). However, in reality, the signal can reflect many times inside the subsurface before being recorded at the surface. When treating the arrivals that have reflected many times as arrivals that have only reflected once, we incorrectly image the subsurface and create ghost reflectors that do not exist. This problem is particularly strong in geological settings that have a complex structure with many strong velocity and/or density contrasts above an area of interest. This may happen, for example, when there is a reservoir of oil below a thick stratified salt layer. In such cases, the image of the area of interest is unreliable due to the presence of many ghost reflectors. Therefore, we have to use knowledge of wave propagation to predict and attenuate the internal multiples in the data prior to imaging.

In this thesis, I further develop the data-driven and wave-equation-based Marchenko method to make it suitable for the attenuation of internal multiples in seismic field data. In addition, I evaluate the performance of suitable methods by applying them to field datasets recorded in different geological settings. I start this evaluation by demonstrating that what we call the conventional Marchenko method is perhaps not the most suitable Marchenko method for the application to field data. I develop an alternative Marchenko method instead: the adaptive double-focusing method. I show that this method indeed produces improved results compared to the conventional Marchenko method when applying it to a line of 2D data of the Santos Basin, Brazil.

Since the 2D results show promise, I continue with the extension to 3D applications. I first identify the key acquisition parameters that affect the result of our Marchenko method on 3D synthetic data and conclude that the limited crossline aperture and the coarse sail line spacing have the strongest effect on the quality of the result. Based on this evaluation, I interpolate the sail line spacing on 3D field data acquired in the Santos Basin and use the adaptive double-focusing method to predict and subtract internal multiples. I conclude that 3D Marchenko internal multiple attenuation seems to be sufficiently robust for the application to narrow azimuth streamer data in a deep marine setting, provided that there is sufficient aperture in the crossline direction and that the sail lines are interpolated. In addition, the adaptive double-focusing method is suitable for the attenuation of internal multiples generated by a complex overburden and for simultaneously redatuming to a level below this overburden.

Next, I modify the adaptive double-focusing method to obtain an adaptive double dereverberation method that is suitable when only aiming to attenuate internal multiples generated in an overburden without redatuming. Moreover, this method does not require a velocity model. I apply this method to a 2D line of data acquired in the very shallow Arabian Gulf. Also, I assess how to meet the data requirements for the Marchenko method in shallow water environments (e.g., the removal of surface-related multiples, the deconvolution of the source signature) and demonstrate that the state-of-the-art Robust Estimation of Primaries by Sparse Inversion (R-EPSI) method is capable of producing the correct input data for the Marchenko method in such settings.

Subsequently, I discuss the role of the adaptive filter in the application of the Marchenko method to field data. I argue that developments in seismic data processing allow us to predict internal multiples with more accuracy, such that only a conservative adaptive filter is needed to correct for the unavoidable minor amplitude and phase discrepancies between the internal multiples in the data and the predicted internal multiples. I demonstrate this by using a conservative adaptive filter to subtract internal multiples that were predicted by applying an adaptive Marchenko multiple elimination method to a 2D line of field data acquired in the Norwegian North Sea.

Finally, based on the results presented in this thesis, I conclude that the Marchenko method is an effective, data-driven and robust method for the prediction of internal multiples in marine seismic data. Different Marchenko methods are suitable for different purposes. There are two key elements for the successful application of a Marchenko method to field data: 1) the acquisition geometry needs to be sufficiently dense and 2) a careful processing workflow needs to be constructed that accounts for the specifics of the geological setting at hand, with significant emphasis on amplitude and phase preservation.

Mudstones play an important role in hydrocarbon exploration and production, carbon capture and storage, and nuclear waste disposal. The high concentration of clay minerals contribute to the high intrinsic anisotropy (e.g., velocity, strength, permeability, and resistivity changes with direction) of mudstones. This high anisotropy complicates, among other things, seismic interpretation for hydrocarbon exploration and production, as well as predictions on the mechanical behaviour of these clayrich rocks. Mudstones are also characterized by a low-permeability matrix, which makes it difficult for fluids to flow through the rock. This impermeable character of mudstones makes them a potential natural seal for long-term CO2 storage and a potential host rock for nuclear waste disposal. For hydrocarbon production, open fractures are needed to enhance the productivity of oil and gas reservoirs, whereas the presence of such fractures can result in unwanted leakage of CO2 or nuclear waste in the subsurface. Fracture formation depends on, among other things, the mechanical properties of the mudstone. It is thus important to understand the elastic anisotropy and mechanical properties of mudstones for successful hydrocarbon exploration and production, and to safely store CO2 and radioactive waste in the subsurface. Although mudstones are important in the energy sector, the understanding of their elastic anisotropy and deformation behaviour under various physical conditions is limited, due to their complex character and the lack of laboratory experiments performed on well-preserved samples.

The quality and business aspects are both of particular importance in determining the type of seismic acquisition. Usually, a strong emphasis on cost reduction is inevitable. On the other hand, there is an increasing demand for the acquisition of high-quality seismic data that can contribute to the various stages in the field development profile. These conflicting desires eventually make conventional seismic surveys an inadequate option. The application of blended acquisition along with efficient detector and source geometries is capable of providing high-quality seismic data in a cost-effective and productive manner. This way of data acquisition also contributes to minimizing health, safety and environment exposure in the field. Blended acquisition allows multiple source-wavefields to be overlapped in time, space, and temporal and spatial frequency, causing blending interference. The acquisition of less data via sparse detector and source geometries likely violates the Nyquist sampling criterion. Therefore, to make the aforementioned approach technically justifiable, deficiencies in recorded data have to be dealt with through the course of subsequent processing steps. One way to encourage this technique is to minimize any imperfection in processing algorithms. In addition, one may derive survey parameters that enable a further improvement in these processes, which is the primary focus in this thesis.

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.

This thesis investigates the potential of passive seismic methods that make use of body waves, and especially the passive reflection method, as cost-effective applications for multiscale subsurface imaging and characterization. For this purpose, we develop several seismic techniques for different scales: basin, crustal, and lithospheric. For the basin scale, we developed horizontal- and vertical-components spectral ratio of global earthquake phases to estimate the basin depth. We also used the Sp-wave method and analysis of the frequency-dependent quality factor to characterize the basin’s heterogeneities. The results show good agreement with active-seismic profiles. At the crustal scale, we investigated the application of seismic interferometry (SI). Comparison among different SI methodologies suggests that multidimensional deconvolution based on the truncated singular-value decomposition gives better structural imaging than do the conventional crosscorrelation or crosscoherence approaches, but also better than multidimensional deconvolution based on the damped least-squares scheme. This crustal-scale SI could be useful, for example, as a prescreening-exploration tool for deep geothermal reservoirs whose targets can be as deep as 10 km. At the lithospheric scale we studied not only the Earth, but also the Moon. For the Earth, we applied SI with global phases to obtain detailed images of aseismic parts of a subduction slab. Although the interpretation of the imaging results of the aseismic parts is not sufficiently decisive, the results suggest that the applied method is helpful for imaging aseismic parts of slabs. Furthermore, the radiation efficiency of intermediate-depth earthquakes is estimated to understand the source mechanism as a function of focal depth. The results indicate that there is a larger amount of non-radiated energy for intermediate-depth earthquakes. This implies one of the mechanisms for the slabs to be aseismic at certain depths. For the Moon, we applied SI to deep moonquakes to obtain reflection imaging of the lunar subsurface. With this application, the lunar Moho is interpreted to be around 50 km depth, indicating the potential usefulness of SI for other celestial bodies. Following the results obtained in this thesis, we conclude that the passive seismic methods with natural quakes have excellent potential usage in both the resource industry and academia.

In this thesis, I study coupled poroelastic waves and electromagnetic fields in layered media. The focus is two-fold:
1. Increase the theoretical and physical understanding of the seismo-electromagnetic phenomenon by analytically-based numerical modeling.
2. Investigate the potential of seismo-electromagnetic interferometry.

After presenting the governing equations that form the basis of the theoretical framework, I capture this system into a matrix-vector representation of the wave equation. I first use literature eigenvector sets, which I normalize with respect to power-flux. I then derive new, alternative power-flux normalized eigenvector sets that I prove to be numerically more stable and accurate. The eigenvector sets form the basis of the analytically-based numerical modeling code `ESSEMOD' that I developed to model seismo-electromagnetic wave/field propagation/diffusion in layered-Earth media. The alternative eigenvector set models scenarios with no seismo-electromagnetic coupling correctly, where the literature eigenvector sets fail. In addition, the alternative set properly deals with scenarios where both small amplitude signals and large amplitude signals occur in the record, whereas the literature eigenvector sets result in noise levels masking the small events. The same holds for scenarios with a small seismo-electromagnetic coupling coefficient. I design an effective global reflection scheme that properly describes the primary and multiple reflections in the models. I implement the correct boundary conditions to account for scenarios with a free-surface, and also for scenarios containing fluid/porous medium/fluid transitions.
To transform all the seismo-electromagnetic source-receiver combinations in a numerically effective way back from the horizontal wavenumber-frequency domain to the space-frequency domain, I derive and implement explicit Fourier-Bessel transformations.

I then validate the developed modeling code in numerous ways. First of all, I compare the results of seismo-EM layer-code modeling in a homogeneous medium with explicit homogeneous space Green's function expressions. This comparison provides a clear validation that the layer-code models the dynamic responses in homogeneous scenarios correctly. Next, I check numerical consistency by carrying out reciprocity checks. I study homogeneous space models, models containing a free-surface and models with interfaces.
As a next step, I validate the modeling results of seismo-EM layer-code modeling for typical seismo-electromagnetic laboratory configurations, i.e. models containing fluid/porous medium/fluid transitions. I first compare the purely electromagnetic part of the seismo-EM layer-code with an independently developed purely electromagnetic layered-Earth code. The results match perfectly in both phase and amplitude for full transmission and pure reflection experiments, as well as for a combination of both. I then carry out a seismo-electromagnetic reciprocity test for a fluid halfspace overlying a porous medium halfspace, proving that the coupled poroelastic and electromagnetic fields are modeled consistently and yield the expected results.
As a final validation step, I compare ESSEMOD with an independently developed seismo-electromagnetic layered-Earth modeling code. The results display an almost perfect match in both phase and relative amplitudes, and a constant amplitude correction factor of 4 needs to be applied to let the absolute amplitudes match.

I then carry out a small feasibility test to study the potential of the seismo-electromagnetic effect for exploration purposes. I investigate different source-receiver combinations for the same model, and focus on the signal strength recorded at different distances from the target depth level. I conclude that for the source-receiver combinations studied, the electric field due to a volume injection monopole source, as well as the magnetic field due to a seismic bulk force source, yield the strongest converted signals. The receiver-distance from the target of interest plays an important role in the signal measurability. The closer the receivers to the target, the higher the signal strengths. However, when the receivers are located too close to the target, the coseismic reflected fields can mask the interface response fields that we are mainly interested in.

Next, I study if nature itself can help us to overcome the very low signal-to-noise ratio of seismo-electromagnetic converted fields, by investigating the effects of thin-bed geological structures on the seismo-electromagnetic signal.
To investigate the effects of bed-thinning on the seismo-electromagnetic interference patterns, I numerically
simulate seismo-electromagnetic wave propagation through horizontally layered media with different amounts and thicknesses of thin-beds. I demonstrate seismo-electromagnetic sensitivity to changes in medium parameters on a spatial scale much smaller than the seismic resolution. By simulating moving oil/water contacts during
production, where the oil layer is gradually being thinned, seismo-electromagnetic signals are proven very sensitive to oil/water contacts.

I now explore the application of interferometric techniques to the seismo-electromagnetic system, which might eventually lead to an improved signal-to-noise ratio of the weak converted fields.
I derive the theory for interferometric retrieval of 2D SH-TE seismo-electromagnetic Green's functions.
Using both a circular source configuration and a line source configuration, I show that it is possible to correctly retrieve the dynamic seismo-electromagnetic 2D SH-TE response in a homogeneous medium, using seismic boundary sources only. Using seismo-EM layer-code data, I then show that it is also possible to correctly retrieve the direct shear wave-related causal coseismic field in a homogeneous medium, in both phase and amplitude. To obtain a perfect match in absolute amplitudes, I apply a single linear scaling factor. I finally carry out interferometric experiments in a model containing a single interface at 800 m depth, proving that it is possible to correctly
retrieve all 2D SH-TE causal seismic-related direct and reflected coseismic fields, as well as interface response fields, by cross-correlation interferometry, using seismic boundary sources only.
These results are promising for the application of 3D seismo-electromagnetic interferometry using seismo-EM layer-code modeling, and later on, in the field.

Next, I present an alternative way to effectively decompose fields into their up- and downgoing components and different field types, using recordings at multiple depth levels. I present the theory of this MDL decomposition scheme, followed by successful decomposition of synthetic elastodynamic data sets. I additionally study the implications of laterally-varying media on the horizontal wavenumber-frequency domain MDL decomposition scheme.
I demonstrate successful decomposition, using an acoustic approximation and applying a combined multi-component / MDL decomposition approach, of a field data set recorded in Annerveen, in the North of the Netherlands. I address how to effectively use the MDL decomposition scheme in a unified fashion, applied to all wave phenomena including seismo-electromagnetic phenomena.

I then make a step towards seismo-electromagnetic inversion, presenting an effective way to carry out a seismo-electromagnetic sensitivity analysis using resolution functions. I start by explaining the theory of resolution functions using a seismo-electromagnetic example. I define the seismo-electromagnetic resolution function for inversion for a bulk density perturbation. I demonstrate the effectiveness of this method by first carrying out a purely electromagnetic sensitivity analysis for a point perturbation in conductivity, located in an isotropic homogeneous half-space. These results are compared with literature results based on analytical homogeneous space Green's function expressions. The result using the seismo-EM layer-code is nearly identical to the literature result. The position of the scatterer is correctly resolved. At the end of this section, I present the results of the fully-coupled seismo-electromagnetic senstivity analysis for a bulk density contrast for a specific source-receiver combination, using single-frequency multi-component line data. I show that the coupled seismo-electromagnetic system is sensitive to a perturbation in bulk density and that the position of the perturbation can be correctly recovered.

I finalize this thesis by discussing potential seismo-electromagnetic applications, as well as by providing a brief outlook for future research.