Natural gas hydrate, a type of inclusion compound or clathrate, are composed of gas molecules trapped within a cage of water molecules. The presence of gas hydrate has been confirmed by core samples recovered from boreholes. Interests in the distribution of natural gas hydrate stem from its potential as a future energy source, geohazard to drilling activities and their possible impact on climate change. However the current geophysical investigations of gas hydrate reservoirs are still too limited to fully resolve the location and the total amount of gas hydrate due to its complex nature of distribution. The goal of this thesis is twofold, i.e., to model (1) the heterogeneous gas hydrate reservoirs and (2) seismic wave propagation in the presence of heterogeneities in order to address the fundamental questions: where are the location and occurrence of gas hydrate and how much is stored in the sediments. Seismic scattering studies predict that certain heterogeneity scales and velocity contrasts will generate strong scattering and wave mode conversion. Vertical Seismic Profile (VSP) techniques can be used to calibrate seismic characterization of gas hydrate expressions on surface seismograms. To further explore the potential of VSP in detecting the heterogeneities, a wave equation based approach for P- and S-wave separation is developed. Tests on synthetic data as well as applications to field data suggest alternative acquisition geometries for VSP to enable wave mode separation. A new reservoir modeling technique based on random medium theory is developed to construct heterogeneous multi-variable models that mimic heterogeneities of hydrate-bearing sediments at the level of detail provided by borehole logging data. Using this new technique, I modeled the density, and P- and S-wave velocities in combination with a modified Biot-Gassmann theory and provided a first order estimate of the in situ volume of gas hydrate near the Mallik 5L-38 borehole. Our results suggest a range of 528 to 768×10^6 m^3/km^2 of natural gas trapped within hydrate, nearly an order of magnitude lower than earlier estimates which excluded effects of small-scale heterogeneities. Further, the petrophysical models are combined with a 3-D Finite Difference method to study seismic attenuation. Thus a framework is built to further tune the models of gas hydrate reservoirs with constraints from well logs other disciplinary data.

Natural gas hydrate, a type of inclusion compound or clathrate, are composed of gas molecules trapped within a cage of water molecules. The presence of gas hydrate has been confirmed by core samples recovered from boreholes. Interests in the distribution of natural gas hydrate stem from its potential as a future energy source, geohazard to drilling activities and their possible impact on climate change. However the current geophysical investigations of gas hydrate reservoirs are still too limited to fully resolve the location and the total amount of gas hydrate due to its complex nature of distribution. The goal of this thesis is twofold, i.e., to model (1) the heterogeneous gas hydrate reservoirs and (2) seismic wave propagation in the presence of heterogeneities in order to address the fundamental questions: where are the location and occurrence of gas hydrate and how much is stored in the sediments. Seismic scattering studies predict that certain heterogeneity scales and velocity contrasts will generate strong scattering and wave mode conversion. Vertical Seismic Profile (VSP) techniques can be used to calibrate seismic characterization of gas hydrate expressions on surface seismograms. To further explore the potential of VSP in detecting the heterogeneities, a wave equation based approach for P- and S-wave separation is developed. Tests on synthetic data as well as applications to field data suggest alternative acquisition geometries for VSP to enable wave mode separation. A new reservoir modeling technique based on random medium theory is developed to construct heterogeneous multi-variable models that mimic heterogeneities of hydrate-bearing sediments at the level of detail provided by borehole logging data. Using this new technique, I modeled the density, and P- and S-wave velocities in combination with a modified Biot-Gassmann theory and provided a first order estimate of the in situ volume of gas hydrate near the Mallik 5L-38 borehole. Our results suggest a range of 528 to 768×10^6 m^3/km^2 of natural gas trapped within hydrate, nearly an order of magnitude lower than earlier estimates which excluded effects of small-scale heterogeneities. Further, the petrophysical models are combined with a 3-D Finite Difference method to study seismic attenuation. Thus a framework is built to further tune the models of gas hydrate reservoirs with constraints from well logs other disciplinary data.

