Joint analysis of the surface and borehole gravity and gravity gradiometry data for hydrocabrbon (HC) reservoir monitoring and carbon dioxide (CO2) sequestration

Gravity and gravity gradiometry surveys have been widely used in mining and petroleum exploration. Interest in borehole gravity measurements has grown because they can help to detect deep targets for a relatively low cost. The best way to obtain a 3D density distribution is by a joint interpretation of the surface and borehole gravity data, which is a very challenging problem. The 3D inversion would be a method of choice for the quantitative interpretation of gravity and gravity gradiometry data. However, the rigorous inversion is a complicated and very time-consuming procedure, ant it is dependent on the a priori model and constraints used. In this dissertation, I show that a joint iterative migration of the surface and borehole gravity and gravity gradiometry data can be an effective instrument of imaging the subsurface density distribution. I present the results of the new development in the joint migration of the surface and borehole gravity and gravity gradiometry data. In this dissertation, I also consider an application of the developed approach to geophysical monitoring of carbon dioxide (CO2) injections in a deep reservoir, which has become an important component of carbon capture and storage (CCS) projects. Until recently, the seismic method was the dominant technique used for reservoir monitoring. The cost of seismic surveys, however, makes this method prohibitive in monitoring sequestration projects where there is not any direct profit available. Moreover, some environments present challenges for seismic acquisition, e.g., in urban areas. In this dissertation, I present a feasibility study of gravity gradiometry monitoring of CO2 iv sequestration in a deep reservoir using a novel approach involving both borehole and surface measurements. The interpretation of the observed data is based on joint iterative migration of the surface and borehole data. The advantage of this method is that the surface data provide a good estimate of the horizontal extent of the injection zone, while the borehole data control the depth of the target, which increases the sensitivity and resolution of the method. We illustrate the effectiveness of the gravity gradiometry method by computer simulating CO2 injection monitoring in the Kevin Dome sequestration site in Montana, USA.