Seismic imaging is a technique that uses seismic echoes to map and detect underground geological structures. The conventional seismic image has the resolution limit of λ/2, where λ is the wavelength associated with the seismic waves propagating in the subsurface. To exceed this resolution limit, this thesis develops a new imaging method using resonant multiples, which produces superresolution images with twice or even more the spatial resolution compared to the conventional primary reflection image. A resonant multiple is defined as a seismic reflection that revisits the same subsurface location along coincident reflection raypath. This reverberated raypath is the reason for superresolution imaging because it increases the differences in reflection times associated with subtle changes in the spatial location of the reflector. For the practical implementation of superresolution imaging, I develop a post-stack migration technique that first enhances the signal-to-noise ratios (SNRs) of resonant multiples by a moveout-correction stacking method, and then migrates the post-stacked resonant multiples with the associated Kirchhoff or wave-equation migration formula. I show with synthetic and field data examples that the first-order resonant multiple image has about twice the spatial resolution compared to the primary reflection image. Besides resolution, the correct estimate of the subsurface velocity is crucial for determining the correct depth of reflectors. Towards this goal, wave-equation migration velocity analysis (WEMVA) is an image-domain method which inverts for the velocity model that maximizes the similarity of common image gathers (CIGs). Conventional WEMVA based on subsurface-offset, angle domain or time-lag CIGs requires significant computational and memory resources because it computes higher dimensional migration images in the extended image domain. To mitigate this problem, I present a new WEMVA method using plane-wave CIGs. Plane-wave CIGs reduce the computational cost and memory storage because they are directly calculated from prestack plane-wave migration, and the number of plane waves is often much smaller than the number of shots. In the case of an inaccurate migration velocity, the moveout of plane-wave CIGs is automatically picked by a semblance analysis method, which is then linked to the migration velocity update by a connective function. Numerical tests on synthetic and field datasets validate the efficiency and effectiveness of this method.

This thesis focuses on the concept of Single Well Imaging (SWI) in which a seismic source and receivers are deployed in a borehole to investigate the surrounding geology. The Uniwell project (1997-1999) was the first attempt to develop the SWI method; it used a fluid-coupled downhole source, which unfortunately generated high amplitude guided waves in the borehole which obscured all other useful information. Initial research work detailed in this thesis focused on removing the high amplitude guided waves, known as tube waves. Two-step source signature deconvolution using first the recorded source signature, and then the tube-wave reflected from the bottom of the well, succeeded in compressing the tube wave. The results were not consistent across all receivers, but there is enough correlation to identify a P-wave. Further work concentrates on using a new technique called Empirical Mode Decomposition to separate the tube-wave mode from the data. This identifies three dominant modes and a possible body wave arrival, but the results are ambiguous due to the inability of the decomposition to focus on the narrow bandwidth of interest. The source signature deconvolution technique can also be used to process real-time vertical seismic profiling (VSP) data down-hole, during pauses in drilling, in what is referred to as a Seismic-While-Drilling (SWD) setup. Results show that the technique is versatile and robust, giving 1 ms precision on first-break picking even in very noisy data. I also apply the technique to normal VSP data to improve both the resolution and the signal-to-noise ratio. A major effort in this thesis is to consider the effect of a clamped downhole source to overcome the tube-wave problem, using a magnetostrictive source. Earlier work established that the use of a reaction mass tended to excite resonances in the tool which caused the transducer to break. A new design for the source was developed in cooperation with colleagues which utilises a hydraulic amplifier design and a low power coded waveform driving method exploiting the time-bandwidth product to extract the signal from the noise. My results show that as the run time increases the resolution improves. With a run length of 80s it is possible to resolve the signal transmitted 50 cm through a granite formation. This analysis led to a revised design of the source to improve its efficiency. I have used finite difference modelling, with a variable grid technique, to compare an ideal explosive source with an ideal clamped source. The fluid-coupled source emits high amplitude tube waves which virtually obscure the body wave, whereas the clamped source emits a clearly identifiable P-wave along with lower amplitude tube waves. This clearly illustrates the advantage of an ideal clamped source. To model the source more accurately the idealwavelet is replaced by the respective recorded source signatures, and the data is then processed by cross correlation with the appropriate signature. The results show that the coded waveform approaches the resolution of the ideal wavelet very well, with all major events being visible. However, the fluid-coupled source performs very poorly with only the highest amplitude tube-wave visible. This work illustrates that by replacing a fluid-coupled source by a clamped source driven by a coded waveform, and by processing the data using cross correlation or signature deconvolution, it is possible to minimise or eliminate tube-wave noise from a SWI survey. It is hoped that the results outlined here will provide the basis for a new SWI method than can be used to prolong the supply of North Sea oil.

Modern oil and gas exploration depends on a variety of geophysical prospect tools. One of them is reflection seismology that allows to obtain interwell information of sufficient resolution economically. This exploration method collects reflection seismic data on the surface of an area of prospect interest and then uses them to build seismic images of the subsurface. All imaging approaches can be divided into two groups: wave equation-based methods and integral schemes. Kirchhoff migration, which belongs to the second group, is an indispensable tool in seismic imaging due to its flexibility and relatively low computational cost. Unfortunately, the classic formulation of this method images only a part of the surface data, if so-called multipathing is present in it. That phenomenon occurs in complex geologic settings, such as subsalt areas, when seismic waves travel between a subsurface point and a surface location through more than one path. The quality of imaging with Kirchhoff migration in complex geological areas can be improved if multiple paths of ray propagation are included in the integral. Multiple arrivals can be naturally incorporated into the imaging operator if it is expressed as an integral over subsurface take-off angles. In this form, the migration operator involves escape functions that connect subsurface locations with surface seismic data values through escape traveltime and escape positions. These escape quantities are functions of phase space coordinates that are simply related to the subsurface reflection system. The angle-domain integral operator produces output scattering- and dip-angle image gathers, which represent a convenient domain for subsurface analysis. Escape functions for angle-domain imaging can be simply computed with initial-value ray tracing, a Lagrangian computational technique. However, the computational cost of such a bottom-up approach can be prohibitive in practice. The goal of this work was to construct a computationally efficient phase space imaging framework. I designed several approaches to computing escape functions directly in phase space for mapping surface seismic reflection data to the subsurface angle domain. Escape equations have been introduced previously to describe distribution of escape functions in the phase space. Initially, I employed these equations as a basis for building an Eulerian numerical scheme using finite-difference method in the 2-D case. I show its accuracy constraints and suggest a modification of the algorithm to overcome them. Next, I formulate a semi-Lagrangian approach to computing escape functions in 3-D. The second method relies on the fundamental property of continuity of these functions in the phase space. I define locally constrained escape functions and show that a global escape solution can be reconstructed from local solutions iteratively. I validate the accuracy of the proposed methods by imaging synthetic seismic data in several complex 2-D and 3-D models. I draw conclusions about efficiency by comparing the compute time of the imaging tests with the compute time of a well-optimized conventional initial-value ray tracing. / text

Seismic tomography is a promising tool to help understand and evaluate the stability of a rock mass in mining excavations. Lab measurements give evidence that velocities of seismic wave propagations increase in high stress areas of rock samples. It is well known that closing effects of cracks under compressive pressures tend to increase the effective elastic moduli of rocks. Tomography can map stress transfer and redistribution and further forecast rock burst potential and other seismic hazards, which are influenced by mining. Recorded by seismic networks in multiple underground mines, arrival time of seismic waves and locations of seismic events are used as sources of tomographic imaging survey. An initial velocity model is established according to properties of a rock mass, then velocity structure is reconstructed by velocity inversion to reflect the anomalies of the rock mass. Mining-induced seismicity and double-difference tomographic images of rock mass in mining areas are coupled to show how stress changes with microseismic activities. Especially, comparisons between velocity structures of different periods (before and after rock burst) are performed to analyze effects of rock burst on stress distribution. Tomographic results show that high velocity anomalies form in the vicinity of rock burst before the occurrence, and velocity subsequently experiences a significant drop after the occurrence of rock burst. In addition, regression analysis of travel time and distance indicates that the average velocity of all the monitored region appears to increase before rock burst and reduce after them. A reasonable explanation is that rock bursts tend to be triggered in highly stressed rock masses. After the energy release of rock bursts, stress relief is expected to exhibit within rock mass. Average velocity significantly decreases because of lower stresses and as a result of fractures in the rock mass that are generated by shaking-induced damage from nearby rock burst zones. Mining-induced microseismic rate is positively correlated with stress level. The fact that highly concentrated seismicity is more likely to be located in margins between high-velocity and low-velocity regions manifests that high seismic rates appear to be along with high stress in rock masses. Statistical analyses were performed on the aftershock sequence in order to generate an aftershock decay model to detect potential hazards and evaluate stability of aftershock activities. / Ph. D.