JOINT ANALYSIS OF THE SURFACE AND BOREHOLE GRAVITY AND GRAVITY GRADIOMETRY DATA FOR HYDROCARBON (HC) RESERVOIR MONITORING AND CARBON DIOXIDE (CO2) SEQUESTRATION by Muran Han A dissertation submitted to the faculty of The University of Utah in partial fulfillment of the requirements for the degree of Doctor of Philosophy in Geophysics Department of Geology and Geophysics The University of Utah May 2020 Copyright © Muran Han 2020 All Rights Reserved The University of Utah Graduate School STATEMENT OF DISSERTATION APPROVAL The dissertation of Muran Han has been approved by the following supervisory committee members: Michael S. Zhdanov , Chair 12/1/2019 Erich U. Petersen , Member 12/1/2019 Keith Koper , Member 12/1/2019 Le Wan , Member 12/1/2019 Martin Čuma , Member 12/1/2019 and by the Department of Geology and Geophysics and by David B. Kieda, Dean of The Graduate School. ABSTRACT Gravity and gravity gradiometry surveys have been widely used in mining and petroleum exploration. Interest in borehole gravity measurements has grown because they can help to detect deep targets for a relatively low cost. The best way to obtain a 3D density distribution is by a joint interpretation of the surface and borehole gravity data, which is a very challenging problem. The 3D inversion would be a method of choice for the quantitative interpretation of gravity and gravity gradiometry data. However, the rigorous inversion is a complicated and very time-consuming procedure, ant it is dependent on the a priori model and constraints used. In this dissertation, I show that a joint iterative migration of the surface and borehole gravity and gravity gradiometry data can be an effective instrument of imaging the subsurface density distribution. I present the results of the new development in the joint migration of the surface and borehole gravity and gravity gradiometry data. In this dissertation, I also consider an application of the developed approach to geophysical monitoring of carbon dioxide (CO2) injections in a deep reservoir, which has become an important component of carbon capture and storage (CCS) projects. Until recently, the seismic method was the dominant technique used for reservoir monitoring. The cost of seismic surveys, however, makes this method prohibitive in monitoring sequestration projects where there is not any direct profit available. Moreover, some environments present challenges for seismic acquisition, e.g., in urban areas. In this dissertation, I present a feasibility study of gravity gradiometry monitoring of CO2 sequestration in a deep reservoir using a novel approach involving both borehole and surface measurements. The interpretation of the observed data is based on joint iterative migration of the surface and borehole data. The advantage of this method is that the surface data provide a good estimate of the horizontal extent of the injection zone, while the borehole data control the depth of the target, which increases the sensitivity and resolution of the method. We illustrate the effectiveness of the gravity gradiometry method by computer simulating CO2 injection monitoring in the Kevin Dome sequestration site in Montana, USA. iv To my parents and my wife, Xiao 25 TABLE OF CONTENTS ABSTRACT ....................................................................................................................... iii ACKNOWLEDGMENTS ............................................................................................... viii Chapters 1. INTRODUCTION .......................................................................................................... 1 2. HYDROCARBON RESERVOIR MONITORING ........................................................ 7 2.1 Geophysical monitoring .......................................................................................... 10 2.2 HC reservoir properties and monitored quantities ................................................... 13 2.3 Overview of the geophysical monitoring techniques .............................................. 14 2.4 Conclusions ............................................................................................................. 17 3. CARBON DIOXIDE (CO2) SEQUESTRATION ........................................................ 18 3.1 Carbon capture and storage (CCS) .......................................................................... 20 3.2 Geological reservoirs for CO2 sequestration ........................................................... 20 3.3 Trapping mechanism ............................................................................................... 23 3.4 Conclusions ............................................................................................................. 29 4. MODELING AND SENSITIVITY ANALYSIS OF THE SURFACE AND BOREHOLE GRAVITY AND GRAVITY GRADIOMETRY DATA ........................... 30 4.1 Forward modeling of gravity and gravity gradiometry data.................................... 30 4.2 Description of the synthetic model .......................................................................... 34 4.3 Forward modeling of observed data and comparison between the results with and without the noise...................................................................................................... 34 4.4 Conclusions ............................................................................................................. 48 5. GRAVITY AND GRAVITY GRADIOMETRY INVERSE PROBLEM.................... 49 5.1 Inversion of gravity and gravity gradiometry data .................................................. 49 5.1.1 Tikhonov regularization theory ...................................................................... 50 5.1.2 Regularization parameter ................................................................................ 51 5.1.3 Stabilizing functionals .................................................................................... 52 5.1.4 Reweighted regularized conjugate gradient (RRCG) method ........................ 57 5.2 Conclusions ............................................................................................................. 59 6. JOINT MIGRATION OF SURFACE AND BOREHOLE GRAVITY AND GRAVITY GRADIOMETRY DATA .............................................................................. 60 6.1 Gravity gradiometry data ......................................................................................... 60 6.2 Adjoint operators for gravity gradiometry inversion ............................................... 62 6.3 Migration of the surface gravity and gravity tensor fields and 3D density imaging .................................................................................................................... 64 6.4 Migration of the borehole gravity and gravity tensor fields and 3D density imaging .................................................................................................................... 67 6.5 Joint migration ......................................................................................................... 71 6.6 Iterative migration ................................................................................................... 72 6.7 Summary.................................................................................................................. 73 7. MODEL STUDY .......................................................................................................... 74 7.1 Model 1 .................................................................................................................... 74 7.1.1 Migration results for individual surface data .................................................. 75 7.1.2 Migration results of the borehole gravity data................................................ 83 7.1.3 Joint migration results of surface and borehole data ...................................... 92 7.1.4 Joint migration using multiple borehole data ................................................. 92 7.2 Model 2 .................................................................................................................... 96 7.3 Joint migration using multiple borehole data ........................................................ 104 7.4 Conclusions ........................................................................................................... 111 8. CASE STUDY ............................................................................................................ 112 8.1 Big sky carbon sequestration partnership (BSCSP) .............................................. 113 8.2 Kevin Dome project............................................................................................... 115 8.3 Kevin Dome geology ............................................................................................. 116 8.4 Kevin Dome model study ...................................................................................... 120 8.4.1 Kevin Dome no leakage model .................................................................... 122 8.4.2 Kevin Dome leakage model ......................................................................... 133 8.5 Conclusions ........................................................................................................... 139 9. CONCLUSIONS......................................................................................................... 144 REFERENCES ............................................................................................................... 146 vii ACKNOWLEDGMENTS Before I joined the Consortium for Electromagnetic Modeling and Inversion (CEMI) at the Department of Geology and Geophysics, University of Utah, I did not have much knowledge about electromagnetic methods, especially of forward modeling and inversion, which requires good math skills. There are too many people to thank. Therefore, I will be lengthy here and brief elsewhere. I am deeply indebted to my advisor and committee chair, Dr. Michael Zhdanov, who made this work possible and shared his vast knowledge of EM and inversion, for giving me a chance to work under his guidance and for ultimately leading me along an academic path. Dr. Le Wan not only sacrificed a lot of time to explain the problems of EM theory to me, but also made many invaluable insights into the application of my research in the real world and forced me to look beyond my work. Dr. Le Wan, who is an EM expert in my mind, provided helpful guidance for this work and answered numerous questions about programming and inversion techniques. I would also like to thank my committee member, Dr. Erich Petersen, who guided me to connect geophysics with geology, especially when it comes to minerals and rocks. I thank my committee member, Dr. Martin Čuma, who advised me to correct the parts of my dissertation. Finally, my family provided empathy and encouragement along the way. I especially thank my parents and parents-in-law. Without their encouragement and help, I could not viii have focused on my work completely. Most of all, I would like to thank my wife, Xiao, for supporting my study, which was a long journey from beginning to the end. ix 1 CHAPTER 1 INTRODUCTION Gravity gradiometry has become widely used in geophysical exploration since it can provide an independent measure of subsurface density distribution (Zhdanov et al., 2010b). The advantage of gravity gradiometry over other gravity methods is that the data are extremely sensitive to local density anomalies within regional geological formations (Wan & Zhdanov, 2008). High quality data can be acquired from either airborne or marine platforms over very large areas for a relatively low cost. The sensitivity of the gravity field is inversely proportional to the square of the distance. Making use of the borehole gravity measurement can significantly improve the inversion or migration results because the sensors in the borehole are located close to the target. The goal of this dissertation is to develop a method based on integration of the surface and borehole gravity and gravity gradient data to obtain better migration images of the subsurface. The borehole gravity method was pioneered by Smith in 1950 (Goodell & Fay, 1964; Smith, 1950) and then was applied to problems of reservoir evaluation by (McCulloh et al., 1968). Unlike the shallower-sensing density log, the borehole gravimeter is insensitive to wellbore conditions such as rugosity and the presence of casing. 2 The advantages of measuring gravity gradients rather than the gravity field have also been recognized (e.g., Nekut, 1989). The borehole gravity log results have also been reported in many papers. Gravitec (MacQueen, 2007) and Lockheed Martin have since developed the prototype of borehole gravity gradiometers. The first generation of borehole gravimeters (BHGMs) were limited to largediameter, near-vertical boreholes, and they were deployed in hydrocarbon wells. A secondgeneration BHGMs have been developed for mining and geotechnical applications (Seigel et al., 2009). This BHGMs probe can be operated on a standard 4- or 7- conductor wireline in NQ (drill rod that has 58 mm diameter) mining drill rods with an inclination from -300 to vertical and with ambient conditions limited to borehole temperatures and pressures typical in mining exploration. Sensitivity is better than 5 µGal (Gasperikova & Hoversten, 2006). Many researchers have tried to combine the borehole and the surface geophysical data to improve the results of interpretation. For example, Cao inverted the surface and borehole seismic data jointly to overcome the narrower frequency bandwidth defects of the surface seismic data (Cao et al., 2013). Li and Oldenburg inverted the surface and borehole magnetic data jointly to better define the deeper target (Li & Oldenburg, 2000). Krahenbuhl and Li and Sun and Li published the results of the joint inversion of surface and borehole gravity data (Krahenbuhl & Li, 2008; Sun & Li, 2010). Also, Rim and Li and Liu and Zhdanov conducted the research for a single borehole data imaging (Rim & Li, 2010; Zhdanov et al., 2011). Migration of gravity fields is a fast imaging tool to locate a target using a transformation of the observed data into a 3D density image. The results are the same as 3 those of 3D inversion; however, the numerical implementation and physical interpretation are different. Potential field migration theory was originally introduced by Zhdanov (Zhdanov, 2002), and subsequently developed for gravity gradiometry (Zhdanov et al., 2010a). It has been used in the interpretation of practical gravity and gravity gradient data as a fast imaging technique (Zhdanov et al., 2011). The migration can be done iteratively and can be applied in the same way as inversion for the gravity and gravity gradiometry data (Wan et al., 2016; Wan & Zhdanov, 2013). This dissertation demonstrates that the joint iterative migration of the surface and borehole gravity gradiometry data can provide better representation of the subsurface density distribution and improve the imaging of a deep target than the surface measurements alone. This dissertation also examines an application of the developed joint migration of the surface and borehole gravity gradiometry data for hydrocarbon (HC) reservoir monitoring and CO2 sequestration. With the industrial revolution, the average CO2 concentration level has increased from 280 ppm to reach its current level of more than 400 ppm (Huang et al., 2017).The increase in concentration level of CO2 and other greenhouse gases is the main reason for global warming (Bruant et al., 2002). A growing consensus that the global climate is changing has generated significant efforts in developing effective methods for carbon capture and storage (CCS). Many international research programs have been established in order to address this problem, e.g., the Australian government sponsored Cooperative Research Centre for Greenhouse Gas Technologies (CO2CRC), the Canadian and Saskatchewan governments sponsored Aquistore Program, and industry sponsored CO2 Capture Project (CCP). These programs intend to advance technologies that will underpin 4 the deployment of industrial-scale CCS. There are different types of suitable geological structures for CCS, like depleted hydrocarbon reservoirs, unminable coal seams, basalt formations, and deep saline aquifers (NETL, 2010). Most approaches currently proposed for CCS rely on storing CO2 in a supercritical state in deep saline reservoirs where buoyancy forces drive the injected CO2 upward in the aquifer until a seal is reached. The CO2 is stratigraphically and structurally trapped below an impermeable rock layer. The secondary mechanisms include the residual trapping of small amounts of CO2 in pore spaces as the supercritical fluid moves through the formation, whereby CO2 dissolves in existing formation fluids, becoming denser and sinking in the formation over time. The maximum storage security occurs through mineral trapping. CO2 dissolves in brine, forming a weak carbonic acid. Over time, this compound interacts with the minerals in the surrounding rock or with the minerals in the formation fluid to form solid carbonate minerals. The permanence of this type of sequestration depends entirely on the long-term geological integrity of the seal. The Big Sky Carbon Sequestration Partnership (BSCSP) is one of the seven partnerships initiated by the U.S. Department of Energy. BSCSP along with the other partnerships were created to develop the technologies, infrastructures, and regulations required to implement large-scale carbon dioxide (CO2) capture and storage within the nation. The BSCSP has chosen the Kevin Dome in north central Montana, among many potential CO2 storage sites. The flanks of the Kevin Dome are saline-saturated and have the potential to store more than 1.5 billion tons of CO2 (NETL, 2010). Government regulations require continuous monitoring of CO2 sequestration sites to ensure seal integrity. There is always a slight possibility of leakage even with a good 5 seal characterization. A leakage of CO2 to shallow water aquifers will alter the geochemistry, water quality, and ecosystem health. In addition, monitoring fluid movements within the reservoir will lead to informed reservoir management decisions (Bruant et al., 2002). In time-lapse reservoir monitoring, a baseline survey is conducted before the start of the project. Then, several time-lapse surveys are conducted periodically. The difference between the baseline survey and each future survey is calculated, which is associated with changes in the CO2 reservoir. Reservoir monitoring is dominated by seismic methods, which can provide unmatched resolution and accuracy. However, seismic monitoring is usually expensive and sometimes difficult or prohibitive to acquire, such as in urban or industrial areas. Therefore, an alternative or complimentary geophysical method is needed. The advancement in gravimeters and gravity gradiometers technologies has made these tools capable of monitoring the subtle changes in the surface gravitational field due to subsurface change in the fluid contents in the deep geological structures. Bradly et al. (2008) successfully monitored the movements associated with gas cap water injection (GCWI) at Prudhoe Bay using time-lapse gravity. AlJanobi (2017) examined a possibility of using the time-lapse airborne gravity and gravity gradiometry data for monitoring CO2 sequestration (AlJanobi, 2017). In this dissertation, we present feasibility studies of time-lapse gravity and gravity gradiometry monitoring of CO2 sequestration in the deep Duperow formation in the Kevin Dome project using the surface and borehole data. The interpretation is based on joint iterative migration imaging of the surface and borehole data (Zhdanov, 2002; Zhdanov et 6 al., 2010a). The advantage of this method is that the surface data can be used to estimate the horizontal extent of the injection zone, while the borehole data control the depth of the target and thereby increase the sensitivity and resolution of the method. 25 CHAPTER 2 HYDROCARBON RESERVOIR MONITORING Hydrocarbon (HC) reservoirs are composed of rocks containing various minerals, generally quartz, feldspar, rock fragments and clay minerals, fluids (water), and hydrocarbons (oil or natural gas). Natural reserves can be trapped in geological formations at various depths depending on the geological conditions in the reservoir. Each reservoir is characterized by its physical parameters, such as faulting and folding that contributed to the formation of an existing network. In unconventional reservoirs, hydrocarbons might be distributed throughout the reservoir at the basin scale and trapped in low permeability formations such as, “pitchouts,” where buoyant forces are not significant. Therefore, rock-fracturing techniques such as hydraulic stimulation (Hydraulic Fracturing) are needed to optimize the HC flow recovery. Typical unconventional reservoirs are formed by tight sand, shale, and sandstone. Hydraulic-Fracturing (HF) is a sophisticated technique, which is widely applied in low-permeability geological formations to enhance the production. It consists of breaking the rock by high-pressure fluid injection, creating an extensive network of fractures that allow for increased HC flow. HF is a component of the hydrocarbon production, along with the usual practice in oil-gas industry of injection of coproduced water (“wastewater disposal”) that can be reused, for instance, in reservoir depletion. 8 Hydrocarbon extraction in massive usage represents the most economically efficient means of energy production. However, despite the considerable economic benefit, a plausible, yet concerning, relationship between oil-gas industrial activity and serious geohazards is recognized, such as an increase in seismicity. In addition, episodes of freshwater contamination due to stray gas migration propagating in abandoned or improperly cased wells are reported in the literature. In principle, similar HF techniques have been applied also in Europe for a long time, yet in conventional reservoirs, and are likely to be intensified in future. When HF is used, especially in the form of the much-discussed “fracking,” knowledge of the state of the reservoir becomes important, both for optimizing operations, and to safeguard against potential environmental hazards. This suggests an increasing demand for technological development, including updating and adapting existing techniques in applied geophysics (Caffagni & Bokelmann, 2016). The first attempts of tracking temporal changes (monitoring) in hydrocarbon reservoirs are dated 40-50 years ago. These surveys were mainly focused on characterizing the subsurface using probes, such as seismic and electric (Figure 1), deployed in boreholes or at the surface. Reservoirs are nowadays monitored during exploration and exploitation of the hydrocarbons by using different monitoring techniques, according to geological conditions. Different techniques have been developed and are still in usage in offshore plays, such as seismic Permanent Reservoir Monitoring (PRM), with a different design from the onshore case. Some of these techniques aim at constructing images of the reservoir compartments; others can estimate important parameters directly in situ. 9 Figure 1. Detection geometries in geophysical monitoring of a hydrocarbon reservoir; the blue spot identifies the target zone(Caffagni & Bokelmann, 2016). 10 2.1 Geophysical monitoring Hydrocarbons are extracted from unconventional reservoirs by artificial stimulation, such as HF, which can dramatically affect the geomechanical response of the reservoir. Stress alteration and sudden increase in pressure can generate significant changes of the reservoir characteristics. Geophysical monitoring implies keeping track of temporal changes of characteristics of a reservoir that are imposed by external or internal perturbation of the reservoir state. Such changes may not only allow inferences to be made on the reaction of the reservoir to the perturbation (e.g., in terms of deformation, fluid flow, or temperature changes), but also give information on reservoir properties that were previously unmeasured. Monitoring can be achieved by implementing probes in situ (directly into the reservoir), or outside of the reservoir. The noninvasive nature of the method prevents any further perturbation of the reservoir. Monitoring can be achieved by implementing probes in-situ (directly into the reservoir), or outside of the reservoir, where the noninvasive nature of the method prevents any further perturbation of the reservoir. Consider a hydrocarbon reservoir analyzed in situ at time t0. If no external stimulation is applied, the reservoir remains unaltered, and retains its initial state. Natural variations can be due to the background stress in the underground rocks, forces of tectonic, volcanic, or tidal origin, presence of faulting/folding, fracturing, or presence of fluids, such as water and hydrocarbons. At time t1 an artificial stimulation starts, and it is terminated at time t2. The time interval between t1 and t2 refers to the stimulation; much of its effects should occur within this time interval, in addition to any occurring natural variations. After t2, the stimulation is stopped, however, significant effects may still occur, and the reservoir 11 will unlikely return to its exact initial condition at t0. After a certain time, changes may not be detected anymore; however, the reservoir characteristics will usually differ from those at t1. Monitoring refers to tracking the temporal changes of the reservoir, starting from its initial state at time t0, to a possibly long time after t2. A baseline study of the hydrocarbon reservoir should be conducted before the stimulation (e.g., fracking) to obtain the background state, which provides a reference point of comparison for the reservoir state as measured during and after the stimulation. In this way, both the natural and artificial components of the deformation processes can be tracked. Where possible, these two contributions should be distinguished, even though the natural and anthropic effects can be tightly bound. Geophysical monitoring requires a detection of the geometry, which can be adapted to the type of reservoir. Different detection geometries are shown in Figure. 2. The specific design may include sensors suspended in air or water, such as in the monitoring of offshore reservoirs, where instrumentation can be deployed in the water and/or towed by a ship, or in airborne surveys. Limiting factors in reservoir monitoring are specific for the different types of reservoir, and may include available capital, time, current technology level, geometry, feasibility, etc. Geophysical methods can be used to create images of a reservoir at different times, constituting snapshots at the time of data acquisition. Such imaging is one of the most powerful tools and can be used to visualize temporal changes in a reservoir. Changes can be tracked by inspecting the variations of a measured variable after time t0 on a Cartesian 12 Figure 2. InSAR-observed subsidence north of Bakersfield, California, associated with hydrocarbon extraction between August 1997 and July 1998. 13 plot, where one of the axes is time or a different variable to monitor. 2.2 HC reservoir properties and monitored quantities In this section, the term property will refer to a specific feature that characterizes a certain object. Applying this concept to a hydrocarbon reservoir, a reservoir property refers to a specific characteristic of the reservoir. In the following, we assume that the property can be quantified. The property will be associated with one or more quantities, each of which carries a physical unit, at least in principle. A monitored quantity can be measured by using a specific measurement tool or monitoring technique. We do not make a distinction between directly measured quantities and inferred ones (e.g., by some inversion procedures such as seismic tomography, or derived by empirical relationships). Tracking changes of a property means measuring temporal changes of its associated monitored quantities. We deal in principle with three overlapping groups of characteristics, which can be identified as the following: 1) Quantifiable in principle (these are the “properties” according to our definition), 2) Measurable by some specific techniques (these are the Monitored Quantities (MQs)), and 3) Monitor-able (implying that repeated measurements are possible and useful for a reservoir monitoring). For the purpose of our study (“monitoring”), we are mostly interested in the third group, which is a subset of the second group. Geophysical Properties (GPs) can be identified and distinguish them from NonGeophysical Properties (NGPs), which we do not cover in this work. According to the different fields of reservoir analysis, we associate to each GP one or more specific MQs. In some cases, the same MQ can be associated to different GPs; for 14 instance the identified MQ “strains and stresses” can be associated to GPs, such as “stress regime in situ” and “rock cohesion” or “fracturing.” 2.3 Overview of the geophysical monitoring techniques Each Geophysical Monitoring Technique (GMT) will be presented relating to its associated MQs. According to the methodologies currently used in monitoring, the GMTs can be divided mainly into eight classes (the word techniques in each item is omitted): Magnetic (MA), Electrical (EL), Electromagnetic (ELM), Borehole Logging (LO), Nuclear (NU), Static (ST), Seismic (SE), and Geodetic (GD). In our analysis, we include the oldest GMTs, namely those in usage 40-50 years ago in applied geophysics, such as gravimeters, hydrometers, tiltmeters, strain meters, etc. and the techniques applied in water contamination studies such as the spontaneous potential, nuclear logs. Other techniques such as time-lapse techniques are considered, since they provide an efficient way of imaging reservoir changes, based on comparison among snapshots of the reservoir taken at various time. We include recent developments in geodetic techniques (remote sensing) such as the GPS satellite and the InSAR interferometry that have practical usage in the estimation of the subsurface deformation. InSAR measurements can be used to monitor the fluid flow in the subsurface and the sealing properties of faults. Figure 2 shows the subsidence observed with InSAR and associated with hydrocarbon extraction between August 1997 and July 1998; differences in color scale indicate the phase-wrapped vertical displacement. In the context of the hydrocarbon reservoir monitoring, bulk changes of a reservoir 15 can be mapped and tracked. During HF treatments, borehole/surface arrays can detect the induced microseismicity, inferring some important rock properties such as rock cohesion in presence of fracturing or faulting. Some other seismic techniques are included, as well, such as noise cross-correlation techniques that, by using the cross-correlation functions, have the capability to detect small variations in the correlations. The basic idea is to associate these small variations to perturbations in structures and velocity. This feature holds considerable promise for monitoring a hydrocarbon reservoir. With these techniques, it is possible to measure surface wave group and phase velocities, changes in velocities of body wave arrival, conversion P-to-S and S-to-P. In addition, the anisotropy can be inferred in structure changes, combined with shear-wave splitting, due for example to the preferred orientation or presence of faulting/folding and fracturing. Examples from the literature show applications to reservoir monitoring by ambient noise techniques (Lehujeur et al., 2014; Lehujeur et al., 2015; Obermann et al., 2015). Lehujeur et al. (2014) investigated the reliability of daily reservoir-scale near-surface continuous monitoring of the subsurface by ambient noise techniques, concluding that it may be useful for early detection of shorttime-scale hazards (days to weeks) such as migrating gases and fluids. By analyzing ambient seismic noisy cross correlations, they observed a significant loss of waveform coherence, horizontally and vertically constrained to the injection location of the fluid in a geothermal reservoir. The loss of waveform coherence (σ) has been interpreted as a local perturbation of the medium, allowing for causal relationships with the well activities (Figure 3). 16 Figure 3. Observed changes around the injection well using ambient seismic noise techniques, indicating a causal relationship with the activities at the well. 17 2.4 Conclusions A continuous massive usage of hydraulic stimulation in unconventional reservoirs requires an increasing demand for technology development in monitoring systems. This indicates the importance of adapting and updating the existing monitoring techniques. Hydrocarbon reservoirs appear in their initial state (unaltered) mainly under the regime of natural variations. However, significant deformations may happen after hydraulicfracturing treatments in reservoirs; these can generate significant changes in the rock properties, temperature, or fluid flow. Keeping track of such changes is useful on one hand for an enhanced hydrocarbon recovery, and to verify the safe and proper mode of operation. We describe the basic elements of geophysical monitoring, identifying properties, and monitored quantities in a hydrocarbon reservoir. The current available monitoring techniques are briefly reviewed according to the different field of analysis in reservoir. 18 CHAPTER 3 CARBON DIOXIDE (CO2) SEQUESTRATION In 2013, the global average atmospheric CO2 concentration surpassed 400 ppm (parts per million) for the first time since CO2 measurements began in 1958. Atmospheric global average of CO2 concentrations has fluctuated between 180 and 280 ppm over the past 400,000 years, as shown in Figure 4 (Huang et al., 2017). However, CO2 concentrations have increased since the industrial revolution from around 280 ppm to reach the current level of over 400 ppm. This high concentration is believed to be the result of anthropogenic release of CO2 from hydrocarbon production and refining, power generation, and the manufacture of cement, iron, and steel. The increase in CO2 and other greenhouse gases is believed to be the main reason for global warming. The Intergovernmental Panel on Climate Change has estimated the average atmospheric CO2 concentration to reach 750 ppm by 2100 in the absence of new polices or supply constraints. There are several solutions to alleviate CO2 emission: improve energy conversion and efficiency of fossil fuels, use of low-carbon or noncarbon energy sources, enhance uptake by biological sinks, or CO2 capture and storage. The most promising solution of them is the carbon capture and storage (CCS) that could mitigate the release of CO2 to the atmosphere while helping to meet present and future energy demands (Bruant et al., 2002). 19 Figure 4. Atmospheric global average level of carbon dioxide (CO2) over the past 400,000 years. 20 This chapter presents the carbon capture and storage process. Then, I will discuss four main suitable CO2 geological storage reservoirs. Finally, different types of CO2 trapping mechanism are shown. 3.1 Carbon capture and storage (CCS) The carbon capture and storage strategy involves five main steps. The first step is separating and capturing the CO2 from large sources of CO2 such as a fossil fuel power plant, as shown in Figure 5. The second step is compressing the captured CO2 to supercritical state. CO2 in its supercritical state is denser than in its gaseous state, leading to larger CO2 quantities to be stored in unit volume compared to the CO2 in gaseous state. In the third step, the compressed CO2 will be transported to the injection site where the most efficient means of transporting large quantities of CO2 is through pipelines. The fourth step is injecting the compressed CO2 into subsurface geological structure for permanent storage like depleted oil and gas fields. Finally, monitoring the sequestered CO2 to insure permanent storage (Rubin & De Coninck, 2005) is the fifth step. 3.2 Geological reservoirs for CO2 sequestration Large geological structures are favorable for CO2 sequestration projects like depleted oil and gas reservoirs, unmineable coal seams, basalt formations, and deep saline aquifers (NETL 2010). The hydrocarbon reservoirs have trapped oil and gas for millions of years, which suggest good seal for potential CO2 sequestration sites. These reservoirs have an existing infrastructure and have been well characterized and mapped. CO2 is injected to enhance oil 21 Figure 5. Carbon capture and storage (CCS) main process. 22 and gas recovery during production, which will offset the cost of CCS. Then, these reservoirs could be used as a permanent CO2 storage after oil and gas have been depleted. However, these reservoirs are considered challenging sites for geophysical monitoring tools since they contain multiple in situ fluids like oil, gas, brine, and CO2 (Gasperikova & Hoversten, 2008). The unmineable coal seam is a coal that is located at depths greater than 300 m. The injected CO2 will displace the methane within the coal, then the methane will be captured in a process known as Enhance Coal Bed Methane (ECBM). The recovered methane will contribute to offset the cost of the CCS (NETL 2010). Another advantage of the coal formations is the shallow depth compared to the other geological storage reservoirs. However, the coal formation is thin (3-10 m), which will be a challenging site for geophysical monitoring tools. Also, they could only sequester small volumes of CO2 (Gasperikova & Hoversten, 2008). The basalt formations are globally distributed and could store 10,000 years of the world’s emission of CO2. However, they are the least studied among the other geological storage reservoirs. Laboratory results showed that CO2 reacts with basalt minerals and water to form limestone and calcium carbonates. This mineral trapping mechanism is considered one of the securest trapping mechanisms. It transforms the CO2 from supercritical state to solid state, which isolates it from the atmosphere in a short period of time (NETL, 2010). The large storage capacity of deep saline reservoirs, along with the worldwide distribution and proximity to emission sources, have led to consider these reservoirs as ideal sequestration sites. The worldwide deep saline reservoirs capacity of CO2 storage has 23 been estimated to range from 100 to 10,000 GtCO2 (1 GtCO2 = 109 metric tons of CO2). Moreover, the large depth of these reservoirs will allow the CO2 to be injected and remains in its supercritical state leading to larger CO2 quantities to be stored in unit volume compared to the gaseous state in shallower depths. Using the geothermal and pressure gradients of 300 C/km and 10.5 MPa/km, respectively, engineers have found that CO2 at depths greater than 800 m could be stored as supercritical fluid with a density of 260 kg/m3. Due to the high salinity concentrations and depth of these reservoirs, they have little value as drinking or agricultural water usage. However, these sites need to be characterized and infrastructure would have to be built (Bruant et al., 2002). Figure 6 shows three types of water wells used in the field. 3.3 Trapping mechanism The physical, chemical, and hydrostatic state of the reservoir will determine the trapping mechanism of the CO2. These mechanisms are structural and stratigraphic trapping, residual trapping, solubility trapping, and mineral trapping. Structural and stratigraphic trapping refers to the physical accumulation of CO2 under low permeable layer. The change in densities between the injected CO2 and the brine will create a bouncy force that migrates the CO2 upward until it encounters a low permeable layer. Then, CO2 will move laterally to fill that structure until spill point. There are numerous scenarios of structural or stratigraphic trapping or a combination of both like anticline folds or sealed fault blocks in Figure 7 and other examples for CO2 trapping geology structure (Figure 8 and Figure 9). The efficiency of this trapping mechanism is measured in the strength of the cap rock sealing and the storage capacity of the reservoir, 24 Figure 6. Three types of water wells used in the field. 25 Figure 7. CO2 traps examples (a) structural and (b) stratigraphic traps (Zhang & Song, 2013). 26 Figure 8. Oil generation, migration, and accumulation. (http://large.stanford.edu/courses/2013/ph240/malyshev2/). 27 Figure 9. A typical anticline oil and gas reservoir. Oil is trapped by an impermeable cap rock and rests within a porous reservoir rock. 28 which is a function of porosity and permeability of the reservoir. This mechanism accumulates CO2 for a very long period, long enough for other trapping mechanism to come into effect (Zhu et al., 2013). Residual trapping or capillary refers to the immobile amount of CO2 that is trapped when the injection phase ends. The injected CO2 will force the brine away due to the injection current. As the injection stops, the brine pushes back the CO2 by the buoyancy force. During that process, part of the CO2 will be separated from the main CO2 plume to accumulate as immobile fluids (Zhu et al., 2013). Solubility trapping occurs when CO2 dissolves in the formation brine. The dissolution rate is a function of salinity, pressure, and temperature of the brine. The CO2 at the interface with the brine will dissolve into the brine by molecular diffusion. Due to the small molecular diffusion coefficient, CO2 will take thousands of years to completely dissolve in the brine. Even though this is a slow process, it enhances the structural and stratigraphic trapping efficiency in the long term. The CO2 saturated brine will be denser than the unsaturated brine by 1%, resulting in a downward movement of the CO2 saturated brine to the bottom of the reservoir. This process will enhance mixing, leading to more CO2 brine contact and therefore more CO2 dissolution. As a result of that, the upward migration of CO2 will decrease, which in turn will enhance the efficiency of the structural and stratigraphic trapping and also increase storage capacity (Zhu et al., 2013). Mineral trapping occurs when the dissolved CO2 forms weak carbonic acid that in turns reacts with the surrounding rock minerals to form carbonate precipitates. This is a very stable storing mechanism but very slow. Mineral trapping is only significant at geological time scale. 29 According to the Carbon Dioxide Capture and Storage report by the Intergovernmental Panel on Climate Change (IPCC), for well-selected, designed and managed geological storage sites, the vast majority of the CO2 will gradually be immobilized by various trapping mechanisms and, in that case, could be retained for up to millions of years. Because of these mechanisms, storage could become more secure over longer timeframes. 3.4 Conclusions In this chapter, I have briefly reviewed the carbon capture and storage (CCS) strategy. Also, I have discussed different potential geological storages for CCS. I have shown that deep saline reservoirs are ideal sequestration sites since they are widely distributed and have the potential to store large quantities of CO2. Then, I have presented four main trapping mechanisms of the CO2 in the subsurface geological structures. I have shown that with time CO2 trapping integrity increases due to solubility and mineral trapping mechanisms. However, under certain circumstances the caprock can fail such as when leakage occurs through fractures, faults, or pore spaces when capillary pressure has been exceeded. Therefore, continuous monitoring of CO2 sequestration sites is mandatory to detect any possible leakage of CO2 and to ensure the safety of the environment. Monitoring of CO2 sequestration sites will be discussed in depth in the next chapter. 30 CHAPTER 4 MODELING AND SENSITIVITY ANALYSIS OF THE SURFACE AND BOREHOLE GRAVITY AND GRAVITY GRADIOMETRY DATA In this chapter, forward modeling results are shown, which consist of modeling the surface data and borehole data. First, I show results of the forward modeling of gravity and gravity gradiometry data. The sensitivity tests of the gravity and gravity gradiometry instruments to different configurations of the model parameters follow. These tests show the sensitivity of the observed data to different model depth, thickness, size, and density values. These tests demonstrate how early in the sequestration process it is plausible to conduct the gravity and gravity gradiometry surveys. Finally, we show how the observed data change by moving the borehole to different positions. 4.1 Forward modeling of gravity and gravity gradiometry data A subsurface density contrast between anomalous body 𝜌(𝒓) and a homogenous background 𝜌𝑏 will generate a gravity anomaly. Studying the gravity anomaly is the basis of gravity prospecting. The gravity field is expressed in terms of the gradient of the gravity potential U: 𝒈 = ∇𝑈 where the gravity potential U satisfies the Poisson equation: (4.1) 31 𝛁𝟐 𝑈 = −4𝜋𝛾𝜌 (4.2) where 𝛁𝟐 is Laplace operator, 𝛾 is the universal gravitational constant 6.67x10-11 m3/kg.s2, 𝜌 is the body mass density. It is well known that the solution of the Poisson equation can be written in the volume integral form: 1 𝑈(𝒓′ ) = 𝛾 ∭𝐷 𝜌(𝒓) |𝒓−𝒓́ | 𝑑𝑣 (4.3) where D is the domain of mass concentration. By substituting Equation (4.3) in (4.1) and using the following identity: 1 𝒓−𝒓́ ∇′ |𝒓−𝒓́ | = |𝒓−𝒓́ |𝟑 (4.4) (where the prime at the gradient operator (∇′ ) means differentiation with respect to the variable r´), the gravity field of a mass distribution will be governed by the following: 𝒓−𝒓́ 𝒈(𝒓́ ) = 𝐴 𝑔 (𝜌) = 𝛾 ∭𝐷 𝜌(𝒓) |𝒓−𝒓́ |3 𝑑𝑣 (4.5) where 𝐴 𝑔 is the forward gravity operator. Therefore, the gravity anomaly generated by the anomalous density distribution ∆𝜌(𝒓) = 𝜌(𝒓) − 𝜌𝑏 is equal to the following: 𝒓−𝒓́ ∆𝒈(𝒓́ ) = 𝐴 𝑔 (∆𝜌) = 𝛾 ∭𝐷 ∆𝜌(𝒓) |𝒓−𝒓́ |3 𝑑𝑣 (4.6) The gravimeter instruments measure the vertical component of the gravity field gz, which is given as follows: 𝑧−𝑧́ 𝑔𝑧 (𝒓′ ) = 𝛾 ∭𝐷 𝜌(𝒓) |𝒓−𝒓́ |3 𝑑𝑣 (4.7) whereas the gravity gradiometry instruments measure the second spatial derivatives of the gravity potential U(r): 𝑔𝛼𝛽 (𝒓) = 𝜕2 𝜕𝛼𝜕𝛽 𝑈(𝒓) which form symmetric gravity gradiometry tensor: 𝛼, 𝛽 = 𝑥, 𝑦, 𝑧 (4.8) 32 𝑔𝑥𝑥 𝑔̂ = [𝑔𝑦𝑥 𝑔𝑧𝑥 𝑔𝑥𝑦 𝑔𝑦𝑦 𝑔𝑧𝑦 𝑔𝑥𝑧 𝑔𝑦𝑧 ] 𝑔𝑧𝑧 (4.9) where 𝑔𝛼𝛽 (𝒓) = 𝜕𝑔𝛼 𝛼, 𝛽 = 𝑥, 𝑦, 𝑧 𝜕𝛽 (4.10) Substituting Equation (4.3) in (4.8), the gravity tensor components can be calculated as follows: 𝑔𝛼𝛽 (𝒓) = 𝛾 ∭𝐷 𝜌(𝒓́ ) |𝒓́ 1 −𝒓|𝟑 𝐾𝛼𝛽 (𝒓́ − 𝒓)𝑑𝑣́ (4.11) where kernels 𝐾𝛼𝛽 are equal to 𝐾𝛼𝛽 (𝒓́ − 𝒓) = { 3 3 (𝛼−𝛼́ )(𝛽−𝛽́ ) |𝒓́ −𝒓|𝟐 (𝛼−𝛼́ )2 |𝒓́ −𝒓|𝟐 , 𝑖𝑓 𝛼 ≠ 𝛽 , 𝛼, 𝛽 = 𝑥, 𝑦, 𝑧 − 1, (4.12) 𝑖𝑓 𝛼 = 𝛽 The gravity gradiometers (FALCON) measures the horizontal curvature: (𝑔𝑥𝑥 −𝑔𝑦𝑦 ) 𝑔∆ = (4.13) 2 which can be expressed as 𝑔∆ (𝒓) = 𝛾 ∭𝐷 𝜌(𝒓́ ) |𝒓́ 1 −𝒓|𝟑 𝐾∆ (𝒓́ − 𝒓)𝑑𝑣́ (4.14) where 𝐾∆ (𝒓́ − 𝒓) = 3 (𝑥́ −𝑥)2 −(𝑦́ −𝑦)2 |𝒓́ −𝒓|𝟐 2 (4.15) One of the methods to solve the 3D forward and inverse gravity and gravity gradiometry problem is by using the point-mass approximation in discrete form. Thus, I 𝑁 𝑚 divided the domain D into Nm small rectangular cells Dk (𝐷 =∪𝑘=1 𝐷𝑘 ) with constant density for each cell 𝜌(𝒓) = 𝜌𝑘 . Therefore, Equation (4.7) can be rewritten as 𝑁 𝑚 𝑔𝑧 (𝒓́ ) = 𝛾 ∑𝑘=1 𝜌𝑘 ∭𝐷 𝑧−𝑧́ 𝑑𝑣 |𝒓−𝒓́ |3 𝑘 (4.16) 33 Following Zhdanov (2015), I have denoted the coordinates of the center of the cells as 𝒓𝑘 = (𝑥𝑘 , 𝑦𝑘 , 𝑧𝑘 ), 𝑤ℎ𝑒𝑟𝑒 𝑘 = 1, 2, … 𝑁𝑚 and the cell sides as dx, dy, and dz. Also, I have a discrete number (Nd) of observation points 𝒓𝑛́ = (𝑥𝑛́ , 𝑦𝑛́ , 0), 𝑤ℎ𝑒𝑟𝑒 𝑛 = 1, 2, … 𝑁𝑑 . Therefore, the discrete forward gravity modeling could be represented as follows: 𝑁 𝑔 𝑚 𝑔𝑧 (𝒓′𝒏 ) ≈ ∑𝑘=1 𝐴𝑛𝑘 𝜌𝑘 𝑤ℎ𝑒𝑟𝑒 𝑛 = 1, … . 𝑁𝑑 (4.17) 𝑔 where the gravity field kernel 𝐴𝑛𝑘 is 𝑔 𝐴𝑛𝑘 = 𝛾 𝑧𝑘 𝑑𝑥𝑑𝑦𝑑𝑧 (4.18) 3 𝑟𝑛𝑘 𝑟𝑛𝑘 = √(𝑥𝑘 − 𝑥𝑛′ )2 + (𝑦𝑛 − 𝑦𝑛′ )2 + 𝑧𝑘2 (4.19) By considering each cell as a point mass, Equation (4.11) and (4.12) of the gravity tensor components can be rewritten in discrete form as 𝑁 𝛼𝛽 𝑚 𝑔𝛼𝛽 (𝒓) ≈ ∑𝑘=1 𝐴𝑛𝑘 𝜌𝑘 𝑤ℎ𝑒𝑟𝑒 𝑛 = 1, … . 𝑁𝑑 (4.20) 𝛼𝛽 where 𝐴𝑛𝑘 is the forward operator of the gradiometry tensor components defined as 1 𝛼𝛽 𝐴𝑛𝑘 (𝒓) = 𝛾 ∑𝐷 𝜌(𝒓′ ) |𝒓́ −𝒓|3 𝐾𝛼𝛽 (𝒓́ − 𝒓)∆𝑥 ′ ∆𝑦 ′ ∆𝑧 ′ (4.21) where 𝐾𝛼𝛽 (𝒓́ − 𝒓) is the kernel defined as 𝐾𝛼𝛽 (𝒓́ − 𝒓) = { 3 3 (𝛼́ −𝛼)(𝛽́−𝛽) |𝒓́ −𝒓|𝟐 (𝛼́ −𝛼)2 |𝒓́ −𝒓|𝟐 , 𝑖𝑓 𝛼 ≠ 𝛽 , 𝛼, 𝛽 = 𝑥, 𝑦, 𝑧 − 1, (4.22) 𝑖𝑓 𝛼 = 𝛽 For example, 𝑔𝑧𝑧 (𝒓) = 𝛾 ∑𝐷 𝜌(𝒓′ ) |𝒓́ (𝑧́ −𝑧)2 1 −𝒓| 𝑔𝑥𝑦 (𝒓) = 3𝛾 ∑𝐷 𝜌(𝒓′ ) (3 3 (3 |𝒓́ −𝒓|2 (𝑥́ −𝑥)(𝑦́ −𝑦) |𝒓́ −𝒓|5 − 1) ∆𝑥 ′ ∆𝑦 ′ ∆𝑧 ′ (4.23) ) ∆𝑥 ′ ∆𝑦 ′ ∆𝑧 ′ (4.24) and the horizontal curvature component will have the form of the following: 3 (𝑥́ −𝑥)2 −(𝑦́ −𝑦)2 2 |𝒓́ −𝒓|𝟓 𝑔∆ (𝒓) = 𝛾 ∑𝐷 𝜌(𝒓́ ) ∆𝑥 ′ ∆𝑦 ′ ∆𝑧 ′ (4.25) 34 A paper by Jessop and Zhdanov (2005) demonstrated that the point-mass approximation yields very fast and accurate results when the depth to the center of the cell exceeds twice the dimension of the cell. They also presented quantitative comparison between the point-mass approximation and the exact prism body methods. The comparison showed that depending on the ratio of observation points to cells, the point-mass approximation method is up to 10 times faster and requires one-tenth the amount of memory required by the exact prism body method (Jessop & Zhdanov, 2005). 4.2 Description of the synthetic model In our numerical model study, we consider an injection borehole and a production borehole. The typical petroleum reservoir layers have an anticline shape, and they are formed of multiple layers. Based on this characteristic of reservoir, we have constructed an anticline type model for the synthetic study. Our anticline-type reservoir model extends from 800 m to 1000 m, as shown in Figure 10. The anomalous density is -1 g/cm3. The blue dots indicate the locations of the surface receiver, while the black dots represent the borehole receivers. In the numerical model study, we will make a comparison between the results with and without the noise in the data. 4.3 Forward modeling of observed data and comparison between the results with and without the noise After calculation of the forward modeling response, the observed data were plotted as Figures 11 and 12 for surface and borehole data, respectively. We can compare the results with and without the noise in the data. In addition, different configurations with 35 Figure 10. 3D view of the synthetic model of a reservoir. 36 Figure 11. Maps of the observed data without noise. 37 Figure 12. Plots of the observed data in the borehole without noise. 38 different borehole-target distances will be discussed in this section. From the result Figure 11, we can find from the surface field, that our target is symmetric in horizontal direction, and from the borehole observed data of Figure 12, we can easily find out the z direction location of our target. When we add 5% random noise to the field, the results are shown as Figure 13 for surface observed data, and Figure 14 for borehole data. The borehole position is along the y direction and in the middle of the x direction. That is why the responses from Gxy and Gxz components are very weak, almost equal to zero, the response at each side of borehole should counteract. Therefore, for the joint combination of components migration we will not consider Gxy and Gxz. After trying several different combinations, finally, the combination of Gzz and Gyz components is chosen. As it can be seen from the borehole observed data (Figure 13), the Gyz has either positive value or negative value. In order to remove the negative part, Gzz component has to be added. Finally, a horizontal sensitivity test in the borehole is made. The configuration figure, in which the only difference parameter is the distance between the well and our target, is shown as Figure 15. In different surveys, the response of borehole observed data are shown in the figures. In Figure 16-20, the observed data in the borehole located in different position away from the target are plotted. And Figure 21 is the comparison the response from these five different surveys for only Gzz and Gyz components, which are used in the following migration method. This figure indicates the horizontal sensitivity in the borehole more 39 Figure 13. Maps of the observed data with 5% random noise. 40 Figure 14. Plots of the observed data in the borehole with 5% random noise. 41 Figure 15. The wells with the receivers are located at different distances from the target. 42 Figure 16. Plots of the observed data in the borehole located 200 m away from the target. 43 Figure 17. Plots of the observed data in the borehole located 400 m away from the target. 44 Figure 18. Plots of the observed data in the borehole located 600 m away from the target. 45 Figure 19. Plots of the observed data in the borehole located 800 m away from the target. 46 Figure 20. Plots of the observed data in the borehole located 1000 m away from the target. 47 Figure 21. Comparison the different borehole survey results for Gzz and Gyz data. 48 efficiently. 4.4 Conclusions Even after adding 5% random noise to the data, the migration method can produce a similar result as compared to the noise free data. The results show symmetrical images in the plane view obtained from the surface data, and a correct depth location of the target in the images of the borehole data migration. The following chapter will demonstrate that migration and iterative migration methods are robust with respect to the noise in the data. The modeling results for different distances between the borehole and a target show that with the distance increase the response decays nonlinearly. The distance of 1000 m is the maximum distance at which the target still can be detect. 49 CHAPTER 5 GRAVITY AND GRAVITY GRADIOMETRY INVERSE PROBLEM The gravity and gravity gradiometry, inversion is the process of reconstructing the density model from the observed data set. It is well known that this is an ill-posed problem; therefore, regularization must be applied to produce a stable solution. In this chapter, the reweighted regularized conjugate gradient (RRCG) method is reviewed, which is used throughout the dissertation. The RRCG method is an iterative method used to minimize the Tikhonov parametric functional to obtain a stable solution of the gravity and gravity gradiometry inverse problem. 5.1 Inversion of gravity and gravity gradiometry data The solution of the gravity and gravity gradiometry inverse problems in discrete form is a linear inverse problem. Therefore, the discrete forward modeling operator for gravity anomaly can be expressed in general matrix notations as 𝑨𝒎 = 𝒅 (5.1) where (m) is the vector column of model parameter (density blocks 𝜌𝑖 ) of length Nm, (d) is the vector column of data like observed vertical gravity data gz of length Nd, and (A) is a rectangular matrix of size Nd x Nm formed by the gravity kernel, which is a system of equations that model the gravity response. Equation (5.1) is formulated as 50 𝑁 𝑚 𝑔11 𝑔12 ⋯ 𝑔1 𝑁 𝑔21 𝑔22 𝑔2 𝑚 ⋯ ⋮ ⋮ ⋮ ⋯ 1 2 𝑁 𝑔𝑁𝑑𝑚 ] [𝑔𝑁𝑑 𝑔𝑁𝑑 𝜌1 𝜌2 [ ⋮ ] = 𝜌𝑁𝑚 𝑑10 𝑑20 ⋮ 0 [𝑑𝑁𝑑 ] (5.2) This section will discuss how to reconstruct the density model from the observed gravity and gravity gradiometry data. The inversion of these data is an ill-posed problem, which means that the solution may not exist, may not be unique, or may not be stable. The solution does not exist if there is not an adequate numerical model from the given model set that could fit the observed data. The solution is called nonunique if there are more than one model that could fit the data. Finally, the solution is called not stable if a small perturbation of the data will produce a significant perturbation in the solution. It was believed that an ill-posed problem was not physically or mathematically meaningful. In 1977, Andrei Tikhonov a Russian mathematician laid the foundation of the regularization theory of the solution of ill-posed problems (Zhdanov, 2015). 5.1.1 Tikhonov regularization theory It was discussed earlier that the solution of the gravity and gravity gradiometry inversion is an ill-posed problem. Therefore, regularization is essential to obtain a stable solution. Tikhonov regularization is based on minimizing the Tikhonov parametric functional: 𝑃𝛼 (𝝆) = 𝜑(𝝆) + 𝛼𝑠(𝝆) = 𝑚𝑖𝑛 (5.3) where 𝑃𝛼 is the Tikhonov parametric functional, 𝜑(𝜌) is the misfit functional between the theoretical predicted data 𝐴(𝜌) and the observed data d, 𝑠(𝜌) is the stabilizer functional, and 𝛼 is the regularization parameter. Therefore, 51 2 𝑃𝛼 (𝝆) = ‖𝑊𝑑 (𝐴(𝝆) − 𝒅)‖2 + 𝛼‖𝑊𝑚 (𝝆 − 𝝆𝑎𝑝𝑟 )‖ = 𝑚𝑖𝑛 (5.4) where 𝝆𝑎𝑝𝑟 is a priori model, 𝑾𝑑 is the diagonal data weighting matrix, which is equal to the norms of each rows of matrix A 1⁄ 2 𝑇 2 𝑾𝑑 = 𝒅𝒊𝒂𝒈√∑𝑘(𝐴𝑤 𝑖𝑘 ) = 𝒅𝒊𝒂𝒈(𝑨𝑨 ) (5.5) and 𝑾𝑚 is the model weighting matrix which is equal to the diagonal integrated sensitivity matrix (Zhdanov, 2015): 𝑾𝑚 = [𝑊𝑗 ] = [𝑠𝑗 ] = 𝑺 𝑇 2 𝑆𝑘𝑤 = 𝒅𝒊𝒂𝒈√∑𝑖(𝐴𝑤 𝑖𝑘 ) = 𝒅𝒊𝒂𝒈(𝑨 𝑨) (5.6) 1⁄ 2 (5.7) There are many approaches to solve the inverse problem by minimization of the parametric functional. Also, there are many choices of stabilizer functionals and different methods of calculating the regularization parameter. Each of these subjects will be discussed in the following sections. 5.1.2 Regularization parameter The regularization parameter 𝛼 is a tradeoff between the best fitting and most reasonable stabilization. If 𝛼 is selected to be very large, then the minimization of the parametric functional 𝑃𝛼 (𝝆) is equivalent to the minimization of the stabilizing functional 𝑠(𝝆), which will lead to a solution that is very similar to the priori model. If 𝛼 is selected to be very small (almost no regularization is applied), then the minimization of the parametric functional 𝑃𝛼 (𝝆) is equivalent to the minimization of the misfit functional which could lead to unstable solution. Therefore, the choice of the optimal regularization parameter is essential for stable realistic solution. There are many methods of calculating the optimal regularization parameter. 52 The L-curve method is a simple graphical tool for quasi optimal selection of the regularization parameter. It is based on plotting all possible 𝛼 of the misfit functional versus stabilizing functional in a log-log scale. The produced curve will have L-shape appearance, as shown in Figure 22. The quasi optimal value of the regularization parameter is chosen at the corner of the L-curve (Hansen, 2005). Adaptive regularization, following Zhdanov, is another method for determining the regularization parameter empirically (Zhdanov, 2015). In an iterative inversion method, 𝛼0 = 0 for the first iteration that is running the inversion code without regularization. After the first iteration, 𝛼1 is calculated as the ratio of the misfit functional and the stabilizer functional for the initial model as follows 𝛼1 = ‖𝑾𝑑 𝑨(𝒎1 )−𝑾𝑑 𝒅‖2 2 ‖𝑾𝑚 𝒎1 −𝑾𝑚 𝒎𝑎𝑝𝑟 ‖ (5.8) Then, for each subsequent iteration 𝛼 is reduced as follows 𝛼𝑘 = 𝛼1 𝑞 𝑘−1 𝑘 = 1,2,3, … … . . 𝑞 > 0 (5.9) where q is an empirical coefficient between (0.5 to 0.9) based on numerical experiment. 5.1.3 Stabilizing functionals Stabilizing functional is incorporated in the inversion algorithm to select the appropriate class of models for the inversion solution. The traditional choices like minimum norm and maximum smoothness stabilizers are considered smooth stabilizing functionals since they produce smooth solutions. In environments with sharp geological boundaries, focusing stabilizers are applied like minimum support, minimum gradient support, and minimum vertical support stabilizers. These stabilizing functionals produce sharp and compact model parameter distributions (Zhdanov, 2015). 53 Figure 22. Shows the L-curve representation of selecting the quasi optimal regularization parameter by plotting the misfit functional on the x-axis versus the stabilizing functional in the y-axis in log-log scale for all possible 𝛼. 54 The minimum norm stabilizer is defined as the square L2 norm of the difference between the calculated model and a priori model: 2 𝑠𝐿2𝑎𝑝𝑟 (𝑚) = ‖𝒎 − 𝒎𝑎𝑝𝑟 ‖𝐿 = 𝑚𝑖𝑛 (5.10) 2 Another smooth stabilizing functional is the maximum smoothness stabilizer, which is defined as the L2 norm square of the gradient of the model parameters: 𝑠max 𝑠𝑚 (𝑚) = ‖𝛁𝒎‖2𝐿2 = 𝑚𝑖𝑛 (5.11) The minimum support stabilizer is a focusing stabilizer, which is proportional to the volume (support) of the nonzero values of the difference between the calculated model and a priori model. It provides a model with a minimum area of the distribution of anomalous parameters: 𝑠𝑀𝑆 (𝑚) = ∫𝑉 (𝑚−𝑚𝑎𝑝𝑟 ) 2 2 (𝑚−𝑚𝑎𝑝𝑟 ) +𝛽2 𝑑𝑣, (5.12) where 𝛽 is the focusing parameter that influences the sharpness and stability of the inversion. A very small 𝛽 may lead to unstable solution, whereas very large 𝛽 may lead to smooth inversion. Another focusing stabilizer is the minimum gradient support stabilizer, which is defined as the gradient of the minimum support functional integrand: 𝑠𝑀𝐺𝑆 (𝑚) = ∫𝑉 ∇ (𝑚−𝑚𝑎𝑝𝑟 ) 2 2 (𝑚−𝑚𝑎𝑝𝑟 ) +𝛽2 𝑑𝑣 (5.13) The vertical minimum support stabilizer is also a focusing stabilizing functional, which is defined in analogy to the minimum support stabilizer as (Zhdanov, Gribenko & Čuma, 2007): 2 𝑠𝑉𝑀𝑆 (𝑚) = ∫𝑉 (𝑚−𝑚𝑎𝑝𝑟 ) 2 ∬𝑠 (𝑚−𝑚𝑎𝑝𝑟 ) 𝑑𝑥𝑑𝑦+𝛽2 where S is a horizontal section of the rectangular domain V. 𝑑𝑣 (5.14) 55 Figure 23 shows a comparison between smooth and focusing inversions of a rectangular model. The smooth inversion produced a smooth image while the focusing inversion using the minimum gradient support stabilizer produced a sharp image of the true rectangular model. Thus, stabilizer controls the selection of class of the inversion solution. These stabilizers could be expressed in a form of a pseudo quadratic functional as squared L2 norm of the weighted model parameters: 𝑠(𝒎) = (𝑊𝑒 (𝒎 − 𝒎𝑎𝑝𝑟 ), 𝑊𝑒 (𝒎 − 𝒎𝑎𝑝𝑟 )) 𝐿2 2 = ∫𝑉 |𝑤𝑒 (𝒓) (𝑚(𝒓) − 𝑚𝑎𝑝𝑟 (𝒓))| 𝑑𝑣, (5.15) where We is a linear operator resulted from multiplication the model parameters function m(r) by the function 𝑤𝑒 (𝒓), which depends on m. Using stabilizer functionals in pseudo quadratic form simplifies the regularization problem solution and gives a unified approach to regularization using different stabilizers. In matrix notations and assuming mapr = 0, these stabilizers will be expressed as follows: 𝑠𝑚 𝑾max = 𝒅𝒊𝒂𝒈 [ 𝑒 𝑾𝑀𝑆 𝑒 = 𝒅𝒊𝒂𝒈 [ ∇𝑚 1 (𝑚2 +𝑒 2 ) ⁄2 1 1 (𝑚2 +𝛽2 ) ⁄2 𝑾𝑀𝐺𝑆 = 𝒅𝒊𝒂𝒈 [ 𝑒 ] (5.16) ] (5.17) ∇𝑚 1 1 ] ((∇𝑚)2 +𝛽2 ) ⁄2 (𝑚2 +𝑒 2 ) ⁄2 𝑾𝑉𝑀𝐺 = 𝒅𝒊𝒂𝒈 [ 𝑒 1 2 +𝛽 2 ) (∑𝑖𝑗 𝑚𝑖𝑗𝑘 1⁄ 2 ] , (5.18) (5.19) where e is included in the denominator to avoid singularity when m = 0. Therefore, the corresponding parametric functional could be written as follows: 56 Figure 23. Example of focusing stabilizer application to the solution of 1D inverse problem. Minimum gradient support stabilizer provides a sharp boundary solution by minimizing the area where strong model parameter variations and discontinuity occur (the gradient support is minimum). 57 2 𝑃𝛼 (𝒎, 𝒅) = ‖𝑾𝑑 (𝑨(𝒎) − 𝒅)‖2 + 𝛼‖𝑾𝑒 𝑾𝑚 (𝒎 − 𝒎𝑎𝑝𝑟 )‖ = 𝑚𝑖𝑛, (5.20) where 𝑾𝑒 is a variable matrix, which depends on m, and 𝑾𝑚 is the conventional fixed diagonal matrix of the model weighting parameter (Zhdanov, 2015). 5.1.4 Reweighted regularized conjugate gradient (RRCG) method There are various methods to solve the inversion by minimizing the parametric functional. One of the simplest methods is the least square method. However, the least square method is not practical if the forward modeling operator is huge, which makes finding the inverse matrix is so difficult. Therefore, iterative methods are more practical in dealing with large size forward operators. The reweighted regularized conjugate gradient (RRCG) method that is a form of iterative gradient-type methods is used to developed by (Portniaguine & Zhdanov, 1999). In this method, the variable weighting matrix 𝑾𝑒 is precomputed before each iteration using the computed model from the previous iteration. Therefore, it is considered a fixed matrix during each iteration. In RRCG method, the regularized steepest ascent direction (𝒍𝛼 ) is simplified as follows: 𝒍𝛼 (𝒎𝑛 ) = 𝑭𝑇𝑚𝑛 𝑾2𝑑 (𝑨(𝒎𝑛 ) − 𝒅) + 𝛼𝑾2𝑒 𝑾2𝑚 (𝒎𝑛 − 𝒎𝑎𝑝𝑟 ) (5.21) where 𝑭𝑚𝑛 is the Fréchet derivative matrix. The direction of ascent of the conjugate gradient 𝒍̃𝛼 (𝒎0 ) in the first iteration is equal to the direction of the regularized steepest ascent of the initial model 𝒎0 . For the following iteration, the direction of ascent of the conjugate gradient 𝒍̃𝛼 (𝒎𝑛 ) is a linear combination of the regularized steepest ascent direction on the current iteration and the direction of ascent of the conjugate gradient from the previous iteration (Zhdanov, 2015): 58 𝒍̃𝛼 (𝒎0 ) = 𝒍𝛼 (𝒎0 ) (5.22) 𝒍̃𝛼 (𝒎1 ) = 𝒍𝛼 (𝒎1 ) + 𝛽1𝛼 𝒍̃𝛼 (𝒎0 ) (5.23) 𝛼 The length of the iteration step using the linear line search 𝑘̃𝑛 𝑛 is expressed as 𝛼 𝑘̃𝑛 𝑛 = 𝛼 𝛼 (𝒍̃𝑛𝑛 ,𝒍𝑛𝑛 ) 2 2 𝛼 𝛼 ‖𝑾𝑑 𝑭𝑚𝑛 𝒍̃𝑛𝑛 ‖ +𝛼‖𝑾𝑚 𝒍̃𝑛𝑛 ‖ (5.24) 𝛼 whereas the 𝛽𝑛 𝑛 coefficient is defined as ‖𝒍𝛼 (𝒎𝑛 )‖2 𝛽𝑛𝛼 = ‖𝒍𝛼(𝒎 (5.25) 2 𝑛−1 )‖ The stabilizer functional can increase with iterations due to reweighting. Therefore, the gamma term is introduced to quantify the change in the stabilizer function from current to the previous iterations. 𝑠(𝒎𝑛 ) 𝛾𝑛 = (5.26) 𝑠(𝒎𝑛−1 ) To ensure the convergence of the parametric functional to global minimum, adaptive regularization (discussed in Section 4.2.2) was used and decreased if 𝛾 > 1: 𝛼𝑛 𝛼𝑛 = {𝛼𝑛 ⁄𝛾 𝑖𝑓 𝛾 ≤ 1 (5.27) 𝑖𝑓 𝛾 > 1 The numerical schemes for the reweighted conjugate gradient method is summarized as follows (Zhdanov, 2015): 𝒓𝑛 = 𝐴(𝒎𝑛 ) − 𝒅 𝒍𝜶𝑛 = 𝑭𝑇𝑚𝑛 𝑾2𝑑 𝒓𝑛 + 𝛼𝑛 𝑾2𝑒𝑛 𝑾2𝑚 (𝒎𝑛 − 𝒎𝑎𝑝𝑟 ) 𝜶 𝜶 𝒍̃𝟎 𝟎 = 𝒍𝟎 𝟎 𝛼 𝛼 𝛼 𝒍̃𝑛𝑛 = 𝒍𝑛𝑛 + 𝛽𝑛 𝑛 𝒍̃𝛼−1 𝑛−1 𝛼 𝛼 𝛽𝑛 𝑛 = 2 ‖𝒍𝑛𝑛 ‖ 𝛼 2 𝑛−1 ‖ ‖𝒍𝑛−1 𝛼 𝑘̃𝑛 𝑛 = 𝛼 𝛼 (𝒍̃𝑛𝑛 ,𝒍𝑛𝑛 ) 𝛼 2 𝛼 2 ‖𝑾𝑑 𝑭𝑚𝑛 𝒍̃𝑛𝑛 ‖ +𝛼‖𝑾𝑚 𝒍̃𝑛𝑛 ‖ 𝛼 𝛼 𝒎𝑛+1 = 𝒎𝑛 − 𝑘̃𝑛 𝑛 𝒍̃𝑛𝑛 𝑠(𝒎𝑛+1 ) 𝛾= 𝑠(𝒎𝑛 ) 59 𝛼𝑛 𝑖𝑓 𝛾 ≤ 1 𝛼𝑛+1 = {𝛼𝑛 ⁄𝛾 𝑖𝑓 𝛾 > 1 The iterative process is terminated if the misfit reaches a predefined noise level: 𝜙(𝒎𝑛 ) = ‖𝒓𝑁 ‖2 ≤ 𝜀0 (5.28) (5.29) 5.2 Conclusions Most of geophysical data inversions are ill-posed problems, where regularization must be applied to obtain a stable solution. Tikhonov regularization provides a stable solution of an ill-posed problem through minimization of the parametric functional. The parametric functional is a combination of the misfit functional and a product of regularization parameter and a stabilizer functional. In this chapter, the discrete form of 3D gravity and gravity gradiometry inversion method, in particular, reweighted regularized conjugate gradient (RRCG) method, has been described as a solution of the inversion process. The RRCG method is an iterative gradienttype method where the variable weighting matrix We is precomputed before each iteration. Also, I have discussed different choices of stabilizer functionals that determine the class of the inversion solution. Unlike smooth stabilizers, which provide a smooth distribution of the model parameter, the focusing stabilizers provide a compact solution making them suitable for environments with sharp model parameter contrast and for thin layers as well. 27 CHAPTER 6 JOINT MIGRATION OF SURFACE AND BOREHOLE GRAVITY AND GRAVITY GRADIOMETRY DATA Migration is an integral transform of the observed gravity gradients to a 3D density model. It is based on the adjoint operator, which is the key point of the migration method. Unlike the inverse operator, the adjoint operator can generate a well-behaved analytical function, which means the process of migration method is more stable and faster than the inversion. Besides, it is a linear problem for potential field, described in close form as an integration, which makes it robust to the noise level. In this chapter, I will demonstrate the algorithm of migration and joint migration method. 6.1 Gravity gradiometry data The gravity field, g, satisfies the following equations: ∇ ∙ 𝐠 = −4πρ, ∇ × 𝐠 = 0, (6.1) where γ is the universal gravitational constant (γ = 6.67384 × 10−11 m3 /kgs), and ρ is the anomalous density distribution within a domain 𝐷. The solution of these equations is given by the following well-known integral formula: 61 r′ −r 𝐠(𝐫) = γ ∭𝐷 ρ(r ′ ) |r′ −r|3 𝑑𝑣 ′ , (6.2) where integration is conducted over the variable 𝐫 ′ . The gravity field can be expressed by the gravity potential U(r) as 𝐠(𝐫) = ∇U(𝐫), (6.3) where 𝐫 ′ −𝐫 U(𝐫) = γ ∭𝐷 ρ(𝐫 ′ ) |𝐫 ′ −𝐫| 𝑑𝑣 ′ . (6.4) The second spatial derivatives of the gravity potential U(r), 𝐠 𝛼𝛽 = 𝜕2 𝜕𝛼𝜕𝛽 U(𝐫), α, β = x, y, z, (6.5) 𝐠 𝑥𝑦 𝐠 𝑦𝑦 𝐠 𝑧𝑦 (6.6) form a symmetric gravity tensor: 𝐠 𝑥𝑥 𝐠 𝐠̂ = [ 𝑦𝑥 𝐠 𝑧𝑥 𝐠 𝑥𝑧 𝐠 𝑦𝑧 ], 𝐠 𝑧𝑧 where 𝐠 𝛼𝛽 = 𝜕g𝛼 𝜕𝛽 , α, β = x, y, z, (6.7) The expressions for the gravity tensor components can be calculated based on Equations (6.4) and (6.5): 𝐠 𝛼𝛽 = γ ∭𝐷 ρ(𝐫 ′ ) |𝐫 ′ −𝐫|3 𝐾𝛼𝛽 (𝐫 ′ − 𝐫)𝑑𝑣 ′ , (6.8) where the kernels, 𝐾𝛼𝛽 , are equal to (𝛼−𝛼 ′ )(𝛽−𝛽′ ) 3 |𝐫 ′ 2 , 𝛼 ≠ 𝛽 −𝐫| ′ 𝐾𝛼𝛽 (𝐫 − 𝐫) = { , α, β = x, y, z. 2 ′ (𝛼−𝛼 ) 3 |𝐫 ′ 2 − 1, 𝛼 = 𝛽 (6.9) −𝐫| In addition to the gravity tensor components described by Equations (6.8) and (6.9), the gravity gradiometers also measure the difference between the gradients, 62 1 𝐠 ∆ = (𝐠 𝑥𝑥 − 𝐠 𝑦𝑦 ), (6.10) 2 which can be expressed as 𝐠 ∆ = γ ∭𝐷 ρ(𝐫 ′ ) |𝐫 ′ −𝐫|3 𝐾∆ (𝐫 ′ − 𝐫)𝑑𝑣 ′ , (6.11) where 2 𝐾∆ (𝐫 ′ − 𝐫) = 2 3 (𝑥 ′ −𝑥) −(𝑦 ′ −𝑥) |𝐫 ′ −𝐫| 2 . (6.12) 6.2 Adjoint operators for gravity gradiometry inversion Consider a problem of the inversion of gravity and gravity gradiometry data using gradient-type methods. According to Equation (6.2), we have the following expression for the gravity gradiometry field: 𝐠 𝛼𝛽 = γ ∭𝐷 ρ(𝐫 ′ ) |𝐫 ′ −𝐫|3 𝐾𝛼𝛽 (𝐫 ′ − 𝐫)𝑑𝑣 ′ , (6.13) 𝑔 where 𝐠 𝛼 (r) is the predicted gravity field on the observation surface; 𝐴𝛼 (𝜌), 𝛼 = 𝑥, 𝑦, 𝑧, denotes the forward modeling operator for different gravity field components; and the kernel 𝐾𝛼 (r ′ − r) is equal to: 𝐾𝛼 (r ′ − r) = α′ − α, α = x, y, z. (6.14) Assume that we have observed some gravity field 𝐠 𝑜𝑏𝑠 𝛼 (r) on the observational surface S, and domain D is in the lower half-space. The problem is to determine the density distribution, 𝛒(r ′ ). For simplicity, we first ignore the ill-posedness of gravity inversion and reduce the inversion problem to a minimization of the misfit functional between the observed and predicted data: ∅(𝜌) = ‖𝐠 𝛼 − 𝐠 obs 𝛼 ‖ = min. (6.15) 63 To solve the minimization problem (6.15), the direction of steepest ascent at the point 𝜌 of the model space M can be expressed as g ∗ obs 𝒍(𝜌) = 𝐴∗ 𝑟 = 𝐴∗ (𝐴𝛼 (𝜌) − 𝐠 obs 𝛼 ) = 𝐴 (𝐠 𝛼 − 𝐠 𝛼 ), (6.16) where the star ∗ denotes the adjoint operator. The adjoint operator 𝐴∗ for the gravity problem is equal to (Zhdanov et al., 2011) 𝐴∗𝛼 (𝑓) = γ ∬𝑆 𝑓(r) 𝐾 (r ′ |r′ −r|3 𝛼 − r)𝑑𝑠. (6.17) Therefore, according to equation (6.16), the direction of steepest ascent is equal to 𝒍(𝜌) = γ ∬𝑆 𝐠 𝛼 (r)−𝐠 obs 𝛼 (r) 𝐾𝛼 (r ′ |𝐫 ′ −𝐫|3 − r)𝑑𝑠, (6.18) where 𝐠 𝛼 (r) is the predicted gravity field on the observation surface. Let us now consider gravity gradiometry inversion. We assume that we have observed some gravity gradients 𝐠 obs 𝛼𝛽 (r) on the observational surface S, and the domain D is in the lower half space. Again, the problem is to determine the density distribution, ρ(r ′ ). The expression for the gravity tensor field is 𝐠 𝛼𝛽 = 𝐴𝛼𝛽 (𝜌) = γ ∭𝐷 𝛒(r′ ) |𝐫 ′ −𝐫|3 𝐾𝛼𝛽 (𝐫 ′ − 𝐫)𝑑𝑣 ′ , r ∉ 𝐷. (6.19) The adjoint operator for gravity gradients is given by the following equation (Zhdanov et al., 2011): 𝐴∗𝛼𝛽 (𝑓) = γ ∬𝑆 𝑓(r) 𝐾 (r ′ |r′ −r|3 𝛼𝛽 − r)𝑑𝑠. (6.20) According to Equation (1.16), the direction of the steepest ascent is equal to 𝒍(𝜌) = γ ∬𝑆 𝐠 𝛼𝛽 (r)−𝐠 obs 𝛼𝛽 (r) |𝐫 ′ −𝐫|3 𝐾𝛼𝛽 (𝐫 ′ − 𝐫)𝑑𝑠, where 𝐠 𝛼𝛽 (r) is the predicted gravity gradient field on the observation surface. (6.21) 64 For the 𝐠 ∆ component, the adjoint operator is given by the following equation (Zhdanov et al., 2011): 𝐴∗∆ (𝑓) = γ ∬𝑆 𝑓(r) 𝐾 (𝐫 ′ |𝐫 ′ −𝐫|3 ∆ − 𝐫)𝑑𝑠. (6.22) 6.3 Migration of the surface gravity and gravity tensor fields and 3D density imaging Let us assume that we have observed some component of the surface gravity field 𝐠 S𝛼 (r) and/or some surface gravity gradients 𝐠 S𝛼𝛽 (r) over an observational surface S, located in the air or on the ground. The problem is to determine the 3D density distribution, ρ(r ′ ) under the ground. Following Zhdanov (2002), the surface migration gravity field, 𝐠 Sm 𝛼 (r) is introduced as a result of application of the adjoint gravity operator, 𝐴𝛼𝑆∗ to the observed component of the surface gravity field 𝐠 S𝛼 : 𝑆∗ S 𝐠 Sm 𝛼 (r) = 𝐴𝛼 𝐠 𝛼 , (6.23) where the adjoint operator 𝐴𝛼𝑆∗ for the gravity problem (Zhdanov et al., 2011) is equal to 𝐴𝛼𝑆∗ (𝑓) = ∬𝑆 𝑓(𝐫) 𝐾 (𝐫 ′ |𝐫 ′ −𝐫|3 𝛼 − 𝐫)𝑑𝑠. (6.24) From the physical point of view, the migration field is obtained by moving the sources of the observed gravity field above the observational surface. Nevertheless, the migration field contains some remnant information about the original sources of the gravity anomaly. That is why it can be used in imaging the sources of the gravity field. In a similar way, we can introduce a surface migration gravity tensor field 𝐠 Sm 𝛼𝛽 (r) and use the following notations for the components of this tensor field: 65 𝑆∗ S 𝐠 Sm 𝛼𝛽 (r) = 𝐴𝛼𝛽 𝐠 𝛼𝛽 , (6.25) 𝑆∗ where the adjoint operators, 𝐴𝛼𝛽 applied to some function 𝑓(r), are given by the formulas: 𝑆∗ (𝑓) = ∬𝑆 𝐴𝛼𝛽 𝑓(𝐫) 𝐾 (𝐫 ′ |𝐫 ′ −𝐫|3 𝛼𝛽 − 𝐫)𝑑𝑠. (6.26) We should note, however, that the direct migration of the observed gravity and/or gravity tensor fields does not produce an adequate image of the subsurface density distribution because the migration fields rapidly attenuate with the depth, as one can see from expressions (6.24) and (6.26). In order to image the sources of the gravity fields at their correct location, one should apply an appropriate spatial weighting operator to the migration fields. This weighting operator is constructed based on the integrated sensitivity of the data to the density. We can find a distribution of the density of the gravity field sources, described by the following expression: 𝜌𝛼𝑆𝑚 (𝐫) = k 𝛼𝑆 (𝑤𝛼𝑆 )−2 (𝑧)𝐠 Sm 𝛼 (𝐫), (6.27) where unknown coefficient k 𝛼𝑆 can be determined by a linear line search (Zhdanov, 2002) according to the following expressions: 2 k 𝛼𝑆 = S ‖𝐴𝑤∗ 𝛼 𝐠 𝛼 ‖𝑀 2 𝑤∗ S ‖𝐴𝑤 𝛼 𝐴𝛼 g𝛼 ‖𝐷 𝑆 −1 𝐴𝑤 𝛼 = 𝐴𝛼 𝑊𝛼 , , (6.28) (6.29) and the linear weighting operator 𝑊𝑚 = 𝑊𝛼 is selected as a linear operator of multiplication of the density 𝜌 by a function, 𝑤𝛼 , equal to the square root of the integrated sensitivity of the gravity field, 𝑆𝛼 : 𝑤𝛼𝑆 = √𝑆𝛼𝑆 . (6.30) The integrated sensitivity, 𝑆𝛼𝑆 , can be computed using the following formula: 66 1 𝑆𝛼𝑆 = 𝑐𝛼 |𝒁| , 𝑧 < 0, 𝛼 = 𝑥, 𝑦, 𝑧, (6.31) where coefficients, 𝑐𝛼 , are the corresponding constants for the different components equal to 𝜋 𝑐𝑥 = 𝑐𝑦 = 𝛾√ , 𝑐𝛼 = 𝛾√𝜋. 2 (6.32) Equation (6.27) is called a migration density 𝝆𝛼𝑆𝑚 (𝜁): 𝝆𝛼𝑆𝑚 (r) = k 𝛼𝑆 (𝑤𝛼𝑆 )−2 𝐠 𝐒𝐦 𝜶 (r), (6.33) and is proportional to the weighted migration field, 𝐠 Sm 𝛼 : 𝝆𝛼𝑆𝑚 (r) = k𝑆 𝛼 𝑐𝛼 |𝑧|𝐠 Sm 𝛼 (r), (6.34) where 𝐠 Sm 𝛼 (r) = ∬𝑆 ′ 𝐠S 𝛼 (𝐫 ) |𝐫−𝐫 ′ |3 𝐾𝛼 (𝒓 − 𝐫 ′ )𝑑𝑠 ′ . (6.35) Thus, the migration transform with spatial weighting provides a stable algorithm for calculating the migration density. In a similar way, we can introduce a migration density based on the gravity tensor migration: −2 𝑆𝑚 𝑆 𝑆 (r) = k 𝛼𝛽 𝝆𝛼𝛽 (𝑤𝛼𝛽 ) 𝐠 Sm 𝛼𝛽 (r), (6.36) where 2 𝑆 k 𝛼𝛽 = S ‖𝐴𝑤∗ 𝛼𝛽 𝐠 𝛼𝛽 ‖ 𝑀 2 𝑤∗ S ‖𝐴𝑤 𝛼𝛽 𝐴𝛼𝛽 𝐠 𝛼𝛽 ‖ . (6.37) 𝐷 𝑆 Functions 𝑤𝛼𝛽 are equal to the square root of the integrated sensitivity of the 𝑆 gravity tensor fields, 𝑆𝛼𝛽 , respectively: 𝑆 𝑆 𝑤𝛼𝛽 = √𝑆𝛼𝛽 . (6.38) 67 The integrated sensitivity of the gravity tensor field is calculated as follows: 𝑆 𝑆𝛼𝛽 = 𝑐𝛼𝛽 1 𝒛2 , (6.39) where 𝑐𝛼𝛽 are the corresponding constants for different tensor components equal to 𝑐𝑧𝑧 = 𝑐𝑧𝑥 = 𝑐𝑧𝑦 = 𝛾 √3𝜋 , 𝑐𝑥𝑥 2 = 𝑐𝑦𝑦 = 𝛾 3√𝜋 . 4 (6.40) Expression (6.36) is called a tensor field migration density. It is proportional to the magnitude of the weighted tensor migration field 𝐠 Sm 𝛼𝛽 . Thus, migration transformation provides a stable algorithm for calculating migration density. Substituting Equation (6.39) for the weighting function 𝑤 𝑆 back into Equations (6.38) and (6.36), find that 𝑆𝑚 (r) = 𝝆𝛼𝛽 k𝑆 𝛼𝛽 𝑐𝛼𝛽 𝒛2 𝐠 Sm 𝛼𝛽 (r), (6.41) where 𝐠 Sm 𝛼𝛽 (r) = ∬𝑆 ′ 𝐠S 𝛼𝛽 (𝐫 ) |𝐫 ′ −𝐫|3 𝐾𝛼𝛽 (𝐫 ′ − 𝐫)𝑑𝑠 ′ . (6.42) 6.4 Migration of the borehole gravity and gravity tensor fields and 3D density imaging Consider that we have observed the components of the borehole gravity field 𝐠 B𝛼 (𝐫) and borehole gravity gradients 𝐠 B𝛼𝛽 (𝐫) over an observational line L, associated with a given borehole. The problem is to determine the 3D density distribution, 𝛒(𝐫 ′ ), around the borehole. Following (Zhdanov, 2002; Zhdanov et al., 2010b; Zhdanov et al., 2011), the borehole migration gravity field, 𝐠 Bm 𝛼 (r) is introduced as a result of application of the adjoint gravity operator, 𝐴B∗ 𝛼 to the observed component of the borehole gravity field: 68 B∗ B 𝐠 Bm 𝛼 (𝐫) = 𝐴𝛼 𝐠 𝛼 , (6.43) where the adjoint operator 𝐴B∗ for the gravity problem is equal to 𝛼 𝐴B∗ 𝛼 (𝑓) = ∫𝐿 𝑓(𝐫) 𝐾 (𝐫 ′ |𝐫 ′ −𝐫|3 𝛼 − 𝐫)𝑑𝑙 . (6.44) In a similar way, we can introduce a borehole migration gravity tensor field 𝐠 Bm 𝛼𝛽 (𝐫) along the borehole L, and use the following notations for the components of this tensor field: B∗ B 𝐠 Bm 𝛼𝛽 (𝐫) = 𝐴𝛼𝛽 𝐠 𝛼𝛽 , (6.45) where the corresponding adjoint operators, 𝐴B∗ 𝛼𝛽 applied to some function 𝑓(r), are given by the formulas: 𝐴B∗ 𝛼𝛽 (𝑓) = ∫𝐿 𝑓(𝐫) 𝐾 (𝐫 ′ |𝐫 ′ −𝐫|3 𝛼𝛽 − 𝐫)𝑑𝑙 . (6.46) Thus, the migration field can be calculated everywhere around the borehole for the given values of the gravity and/or gravity gradient field, measured along the borehole. The direct migration of the observed gravity and/or gravity tensor fields does not produce an adequate image of the subsurface density distribution because the migration fields rapidly attenuate with the depth, as shown by expressions (6.44) and (6.46). In order to image the sources of the gravity fields at their correct location, one should apply an appropriate spatial weighting operator to the migration fields. This weighting operator is constructed based on the integrated sensitivity of the data to the density. Formula (6.43) serves as the basis of the migration. Considering Equation (6.43) and the direction of the steepest ascent, one can find an approximation to the distribution of the density as follows: 69 B B −2 Bm 𝝆B𝑚 𝛼 (𝐫) = k 𝛼 (𝑤𝛼 ) (𝑧)𝐠 𝛼 (𝐫), (6.47) where unknown coefficient k B𝛼 can be determined by a linear line search (Zhdanov, 2002) according to the following equations: 2 k B𝛼 = B ‖𝐴𝑤∗ 𝛼 𝐠 𝛼 ‖𝑀 2 𝑤∗ B ‖𝐴𝑤 𝛼 𝐴𝛼 𝐠 𝛼 ‖𝐷 , B −1 𝐴𝑤 𝛼 = 𝐴𝛼 𝑊𝛼 , (6.48) (6.49) and the linear weighting operator 𝑊𝛼 is selected as a linear operator of multiplication of the density 𝜌 by a function, 𝑤𝛼B , equal to the square root of the integrated sensitivity of the integrated sensitivity of the complex intensity of the gravity field, 𝑆𝛼B : 𝑤𝛼B = √𝑆𝛼B . (6.50) The integrated sensitivity, 𝑆𝛼B , in accordance with the definition, is calculated as follows (Zhdanov, 2002). The integrated sensitivity of the gravity field measured in the borehole decreases with the distance from the borehole, 𝑅, as inverse 𝑅1.5 : 𝑆𝛼B ~ 1 𝑅1.5 . (6.51) Equation (6.47) is called a migration density 𝜌𝛼B𝑚 (𝜁): B B −2 Bm 𝝆B𝑚 𝛼 (r) = k 𝛼 (𝑤𝛼 ) 𝐠 𝛼 (𝐫), (6.52) and is proportional to the weighted migration field, 𝐠 Bm 𝛼 : B 1.5 Bm 𝝆B𝑚 𝛼 (𝐫) = k 𝛼 𝑅 𝐠 𝛼 (𝐫), (6.53) where 𝐠 Bm 𝛼 (𝐫) = ∫𝐿 ′ 𝐠B 𝛼 (𝐫 ) |𝐫−𝐫 ′ |3 𝐾𝛼 (𝐫 − 𝐫 ′ ) 𝑑𝑙′ . (6.54) Thus, the migration transformation with spatial weighting provides a stable algorithm for calculating the migration density. 70 In a similar way, we can introduce a migration density based on the gravity tensor migration: −2 B B Bm 𝝆B𝑚 𝛼𝛽 (𝐫) = k 𝛼𝛽 (𝑤𝛼𝛽 ) 𝐠 𝛼𝛽 (𝐫), (6.55) where 2 𝑆 k 𝛼𝛽 = B ‖𝐴𝑤∗ 𝛼𝛽 𝐠 𝛼𝛽 ‖ 𝑀 2 𝑤∗ B ‖𝐴𝑤 𝛼𝛽 𝐴𝛼𝛽 𝐠 𝛼𝛽 ‖ , (6.56) 𝐷 B And the weighting function 𝑤𝛼𝛽 is equal to the square root of the integrated sensitivity of B the gravity tensor fields, 𝑆𝛼𝛽 : B B 𝑤𝛼𝛽 = √𝑆𝛼𝛽 . (6.57) The integrated sensitivity of the gravity tensor field measured in the borehole is proportional to inverse 𝑅2.5 , where 𝑅 is a distance from the borehole to the target: B 𝑆𝛼𝛽 ~ 1 𝑅2.5 . (6.58) Expression (6.55) is called a tensor field migration density, and it is proportional to the magnitude of the weighted tensor migration field 𝐠 Bm 𝛼𝛽 . Thus, migration transformation provides a stable algorithm for calculating migration density. Substituting Equation (6.58) B for the weighting function 𝑤𝛼𝛽 back into Equations (6.55), we find that B 2.5 Bm 𝝆B𝑚 𝛼𝛽 (𝐫) = k 𝛼𝛽 𝑅 𝐠 𝛼𝛽 (𝐫), (6.59) where 𝐠 Bm 𝛼𝛽 (𝐫) = ∫𝐿 ′ 𝐠B 𝛼𝛽 (𝐫 ) |𝒓−𝐫 ′ |3 𝐾𝛼𝛽 (𝐫 − 𝐫 ′ )𝑑𝑙′ . (6.60) The following expressions describe the migration for the migration densities for the different gravity gradients as measure from a single borehole: 71 B 2.5 𝛒Bm ∫𝐿 𝛼𝛽 (𝐫) = k 𝛼𝛽 𝑅 ′ 𝐠B 𝛼𝛽 (𝐫 ) |𝐫−𝐫 ′ |3 𝐾𝛼𝛽 (𝐫 − 𝐫 ′ )𝑑𝑙′ , (6.61) where 𝑅 is the horizontal distance from the center of the borehole. 6.5 Joint migration The goal is to jointly migrate the surface and borehole gravity fields to make a clear image of a deep target. We consider a joint migration of the multiple components of the surface and borehole gravity and gravity tensor fields according to the following formula: 𝑆 𝑆 B B (𝐫) + 𝑐𝛼B 𝛒B𝛼 (𝐫) + 𝑐𝛼𝛽 𝛒𝑚 (r) = 𝑐𝛼𝑆 𝛒𝛼𝑆 (𝐫) + ∑ 𝑐𝛼𝛽 𝛒𝛼𝛽 𝛒𝛼𝛽 (𝐫), (6.62) 𝑆 B where 𝑐𝛼𝑆 , 𝑐𝛼𝛽 , 𝑐𝛼B and 𝑐𝛼𝛽 can be treated as the weights of the corresponding migration fields in the density model, which can be empirically determined from the results of the model studies. We use the following expressions for joint migration, which provides an averaging for both the surface and borehole data 𝑆 (𝐫))⁄(𝑁𝑠 + 1) + 𝑐2 (𝝆B𝛼 (𝐫) + ∑ 𝝆B𝛼𝛽 (𝐫))⁄(𝑁𝑏 + 1) (6.63) 𝝆𝑚 (𝐫) = 𝑐1 (𝝆𝛼𝑆 (𝐫) + ∑ 𝝆𝛼𝛽 if the gravity field is used in the migration, or 𝑆 (𝐫))⁄𝑁𝑠 + 𝑐2 (∑ 𝝆B𝛼𝛽 (𝐫))⁄𝑁𝑏 . 𝝆𝑚 (𝐫) = 𝑐1 (∑ 𝝆𝛼𝛽 (6.64) if the gravity field is not used. In the last formulas, 𝑁𝑠 is the number of the surface gravity gradient components, the 𝑁𝑏 is the number of borehole gravity gradient components, and 𝑐1 = 0.5 and 𝑐2 = 0.5. 72 6.6 Iterative migration Equation (6.62) can produce a migration image of the density distribution in the lower half-space. A better quality migration image can be produced by repeating the migration process iteratively (Wan & Zhdanov, 2013). We begin with the migration of the gravity and/or gravity tensor field data observed on the surface and/or in the borehole. In order to check the accuracy of our migration imaging, we apply the forward modeling and compute a residual between the observed and predicted data for the given density model: 𝑟1 = 𝐠 𝑝𝑟𝑒 − 𝐠 𝑜𝑏𝑠 , (6.65) where 𝐠 𝑜𝑏𝑠 is the observed gravity or gravity gradient component; 𝐠 𝑝𝑟𝑒 is the predicted gravity or gravity gradient component calculated according to formulas (6.2) or (6.8) with the density 𝝆1𝑚 obtained from Equation (6.64). If the residual is smaller than the prescribed accuracy level, we use the migration image as a final density model. In a case where the residual is not small enough, we migrate the residual field and produce the density correction, 𝛿𝝆1𝑚 , to the original density model using the same analysis, which we have applied for the original migration: 𝛿𝝆1𝑚 = 𝑘𝑤 −2 𝒓𝟏 , (6.66) 𝑚 𝑚 𝝆𝑚 2 = 𝝆1 − 𝛿𝝆1 , (6.67) where 𝛿𝝆1𝑚 stands for the migration image obtained by the residual field migration, Equation (6.64). A general scheme of the iterative migration can be described by the following formula: 𝑚 𝑚 𝝆𝑚 𝑛+1 = 𝝆𝑛 − 𝛿𝝆𝑛 . (6.68) The iterative migration is terminated when the residual field becomes smaller than the required accuracy level of the data fitting. 73 6.7 Summary In this chapter, the migration method for the surface and borehole gravity and gravity gradiometry data has been introduced. Migration is an integral transformation of the observed gravity gradients to a 3D density model, and it is based on the adjoint operator, which is the key point of the migration method. Unlike the inverse, the adjoint operator generates a well-behaved analytical function, which means the process of migration method is more stable, and faster than the inversion. A joint migration method has also been introduced, in which the trade-off coefficient 𝑐, relating different components of the migration field, is important. Finally, the iterative migration can be combined with the regularization method. This also allows us to apply the smoothing or focusing stabilizers to produce a more focused image of the target (Wan & Zhdanov, 2013). 74 CHAPTER 7 MODEL STUDY In this chapter, I present a few examples of 3D joint migration of the surface and borehole gravity gradient field data. Firstly, based on the characteristic of hydrocarbon reservoir, the anticline-type model is built for a synthetic model test. Then, the multiple-layer model is introduced since in real case, the reservoir can consist of multiple layers. In the previous chapter, we have discussed the availability of joint iterative migration of the surface and borehole gravity gradiometry data for reservoir monitoring. Here we apply this technology to synthetic monitoring model to demonstrate its utility for this kind of application. 7.1 Model 1 In this section, we consider Model 1 that contains one reservoir with an anticline shape, having the size of 1000 m x 1000 m x 200 m (L x W x H). It is known that the density of sandstone is between 2.2~2.8 g/cm3, and the density of shale is between 2.4~2.8 g/cm3, the density of petroleum is 0.64 g/cm3, the density of saltwater is about 1.02 g/cm3. In this synthetic model, consider that the reservoir is filled with gas, which has lower density than the surrounding rock, so the anomalous density is set to -1 g/cm3. The reservoir 75 is located at a depth of 0.9 km below the surface (see Figure 24). In Figure 24, the model is anticline type anomalous body, which is located in the center of the domain at a depth from 800 m to 1000 m. The blue dots indicate the surface receivers’ locations, while the black dots indicate the borehole location. The surface gravity sensors are distributed within the range of 6 km in the x direction and 5.6 km in the y direction with 100 m separation in the x and in y directions. A vertical borehole has the horizontal coordinates of x= 3000 m and y= 2200 m. The interval between the receivers in the borehole is 5 m in the z direction. The observed data contained 5% random noise (Figures 25 and 26). From Figure 25, the response is moreless symmetric in the horizontal direction, while in Figure 26, the depth location of the target is indicated by the borehole data. 7.1.1 Migration results for individual surface data Figures 27 and 28 present the result of iterative migration of the surface g 𝑧𝑧 component only at the cross section of x = 3000 m and y = 2800 m (bottom panel). The blue line shows the observed data, and the red line shows the predicted data at the surface (top panel). The data fitting is very good, and the migration was run for six iterations only. The inverse density image outlines the target well. However, it is extended far beyond the HC reservoir in the vertical direction. The white line shows the target location. Figures 29 and 30 present the result of iterative migration of the surface g 𝑥𝑧 component only at the cross section of x = 3000 m and y = 2800 m. It took just six iterations to reach the noise level in the misfit. The blue line shows the observed data, and the red line presents the predicted data at the surface. As in the previous case, the migration image 76 Figure 24. 3D view of the synthetic Model 1. 77 Figure 25. Maps of the observed data with 5% random noise. 78 Figure 26. Plots of the observed data in the borehole with 5% random noise. 79 Figure 27. Results of migration of the surface g 𝑧𝑧 data with 5% random noise at x = 3000 m. 80 Figure 28. Results of migration of the surface g 𝑧𝑧 data with 5% random noise at y = 2800 m. 81 Figure 29. Results of migration of the surface g 𝑥𝑧 data with 5% random noise at x = 3000 m. 82 Figure 30. Results of migration of the surface g 𝑥𝑧 data with 5% random noise at y = 2800 m. 83 shows the location of the target. However, the image is much thicker than the actual reservoir. Figures 31 and 32 present the result of iterative migration of the surface g 𝑦𝑧 component only at the cross section of x = 3000 m and y = 2800 m. We observed a similar result in this case to those in the previous gravity gradient components. Several other components were also migrated. The results show very good data fitting after even few iterations (e.g., six iterations), and the misfit decreases to a very low level, even with added 5% noise. Since the migration of gravity and gravity gradiometry data is the integration process, the noise does not affect the results. We can also see the typical characteristic of migration or inversion for the potential field data -- the results are diffuse and smeared-out compared to the reference model, while the maximum density of the migration results is significantly less than that of the true model. Three different components, g 𝑧𝑧 , g 𝑦𝑧 , g 𝑥𝑧 , are shown here, while in the combined migration, only g 𝑧𝑧 and g 𝑦𝑧 will be used. Besides these three components, I have tried other tensor components, but all the results show the diffusion problem like one discussed above. From those results, we can conclude that the individual surface data are not enough to recover the correct location of the reservoir, especially in the z direction. 7.1.2 Migration results of the borehole gravity data In this section, the 3D migration of the borehole gravity data is considered. Figures 33 and 34 present the result of iterative migration of the borehole g 𝑧𝑧 component only at the cross section of x = 3000 m and y = 2800 m (left panel). The blue 84 Figure 31. Results of migration of the surface g 𝑦𝑧 data w/ 5% random noise at x = 3000 m. 85 Figure 32. Results of migration of the surface g 𝑦𝑧 data w/ 5% random noise at y = 2800 m. 86 Figure 33. Results of migration of the borehole g 𝑧𝑧 data with 5% random noise at x = 3000 m. 87 Figure 34. Results of migration of the borehole g 𝑧𝑧 data with 5% random noise at y = 2800 m. 88 line shows the observed data and the red line shows the predicted data at the surface (right panel). The data fitting is very good, and the migration ran 20 iterations only. The inverse density image outlines the target well (in the z direction). Figures 35 and 36 present the result of iterative migration of borehole g 𝑦𝑧 component only at the cross section of x = 3000 m and y = 2800 m. It took just 20 iterations to reach the noise level in the misfit. The blue line shows the observed data, and the red line presents the predicted data at the surface. As in the previous case, the migration image shows the location (in the z direction) of the target. The response from the target can be seen everywhere in the horizontal direction of the response. However, another opposite response can be seen in the other side of the well. The white line shows the target location. Finally, Figure 37 presents the result of iterative migration of the borehole g ∆ component only at the cross section of x = 3000 m. We observed a similar result in this case to those in the previous studies. All other components were also migrated but the results are not present here. They show a very good data fit after 20 iterations, and the misfit decreases to a very low level, even with added the 5% noise. The results show the typical characteristic of borehole migration or inversion for the potential field data, and they recover the true model exactly at the right position in the z direction. However, the adjoint operator in the borehole migration algorithm is based on the line integration, which provides the correct z position of the reservoir only. In other words, the borehole migration does not have a horizontal sensitivity. From these results, we conclude that the migration of individual borehole data is not enough to determine the correct anomaly location, especially in the horizontal 89 Figure 35. Results of migration of the borehole g 𝑦𝑧 data with 5% random noise at x = 3000 m. 90 Figure 36. Results of migration of the borehole g 𝑦𝑧 data with 5% random noise at y = 2800 m. 91 Figure 37. Results of migration of the borehole g ∆ data with 5% random noise at x = 3000 m 92 direction. A better approach is to combine the surface and borehole data together to do the migration, which can present the correct location of the target in both the horizontal and vertical directions. 7.1.3 Joint migration results of surface and borehole data Figure 38 presents the results of joint iterative migration of the surface and borehole 𝑔𝑦𝑧 data. The blue line shows the observed data and the red line shows the predicted data at the surface (top panel) and in the borehole (right panel). Compared with Figure 31, the density image is much smaller than the result of migration of the surface 𝑔𝑦𝑧 data only, but the artifacts on the left side are still there. Figure 39 shows the results of joint iterative migration of the surface and borehole 𝑔𝑧𝑧 + 𝑔𝑦𝑧 data. One can see that the artifacts shown in Figure 38, have disappeared. We can also notice that the joint migration of the surface and borehole 𝑔𝑧𝑧 + 𝑔𝑦𝑧 components produce a more focused image of the target. Comparing with the results from the individual migration method, the resolution improves significantly. In particularly, the vertical resolution is much better than the individual surface migration results, while the horizontal resolution is much better than the individual borehole migration results. 7.1.4 Joint migration using multiple borehole data Although the multiple component joint migration provides a good resolution, we have noticed that the resolution improves further when the anomaly is closer to the borehole. For this reason, I have studied the case where an additional borehole is located on the opposite side of the target to gather more borehole data (Figure 40). 93 Figure 38. The result of joint iterative migration of surface and borehole 𝑔𝑦𝑧 component only at the cross section of x = 3000 m (bottom left panel). The blue line shows the observed data, and the red line shows the predicted data at the surface (top panel) and in the borehole (right panel). 94 Figure 39. The result of joint iterative migration of surface and borehole 𝑔𝑧𝑧 + 𝑔𝑦𝑧 data only at the cross section of x = 3000 m (bottom left panel). The blue line shows the observed data, and the red line shows the predicted data at the surface (top panel) and in the borehole (right panels). 95 Figure 40. 3D view of the synthetic model for two boreholes, shown by the solid vertical black lines. 96 The results are shown in Figures 41 and 42 for the joint migration of g 𝑧𝑧 and g 𝑦𝑧 components. In this case, we can see from the xz view that the anomalous response is larger than in the previous result based on one borehole data. Also, in the yz view, not only the left part of the anomaly is resolved, but also its right part is shown very well. 7.2 Model 2 In the next study, we consider a model, which contains two reservoirs, one above the other. Both reservoirs have the same size and shape as in Model 1. The upper reservoir is located at a depth of 0.9 km below the surface, and the lower reservoir is located at 1.9 km below the surface (see Figure 43). The observation system is the same as in Model 1. Figure 44 shows the results of iterative migration of the surface 𝑔𝑧𝑧 component only at the cross section of x = 3000 m. The blue line shows the observed data, and the red line represents the predicted data at the surface (top panel). One can see that the migration density image shows one target only, the lower reservoir cannot be seen at all from the surface data, even though the data fitting is very good, as shown in Figure 45(top panel). Figure 45 presents the results of iterative migration of the surface 𝑔𝑥𝑧 component only at the cross section of x = 3000 m. We use the same notations in this figure as above. Once again, the migration density image outlines one target only confirming the fact that the surface data do not have enough vertical resolution to reveal two stacked reservoirs. Figure 46 shows similar results of the iterative migration of the surface 𝑔𝑦𝑧 component only at the cross section of x = 3000 m. Now, we consider the migration of the borehole data. Figure 47 presents the results of iterative migration of the borehole 𝑔𝑧𝑧 data. The blue line shows the observed data and 97 Figure 41. Joint migration results of surface and borehole g 𝑧𝑧 and g 𝑦𝑧 data at x = 3000 m. 98 Figure 42. Joint migration results of surface and borehole g 𝑧𝑧 and g 𝑦𝑧 data at y = 2800 m. 99 Figure 43. Model of two HC reservoirs. The blue dots show the observation stations at the surface. The vertical black dots denote the position of the borehole. 100 Figure 44. The result of iterative migration of surface 𝑔𝑧𝑧 component only at the cross section of x = 3000 m (bottom panel). The blue line shows the observed data, and the red line shows the predicted data at the surface (top panel). 101 Figure 45. The result of iterative migration of surface 𝑔𝑥𝑧 component only at the cross section of x = 3000 m (bottom panel). The blue line shows the observed data, and the red line shows the predicted data at the surface (top panel). 102 Figure 46. The result of iterative migration of surface 𝑔𝑦𝑧 component only at the cross section of x = 3000 m (bottom panel). The blue line shows the observed data, and the red line shows the predicted data at the surface (top panel). 103 Figure 47. The result of iterative migration of surface 𝑔𝑧𝑧 component only at the cross section of x = 3000 m (left panel). The blue line shows the observed data, and the red line shows the predicted data at the surface (right panel). 104 the red line represents the predicted data in the borehole (right panel). The migration density image shows that it is possible to resolve two reservoirs from the borehole data. However, the artifacts on the opposite side from the borehole complicate the image. We have a similar result shown in Figure 48 for iterative migration of the borehole 𝑔𝑦𝑧 data. The density image resolves two reservoirs well, but the artifacts appear in the opposite direction from the borehole A similar picture is shown in Figure 49, where iterative migration was applied to the borehole 𝑔Δ data In the next step of our numerical study, we consider a joint migration of the surface and borehole gravity data. Figure 50 shows the result of the iterative migration of the surface 𝑔𝑦𝑧 component and borehole 𝑔𝑦𝑧 component jointly. The migration density image reconstructs two reservoirs clearly enough in this figure. Finally, Figure 51 presents the results of joint iterative migration of the surface and borehole 𝑔𝑧𝑧 + 𝑔𝑦𝑧 data. The migration transformation in this case provides a clear image of the two targets without any artifacts. 7.3 Joint migration using multiple borehole data The goal of the final model study is to demonstrate that by using more boreholes one can increase the resolution of migration imaging. In this case, we have simulated the gravity data for Model 2 with two boreholes (Figure 52). Figure 53 presents the results of joint iterative migration of the surface 𝑔𝑧𝑧 + 𝑔𝑦𝑧 data and the same data collected inside two boreholes shown in Figure 52. One can see that the migration density image is improved in this figure in comparison with the previous 105 Figure 48. The result of iterative migration of surface 𝑔𝑦𝑧 component only at the cross section of x = 3000 m (left panel). The blue line shows the observed data, and the red line shows the predicted data at the surface (right panel). 106 Figure 49. The result of iterative migration of surface 𝑔Δ component only at the cross section of x = 3000 m (left panel). The blue line shows the observed data, and the red line shows the predicted data at the surface (right panel). 107 Figure 50. The result of joint iterative migration of surface and borehole 𝑔𝑦𝑧 component only at the cross section of x = 3000 m (bottom left panel). The blue line shows the observed data, and the red line shows the predicted data at the surface (top panel) and in the borehole (right panel). 108 Figure 51. The result of joint iterative migration of surface and borehole 𝑔𝑧𝑧 + 𝑔𝑦𝑧 data only at the cross section of x = 3000 m (bottom left panel). The blue line shows the observed data, and the red line shows the predicted data at the surface (top panel) and in the borehole (right panels). 109 Figure 52. Model of two HC reservoirs and two boreholes. The blue dots show the observation stations at the surface. The vertical black dots denote the position of the boreholes. 110 Figure 53. The result of the joint iterative migration of the surface and borehole 𝑔𝑧𝑧 + 𝑔𝑦𝑧 data at the cross section of x = 3000 m. 111 model study where we used the gravity data from one borehole only. The study shows that by using multiple components of the surface and borehole gravity gradient data one can image the deep targets very well. 7.4 Conclusions We have developed a novel approach to interpretation of the gravity gradiometry data by considering the joint iterative migration of the data observed both on the surface and in the borehole. According to the results from two different models, joint migration of the surface and borehole data can improve the resolution, especially in the z direction, which can help us to locate the lower anomalous body. Furthermore, using multiple components gzz and gyz, can help to represent more accurately the horizontal location of the target as well. And finally, the use of the data from additional borehole improves the resolution of the migration image significantly. Our results also show that by migrating the borehole data jointly with the surface data, we obtain a vertical resolution of the inversion, which would be otherwise impossible to achieve with surface observations only. Joint iterative migration can significantly help to locate the lower target, which has been covered by the upper layers. This research shows the importance of developing gravity gradiometry systems capable of measuring the gravity tensor field in the borehole. 112 CHAPTER 8 CASE STUDY Geophysical monitoring of carbon dioxide (CO2) injections in a deep reservoir has become an important component of carbon capture and storage (CCS) projects. Until recently, the seismic method was the dominant technique used for reservoir monitoring. However, the cost of seismic surveys makes this method prohibitive in monitoring sequestration projects where there is not a direct profit. Moreover, some environments present challenges for seismic acquisition, such as in the urban areas. In this chapter, I present a feasibility study of permanent gravity gradiometry monitoring of CO2 sequestration in a deep reservoir using a novel approach involving both borehole and surface measurements. The interpretation is based on joint iterative migration imaging of the surface and borehole data. The advantage of this method is that the surface data provide a good estimate of the horizontal extent of the injection zone, while the borehole data control the depth of the target, which increases the sensitivity and resolution of the method. It illustrates the effectiveness of the gravity gradiometry method by computer simulating CO2 injection monitoring in the Kevin Dome sequestration site in Montana, USA. 113 8.1 Big sky carbon sequestration partnership (BSCSP) The Big Sky Carbon Sequestration Partnership (BSCSP) is one of the seven carbon sequestration programs initiated by the U.S. Department of Energy. The BSCSP is based at Montana State University’s Energy Research Institute and encompasses six regions Montana, Wyoming, Idaho, South Dakota, eastern Oregon, and Washington. Figure 54 shows the distribution the CO2 emission in BSCSP region as a function of the amount and source types of CO2 emission. The region produces more than 120 million tons of CO2 annually. It is estimated that 82% of the CO2 emission in the region is the result of electricity generation due to high dependence on fossil fuel. Two-third of the region’s emission of CO2 is from Montana and Wyoming because they rely on coal to generate electricity, unlike Oregon, Idaho, and Washington, which rely on hydroelectric power generations (NETL 2010). The BSCSP region has four main CO2 geological storage targets: depleted oil and gas reservoirs, unminable coal seams, basalt formations, and deep saline aquifers. The depleted oil and gas reservoirs are mainly located in Montana and Wyoming. It is estimated that these reservoirs have the potential of more than 1.6 billion tons of CO2. The unminable coal seams in the Powder River, Green River, and Hanna Basins could hold up to 12 billion tons of CO2. The Columbia River Basalt Group in Idaho, Oregon, and Washington has the potential of long-term storage ranging from 36 to 148 billion tons. The BSCSP has estimated CO2 storage capability of 220 billion tons in the saline formations in Montana and Wyoming (NETL 2010). The BSCSP has chosen the saline formation for CO2 storage. The saline formations in the region are characterized by porous and permeable layers of sandstone and limestone 114 Figure 54. BSCSP region with the amount and source types of CO2 emission. 115 with poor water quality of more than 10,000 ppm of total dissolved solids (TDS). The extensive saline formation and proximity to CO2 stationary sources and to existing infrastructure in the region made these formations favorable for carbon storage in the BSCSP region. Among many potential CO2 geological structures, the BSCSP has chosen the Kevin Dome. The Kevin Dome is a unique structure that trapped CO2 for tens of millions of years. The naturally trapped CO2 in the Kevin Dome insures good seal and formation rocks compatibility with CO2. Also, the Kevin Dome is geologically like other structures in Montana and Wyoming like Poplar Dome and Cedar Creek Anticline. Therefore, choosing the Kevin Dome as a CO2 sequestration site will add significant knowledge about carbon storage potential (NETL 2010). 8.2 Kevin Dome project The Kevin Dome is a large underground geological feature of more than 750 square miles (1900 km2) located in the north central part of Montana. The Kevin Dome consists of layers of porous rocks underneath thick impermeable rocks that prevent fluid from migrating to surface. Kevin Dome geology has naturally stored CO2 for millions of years in the upper Devonian carbonate Duperow Formation. The top of the dome contains the CO2 that is above the spill point. On the other hand, the flanks are saline saturated that do not contain CO2. It has estimated that the Kevin Dome has the potential to store 1.51 billion tons of CO2 at the flanks of the dome (NETL, 2010). In 2014, BSCSP drilled, cored, and logged two wells and conducted and analyzed 3D seismic surveys in the project area. Using these valuable geological and geophysical information, 3D static geologic model was constructed and will be updated as more data 116 becomes available. The geostatic modelling will be used to simulate CO2 flow and the effects of CO2 injection to reservoir rocks such as dissolution or precipitation of minerals. Also, it will guide infrastructure design and risk evaluation and management. The BSCSP program, as part of a long-term CO2 monitoring, collected CO2 measurements in soil and atmosphere, water samples, and advance imaginary techniques to establish baseline assessments that will help researchers to assess the impact of CO2 injection on the environment and geology of the area (Spangler, 2016). The BSCSP is planning to produce 1 million tons of naturally occurring CO2 in the Duperow Formation in the Kevin Dome. Then, transport it north of the dome in an underground pipeline of 2-inch diameter and approximately 6 miles in length to the injection site. The CO2 will then be reinjected back to the Duperow formation at the edge of the dome (Figure 55). The Nisku and Souris River formations above and below the Duperow formation, respectively, will be tested for additional storage during the process. BSCSP is planning to have five production wells and four monitoring wells around the injection well. CO2 will be injected at a constant rate of 250,000 metric tons of CO2 a year over 4 years leading to 1 million metric tons (1.1 million tons) of CO2 storage over the span of the project (Spangler, 2016). 8.3 Kevin Dome geology The Kevin Dome is geologically well known from the results of exploration conducted by the oil and gas industry; however, little is known about the characteristics of the Duperow formation since it is located at larger depth than the hydrocarbon fields. Therefore, the continuity and properties of the reservoir is poorly known. However, from 117 Figure 55. The Kevin Dome project in BSCSP. The figure shows the naturally occurring CO2 in the Duperow formation and the planned injection zones. 118 the Validation Phase, the BSCSP has confirmed no major faults in the injection area from the acquired 3D and 9-component seismic. In addition, BSCSP confirmed good reservoir properties and continuity over large area from the well logs (Litynski et al., 2009; Zhdanov et al., 2013). The Duperow formation is located at depths ranging from 1000 to 1900 m within the Kevin Dome. Figure 56 shows a simplified sketch of the Kevin Dome geology. The upper Duperow formation is about 90 m in thickness of tight carbonates interbedded with anhydrites, which serves as the primary seal for the middle Duperow reservoir. The Potlatch formation plays as a secondary seal, which consist of anhydrites of 50 m in thickness. Both of these layers have very low permeability ranging from 0.001 to 10 mD and low porosity ranging from 1%-10% (Dai et al., 2014). The core test results showed that the density of the anhydrites in the Potlatch formation is 2.5-2.83 g/cm3 close to the theoretical density of anhydrites of 2.97 g/cm3 , indicating nearly pure anhydrites with poor porosity (Spangler, 2016). Additional seals have been proven at shallower depth from oil and gas reservoirs (NETL, 2010). The middle Duperow consists of carbonate rocks with a thickness ranging from 20 to 58 m. It has high porosity and permeability ranging from 5% to 25% and 1 to 210 mD, respectively. The expected depth of the middle Duperow formation at the injection site is 1100 m (Dai et al., 2014). The Nisku formation consists of limestone rocks with 15 to 23 m in thickness. It has small zones of good porosity and permeability, which will be tested during the CO2 sequestration project. The Souris River formation will also be tested as a potential CO2 storage resource since it is mainly carbonates (Spangler, 2016). Dai and others (Dai et al., 2014) estimated the radius of the CO2 plume after 1 year of injection 119 Figure 56. A simplified sketch of the Kevin Dome geology. 120 and at the end of the project using Monte Carlo simulations. Figure 57 shows their finding after one year of injection where CO2 plume will have a radius of 500 m with a probability greater than 75%. In addition, it can be seen from Figure 57 that after 4 years of injection, the CO2 plume will have a radius of 1000 m with a probability greater than 75%, radius of 1500 m with a probability greater than 50%, and radius of 2000 m with a probability of more than 30%. 8.4 Kevin Dome model study Simulation of the synthetic gravity gradiometry data for the Kevin Dome reservoir using the CEMI software is outlined in the accompanying paper (Han et al., 2018). For gravity and gravity gradiometry monitoring, the data will be collected periodically over the same area. A baseline survey will be conducted before the start of the injection project. Then, the difference in the fields between each future survey and the baseline survey will be calculated. This difference is attributed to changes in the reservoir properties. The change in density is one of the main factors controlling the measured gravity field at the surface. We investigate the relationship between the density and the measured field. The density of the reservoir is determined by the properties of the rock and in-situ fluids properties as described in the following equation (Han et al., 2014): 𝜌 = 𝜌𝑚 (1 − ∅) + ∅(𝑆𝑤 𝜌𝑤 + 𝑆𝑔 𝜌𝑔 ), (8.1) where 𝜌 is the density of the reservoir, 𝜌𝑚 is the rock matrix density, ∅ is the reservoir’s porosity, 𝑆𝑤 and 𝜌𝑤 are the saturation and density of the brine, and 𝑆𝑔 and 𝜌𝑔 are the saturation and density of the supercritical CO2. 121 Figure 57. Estimated radius of CO2 after 1 year of injection and at the end of the Kevin Dome project using Monte Carol simulations. 122 We selected the following parameters for the model: 𝜌𝑚 = 2.67 g/cm3 , ∅ = 50%, 𝑆𝑤 = 100%, 𝜌𝑤 = 1.1 g/cm3 , 𝑆𝑔 = 100%, 𝜌𝑔 = 0.2 g/cm3 . As a result, we calculated the anomalous density for the water-filled area as −0.74 g/cm3 and the anomalous density of the area with injected CO2 gas as −1.24 g/cm3 . The density contrast is negative because the injected CO2 is lower in density that the in-situ brine in the reservoir. 8.4.1 Kevin Dome no leakage model A simplified model of the reservoir is shown in Figure 58. We have considered four different stages of CO2 injections in the reservoir (Figure 59). For each of these stages, we computer simulated the surface and borehole gravity gradiometry data. As an example, Figure 60 presents the surface 𝑔𝑧𝑧 and 𝑔𝑦𝑧 components for the initial stage (before injection) for five different CO2 contents at the x = 50 m profile. Figure 61 shows the borehole 𝑔𝑧𝑧 and 𝑔𝑦𝑧 components for the different CO2 contents. We then calculated the differences between the fields observed at the current stage and the reference fields corresponding to the initial stage before the CO2 injections (see Figures 62 and 63). These data were migrated back toward the location of the reservoir. Figures 64 and 65 present the migration images for all four different stages of the CO2 injections. It clearly shows the propagation of CO2 during the different stages of injection. We also calculated the differences of the gravity gradient data between the different stages of CO2 injections and migrated these difference fields as well (Figure 66). Figure 67 presents the difference fields (components 𝑔𝑧𝑧 and 𝑔𝑦𝑧 ) on the surface along the profile 123 Figure 58. A simplified model of the Kevin Dome no leakage reservoir. 124 Figure 59. Schematic representation of four different stages of CO2 injections in the reservoir. 125 Figure 60. The 𝑔𝑧𝑧 and 𝑔𝑦𝑧 components on the surface along x = 50 m profile for different CO2 contents. 126 Figure 61. The 𝑔𝑧𝑧 and 𝑔𝑦𝑧 components in the borehole for different CO2 contents. 127 Figure 62. The difference between the surface fields gzz and gyz along the profile x = 50 m observed at the current stage and the reference fields corresponding to the initial stage before the CO2 injections. 128 Figure 63. The difference between the borehole fields gzz and gyz observed at the current stage and the reference fields corresponding to the initial stage before the CO2 injections. 129 Figure 64. Images produced by the migration of the fields calculated as the differences between the data observed at the current stage of CO2 injection and at the initial stage before the start of the injection. 130 Figure 65. Data plots produced by the migration of the fields calculated as the differences between the data observed at the current stage of CO2 injection and at the initial stage before the start of the injection. 131 Figure 66. Schematic representation of four different cases. 132 Figure 67. The difference fields (components 𝑔𝑧𝑧 and 𝑔𝑦𝑧 ) on the surface along the profile x = 50 m, representing the changes in the fields between the subsequent phases. 133 x = 50 m, while Figure 68 shows the same fields in the borehole for the four cases representing the changes in the fields between the subsequent phases. We jointly migrated the surface and borehole difference fields to produce the images representing the changes within the reservoir for the different phases of the CO2 injection (Figure 68). These images manifest how the front of the injected CO2 moves from the left to the right. A, B, C, D in Figure 69 stand for four different cases with the time from beginning drill to the end. It can be seen that the anomaly negative density model is moving from the most left part from panel A figure to the most right part from D part, which can indicate the density difference (between CO2 and water) is moving from left to the right. 8.4.2 Kevin Dome leakage model One of the important goals of monitoring the CO2 sequestration process is to prevent a leakage of CO2 from a deep reservoir. In order to study the detectability of leakage of CO2 in the Kevin Dome model, we considered a leakage scenario of CO2 from the main CO2 reservoir with a radius of 2000 m. We assumed that the CO2 escaped through the upper, 90 m thick Duperow formation of tight carbonate interbedded with anhydrite and through the Potlatch formation of 50 m thickness of anhydrites. We modelled the escaped CO2 when it reached the dolomitic limestone in the Madison formation at a depth of about 500 m, forming a relatively small gas-filled (leaking gas) structure located above the main gas structure, about 500 m in diameter (Figures 70 and 71). Figure 72 presents the 𝑔𝑧𝑧 and 𝑔𝑦𝑧 data on the surface along the profile x = 50 m, 134 Figure 68. The difference fields (components 𝑔𝑧𝑧 and 𝑔𝑦𝑧 ) in the borehole for the four cases representing the changes in the fields between the subsequent phases. 135 Figure 69. Images produced by the migration of the fields calculated as the differences between the data observed at the subsequent stages of CO2. 136 Figure 70. 3D view of a simplified model of the Kevin Dome leakage reservoir. 137 Figure 71. A cross section of the leakage model. 138 Figure 72. The 𝑔𝑧𝑧 and 𝑔𝑦𝑧 components on the surface along x = 50 m observed at the current stage and the reference fields corresponding to the initial stage before the CO2 injections of leakage model. 139 while Figure 73 shows the same fields in the borehole, as well as Figures 74 and 75. Comparing Figures 60, 61, 62, and 63, one can clearly see the anomaly related to the CO2 leakage in both the surface and borehole gravity data. Figure 76 shows a comparison of the migration images for the models in the "leakage" and "no leakage" situations. One can clearly see that the joint migration of surface and boreholes 𝑔𝑧𝑧 and 𝑔𝑦𝑧 data helps identify the presence of leakage from the CO2- filled gas reservoir. 8.5 Conclusions The most widely considered approach to carbon capture and storage is the one based on storing CO2 in deep, natural saline reservoirs. An important problem arising in this case is monitoring and verification of the injection process and the long-term geological integrity of the reservoir seal. Thus, geophysical methods of reservoir monitoring should play a critical role in the CCS process. In this chapter, a novel approach to monitoring CO2 sequestration has been proposed, which involves both the borehole and surface measurements of the gravity gradiometry data. Gravity gradiometry data, especially collected both on the surface and within the borehole, may represent an effective indicator for monitoring CO2 injection in deep reservoirs. Computer simulation has shown that the gravity gradiometry data provide a clear indication of the location of the CO2 plume in the underground formation and of the movement of the front of the injected CO2. This technique can also be used for controlling the leakage of CO2 from a deep reservoir. 140 Figure 73. The difference between borehole fields 𝑔𝑧𝑧 and 𝑔𝑦𝑧 observed at the current stage and the reference fields corresponding to the initial stage before the CO2 injections of leakage model. 141 Figure 74. The difference fields (components 𝑔𝑧𝑧 and 𝑔𝑦𝑧 ) on the surface along the profile x = 50 m, representing the changes in the fields between the subsequent phases of leakage model. 142 Figure 75. The difference fields (components 𝑔𝑧𝑧 and 𝑔𝑦𝑧 ) in the borehole for the four cases representing the changes in the fields between the subsequent phases of leakage model. 143 Figure 76. Comparison of the migration images for the leakage model at the 25% gasinjected stage: (a) no leakage; (b) with leakage; (c) difference between the two images. 27 CHAPTER 9 CONCLUSIONS Gravity gradiometry is a novel technique used for oil and mineral exploration, which provides detailed information about the density distribution in the subsurface by measuring the rate of change of gravitational acceleration due to underlying rock properties. From this information it is possible to build a picture of subsurface anomalies, which can then be used to more accurately target oil, gas, and mineral deposits. It is also used to image water column density, when locating submerged objects, or determining water depth (bathymetry). Physical scientists use gravimeters to determine the exact size and shape of the earth, and they contribute to the gravity compensations applied to inertial navigation systems. Gravity gradiometry has predominately been used to image subsurface geology to aid hydrocarbon and mineral exploration. Over 2.5 million-line km has now been surveyed using this technique. The surveys highlight gravity anomalies that can be related to geological features such as salt diapirs, fault systems, reef structures, kimberlite pipes, etc. We have developed a novel approach to the joint interpretation of the surface and borehole gravity and gravity gradient data based on the concept of potential field migration. The results of the numerical model study demonstrated that the potential field migration can be used for 3D imaging of deep-seated HC reservoirs. Our results also show that by 145 migrating the borehole data jointly with the surface data, we obtain a vertical resolution of the inversion, which would be otherwise impossible to achieve with surface observations only. Joint iterative migration can significantly help to locate the lower target that has already been covered by upper layers. This research shows the importance of developing gravity gradiometry systems capable of measuring the gravity tensor field in the borehole. The most widely considered approach to carbon capture and storage is the one based on storing CO2 in deep, natural saline reservoirs. An important problem arising in this case is monitoring and verification of the injection process and the long-term geological integrity of the reservoir seal. Thus, geophysical methods of reservoir monitoring should play a critical role in the CCS process. In this dissertation, we have proposed a novel approach to monitoring CO2 sequestration, which involves both the borehole and surface measurements of the gravity gradiometry data. We have demonstrated that gravity gradiometry data, especially collected both on the surface and within the borehole, may represent an effective indicator for monitoring CO2 injection in deep reservoirs. Computer simulation has shown that the gravity gradiometry data provide a clear indication of the location of the CO2 plume in the underground formation and of the movement of the front of the injected CO2. This technique can also be used for controlling the leakage of CO2 from a deep reservoir. 27 REFERENCES Alixant, J.-L., & Mann, E. (1995). In-Situ residual oil saturation to gas from time-lapse borehole gravity. Paper presented at 1995 SPE Annual Technical Conference and Exhibition, Dallas, Texas. AlJanobi, H. A. (2017). Time-lapse airborne gravity and gravity gradiometry Monitoring of CO2 sequestration in the Kevin Dome, (M.S. thesis). Retrieved from ProQuest https://search.proquest.com/docview/2181693572?pq-origsite=gscholar. Salt Lake City, UT: University of Utah. Brady, J., B. Watson, D. Warner, R. North, D. Sommer, J. Colson, R. Kleinberg, & Sezginer, A. (1998). Improved production log interpretation in horizontal wells using a combination of pulsed neutron logs, quantitative temperature log analysis, time lapse LWD resistivity logs and borehole gravity. Paper presented at 1998 SPE Western Regional Meeting, Bakersfield, California. Brady, J. L., Ferguson, J. F., Hare, J. L., Seibert, J. E., Chen, T., Klopping, F., & Niebauer, T. (2008). Results of the world's first 4D microgravity surveillance of a waterflood-Prudhoe Bay, Alaska. SPE Reservoir Evaluation & Engineering, 11, 824-831. https://doi.org/10.2118/101762-PA Bruant, RG. Jr., Guswa, AJ., Celia, MA., & Peters, CA. (2002). Peer reviewed: Safe storage of CO2 in deep saline aquifiers. Environmental Science and Technology, 36(11), 240A-245A. DOI: 10.1021/es0223325 Caffagni, E., & Bokelmann, G. (2016). Geophysical monitoring of a hydrocarbon reservoir. Energy Procedia, 97, 294-301. https://doi.org/10.1016/j.egypro.2016.10.003 Cao, D., Yin, X., Wu, G., & Zhao, X. (2013). Impedance joint inversion of borehole and surface seismic data. Journal of Geophysics and Engineering, 10(4), 045003. https://doi.org/10.1088/1742-2132/10/4/045003 Dai, Z., Stauffer, P.H., Carey, J.W., Middleton, R.S., Lu, Z., Jacobs, J.F., Telleen, K, & Spangler, L.H. (2014). Pre-site characterization risk analysis for commercial-scale carbon sequestration. Environmental Science & Technology, 48, 3908-3915. https://doi.org/10.1021/es405468p 147 Gasperikova, E., & Hoversten, G. M. (2008). Gravity monitoring of CO 2 movement during sequestration: Model studies. Geophysics, 73(6), WA105-WA112. https://doi.org/10.1190/1.2985823 Gasperikova, E., & Hoversten, G. M. (2006). A feasibility study of nonseismic geophysical methods for monitoring geologic CO 2 sequestration. The Leading Edge, 25, 12821288. https://doi.org/10.1190/1.2360621 Goodell, R.R., & Fay, J.H. (1964). Borehole gravity meter and its application. Geophysics, 29, 774-782. https://doi.org/10.1190/1.1439417 Gournay, L. S. (1983). Applications of borehole gravimetric techniques to determine residual oil saturation. Google Patents. New York, NY: Mobil Oil Corporation. Han, M., Alshehri, A. J., Krinis, D., & Lyngra, S. (2014). State-of-the-art of in-depth fluid diversion technology: Enhancing reservoir oil recovery by gel treatments. Paper presented at 2014 SPE Saudi Arabia section technical symposium and exhibition. Society of Petroleum Engineers, Al-Khobar, Saudi Arabia. Han, M., Wan, L., & Zhdanov, M. S. (2018). Joint iterative migration of surface and borehole gravity gradiometry data. SEG Technical Program Expanded Abstracts 2018, (pp. 1399-1403), Anaheim, CA. Hansen, P. C. (2005). Rank-deficient and discrete ill-posed problems: numerical aspects of linear inversion. (Vol. 4) (Rep. ISBN. 0898714036). Lyngby, Denmark: Technical University of Denmark Huang, T., Quach, N., Abercrombie, S. P., Boening, C., Brennan, H. P., Gill, K. M., et al. (2017). NASA sea level change portal-it not just another portal site. Paper presented at 2017 AGUFM, IN33E-05, San Francisco, CA. Jageler, A.H. (1976). Improved hydrocarbon reservoir evaluation through use of boreholegravimeter data. Society of Petroleum Engineers, 28, 709-718. https://doi.org/10.2118/5511-PA Jessop, M. & Zhdanov, M. (2005). Numerical study of gravity gradiometer data for typical kimberlites in the Northwest Territory of Canada. Proceedings of 2005 CEMI Annual Meeting. Salt Lake City, UT: Consortium for Electromagnetic Modeling and Inversion, University of Utah Krahenbuhl, R. A. & Li, Y. (2008). Joint inversion of surface and borehole 4D gravity data for continuous characterization of fluid contact movement. SEG Technical Program Expanded Abstracts 2008, (pp. 726-730), Las Vegas, NV. 148 Krieger, M. H., Smilde, P. L., & Geisler, O. K. (2009). Completing the image with borehole gravity gradients. SEG Technical Program Expanded Abstracts 2009, (pp. 923927), Houston, TX. LaFehr, T. R. (1983). Rock density from borehole gravity surveys. Geophysics, 48, 341356. https://doi.org/10.1190/1.1441472 Lehujeur, M., Vergne, J., Schmittbuhl, J., & Maggi, A. (2014). Investigating a deep geothermal reservoir using ambient noise correlation. Paper presented at 2014 EGU general assembly conference, Vienna, Austria, id.13798. Lehujeur, M., Vergne, J. Schmittbuhl, J., & Maggi, A. (2015). Characterization of ambient seismic noise near a deep geothermal reservoir and implications for interferometric methods: a case study in northern Alsace, France. Geothermal Energy, 3, 3. https://doi.org/10.1186/s40517-014-0020-2 Li, Y., & Oldenburg, D. W. (2000). Joint inversion of surface and three-component borehole magnetic data. Geophysics, 65(2), 540-552. https://doi.org/10.1190/1.1444749. Litynski, J., Plasynski, S., Spangler, L., Finley, R., Steadman, E., Ball, D., Nemeth, K. J., McPherson, B., & Myer. L. (2009). US Department of Energy’s regional carbon sequestration partnership program: Overview. Energy Procedia, 1, 3959-3967. https://doi.org/10.1016/j.egypro.2009.02.200 MacQueen, J. D. (2007). High-resolution density from borehole gravity data. SEG Technical Program Expanded Abstracts 2007, (pp. 741-744), San Antonio, TX. McCulloh, T. H., Kandle, J. R., & Schoellhamer, J. E. (1968). Application Of gravity measurements in wells to problems or reservoir evaluation. Paper presneted at 9th SPWLA 9th Annual Logging Symposium, New Orleans, LA. Nekut, A. G. (1989). Borehole gravity gradiometry. Geophysics, 54, 225-234. https://doi.org/10.1190/1.1442646 NETL, USDOE. (2010). Carbon sequestration atlas of the United States and Canada (3rd Edition). U.S. Department of Energy • Office of Fossil Energy National Energy Technology Laboratory. Obermann, A., Kraft, T., Larose, E., & Wiemer, S. (2015). Potential of ambient seismic noise techniques to monitor the St. Gallen geothermal site (Switzerland). Journal of Geophysical Research: Solid Earth, 120(6), 4301-4316. https://doi.org/10.1002/2014JB011817 Portniaguine, O., & Zhdanov, M. S. (1999). Focusing geophysical inversion images. Geophysics, 64(3), 874-887. https://doi.org/10.1190/1.1444596 149 Rasmussen, N. F. (1975). The successful use of the borehole gravity meter in northern Michigan. The Log Analyst, 16(05). DOI: SPWLA-1975-vXVIn5a1 Rim, H., & Li, Y. (2010). Single-borehole imaging using gravity gradiometer data. SEG Technical Program Expanded Abstracts 2010 (pp. 1137-1141), Denver, CO. Rubin, E., & De Coninck, H. (2005). IPCC special report on carbon dioxide capture and storage. UK: Cambridge University Press. TNO (2004): Cost Curves for CO2 Storage, Part 2, 14. http://www.rite.or.jp/English/lab/geological/geowse/20-31_Rubin.pdf Seigel, H. O., Nind, C. J. M., Milanovic, A., & MacQueen, J. (2009). Results from the initial field trials of a borehole gravity meter for mining and geotechnical applications. Paper presented at 11th SAGA Biennial Technical Meeting and Exhibition, Swaziland, South Africa. Smith, N. J. (1950). The case for gravity data from boreholes. Geophysics, 15(4), 605-636. https://doi.org/10.1190/1.1437623 Spangler, L. (2016). Big Sky regional carbon sequestration partnership–Kevin Dome carbon storage, Montana, U.S.: National Energy Technology Laboratory. Sun, J., & Li, Y. (2010). Inversion of surface and borehole gravity with thresholding and density constraints. SEG Technical Program Expanded Abstracts 2010 (pp. 17981803), Denver, CO. Popta, J. V., Heywood, J. M. T., Adams, S. J., & Bostock, D. R. (1990). Use of borehole gravimetry for reservoir characterisation and fluid saturation monitoring. Paper presented at European Petroleum Conference, The Hague, Netherlands. Wan, L., Han, M., AlJanobi, H. A., & Zhdanov, M. S. (2020). Feasibility study of gravity gradiometry monitoring of CO2 sequestration in deep reservoirs using surface and borehole data. Active Geophysical Monitoring (pp. 123-140). Salt Lake City, UT: Elsevier. Wan, L., Han, M., & Zhdanov, M. (2016). Joint iterative migration of surface and borehole gravity gradiometry data. SEG Technical Program Expanded Abstracts 2016 (pp. 1607-1611), Dallas, TX. Wan, L., & Zhdanov, M. S. (2013). Iterative migration of gravity and gravity gradiometry data. SEG Technical Program Expanded Abstracts 2013 (pp. 1211-1215), Houston, TX. Zhdanov, M. S. (2015). Geophysical inverse theory and regularization problems. Inverse theory and applications in geophysics. Elsevier. 150 Zhdanov, M. S., Endo, M., Black, N., Spangler, L., Fairweather, S., Hibbs, A., et al. (2013). Electromagnetic monitoring of CO2 sequestration in deep reservoirs. First Break, 31(2). Zhdanov, M. S., Gribenko, A., & Čuma, M. (2007). Regularized focusing inversion of marine CSEM data using minimum vertical-support stabilizer. SEG Technical Program Expanded Abstracts 2007 (pp. 579-583), San Antonio, TX. Zhdanov, M. S., Han, M., & Wan, L. (2020). Joint iterative migration of surface and borehole gravity gradiometry data. Active Geophysical Monitoring (pp. 97-121). Salt Lake City, UT: Elsevier. https://doi.org/10.1016/B978-0-08-102684-7.000054 Zhdanov, M. S., Liu, X., & Wilson, G. A. (2010). Rapid imaging of gravity gradiometry data using 2D potential field migration. SEG Technical Program Expanded Abstracts 2010 (pp. 1132-1136), Denver, CO. Zhdanov, M. S., Liu, X., Wilson, G. A., & Wan, L. (2011). Potential field migration for rapid imaging of gravity gradiometry data. Geophysical Prospecting, 59, 10521071. https://doi.org/10.1111/j.1365-2478.2011.01005.x Zhdanov, M. S., Liu, X., & Wilson, G. (2010). Potential field migration for rapid 3D imaging of entire gravity gradiometry surveys. First Break, 28(11), 47-51. https://www.researchgate.net/publication/265493096_Potential_field_migration_f or_rapid_3D_imaging_of_entire_gravity_gradiometry_surveys Zhu, S., Meng, Q., Wang, L., Zhang, J., Song, Y., Jin, H., et al. (2013). Highly photoluminescent carbon dots for multicolor patterning, sensors, and bioimaging. Angewandte Chemie International Edition, 52(14), 3953-3957.