L’atténuation sismique est un paramètre physique très utile pour décrire et imager les propriétés du sous-sol, et tout particulièrement les roches saturées et les nuages de gaz. Les approches classiques analysent l’amplitude du spectre des données ou bien la distorsion de ce spectre, avec des méthodes asymptotiques. L’inversion des formes d’onde (Full Waveform Inversion en anglais, FWI) est une approche alternative qui prend en compte les aspects de fréquences finies. En pratique, à la fois les vitesses et l’atténuation doivent être déterminées. Il est connu que l’inversion multi-paramètre ne conduit pas à un résultat unique.Ce travail se focalise sur la détermination des vitesses et de l’atténuation. La dispersion liée à l’atténuation produit des modèles de vitesse équivalents en termes de cinématique. Je propose une inversion hybride : la « relation cinématique » est un moyen de guider l’inversion des formes d’onde non-linéaire. Elle se décompose en deux étapes. Dans un premier temps, l’information cinématique est remise à jour, et ensuite les vitesses et l’atténuation sont modifiées, pour une cinématique donnée. Différentes approches sont proposées et discutées au travers d’applications sur des données synthétiques 2D, en particulier sur les modèles Midlle-East et Marmousi. / Seismic attenuation is a useful physical parameter to describe and to image the properties of specific geological bodies, e.g., saturated rocks and gas clouds. Classical approaches consist of analyzing seismic spectrum amplitudes or spectrum distortions based on ray methods. Full waveform inversion is an alternative approach that takes into account the finite frequency aspect of seismic waves. In practice, both seismic velocities and attenuation have to be determined. It is known that the multi-parameter inversion suffers from cross-talks.This thesis focuses on retrieving velocity and attenuation. Attenuation dispersion leads to equivalent kinematic velocity models, as different combinations of velocity and attenuation have the same kinematic effects. I propose a hybrid inversion strategy: the kinematic relationship is a way to guide the non-linear full waveform inversion. The hybrid inversion strategy includes two steps. It first updates the kinematic velocity, and then retrieves the velocity and attenuation models for a fixed kinematic velocity. The different approaches are discussed through applications on 2D synthetic data sets, including the Midlle-East and Marmousi models.

Master of Science / Department of Geology / Brice LaCroix / Abdelmoneam Raef / South-central Kansas has experienced an increase in seismic activity within the Proterozoic basement over the past 10 years. In 2009, Oklahoma seismic stations recorded 50 earthquakes statewide, a 200% increase from 2008. Oklahoma Geological Survey (OGS) seismograph stations recorded 1,028 in 2010, an increase of over 2000% from 2009. Between 2000-2012, Kansas experienced only 12 earthquakes statewide. Beginning in September 2013, clusters of seismic events in south-central Kansas began to increase. In 2015 alone, Kansas seismograph stations recorded 448 earthquakes, of which 166 resulted in a magnitude 2.0 or greater. Since 2013, United States Geological Survey (USGS) seismograph stations have recorded over 12,000 earthquakes within Kansas and Oklahoma. Pore fluid pressure increases associated with recent high-rate wastewater injection into the dolomitic Arbuckle disposal zone are hypothesized as cause of reactivation of the faulted study region’s Proterozoic basement. Although the magnitude of fluid-pressure change required for reactivation of these faults is likely low given failure equilibrium conditions in the midcontinent, heterogeneities (i.e. permeability, porosity, fluid pressure) in the basement could allow for a range of fluid pressure changes associated with injection. This research aims to quantify the fluid pressure changes responsible for fault reactivation of the Proterozoic basement. To address this issue, we use 97 earthquake focal mechanisms and over 12,000 seismic events, from the USGS catalog, within an area encompassing ~ 4,000 km². Focal mechanism data was utilized to determine the regional stress field present within the study region. Nodal plane data extracted from the focal mechanisms was crucial to identifying lineaments within the underlying basement complex. A 3D seismic dataset covering the Wellington and Anson Bates Fields in north central Sumner County, Kansas was utilized for enhanced structural delineation of an interpreted faultnetwork affecting the Mississippian and Arbuckle Groups, to investigate whether it impacts the underlying granitic basement and its complex network of potentially interconnected fault planes. Smoothed similarity and spectral whitening analyses were applied to the dataset to improve depth of investigation and uncover fault lineaments masked by seismic attenuation due to increasing depth. An interpreted network of fault planes at depths of 3.5 km was uncovered beneath Wellington Field. The lineaments are well aligned with known structural features present within the Proterozoic basement, the Central Kansas Uplift and the Nemaha Ridge-Humboldt fault zone.

La construction du macro-modèle de vitesse est une étape cruciale de la chaîne d’imagerie sismique afin de produire un modèle adéquat pour la migration ou l'inversion de la formes d’onde complètes (Full Waveform Inversion). Parmi les approches possibles pour construire ce macro-modèle, la tomographie des pentes est fondée sur le pointé d’évènements localement cohérents caractérisés par leur temps de trajet et leurs pentes dans les collections de sismogrammes à source et réflecteur communs. Chaque évènement dans les observables est associé à un petit segment de réflecteur dans le sous-sol caractérisé par sa position et son pendage. Cette thèse propose une reformulation de la tomographie des pentes pour palier aux deux limitations sus-mentionnées. Les temps de trajet sont calculés à l’aide d’un solveur eikonal et la méthode de l’état adjoint est utilisée pour calculer efficacement le gradient de la fonction coût dans le contexte de méthodes d’optimisation locale. Le solveur eikonal est fondé sur une méthode aux différences finies pour la discrétisation des opérateurs différentiels. La méthode est implémentée pour des milieux 2D transverses isotropes avec un axe de symétrie dont l’orientation varie spatialement (milieux TTI). La méthode est évaluée avec différents examples synthétiques comme les cas du modèle isotrope complexe Marmousi et du modèle 2D TTI BP-salt. Dans le cas de milieux TTI, les couplages pouvant exister entre les paramètres de différente nature sont analysés avec un cas synthétique canonique. Cette thèse se conclut par une application à des données de sismique réflexion multitrace large bande (fournie par CGG) afin de reconstruire le modèle de vitesse vertical. / Velocity macro-model building is a crucial step in seismic imaging workflows as it provides the necessary background model for migration or full waveform inversion. Slope tomography, as a reliable alternative for conventional tomography, provides a tool to achieve this purpose where one picks the local coherent events rather than continuous events. These approaches are based on the slopes and traveltimes of the local coherent events which are tied to a reflecting/diffracting point (scatterer) in the subsurface. In this thesis, I introduce an anisotropic slope tomographic approach which aim at macro-model building for subsurface properties in 2D tilted transversely isotropic (TTI) media. In this method, I reformulate the stereotomography, as an slope tomographic tool, such that I replace the ray-based forward engine with a TTI eikonal solver and take advantage of the adjoint state method to calculate the gradients. In result, I can efficiently calculate the traveltimes for complex media and long offset acquisition on a regular grid of the subsurface and formulate a matrix-free framework for the inversion. Different synthetic examples including the isotropic Marmousi and 2D TTI BP-salt model are considered to assess the potential of the method in subsurface parameter estimations. Also, through a simple example, the footprint of parameter cross-talk is investigated for the Thomsen parametrization. As a real data application, the proposed method is applied on a 2D marine BroadSeis data set (provided by CGG) to retrieve the vertical velocity of the subsurface.

Seismic inversion is the leading method to map and quantify changes in time-lapse (4D) seismic datasets, with applications ranging from monitoring hydrocarbon-producing fields to geological CO2 storage. However, the process of inverting seismic data for reservoir properties is a notoriously ill-posed inverse problem due to the band-limited and noisy nature of seismic data. This comes with additional challenges for 4D applications, given the inaccuracies in the repeatability of time-lapse acquisition surveys. Consequently, adding prior information to the inversion process in the form of properly crafted regularization terms is essential to obtain geologically meaningful subsurface models and 4D effects. In this thesis, I propose a joint inversion-segmentation algorithm for 4D seismic inversion, which integrates total variation and segmentation priors as a way to counteract the missing frequencies and noise present in 4D seismic data. I validate the algorithm with synthetic and field seismic datasets and benchmark it against state-of-the-art 4D inversion techniques. The proposed algorithm shows three main advantages: 1. it produces high-resolution baseline and monitor acoustic impedance models, 2. by leveraging similarities between multiple seismic datasets, the proposed algorithm mitigates the non-repeatable noise and better highlights the real seismic time-lapse changes, and 3. it simultaneously provides a volumetric classification of the acoustic impedance 4D difference model based on user-defined classes, i.e., percentages of seismic time-lapse changes. Such advantages may enable more robust stratigraphic/structural and quantitative 4D seismic interpretation and provide more accurate inputs for dynamic reservoir simulations.