Advanced methods for depth-to-basement estimation using gravity, magnetic, and electromagnetic data

There is a strong interest in developing effective methods to estimate the depth-to-basement. Potential field methods have already been widely used in this application by parameterizing the earth's subsurface into 3D cells. I introduce a new method of solving this problem based on the 3D Cauchy-type integral (CTI) method which makes it possible to represent the potential fields as surface integrals and the density or magnetization contrast surface needs to be discretized only for the calculation of the potential fields. Another significant objective is the development of a novel method for inversion of potential field data to recover the depth-to-basement using 3D Cauchy-type integral representation. Numerical studies show that the new method is much faster than the conventional method to compute the potential field. My synthetic model studies also show that the developed inversion algorithm is capable of recovering the geometry and depth of a sedimentary basin effectively with a complex density profile in the vertical direction. By nature, the recovered model from potential field inversion is usually very diffusive. Under these circumstances, one has to consider some other geophysical methods, such as electromagnetic (especially the magnetotelluric) methods, which have higher resolution and acceptable exploration cost. Conventional inversion of magnetotelluric (MT) data is aimed at determining the volumetric conductivity distribution. This dissertation develops a novel approach to 3D MT inversion for the depth-to-basement estimation. The key to this approach is selection of the depth-to-basement being the major unknown parameter. The inversion algorithm recovers both the thickness and the conductivities of a sedimentary basin. The sediment-basement interface is usually characterized by density, magnetization, and electrical conductivity contrasts. This makes realistic the joint inversion of potential field and MT data to recover the depth-to-basement. I have developed a joint inversion algorithm to recover the depth to-basement using MT and gravity data. The inversion can recover the physical properties and electrical conductivity simultaneously. The developed methods are illustrated on several realistic geological models. They have also been applied to the USGS field gravity data and synthetic MT data for the Big Bear Lake area.

ADVANCED METHODS FOR DEPTH-TO-BASEMENT ESTIMATION USING GRAVITY, MAGNETIC, AND ELECTROMAGNETIC DATA by Hongzhu Cai 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 December 2015 Copyright © Hongzhu Cai 2015 All Rights Reserved The Unive r si t y of Utah Graduat e School STATEMENT OF DISSERTATION APPROVAL The dissertation of Hongzhu Cai has been approved by the following supervisory committee members: Michael Zhdanov Chair 10/05/2015 Date Approved Keith Koper Member 08/31/2015 Date Approved Michael Thorne Member 08/19/2015 Date Approved Alexander Gribenko Member 09/14/2015 Date Approved Le Wan Member 08/31/2015 Date Approved and by John Bartley Chair/Dean of the Department/College/School o f _____________Geology and Geophysics and by David B. Kieda, Dean of The Graduate School. ABSTRACT There is a strong interest in developing effective methods to estimate the depth-to-basement. Potential field methods have already been widely used in this application by parameterizing the earth's subsurface into 3D cells. I introduce a new method of solving this problem based on the 3D Cauchy-type integral (CTI) method which makes it possible to represent the potential fields as surface integrals and the density or magnetization contrast surface needs to be discretized only for the calculation of the potential fields. Another significant objective is the development of a novel method for inversion of potential field data to recover the depth-to-basement using 3D Cauchy-type integral representation. Numerical studies show that the new method is much faster than the conventional method to compute the potential field. My synthetic model studies also show that the developed inversion algorithm is capable of recovering the geometry and depth of a sedimentary basin effectively with a complex density profile in the vertical direction. By nature, the recovered model from potential field inversion is usually very diffusive. Under these circumstances, one has to consider some other geophysical methods, such as electromagnetic (especially the magnetotelluric) methods, which have higher resolution and acceptable exploration cost. Conventional inversion of magnetotelluric (MT) data is aimed at determining the volumetric conductivity distribution. This dissertation develops a novel approach to 3D MT inversion for the depth-to-basement estimation. The key to this approach is selection of the depth-to-basement being the major unknown parameter. The inver­sion algorithm recovers both the thickness and the conductivities of a sedimentary basin. The sediment-basement interface is usually characterized by density, magnetiza­tion, and electrical conductivity contrasts. This makes realistic the joint inversion of potential field and MT data to recover the depth-to-basement. I have developed a joint inversion algorithm to recover the depth-to-basement using MT and gravity data. The inversion can recover the physical properties and electrical conductivity simultaneously. The developed methods are illustrated on several realistic geological models. They have also been applied to the USGS field gravity data and synthetic MT data for the Big Bear Lake area. iv To my parents CONTENTS A B S T R A C T ............................................................................................................... iii A CK N OW L E D GM EN T S ...................................................................................... viii 1..... IN T R O D U C T IO N ........................................................................................... 1 2. 3D ANALOGS OF C AU CH Y -TY PE INTEGRALS AND THEIR P R O P E R T IE S .................................................................................................... 6 2.1 Definition........................................................................................................ 6 2.2 Major theorems............................................................................................. 8 3. REPRESENTATION OF POTENTIAL FIELDS USING 3D ANALOGS OF C AU CH Y -TY PE IN T E G R A L S ................................ 11 3.1 Representation of the gravity field using C T I ......................................... 12 3.2 Representation of the magnetic field using C T I ...................................... 14 3.3 CTI representation of the gravity field of the density-contrast surface 15 3.4 CTI representation of the magnetic field of the magnetic-contrast surface 21 4. 3D FORWARD MODELING OF M T D A T A ...................................... 29 5. INVERSION OF POTENTIAL-FIELD AND M T DATA FOR DE P TH-TO-BASEMEN T ANALYSIS .................................................. 32 5.1 Inversion of the potential-field data for the depth-to-basement ........... 32 5.2 Inversion of the magnetotelluric data for the depth-to-basement ......... 35 5.3 Inversion of the geophysical fields for the depth-to-basement and for the physical properties of the sediment .................................................... 36 5.4 Joint inversion of different geophysical data for the depth-to-basement estimation ...................................................................................................... 38 6. MODEL STU D IE S ........................................................................................... 41 6.1 Modeling and inversion of gravity data for the sediment-basement interface model ............................................................................................. 41 6.2 Modeling and inversion of magnetic data for the sediment-basement interface model ............................................................................................. 50 6.3 Modeling and inversion of MT data for the sediment-basement inter­face model ...................................................................................................... 52 6.4 Joint inversion of MT and potential-field data for depth-to-basement estimation ...................................................................................................... 55 7. CASE S TU D IE S ............................................................................................... 96 7.1 Inversion of gravity data at the Big Bear Lake A rea ............................. 96 7.2 Inversion of the MT data for the USGS basin model............................. 101 8. CON C LU S ION S ............................................................................................... 123 A PPEN DIX : FRECHET DERIVATIVE C A LCU LA T IO N .................... 126 REFE R EN C E S.......................................................................................................... 135 vii ACKNOWLEDGMENTS I would like to express my deep gratitude to my graduate advisor, Michael Zh­danov, for his support and guidance in my research and graduate study. I have benefited greatly from his expertise in geophysics and mathematics. Without his help during my study and research, I could not have finished my research. I extend special acknowledgments to Dr. Gribenko and Dr. Cuma for helping with the MT inversion and to Dr. Wan for helping with the gravity inversion during the past years. I am also grateful to my supervisory committee members, Professor Keith Koper and Professor Michael Thorne, for their useful suggestions with my research. I am deeply thankful to the Consortium for Electromagnetic Modeling and Inver­sion (CEMI), which provides an excellent environment for study and research. In the CEMI student office, I met some great friends with whom to work and study. I am also thankful to all the other CEMI researchers and students, especially Daeung Yoon and Yue Zhu, for their help with my research during the past six years. They gave me a lot of suggestions during my study here and always encouraged me to achieve more in my life. Finally, I would like to give my special thanks to my family for their continuous support and encouragement of my study. CHAPTER 1 INTRODUCTION In exploring the oil and gas potential of a region, there is a strong interest in developing effective geophysical methods for the estimation of the depth-to-basement and the physical properties of the sediments. Three main geophysical methods are used for this application: seismic, potential field, and electromagnetic (Afnimar and Nakagawa, 2002; Zevallos et al., 2004; Cai and Zhdanov, 2015a, 2015b). Different geophysical methods operate with different physical properties, such as the velocity, density, magnetic susceptibility, and conductivity, of rocks. Among these methods, potential-field methods (gravity and magnetic) have been widely used for this appli­cation. Due to the process of sedimentation, usually the sediments are less dense and weakly magnetized, whereas the basement rocks are characterized by higher density and stronger magnetization. In this situation, regional gravity and magnetic anomaly can be attributed to variations in the depth-to-basement. As such, the depth to a crystalline basement can be estimated from the gravity and magnetic field data by assuming that the densities of the sediments and basement rock are given or well constrained from some other geophysical information. A large number of research papers have been published over the last decade on this subject (e.g., Barbosa et al., 1997, 1999a, b; Gallardo-Delgado et al., 2002; Martins et al., 2010, 2011a, b; Silva et al., 2001, 2006, 2007, 2010a, b). The conventional approach to solving depth-to-basement potential-field inverse problems is based on parameterization of the earth's subsurface, containing the sedimentary pack, into prismatic cells with known horizontal dimensions and known density/susceptibility contrasts, and on estimation of the thicknesses of the cells. This dissertation presents a novel approach to the solution of this problem based on 3D analogs of Cauchy-type integrals, introduced by Zhdanov (1980, 1984, 1988). These integrals extend to a 3D 2 case all the major properties of the classical Cauchy integrals of the theory of functions of complex variables. In a 2D case, Cauchy integrals can be used to provide an effective representation of the potential field of 2D density/susceptibility distributions and to solve the problems of the upward and downward analytical continuation of the potential-field data. It was demonstrated in the papers by Zhdanov, cited above, that 3D analogs of Cauchy-type integrals make it possible to extend to 3D cases a large body of the research developed for 2D potential fields. For example, in the paper by Zhdanov and Liu (2013), 3D Cauchy-type integrals were applied to solve the problem of terrain correction for gravity and gravity gradiometry data. In this dissertation, I apply the method of 3D Cauchy-type integrals to solving both forward and inverse problems for a sediment-basement interface model. This type of model is used, for example, in the inversion of potential data for the depth-to-basement. In this study, I consider a model formed by two quasi-horizontal layers, the upper layer representing the sediments and the lower layer describing the basement. I assume that the density/susceptibility does not vary in the horizontal direction, but, in a general case, may vary vertically, having a discontinuity at the surface of the basement. The goal is to find the surface of the basement, which is a density-contrast surface in this case, for example. I have developed an inversion scheme to identify the sediment-basement interface using the gravity field, full-tensor gravity-gradiometry data, and magnetic data. The inversion scheme is based on the reweighted regularized conjugate-gradient method (Zhdanov, 2002). Note that the method based on Cauchy-type integrals requires the discretization of the contrast surface only, which reduces dramatically the com­puting resources in comparison with the conventional methods based on the volume discretization of the subsurface into prismatic cells. The conventional methods of solving this problem are based on the spectrum and/or Euler deconvolution analysis of the gravity field. However, these methods have several limitations. In the framework of the spectrum method, the gravity field needs to be analyzed within a moving window, and the size of the window needs to be determined based on an expected depth to the source, which is usually either unavailable or not accurate (Chavez et al., 1999). A complex source structure 3 can complicate the spectrum analysis, which may result in significant errors of the depth estimation (Odegard, 2011). The conventional Euler deconvolution method can be used for fast depth estimation, but it requires input of the source structure index to estimate the depth to the source, which simplifies the source to some specific geometries, such as sphere, cylinder, etc. (Lafehr and Nabighian, 2012). Even though an extended Euler deconvolution method can be used to estimate both the source depth and the structure index simultaneously, it is still difficult to deal with a complex source structure (Lafehr and Nabighian, 2012). Moreover, both the spectrum and Euler deconvolution methods do not provide a direct comparison between the observed and predicted gravity-field data, which makes it difficult to evaluate the correctness of the solution for the depth-to-basement. In comparison to those methods, my method is based on direct evaluation of the misfit between the observed and predicted data. In addition, as we will see below, the a priori information can also be incorporated into the inversion in the framework of the method based on the Cauchy-type integrals. The developed method was tested for inversion of the gravity data computer-simulated for typical contrast-surface models. I have also applied this method to the USGS field gravity data in the Big Bear Lake area in California to recover the depth-to-basement. It is well known that potential-field inversion is inherently diffusive and nonunique, which makes the depth-to-basement estimation using potential-field data a very chal­lenging problem. In order to effectively solve this problem, one can use other geo­physical methods which can provide a higher resolution of the depth-to-basement. It is well known that electromagnetic (EM) data can provide higher resolution for subsurface formation than gravity and magnetic data due to the frequency dependence of the EM field and the depth of investigation (Zhdanov, 2009). The magnetotelluric (MT) method is a passive-source EM method based on the analysis of the impedance tensor defined by the ratio of the electric and magnetic fields (Zhdanov and Keller, 1994; Berdichevsky and Dmitriev, 2008). It has been demonstrated that this method can be used in exploration on local and regional scales (Zhdanov and Keller, 1994). The MT method also provides an effective approach for sedimentary basin analysis such as depth-to-basement estimation based on the conductivity contrast between the sediments and bedrocks (Zevallos et al., 2004; Tournerie and Chouteau, 2005). 4 Conventional inversions of the MT data are usually aimed at determining the vol­umetric distribution of the conductivity within the inversion domain (Berdichevsky and Dmitriev, 2008; Zhdanov, 2002, 2009). By the nature of the MT method, the recovered distribution of the subsurface conductivity is typically diffusive, although it can be focused by adopting more advanced regularization schemes such as focusing stabilizers (Zhdanov, 2002). In the problem of depth-to-basement estimation using geophysical data, the goal is to recover a sharp boundary between a sedimentary basin and a crystalline basement. Therefore, we need to adopt a sharp boundary parameterization of the subsurface for the inversion. This type of parametrization has been widely used for potential-field in­version to recover the depth-to-basement. Gallardo-Delgado et al. (2003) discretized the subsurface with regular columns where the bottom of each column represents the sediment-basement interface. However, these parameterization approaches have not been widely used for MT inversion problems, especially in a 3D case. The main reason is that it is difficult to find a straightforward relationship between electromagnetic data and the thickness of the columns in a sharp boundary parameterization. Chen et al. (2012) introduced a sharp boundary parameterization for a 2D MT inverse problem. In their inversion, the subsurface was discretized into several layers of arbitrary shape. They used a stochastic inversion algorithm to recover the location of the boundaries and also the conductivities within each layer. However, this method cannot be easily extended to a 3D case because of the high computational cost of stochastic inversion. In this dissertation, we suggest using a column parameterization for the MT inversion, same as the discretization used in Gallardo-Delgado et al. (2003) for potential field inversion. Similarly to the formulation of a potential-field problem, we assume that the subsurface comprises a conductive layer of sediments and a resistive bedrock foundation. The interface between the sediments and the bedrock has an arbitrary shape. The sediment packs are discretized into a grid of columns with known horizontal dimensions. The MT response of the geoelectrical model is computed using the integral equation method (the CEMI code PIE3D). In my inversion, I calculate the Frechet derivatives of the data with respect to the thickness of the columns and the 5 conductivity of the sediments using the quasi-Born approximation. Realistic model studies show that the developed method can be used for fast and accurate estimation of the depth-to-basement using MT data. I have also implemented a joint inversion algorithm for the estimation of the depth-to-basement and the physical properties of the sediments using potential-field (gravity) and MT data. In a general case, the density/susceptibility contrast interface corresponds to the electrical conductivity contrast as well. Due to the frequency dependence of the electromagnetic field, we expect that the depth to crystalline basement and the sediment conductivity can be well constrained from the MT data. By joint inversion of MT and gravity-field data, the density of the sediments can also be recovered. The joint inversion algorithm has been tested using several realistic synthetic models and the synthetic MT data computer-generated for the USGS model built from 3D inversion of gravity data based on the Cauchy-type integrals. The numerical studies demonstrate that the developed algorithm can be used for detailed and high-resolution imaging of the sedimentary basin. CHAPTER 2 3D ANALOGS OF CAUCHY-TYPE INTEGRALS AND THEIR PROPERTIES In 2D applications, one can effectively use the theory of the functions of a complex variable to study the potential-field problems. For example, this theory can be used for the upward and downward continuation of potential-field data. Zhdanov et al. (2012) has applied this method to the migration of magnetic vector and magnetic tensor data for the fast imaging of subsurface magnetic susceptibility using the technique of the classical Cauchy integrals. It has been shown by Zhdanov (1988) that most of the properties of the classical 2D Cauchy integrals can be extended to a 3D case in the form of 3D Cauchy-type integral analogs. In this section, I will first introduce the definition of 3D analogs of Cauchy type integrals and discuss the major theorems related to these integrals as well as the basic properties which can be used for solving 3D potential-field problems. The relationship between the 3D analog of a Cauchy-type integral and the classical Cauchy integral formula will also be discussed in this section. 2.1 Definition The 3D analog of a Cauchy-type integral has been introduced and studied by Zhdanov (1988). The detailed derivation of the 3D analog of a Cauchy-type integral is beyond the scope of this research. However, I will present here the basic definitions and derivations of 3D Cauchy-type integrals for completeness. We define a vector function $ (r ) in domain D as the following form: $ (r ) = (C ■ Vh) V P + V P x (Vh x C) (2.1) 7 where h(r) and P(r) are arbitrary functions twice differentiable in domain D (in­cluding its boundary S) and C is an arbitrary constant vector. It is easy to find that V - $ = C- (APVh + A hV P ) , (2.2) $ ■ n = C- [(n -V P ) Vh + (n x V P ) x Vh] = C- [(n -Vh) V P + (n x Vh) x V P ] . (2.3) I apply the Stokes theorem to the vector function $ (r ) and consider the fact that the constant vector C is an arbitrary vector. I arrive at the following two corollaries of the Ostrogradsky-Gauss theorem (Zhdanov, 1988): (APVh + A hV P ) dv = JJ [(n ■ V P ) Vh + (n x V P ) x Vh] ds, (2.4) (APVh + A hV P ) dv = / / [(n ■ V P ) Vh + (n x V P ) x Vh] ds. (2.5) id JJs We consider that P is the fundamental Green's function for the Laplace equation: D P = G(r - r') = - 1 4n|r - r'| ' A P = AG(r - r') = 6 (r - r ' ) , (2.6) (2.7) Ah = 0, r € D, (2.8) where 6 is the Dirac delta function. Then we can obtain the following equation: - 1 4n (n ■ Vh)V 1 + (n x Vh) x V 1 ds Vh (r ) , r 1 D , (2.9) 0, r € CD. V ' The classic Cauchy integral theorem for complex analysis in a 2D application states that if f (Z) is an analytical function defined in domain D bounded by the closed contour L and continous on the contour, there exists the following relationship: f f (Z) dz = 2ni Jl Z - Z' Z f (Z ') ,Z ' € D 0, Z' € CD (2.10) From equation (2.10), we can see that the value of the complex analytical function inside domain D can be calculated from its value on the closed contour which bounds domain D. S 8 It is obvious that the relationship in equation (2.9) is structurally similar to the Cauchy integral formula (2.10). As a result, equation (2.9) enables us to compute the values of the fields, defined by the gradient of a differentiable function h(r), inside domain D from their values on the boundary (S) of the domain (Zhdanov, 1988). Formula (2.9) is called a three-dimensional analog of the Cauchy integral formula. We further introduce a vector function <^(r) defined by the gradient of the scalar function h(r), <^(r) = Vh(r). (2.11) The 3D analog of the Cauchy-type integral for the vector function <^(r) then can be defined as follows: C s ( r > ) = - 1 ,. 4n J js (n ■ (p)V |------- r. + (n x p) x V ds, (2.12) |r - r' | |r - r' | where S is some closed surface bounding a domain D, <p = <^(r) is some vector function defined on the closed surface S, and n is the normal vector to the surface S, pointing outside D. The vector function <p is called the vector density of the Cauchy-type integral. 2.2 Major theorems It was demonstrated in Zhdanov (1988) that everywhere outside of S , the vector function Cs represents the Laplace vector field, which satisfies the following equations: V - Cs = 0, V x C s = 0. (2.13) Thus, the scalar components of vector function Cs are harmonic functions. In a special case where <^(r) stands for the boundary values on S of the gradient of a function harmonic inside domain D, a 3D Cauchy-type integral can be calculated using the following formula: C V , . ) = { $ >e c D D , (214) where CD is a complement of the closed domain D with respect to the whole space. 9 Now, we consider a vector field F(r) which represents the potential field and satisfies the following equations: V ■ F = q, V x F = 0, r eD. (2.15) We now set the following: F = Vh,P = - , 1 „ 4n|r - r'| (2.16) We then reconsider equations (2.4) and (2.5). After repeating a similar deriva­tion for the 3D analog of a Cauchy-type integral, we can obtain the following two relationships: - 1 4n J Js 1 4n (n ■ F)V 1 |r - r'| + (n x F) x V 1 ds + 1 d V -F V |r - r | dv |r - r'| F fr') , r E D 0, r e CD; (2.17) - 1 4n J Js 1 4n (n ■ V,- ^ ) F + (n x V,- ) x F |r - r | |r - r | V - FV - d |r - r | dv ds + F (r j , r E D, 0, r e CD. (2.18) Formulas (2.17) and (2.18) are just 3D analogs of the Pompei formulas (Zhdanov, 1988). These formulas yield a solution of the boundary value problem for potential fields: F (r ) = -1 4n J js (n ■ F)V-- + (n x F) x V 1 1 + "-- 4n |r - r' | / / / d « <r)Vi r i 7 i * '■r' eD |r - r' | ds (2.19) F i r > = T4n1 Jm js (n ■ V-.- )F + (n x V-.- ) x F |r - r'|' | r - r' |' ds + 4n J J J d «<r)V 1 |r - r'| dv, r ED. (2.20) We can see that if the vector field F (r) is a Laplace field in domain D (V ■ F = 0, V x F = 0), formula (2.19) yields a 3D analog of the Cauchy-type integral in formula (2.9). 10 In a compact form, one can formulate a 3D analog of the Pompei formula for the Cauchy-type integral, which is given by the following expression: C "(r , FW ) + 4 - fg(V ■ F')pdv = { ^ C DD , <2-21) where vector field C s is the 3D analog of the Cauchy-type integral for the vector field F(r). Cauchy-type integral formulas can be represented using matrix notations accord­ing to Zhdanov (1988). The matrix form makes it suitable for numerical computation, which is important in practical applications. We take the convention that the z axis is directed upward. In a Cartesian coordinate system {d x, dy, dz}, we can represent the vectors Cs, 7 n and V , Si as follows: ' 7 |r-r | C s = C^da, p = ^ fi dp, n = n7 d7, (2.22) 1 r„ - rv V = - r - n73 dn, (2.23) |r - r | |r - r | where rn = n; ®,P,Y,V = x ,y,z; and we also use the convention that the twice recurring index indicates a summation over the index. Using these notations, we can write the scalar components of the Cauchy-type integral as follows: / _1 p r r - r Ca = - PfT1 -----7733 nYdS a ,^ , l ,n = x , y ,z, (2.24) 4n JJs |r - r | where the four-index A symbol is expressed in terms of the symmetric Kronecker symbol 8ap as follows: A afYn 8af 8^^ + 8«n8Py 8«y8fn; 8«p | 0 ck ^ /3 " (2.25) CHAPTER 3 REPRESENTATION OF POTENTIAL FIELDS USING 3D ANALOGS OF CAUCHY-TYPE INTEGRALS It is well known that the potential field anomaly caused by a 3D body with anomalous density or magnetization can be calculated using a volume integral (Zh­danov, 2002). However, the numerical evaluation of the volume integral for a domain with arbitrary shape can be computationally expensive. It can be demonstrated that the potential field caused by this model can be represented efficiently by 3D analogs of Cauchy-type integrals. In the Cauchy-type integral representation, the calculation of a potential field caused by a 3D body can be reduced to a surface integral over the surface which encloses the volume with anomalous density or magnetization (Zhdanov, 1988). The derivation of the Cauchy-type representation of a potential field caused by a general 3D body filled with masses with anomalous density or magnetization can be found in Zhdanov (1988). For completeness, I will briefly review these derivations in this chapter. I will also extend this formulation to the gravity gradiometry data which have become widely used in oil and mineral exploration recently due to the higher resolution compared to conventional gravity field measurements. One important application of the 3D analog of a Cauchy-type integral is modeling and inversion for a density/susceptibility-contrast model. The representation of a gravity field caused by the density-contrast surface has been elaborated in Zhdanov (1988). This elegant representation can also take into account the variation of density contrast in the vertical direction, which is usually the case in practical application due to the compaction of rocks with depth and pressure increase. In the following 12 sections, I will extend these derivations to the gravity gradiometry data and the magnetic problem. 3.1 Representation of the gravity field using CTI The gravity field g(r) is governed by the following equations: V ■ g = - 4nYpo, V x g = 0, (3.1) where y is the universal gravitational constant and p0 is a constant density inside a 3D domain. It was shown by Zhdanov (1988) that we can construct a vector field F such that 4n F = - Yp0r, and V ■ F = 4nYp0. (3.2) 3 By substituting equation (3.2) into equation (2.21), we find C' ( r - f 7Por)+ Y Uln PoVy -7\dv = { ( S T c D € D ■ (3-3) It is obvious that the second term on the left side of equation (3.3), with the negative sign, is the gravity field g(r) caused by the density distribution within domain D: g (r ) = -Y JJJn PoV jr-7 j dv (3.4) As a result, the gravity field g(r) can be represented using a Cauchy-type integral as follows: { g(r' ) = { - ^YPor' + TYPoCs(r' , r), r' € D (3 5) ) I f 7PoCs(r ', r), r € CD ■ (3.5) In geophysical applications, we usually consider the case where the observation point is located outside of the anomalous body: r € CD. In this case, the gravity field can be written as follows: 4>n g (r' ) = yY P oCV , r ) . (3.6) The previous equation can be rewritten in a matrix notation for the scalar gravity field as follows: I 1 c c r - r 9a = - o YPo A a^lvr^-1 - 7n^i n1ds, = x , y , z . (3.7) 3 JJs jr - r j 13 According to the properties of the Cauchy-type integral, (3.8) equation (3.5) can be written in a uniform form as g(r ) = -3- 7P0 [Cs(r , r) - C s(r ', r ')] = - 7poCs(r ', r - r ). (3.9) One needs to note that equation (3.9) holds for the observation points either within or outside domain D, r' €D or r' €CD. Similarly to equation (3.7), equation (3.9) can be represented in matrix notations as follows: We can use equation (3.10) to calculate the gravity-gradient tensor whose scalar components are equal to the derivatives of the corresponding scalar components of the gravity field with respect to the spatial coordinates: After some algebra, one can express equation (3.11) in a matrix notation as follows: From equations (3.10) and (3.12), we can see that the gravity field caused by a volume D filled by masses with some constant density p0 can be represented as a Cauchy-type integral over the surface S of the volume. Thus, the original formula for calculating the gravity field as a volume integral is reduced to the surface integral. It is also important to point out that the density distribution inside volume D may not necessarily be a constant value. It is shown in Zhdanov (1988) that one can incorporate arbitrary integrable density-depth distribution within the volume in this formula. The advantage is that in applications, we can use this method to simulate the potential field due to the sediment basin, which is usually characterized by a density change with the depth. gav = -TTj, a, n = x,y, z. (3.11) ^ Y'/ I 'I s |r - r | 14 3.2 Representation of the magnetic field using CTI It was shown by Zhdanov (1988) that the magnetic field due to a magnetized 3D body can be represented by a 3D analog of a Cauchy-type integral as follows: H(r ) = C s(r , h(r) - h(r;)-4nIn(r)), (3.13) where a 3D body is bounded by a closed surface S, In(r) is the projection of mag­netization intensity on the normal direction to the surface, and h(r) is an arbitrary solution inside a 3D domain of the following equation: V ■ h = 4nV ■ I, V x h = 0. (3.14) In expression (3.13), term Cs represents a 3D analog of the Cauchy-type integral defined in equation (2.12). For simplicity, we will consider a Laplace distribution of magnetization: V ■ I(r) = 0, V x I(r) = 0, r E D. (3.15) Now we can select a special solution of equation (3.14) as h = 4nI such that equation (3.14) is always valid if magnetization vector I(r) satisfies equation (3.15). In this case, equation (3.13) can be written as follows: H(r ) =4nCs(r ', I(r) - I(r') - In(r)). (3.16) Since I(r) is a Laplacian vector field inside the 3D domain, by virtue of the properties of the Cauchy-type integral, we have 4nCs(r ', I ( r ) )= | ^ ■ (3.17) Also, since I(r;) is independent of r, the following relationship holds: 4nCs(r , I(r')) = | 0<r^.) '^ 0 ^ Comparing formulas (3.17) and (3.18), we arrive at the following identity: C s(r , I(r)) = C s(r ', I(r')). (3.19) 15 Based on this identity, equation (3.16) can be written as follows: H(r ) = - 4nCs(r , I„(r)), (3.20) or in equivalent form as H ( r ) =4nCs (r ', Ir (r) - I(r)), (3.21) where IT(r) is the tangential component of a magnetization vector on surface S. For a more special case where the magnetization vector is a constant inside domain D, we have the following: I(r) = I0, (3.22) and equations (3.20) and (3.21) can be rewritten as follows: H(r ) = - 4nCs(r , £ ( r ) ) , (3.23) H(r ) = 4nCs(r ', I0(r) - I0(r)) = 4nCs(r ', I0(r)), (3.24) where we take into account that C s(r ', I0(r)) = 0, r' eCD. (3.25) Note that, although I(r) = I0 is a constant within domain D, its normal, I^(r), and tangential, I°(r), components on the surface S are some arbitrary functions related to the geometry of S. 3.3 CTI representation of the gravity field of the density-contrast surface In this section, I will show the derivation of the representation of a gravity and gravity-gradiometry field for a density contrast surface model using a 3D Cauchy-type integral. We first assume that the density contrast is a constant which leads to the constant sediment and basement density for a typical sediment-basement interface model. Following this, I will show the derivation for a sediment-basement interface model with density contrast varying with depth. 16 3.3.1 Models with constant density Let us consider a model of sediment-basement interface with a density contrast at some surface r, shown in Figure 3.1. The most idealized sediment-basement interface is just a horizontal plane which is indicated as plane P . In this case, one will not observe any gravity anomaly on the earth's surface. However, the real sediment-basement interface r is a variable surface. We assume that surface r is described by equation z = h (x, y) - H0, and a horizontal plane P is given by equation z = - h0, where H0 > h (x,y) > 0 and h (x, y) - H0 ^ 0 for \Jx2 + y2 ^ ro, where H0 is a constant. Let us draw a sphere OR of radius R with the center in the origin of the Cartesian system of coordinates. We denote by r R and PR the parts of the surfaces r and P , respectively, located within the sphere OR. For the first model, we assume that the real sediment-basement interface r R is located above plane P . We also assume that the sediment layer has a constant density ps and the basement has a constant density pb (pb > ps). We assume also that both rand P extend infinitely in the horizontal direction, and r R ^ P at infinity. The gravity anomaly is caused by the density volume DR which is bounded by a closed surface, formed by r R and PR and the parts of the sphere OR between these two surfaces as shown in Figure 3.1. It is demonstrated in Zhdanov (1988) that the gravity field caused by volume DR is expressed by g (r ') = 4nYP0CrR(r ', (z + H))dz), (3.26) for a case where r R^ ^ ^ P at infinity. As a result, the Cauchy-type integral in equation (3.26) is calculated along an infinitely extended surface r. In equation (3.26), p0 is the density contrast between the sediments and the basement: P0 = Pb - Ps > °. (3.27) In this subsection, we consider for simplicity that pb and ps are constant, which leads to a constant density contrast p0. For the model shown in Figure 3.1, we will always have a positive gravity anomaly. Now we will consider another model presented in Figure 3.2 where the density contrast 17 surface is below the horizontal plane P . In this case, we will have a negative density anomaly within domain Dr. Similarly, the gravity field can be expressed as follows: g(r/) = 4nY(-po)CrR(r' , (z + Ho)dz)• (3.28) We have the following expressions for the scalar components of the normal vector pointing outside domain DR for a model shown in Figure 3.1: dh(x, y) nxds = -------^ - dxdy = bx(x,y)dxdy, dx dh(x, y) ny ds = ------ 7^- dxdy = by (x,y)dxdy, (3.29) nz ds = dxdy = bz (x,y)dxdy• Similarly, for the model shown in Figure 3.2, the scalar components of the normal vector pointing outside domain DR are equal to the following: dh(x, y) nxds = - ^ - dxdy = -bx(x,y)dxdy, dx dh(x, y) ny ds = - ~~~- dxdy = -by (x,y)dxdy, (3.30) nz ds = -dxdy = - bz (x,y)dxdy where bx(x,y) = - dh^x y) ,by(x,y) = - dh^y y) ,bz(x,y) = - 1 h = z + H0• (3.31) It is important to note that, although for the models in Figure 3.1 and Figure 3.2, the equations for the normal vector have different signs, as shown in equations (3.29) and (3.30), respectively, the final expressions for the fields are exactly the same, since the sign for the anomalous density for models 1 and 2 is also different. Thus, in matrix notations, the gravity field caused by the density anomaly for model 1 (Figure 3.1) and model 2 (Figure 3.2) can be expressed using a unified equation as follows: , u ( ( h(x, y)(rn - rn) „ ga(r ) = -YPo Aazjn----- .------- /|3 b7dxdy, a,Y,n = x ,y , z• (3.32) JJs |r - r | Similarly, the gravity gradient for the models in Figure 3.1 and Figure 3.2 can also be unified as follows: ( '\ f f 3Aazin h(x,y) ( ' \( / \7 7 7 gav(r ) = YPo ---- :------- 7:5---- (rv JJS - rv)(rn - rn)bYdxdy |r - r | 18 f A azjnh(xi y) s |r - r' |4 5vn bY dxdy,a,Y,n = x,y,z. (3.33) Equations (3.32) and (3.33) represent the gravity and gravity-gradient fields in the form of Cauchy-type integrals over the density-contrast surface corresponding to the sediment-basement interface. These expressions provide an analytical basis paper by Zhdanov and Liu (2013), both rectangular and triangular discretization of the density contrast surface is introduced. Numerically, rectangular is simpler than triangular discretization. However, triangular discretization is demonstrated to have higher accuracy than rectangular discretization. In my forward modeling part, both of these types of discretization are implemented. In the inversion part, only rectangular discretization is used for simplicity. In this paper, we will only introduce the scheme for rectangular discretization. For these who are interested in triangular discretization, please refer to Zhdanov and Liu (2013). We can discretize the Cauchy-type integral for computing the gravity field and its gradient in equations (3.32) and (3.33) by dividing the horizontal plane X Y into a grid of Nm cells with constant discretization A x and Ay in the x and y directions. As a result, within each cell Pk (k = 1, 2...Nm), the corresponding density-contrast surface can be represented by a flat plane described by the following linear equation: In a discretized form, equations (3.32) and (3.33) can be represented as follows: for a fast method of numerical modeling of gravity and gravity-gradiometry data. Both of these equations need to be discretized to be solved numerically. In the iNm 9 « ( 0 = £ / ^ " ' b f , (3.35) where (3.36) and rXk) = x(k), rjfc) = y(k), r{k) = z(k) = h(k) - H0, (3.37) 19 »(n); x(n); , r (n); = y(n); , r(n); = z(n); As usual, the twice recurring index in equation (3.36) indicates the summation over this index. Similarly, the discretized form of the equation for the gravity-gradient fields can be expressed as follows: Nm gav (O = E FanY?h(k)b f , (3.38) where k=1 Fa(vnYk) YPo A azYn (k) i5 3(rVk) - rVn); )(r(k) - r(™);) - r(k) r; 5vV AxAy (3.39) In the previous discrete formulas, we have approximated the density contrast surface within each cell by a flat plane. In a simple case, we can approximate the density contrast surface within each cell by an element of the horizontal plane (Zhdanov and Liu, 2013): z - h(x,y) - Ho - h(k) - 6ik)(x - xk) - byk)(y - yk) - Ho bik)(x ,y ) - 0,byk) (x ,y ) - 0. In such a special case, equation (3.35) is simplified: Nm g» ( 0 - E /ank)h(k'), where k=1 r(k) r(n); /;<nk) - -YPoia,, n(k) •; ,3 AxAy. | r(k) - r;n| We can obtain a similar formula for the gravity gradient fields: Nm gav ( O -E F<nk) h(k), k=1 where YPo A azYn (k) 5 3(r<k) - r<" ); )(r<k) - r<" );) - r(k) r; ^vn AxAy (3.40) (3.41) (3.42) (3.43) (3.44) We should notice that equations (3.41) through (3.44) may not be accurate enough for forward modeling since the accuracy of approximation by the piece-wise horizon­tal surface may not be sufficient. However, these equations are very effective for calculating the Frechet derivative matrix in the inversion process because of their simplicity. x 2 r 2 r 20 3.3.2 Models with variable density As we mentioned above, in a general case, the density contrast value is a function of depth: A p = f ( z ) . (3.45) In this case, the representation of the gravity field caused by the sediment-basement interface takes the following form (Zhdanov, 1988; Zhdanov and Liu, 2013): g(r ) = 4nGCrR(r , [R(z) - R(-H0)] d^), (3.46) where R(z) = / f (z)dz. (3.47) J-Ho Similar equations can be derived for the gravity gradient component by taking the spatial derivative of the forward operator for the gravity field. Similar to the previous subsection with a constant density contrast, we can write a matrix notation for equation (3.46) as follows: ( ' ) f f a [R(z) - R(-Ho)] (rn - rn) b d d r o4o) ga(r) = -Y / / A aziv-------------,------- TT3----------- b7 dxdy, a,Y,n = x,y,z. (3.48) j j s | r - r | Similarly, the gravity-gradiometry field for a sediment-basement interface model with variable density contrast in the vertical direction can be expressed by a 3D analog of a Cauchy-type integral in a matrix notation as follows: ( ' ) [ f a 3(rv - rV )(rn - r2)[R(z) -R (-Ho)] b d d (349) 9av(r ) = Y A aziv-------------------- ,- -775--------------------b7dxdy - (3.49) -'-'S |r - r | f f a r - r ^vn [R(z)-R (-Ho)] Y A az7n--------------■------- -5--------------bjdxdy, a,Y,n = x, y, z, J J s | r - r | where bY is defined in equation (3.31). One also needs to note that we can write R(z) as R(h(x,y) - H0) in equations (3.48) and (3.49). In a discretized form, equations (3.48) can be represented as follows: Nm / «y -1 ' m 9 0 ( 0 = E f £ f h (k>b<k>, (3.50) k=l where r .. r(k> _ r(n> f f = -Y R(z(k>)-R(-Hofl A „ „ , - 2 - -----L A x A y . (3.51) |r(k> - rn| 21 The gravity gradiometry field can be written in a discretized form as follows: N„ ga, ( 0 = £ Fa;k')h(k)b‘k'), (3.52) k= 1 where 3(r,k) - r,ra)/)(rnk) - - ^ ) R(z(k))-R (-H0) r(r F (nk) = y a - _____ n__________ n ' I '___ '____^ ^ AxAy Fav 7 = /Aaz7 n |r(k) / 15 AxAy A,, R(z(k)) - R(-H0) a a , ^ -YA az7n-------- ^ --------------------A xA y. (3-53) |r(k) - rn| As before, we can approximate the density contrast surface within each cell k by an element of the horizontal plane by considering equation (3.40). In such a special case, equation (3.50) can be represented as follows: Nm 9 a ( 0 = £ /ank)h(k), (3.54) k=1 where -(k) _ -(" )/ /a"k) = y [R(zk) - R(-H0)] Aa,: ,|(k.) _ , 3 AxAy, (3.55) |r rn| Nm is the number of cells, and n is the index of the point of observation, r^. We can obtain a similar formula for the gravity gradient fields: Nm ga, ( 0 = £ /<:k')ft(k), (3.56) where av k=1 f („k) = 3Y(r,k) - rV" )') (r ,k) - -M') [R(zk) - R (-H0)] A., A /av |r(k) - r;,j5 x y Y [R(zk) R( H0)] Av,Aa, 3 AxAy. (3.57) |r(k) - r; | 3 ^ ( ) Again, one needs to note that equations (3.54) and (3.56) are very effective for calculating the Frechet derivative matrix in the inversion process because of their simplicity. 3.4 CTI representation of the magnetic field of the magnetic-contrast surface We consider a typical sediment-basement model shown in Figure 3.1 with plane P as a reference interface between sediments and basement. Surface r represents an 22 actual sediment-basement interface. Domain D is bounded by the surface r and the plane P . We will calculate the magnetic effect of the domain D filled with magnetized material with a magnetization vector I0. To this end, we cut from the domain D a domain DR, bounded by a spherical surface OR with a radius R and the parts r R and PR of the surfaces r and P cut therefrom by the sphere OR (see Figure 3.1). We assume that surface r R tends to plane P at infinity. We consider a nonmagnetized sediment, which is typically the case, and a uniformly magnetized basement with magnetization vector I0. In this case, the magnetic anomaly due to the basement can be calculated as the effect of the bounded domain DR with R ^ ro. According to equation (3.23), the magnetic field caused by the bounded domain DR can be represented as follows: H(r') = -4n {C rR(r', In(r)) + CORUPr(r ' , I ? (r ) ) } , (3.58) where OR is a lateral surface of the sphere OR enclosed between r and P. Since we assume that the surface r R tends to horizontal plane P , the integral over OR tends to zero when R ^ ro. Therefore, equation (3.58) can be reduced to the following form: H(r') = -4n {C rR(r', I?(r)) + C Pr(r ' , I?,(r))} . (3.59) The total magnetization vector I0 can also be decomposed into the vertical and horizontal components: I0 = IV + Ih. (3.60) It is easy to see that both the vertical and the horizontal components are constant within domain D. The Cauchy-type integral over the part of plane PR, CPr(r ', I?(r)), can be evalu­ated in the following way. On the horizontal plane P , the normal component of the magnetization vector is exactly the vertical component, and the tangential component is equal to the horizontal component. Therefore, we have C Pr(r ' , I?(r)) = C Pr(r ', IV), (3.61) where we assume that the positive z direction points downward. Taking into account that IV(r) is a Laplacian vector in the whole space, we have: C^uor upr (r ', 10 ) = 0. (3.62) By ignoring the part of the integral over OR for R ^ x>, we obtain C rRUPR (r ', IV) = 0, (3.63) which leads to the following result: CPR(r ', I0) = -CrR(r ', I0)• (3.64) By substituting equation (3.64) into equation (3.58), we finally find that: H (r ) = -4nCrR ( r , I^(r> - IV )• (3.65) It was shown in Zhdanov (1988) and Zhdanov and Liu (2013) that a Cauchy-type integral can be expressed using the matrix notations as follows: ' _1 p p r _r Ca(r' , v) = - A afYv^^-r- /^3nl ds, a ,P, 4n JJs |r - r | i , n = x , y , z ; (3.66) where the four-index A symbol is expressed in terms of the symmetric Kronecker symbol 8ap as A afjn 8a f + 8an8fy 8ay8fn; 8af ^ 0 CK ^ /? ^ (3.67) We will define the distance between the actual sediment-basement interface and reference plane P (z = H0) as h = z - H0• (3.68) Its derivative in the horizontal dimension is defined as ha = ^h,® = x ,y • (3.69) Using representation (3.65) and formula (3.66), the scalar components of the magnetic field can be written in the following form: Hx(r') = f f s { (x - .t') [I^Wxbx + IJCDyby + I » z ] (3.70) 23 24 +7° [(z' - z) bx - (x' - x)H dxdy, Hy(r') - f f s T^TT'^'T { (y - y') ft(r)xbx + ^Wyby + I j(r)*] (3.71) +7° [(z' - z) by - (y' - y)]j Hz(r') - JJs - { (z - z') [in(r)xbx + I ^ r ) ,b, + I ^ ] (3.72) -7° [(x' - x) bx + (y' - y)b, + (z' - z)]j dxdy, where we use the following notations: bx - - ^ , by - - hy, bz - 1. (3.73) In geophysical exploration, we measure the total magnetic intensity (TMI) field, Ht , which is defined as follows: Ht - l ■ H - ZxHx + lyHy + lzHz, (3.74) where l is a unit vector in the direction of inducing a geomagnetic field: l - (lx,ly ,lz ). (3.75) After taking some algebra, one can find the expression of the TMI field using a Cauchy-type integral as follows: Ht (r') - f f [in(r)xbx + in(r)yby + in(r)z] [/x(x' -x ) + i, (y '-y) + lz(z' - z)]dxdy JJs |r - r'| + rr 7° (z' ~ z) (lxbx + 1„by - L) (3.76) ././s r - r' f f 7° [(x' - x)(lx + lzbx) + (y' - y)(ly + lzby)] ^ ^ y t/t//ss |rr - r'l In the above formula, we have a magnetization vector, which can be calculated from the inducing magnetic intensity H/ and the magnetic susceptibility as follows: I0 - XH/. (3.77) We can discretize the Cauchy-type integral for computing the magnetic field in equations (3.70), (3.71), (3.72), and (3.76) by dividing the horizontal plane X Y 25 into a grid of Nm cells with constant discretization sizes of Ax and Ay in the x and y directions, respectively. As a result, within each cell Pk (k = 1,2...Nm), the corresponding magnetization contrast surface can be represented by a flat plane described by the following linear equation: z = h(x, y) + Ho = h(k> - blf>(x - x - ) - b[k)(y - yk) + H( (3.78) where (xk,yk) denotes the center of the cell Pk, and h(k> = h(xk,yk). In a discretized form, equations (3.70), (3.71), (3.72), and (3.76) can be repre­sented as follows: Nm f f * ( 0 = E f i "k>. (3.79) k=l Nm Hy ( O = e /<"*>, k=l Nm Hz ( O = £ f ' "k>, k=l Nm HT (r; ) = E f ^ 0. k=l where the kernels are defined as follows: f l (nk> 1 - (x- - x j tor-)xb|k> + I;(r-)y byk> + ) |r- - rn| (zn - zk) b|k> - (x; - xk) A xAy (3.80) (3.81) (3.82) (3.83) 1 |r- - rn| tt (y- - y ; ) fI;(rk)xbik> + I;(r-)y byk> + I ^ r - ) + i (z; - z- ) bik> - (y ; - y-^ A xAy , (3.84) fz ( ;k > |r- - rn| tt (y- - yn) [I;(rk)xbik> + I;(r-)y byk> + i ; ^ ) - i 0 (x; - x- ) bik> + (y ;- y- )byk> + (z; - z- ^ A xAy (3.85) AxAy x ^.(rak> _ 1 f T = ] , |rk - rni I; (r- )xbXk>+ I; (r- )ybyk>+ I; (r- )z [ix(x; -x- ) + ^ (y; -y- ) + ^ (z; - z- )] i o + (z; - zk) (1xb|k> + ^ybyk> - ^z) |r- - rn| AxAy (3.86) z z 1 z 26 I 0 ± V (x? - xk )(lx + lz bXk)) + (y? - yk )(ly + lz bik)) AxAy. |rk- rn| In the above derivations, we consider that the actual sediment-basement interface is above the reference plane P. Similar to the derivation for the gravity problem, we can find that the final expression for the magnetic field we derived here will be the same as that for the model in Figure 3.2, where the actual sediment-basement interface is below the reference plane P. 27 Figure 3.1. Illustration of a contrast density model for sediment-basement interface with a positive anomaly, where the plane P is the average depth of sediment-basement interface and r R is the actual sediment-basement interface. 28 Figure 3.2. Illustration of a contrast density model for sediment-basement interface with a negative anomaly, where the plane P is the average depth of sediment-basement interface and r R is the actual sediment-basement interface. CHAPTER 4 3D FORWARD MODELING OF MT DATA In this dissertation, 3D modeling of MT data is mainly based on the integral equation (IE) method. I use the parallelized contraction IE algorithm (Zhdanov et al., 2006), which is capable of modeling large-scale geoelectrical structures. The methodology of a 3D contraction IE method for geophysical electromagnetic modeling was discussed by Zhdanov (2002, 2009). I present here the basic principles of the IE method for completeness. Let us consider a 3D geoelectrical model with background conductivity ab and a domain D with the anomalous conductivity Aa: a(r) = ab(r) + Aa(r), r G D. (4.1) In the frequency domain, the Maxwell's equations for quasi-stationary EM fields considered in the MT method have the following form: V x E = i^^0H, (4.2) V x H = j e + aE, (4.3) where j e is the extraneous currents representing the excitation source. For the MT problem, the excitation source is the vertically propagating plane EM wave. In the framework of the IE formulation, the total electric and magnetic fields are decomposed into their background (Eb, Hb) and anomalous (E" , H" ) parts: E = Eb + E°, (4.4) H = Hb + H" . (4.5) It was shown in Zhdanov (2002, 2009) that the anomalous fields can be expressed as an integral of the excess currents within the anomalous domain as follows: E" (r?) = / / / d G e (rj|r)Aa(r) ■ |Eb(r) + E°(r)j dv, (4.6) 30 H° (rj) = j j j D Gh(r j|r)Aa(r) ■ |Eb(r) + E°(r)j dv, (4.7) where G E(rj|r) and GH(rj|r) are the electric and magnetic Green's tensors defined for a medium with the background conductivity ab. Formula (4.6) becomes an integral equation, when rj e r. The anomalous electric field inside domain D can be obtained by solving the integral equation. Once the electric field in domain D is found, the anomalous field at the receivers outside domain D can be calculated using the above two formulae. In 3D problems, the MT impedance tensor is given as follows: Z 7 7^xy Zy x Zyy (4.8) The MT impedance tensor defines the linear relationships between the horizontal components of the electric and magnetic fields (Cantwell, 1960; Zhdanov and Keller, 1994; Berdichevsky and Dmitriev, 2008): Ex (r) = ZxxHx(r) + Zxy Hy (r), (4.9) Ey (r) = ZyxHx(r) + Zyy Hy (r). (4.10) The background plane-wave MT field can be decomposed into two different po­larizations - TE mode and TM mode (Zhdanov and Keller, 1994; Berdichevsky and Dmitriev, 2008): The background fields for these two different polarizations are given as follows: Eb(1) = {E bx(1), 0, 0} , Hb(1) = { 0,Hby(1), 0} , (4.11) Eb(2) = { 0,Ey(2), 0} , Hb(2) = { Hbx(2), 0, 0} . (4.12) The two polarizations generate the corresponding MT fields on the earth's surface as follows: E(1) = {E(1), E((1), 0} , H(1) = { h ( 1), H(y1), H(1)} , (4.13) E(2) = { E{2),E{2), 0} , H(2) = { H(:2),H(2),H(2)} . (4.14) The full MT impedance tensor can be computed from the field excited by the two polarizations as follows: EW h (2) - e (2)h (1) z = x y x y (4 15) = HHx( 1)HHy(2 ) - HHx( 2)HHy(1) , ( ) 31 z - y y y y ( 4 17) Zyx = - i-(i) zt(2) u(2^(!) ' ( ) 7 _ y x y x Zyy - ^t(1) U(2) tt(2)it(1) . (4.±8) ^7 _ Je-'x( 2)h x( i) _ e x( 1)h x(2 ) Z/| 1 p\ Zxy - H(1) H(2) H(2) H(1) > (4.i0) Hx Hy - Hx Hy e (i)h (2) _ e (2)h (1) y y y y HHx(1 ) HHy(2 ) - HHx(2 ) HHy( e y(2 ) h x(1 ) - e y( 1)h x(2) HHx(1 ) HHy(2 ) - HHx(2 ) HHy( Formulas (4.15) through (4.18) are used for computing the synthetic MT impedance data, required for the inversion algorithm. In this application, I consider that the earth contains a sediment layer and a resistive basement. For simplicity, I assume that the conductivities of the sediment and basement are both constant. The MT anomaly observed on the earth's surface is caused by the variation of the sediment-basement interface. CHAPTER 5 INVERSION OF POTENTIAL-FIELD AND MT DATA FOR DEPTH-TO-BASEMENT ANALYSIS In the previous chapters, I have discussed the forward modeling problem for potential-field data based on a 3D analog of a Cauchy-type integral and the forward modeling problem for MT data based on the integral equation method for a model of the sediment-basement interface. In practical applications, it is important to estimate the location of the sediment-basement interface based on regularized inversion using potential field and electromagnetic field data. In this chapter, I will formulate the inverse problem of the depth-to-basement estimation for potential and MT field data. I first assume that the physical properties such as density and the conductivity of sediment and basement are given. In the inversion, only the depth to basement will be recovered using individual geophysical methods such as gravity or MT. Following this, I will discuss a simultaneous inversion of the depth-to-basement and physical properties of the sediments, such as density and conductivity, using potential field or MT data. In the last section, I will also introduce a joint inversion of potential field and MT data for the depth-to-basement and the physical properties of the sediment. 5.1 Inversion of the potential-field data for the depth-to-basement In forward modeling of the potential-field data, the model parameters are the elevations, h(k) - h(xk,yk), of the density- or magnetization-contrast surface with re­spect to the horizontal plane P , assuming the values of the density and magnetization contrasts are given. As we can see from the forward modeling equations, the forward operator is nonlinear. Correspondingly, the inversion is also a nonlinear problem. 33 The traditional inversion of potential-field data to find the density distribution is a linear problem, and the Frechet derivative can be easily found and it does not change during the iterative inversion. In my inversion, the Frechet derivative is a function of the model parameters and may change from iteration to iteration. Fortunately, in the problem under consideration the Frechet derivative has an analytical form. In Appendices A.1 to A.3, I have derived the explicit expressions for the Frechet derivative for both the gravity and magnetic fields. It is well known that the inversion for potential-field data is an ill-posed problem (Zhdanov, 2002). Theoretically, one can find an infinite number of solutions which can reasonably fit the observed data. In order to obtain a stable and geologically reasonable result, I need to apply regularization to impose some restrictions on the solution. The inversion is based on the minimization of the Tikhonov parametric functional (Tikhonov and Arsenin, 1977): Pa(m, d) = (W d A(m) - W dd)T (W dA(m) - W dd) (5.1) + (W m m - W mmapr )T (W mm - W mmapr), where A is the forward modeling operator for the potential field problem; W d is the data weighting matrix which is determined based on the magnitude of different potential field components (the ratio between the average of different field components is used to scale them to the same level); d is the vector of the observed potential field data; m is the vector of the model parameters, h; mapr is the a priori model based on other known geological information (in most cases, I assume that there is no a priori model); and W m is a diagonal matrix of the model parameter weights based on the integrated sensitivity: Wm = diag(FTF )1/2, (5.2) where F is the Frechet derivative matrix. The minimization problem (5.1) can be reformulated using a space of weighted parameters: mw = W mm. (5.3) In the space of the weighted parameters, the Tikhonov parametric functional is given as follows: 34 P a(mw, d) = (Aw(mw) - d ) T(Aw(mw) - d ) + a(mw-mwpr)T(mw-mwpr), (5.4) where A w is a new forward operator in the space of the weighted parameters, which can be related to the forward operator A in the original space as A w = AWm1. (5.5) There exist several different methods, such as the steepest descent and Newton methods, for solving this minimization problem. In this dissertation, the minimiza­tion of the Tikhonov parametric functional is based on the reweighted regularized conjugate-gradient method, which is generally faster than the steepest descent method whereas it requires less compuation compared to the Newton method. With index n + 1 referring to the iteratively updated model n, the algorithm is given as follows (Zhdanov, 2002): ^au S _w3 + 1 sn+1 sA( d - w)n wwm(( n A wnr - d, (5.6) la3 _ rw + a (mwn_ mwn ) wn wn n ' ^nv11-^ apr/5 (5.7) o r _ i i c j i 7 i i e ( n - 1) ii2 (5.8) l«3 _ l«3 + Q wn wn 1 y an l«ri-1 lao _ n lw(n-1) ' sw0 lao , sw0 , (5.9) )n n 3 nn e> 3 '__ rp __.' l^ (f t F + a Dl^ sws (F wnF wn + anA)lwn , (5.10) mwn _ mn+1 _ - mnw 3 -- knar i lwarni (5.11) mn+1 _ W' T m- 1mnw+711> (5.12) m^ 1 _ Wmmn+1, (5.13) W ++1 - mwpr+1 )s C t _ (mw+1 - mwpr)s (5.14) Y _ iisw++[1 ii2/iisw+1ii2, (5.15) an+1 _ f an, Y < 1 1 an/Y, Y > 1 . (5.16) I solve the problem in the space of the weighted model parameters. In the algorithm given above, rw is a residual vector between the predicted and observed data; lwn is the steepest ascent direction; lwn is the conjugate gradient direction, which is a combination of the current steepest ascent direction and the previous conjugate 35 gradient direction with the coefficient ,5^". We can see from equation 5.9 that the conjugate gradient direction is the same as the steepest ascent direction at the first iteration. The step length is obtained using a linear line search scheme. The regularization parameter an is selected using an adaptive method as shown in my algorithm. Parameter mWpr represents an a priori model, selected based on all known information about the model parameters. During the inversion process, the Frechet derivative matrix changes in every iteration. One of the most expensive parts of inversion is the computation of the Frechet derivative matrix. In order to speed up the inversion, the Frechet derivative can be updated not on every iteration but after every five or ten iterations. 5.2 Inversion of the magnetotelluric data for the depth-to-basement Let us consider the model of the sedimentary basin shown in Figure 5.1. The base­ment has the background conductivity ab, and domain D represents the conductive sediments. I assume for simplicity that the sediments have a uniform conductivity of however, in a general case, the method can be extended to the case of an arbitrary distribution of the conductivity, Os(r) - a6(r) + Aa(r). In the inversion, domain D is discretized into N columns, denoted as subdomains D j , with conductivity as. The horizontal dimension of each subdomain is known and fixed. Contrary to the conventional MT inversion, which recovers a volumetric distribution of the subsurface conductivities, the goal is to find the depth of each column. If the conductivity of the sediments is unknown, the inversion can also recover (r) jointly with the depth-to-basement estimate. We should note that, for IE forward modeling, the columns should be further discretized in the vertical direction. The inversion is based on the regularized conjugate gradient method. One key aspect of nonlinear inversion using the conjugate gradient method is the need for computing the Frechet derivative of the observed MT data with respect to the thick­ness of the column and the sediment conductivity, as. One of the simplest way of 36 solving this problem is using the quasi-Born approximation (Zhdanov, 2009). The quasi-Born approximation has been widely used for deriving the Frechet derivative of the electromagnetic data with respect to the conductivity distributions (Zhdanov, 2009). In this dissertation, I have extended this approach for the calculation of the Frechet derivative of the electromagnetic field with respect to the column thickness, which represents the depth-to-basement in my application. A detailed derivation of the Frechet derivatives for the thickness of the sediments is given in Appendix A.4. As usual, we consider the minimization of the Tikhonov parametric functional (Tikhonov and Arsenin, 1977) for solving the ill-posed inverse problem: Pa(m, d) = (Wd A(m) - Wdd)T (W dA(m) - Wdd) (5.17) + (Wm m - W mmapr )T (W mm - W mmapr), where A is the forward modeling operator for the MT problem; W d is the data weight­ing matrix which is determined based on the magnitude of different MT impedance tensor components (the ratio between the average of different MT impedance compo­nents is used to scale them to the same level); d is the vector of observed MT data; m is the vector of the model parameters, h; and W m is a diagonal matrix of the model parameter weights based on the integrated sensitivity, which is defined in equation 5.2. The minimization of the parametric functional for an MT inverse problem is also solved using the conjugate gradient method, which is described in the previous section. The developed theory and method have been implemented in the computer code, which was tested on several synthetic models, discussed below. 5.3 Inversion of the geophysical fields for the depth-to-basement and for the physical properties of the sediment In the previous section, I have assumed that the physical properties, such as the density, magnetization, and conductivity of rocks, are known from other geophysical or geological data, e.g., well logging and drilling cores. In this scenario, the inversion will only recover the depth-to-basement. As we can see from the previous section, in this case, the model parameter m is represented by the depth-to-basement, h. 37 However, the physical properties may not be well known in real application. In this case, it is desirable to invert for the depth-to-basement and for the physical properties of the sediments. In such a case, the model parameters m can be written as follows: m - [h; Ap] (5.18) for a gravity problem; and m - [h; as] (5.19) for an MT problem. In equations (5.18) and (5.19), the parameters Ap and are the density contrast and the conductivity of the sediments. One needs to note that the conductivity of the basement, which is the background conductivity for the MT problem, is still assumed to be known in this case. This back­ground information can be well determined from 1D inversion or a 3D conventional MT inversion. As we will see below, we can use the conventional MT inversion results in order to create a starting model for the depth-to-basement inversion. In a case where the inversion is aimed at recovering the depth-to-basement and the physical properties of the sediment using gravity field and MT data, one needs to know the Frechet derivative of the gravity data with respect to the density contrast and the Frechet derivative of the MT data with respect to the conductivity of the sediment. Using a Cauchy-type integral representation of the gravity field, one can find the Frechet derivative of the gravity field with respect to the density contrast between the sediment and basement analytically by taking the derivative of the gravity field with respect to the density contrast. A detailed derivation of the Frechet derivative of the gravity field with respect to the density contrast can be found in Appendix A.1 . In the case of the MT inverse problem, I calculate the Frechet derivative with respect to the sediment conductivity using a quasi-Born approximation. A detailed derivation can be found in Appendix A.5. 38 5.4 Joint inversion of different geophysical data for the depth-to-basement estimation In the previous sections, I have discussed the inversion for the depth-to-basement using potential-field or MT data separately. However, there is theoretical and prac­tical nonuniqueness in the solution of these inverse problems (Zhdanov, 2002; Roy et al., 2005). The nonuniqueness is more serious for the potential-field inverse problem due to the existence of equivalent sources. For this reason, the joint inversion and interpretation of different geophysical data have become increasingly popular now. The joint inversion of different geophysical data can produce a geological model that is characterized by different physical properties and benefits from the resolution power of each individual data set (Jegen et al., 2009). The model parameters for different geophysical data sets usually represent dif­ferent physical properties. For example, conventional gravity and electromagnetic inversions deal with different model parameters such as density and conductivity. To construct a geological model from the joint inversion of gravity and electromagnetic data, one needs to impose some constraint on the density and conductivity. The most common way of joint inversion of different geophysical data is based on the minimization of the cross-gradients between the different physical properties such as the density and conductivity (Gallardo and Meju, 2004; Gallardo, 2007; Gallardo et al., 2012). This approach is based on the structural similarity between different model parameters. Zhdanov et al. (2012) proposed a new approach for the joint inversion of different geophysical data based on Gramian constraints. The joint inversion formulation can be simplified if different geophysical data sets (e.g., CSEM and MT) represent the same physical parameters. Gribenko and Zhdanov (2011) have implemented 3D joint inversion of marine MT and CSEM data to image a salt dome and hydrocarbon reservoir. From their result, we can see that the subsurface image can be improved by the joint inversion compared to the inversion using either marine CSEM or MT data alone. For the depth-to-basement inversion, the primary model parameter is the depth-to- basement for both the gravity and MT fields. The joint inverse problem can also be described by the minimization of the Tikhonov parametric functional in equation 39 (5.17). In this case, the data vector d can be written as follows: d dMT; dgravity (5.20) The model vector, m, is the depth-to-basement , h, considering that we know the physical properties of the sediment and basement. In a case of inversion for both depth-to-basement and the physical properties of the sediments, the model vector m becomes m = [h; as;Ap] , (5.21) where h is the depth-to-basement for both the gravity and MT problems, as is the conductivity of the sediment, and Ap is the density contrast between sediment and basement. The Frechet derivative matrix can be calculated for the MT and gravity fields separately and then they can be combined together in the inversion. The Tikhonov parametric functional for the joint inversion case can be written as: Pa(m, d) = (Wd1A 1(m) - Wd1d 1)T(W ^A ^m ) - W ^ ) + (Wd2A 2(m) - Wd2d2)T (Wd2A 2(m) - W d2d2) (5.22) + a (W mm - Wmmapr )T (Wmm - Wmmapr ), where A 1 and A 2 are the forward modeling operators for the MT and potential field problems; d 1 and d2 are the vectors of observed MT and potential field data; m is the vector of the model parameters; W m is a diagonal matrix of the model parameter weights based on the integrated sensitivity; W d1 and W d2 are the data weighting matrix for MT and potential field data. One needs to note that the magnitude of the MT and gravity data are different. Proper data weighting is required to preserve the balance between the two misfits for the MT and gravity data. In this dissertation, I choose the ratio between the average of potential field and MT impedance component to scale them to the same magnitude. In a special case, one can formulate a problem of joint inversion of gravity and gravity gradiometry data for a depth-to-basement estimation. In this case, the joint inversion can be further simplified since the model parameters are always the same for these two different data sets. 40 Figure 5.1. Sketch of a conductivity model of the sedimentary basin. Domain D represents the conductive sediments with conductivity as, discretized into a grid of vertical columns. CHAPTER 6 MODEL STUDIES In this chapter, I will present several model studies for the modeling and inversion of potential field and MT field data for the sediment-basement interface model. In the first section, I will show the modeling and inversion of the gravity field for a sediment-basement interface model using Cauchy-type integral methods. Following this section, I will present the modeling and inversion of magnetic field data for the sediment-basement interface model based on a 3D analog of a Cauchy-type integral. In the third section, I will show the modeling and inversion of MT data for the estimation of depth-to-basement and sediment conductivity. The joint inversion of potential-field and MT data will be discussed in the final section of this chapter. 6.1 Modeling and inversion of gravity data for the sediment-basement interface model As I have mentioned in the previous chapters, the 3D analog of the Cauchy-type integral can be used for the modeling and inversion of the gravity field, especially for a sediment-basement interface model. In this section, I will demonstrate the efficiency of the Cauchy-type integral method in the application of depth-to-basement estimation by several realistic synthetic model studies. I first consider a simple sediment-basement interface model with a symmetric shape. Following this, I will present a more complex sediment-basement interface model with an asymmetric shape. Finally, I will consider a reconstructed sediment-basement interface model for the Ensenada Bay (Gallardo et al., 2003). In the model studies, I consider that the density contrast between sediment and basement can be a constant value or it can vary with depth. For the model studies, I first assume that the density contrast between sediment and basement is given and only the depth-to-basement is inverted. Following this, I consider that the density contrast between sediment and basement 42 is unknown. In the inversion, both the depth-to-basement and the density contrast value will be recovered. 6.1.1 Model 1: Symmetric model In this section, I consider a sediment-basement interface model with a symmetric shape as shown in Figure 6.1 where the black dots indicate the observation stations on the earth's surface with the spacing of 200 m. As per convention, we consider that the positive z directs upward for the potential-field problem in this dissertation. For this model, the maximum depth of the sediment-basement interface is 750 m. In the first subsection, I consider that the density contrast between sediment and basement is a constant value. In the second subsection, I consider a more complex case where the density contrast value changes linearly with depth. In the last subsection, I assume that the density contrast between sediment and basement changed exponentially with depth. For this model, only the depth-to-basement is inverted by assuming that the density contrast profile with depth is well known from other geophysical information such as well logging. 6.1.1.1 Constant density contrast For simplicity, I first assume that the density of both sediment and basement are constant values. Therefore, the density contrast between sediment and basement is also a constant value. I consider the density contrast value between sediment and basement is 0.4 g/cm3, which is a reasonable value. The gravity anomaly caused by this model will be computed using the Cauchy-type integral method. In order to validate the forward modeling result from the Cauchy-type integral method, I also computed the field using the conventional method based on the prism discretization of the sediment-basement interface model. Figure 6.2 shows a prism approximation for this model that will be used for the conventional forward modeling method. Figure 6.3 presents a comparison of the gravity field on the surface calculated from the Cauchy-type integral method and the conventional method. One can see that the fields computed from these two methods are very close to each other. The slight difference can be attributed to the prism approximation of the continuous sediment- 43 basement interface and the center point approximation for the volume integral in each prism. One needs to note that the computation speed of the Cauchy-type integral method is around 30 times faster than the conventional method. The forward modeling result, of the vertical gravity field component, for this model will be used for the inversion to recover the depth-to-basement. I have contaminated the synthetic data by 5% random noise. The inversion process is terminated when the normalized misfit between the observed and predicted data reaches the noise level after five iterations. The starting model for my inversion is a flat surface located at 300 m below the surface. I have tested other starting models at different depth, and all the inversions converge to almost the same models. Figure 6.4 shows a comparison between the true model and the inverted model. One can see that the geometry and depth of the sedimentary basin is well recovered. 6.1.1.2 Linear density contrast In practical application, the density of sediment increases with depth due to compaction caused by pressure. As such, the density contrast between sediment and basement will decrease with depth. As I have mentioned before, the Cauchy-type integral method is capable of modeling the gravity anomaly caused by a sediment-basement interface model with density contrast varying with depth. In this section, I assume that the density of the basement has a constant value of pb - 3 g/cm3. The density of the sediment at the surface is po - 2 g/cm3, and it increases linearly with the depth according to the following formula: Ps - Po + az, (6.1) where p0 - 2 g/cm3, a - 5 x 10-4 kg/m4. (6.2) In practical application the optimized parameter of a is obtained by fitting the density log profile with a linear function. Figure 6.5 shows the density contrast profile with depth. I compared the forward modeling result obtained by the new method of Cauchy-type integrals with that based on the traditional method, using volume integrals over 44 the domain occupied by anomalous masses subdivided into the prismatic cells. The density inside of each prismatic cell was set to be a constant. Figure 6.6 shows the representation of the density contrast distribution in using prismatic cells. Both gravity and vertical gravity-gradient components were computed using both of these two methods. Figure 6.7 shows a comparison of forward modeling results obtained using the Cauchy-type integral and the traditional volume integral methods. We observe a very good fit between these results. A small difference can be attributed to the errors of the prismatic approximation of the volume density distribution in the traditional method and the discretization of the surface for Cauchy-type integral calculation, respectively. 6.1.1.3 Exponential density contrast Previously, I assumed that the density contrast value between sediment and basement decreases linearly with depth. However, the decrease ratio may vary at different depths in a real application. Usually, when the depth increases, the decrease ratio will get smaller. As such, we can construct a more complex density profile with depth using the exponential function Ap - ae-bz + ce-dz, (6.3) where a - 251.5 kg/m3, b - -0.007, c - 197 kg/m3, d - 5.2656 x 10-6. (6.4) One needs to notice that the density profile with depth is usually known from well logging and other geological information. For the application of depth-to-basement inversion using Cauchy-type integral methods, one has to fit the density profile with proper analytical functions. As an example, I use an expontential function in this case. Figure 6.8 shows the density contrast profile with depth for this model. The shape of the sediment-basement interface is still the same as in the previous case. Figure 6.9 shows the representation of the density contrast distribution in Model 1 using prismatic cells. We present the gravity responses computed using the Cauchy-type 45 integral and the traditional volume integral methods in Figure 6.10. One can see that the result produced by the new method practically coincides with that of the traditional method. I applied the inversion algorithm introduced in the previous sections to the inver­sion of the synthetic data simulated for the model with exponential density variation. Figure 6.11 shows the inversion result for the synthetic model with exponential density variation with depth. One can see that the density contrast surface was recovered well by this inversion. The normalized misfit sucessfully converges to 5% which is the noise level of the data. 6.1.2 Model 2: Asymmetric model In this subsection, I consider a more complex sediment-basement interface model with an asymmetric shape. The density contrast between the sediment and the basement is set to a constant value of 0.4 g/cm3. At first, I assume that the density contrast value is given and only the depth-to-basement is inverted. Following this, I invert for both the depth-to-basement and the density contrast value. Both the vertical component of the gravity data and the full-tensor gravity gradiometry will be used for the inversion to compare the resolution of these two different data sets. Figure 6.12 is an illustration of this model where the surface indicates the sediment-basement interface. 6.1.2.1 Inversion for depth-to-basement only For this case, I considered that the density contrast value was given. I used both the vertical component of the gravity field and the full-tensor gravity-gradiometry data to recover the depth-to-basement. The synthetic gravity and gravity-gradiometry data were contaminated by 5% random noise, and the inversion process was to be terminated when the normalized misfit reached the noise level. Similar to the previous model study, I used a flat surface at 300 m below the earth's surface as the starting model for the inversion. Figure 6.13 shows a comparison between the true model and the recovered model, at y=0, from the inversion of the vertical gravity component. One can see that the shape of the sedimentary basin is well recovered. However, the 46 true depth of the sedimentary basin is underestimated due to the sensitivity of the gravity data to the sharp variation of the sediment-basement interface. The difference between the the true model and the recovered model in the left side of this figure is probably caused by the edge effect. Figure 6.14 shows a comparison of the true model and the recovered model, at y=0, from the inversion of the full-tensor gravity-gradiometry data. From this figure, one can clearly see that both the shape of the sedimentary basin and the depth is well recovered. By comparing the inversion result from the vertical gravity component and the full-tensor gravity-gradiometry data, we can see that the full-tensor gravity gradiometry can provide a higher-resolution image of the sedimentary basin structure over the vertical component of the gravity data. However, the full-tensor gravity-gradiometry data are still not widely used for the depth-to-basement estimation in practical application due to certain reasons such as the survey expense. Again, the normalized msifit sucessfully converges to the noise level of 5%. The slight difference between the true model and recovered model in the left side of this figure is probably caused by the edge effect and the nonsmoothness of the inversion result in the right side could be attributed to the over fitting of the observed data. I have also observed that full tensor inversion produces better result than the inversion vertical gravity component and the joint inversion of vertical gravity component and full tensor gravity gradiometry data. 6.1.2.2 Inversion for depth-to-basement and density In practical application, we may not have accurate information about the density contrast value. In this scenario, it is desirable to invert for both the depth-to-basement and the density contrast value. Similar to the previous study, I used both gravity and full-tensor gravity-gradiometry data for this inversion. The initial sediment-basement interface was a flat surface at 300 m below the surface. The initial density contrast value was set to be 0.7 g/cm3 for this inversion, and the selection of the initial density contrast value will be discussed in the following subsection. Figure 6.15 shows the recovered depth-to-basement from the gravity inversion with an unknown density contrast value. One can see that the shape of the sedimentary basin is well recovered 47 but the depth is underestimated. The recovered density contrast value is 0.5 g/cm3 for this inversion. The difference between the the true model and the recovered model in the left side of this figure is probably caused by the edge effect. Figure 6.16 shows the recovered depth-to-basement from full-tensor gravity-gradiometry inversion with an unknown density contrast value. One can see that the shape of the sedimentary basin is well recovered, but the recovered depth shows a certain difference from the true depth. The recovered density contrast value is 0.43 g/cm3 for this inversion. We can see that the inversion of the full-tensor gravity-gradiometry data gives a better estimation for the depth-to-basement and the density contrast value compared to the inversion of the vertical gravity component. I have also observed that full tensor inversion produces better result than the joint inversion of vertical gravity component and full gradiometry data due to the lower resolution of vertical gravity component. The normalized misfit for this inversion also converges to the noise level of 5%. 6.1.3 Model 3: Ensenada Bay model In this section, I will present the model for the sediment-basement interface based on the reconstructed geological section of the Ensenada Bay basement (Gallardo et al., 2003). The Ensenada Bay is located at the Mexico-US border and is characterized by a relatively deep sedimentary basin. Structurally, the basin is controlled by the active South and North Agua Blanca faults (Gallardo et al., 2003). A vertical geological section of the model is shown in Figure 6.17. We have computer-simulated synthetic gravity and gravity-gradiometry data for this model. Several gravity and magnetic surveys were conducted in this area. Gallardo et al. (2003) used the gravity and magnetic data over this area to apply prismatic inversion for the depth-to-basement. In their inversion, a grid of prisms with the width of 1 km in the horizontal direction was used. The thickness and density distribution of the prisms was estimated using gravity and magnetic data. Their inversion started from some initial model with a priori information. The recovered depth corresponded reasonably well to known geological information which validated their approach for depth-to-basement estimation. 48 I set the horizontal reference plane P determining the average depth of sediment-basement interface at a depth of 1000 m under the earth's surface. Figure 6.18 shows a 3D view of the reconstructed Ensenada Bay basement model. The black dots in the figures indicate the gravity stations with the spacing of 200 m. For this model, the left portion of the sediment-basement interface is above the reference plane P, while the right portion is below plane P . The actual interface extends from -1400 m to -600 m in the z direction. I assume that the density contrast between sediment and basement is 0.4 g/cm3. For this model, both the vertical gravity component and the full-tensor gravity-gradiometry data are used in the inversion. The synthetic data are contaminated by 5% random noise. I began my numerical experiments with the inversion of the vertical component of the gravity field; after that, I investigated the inversion of the full tensor gravity gradiometry data. In my inversion process, I used a large domain size of 12 km by 10 km in the x and y directions. However, in the images of the inversion result, I only displayed the area with data coverage. At the beginning, I assumed that the density contrast between sediment and basement was well known. The inversion would only recover the geometry. Then I implemented the inversion to recover both the geometry of the sediment basin and the density contrast between sediment and basement. 6.1.4 Inversion for depth-to-basement with given density contrast I assumed that the density contrast between sediment and basement was given as 0.4 g/cm3, which is the true value for this model. The inversion would only update the location of the sediment-basement interface at each iteration. I ran the inversions for three different data sets: a) vertical component of the gravity field, gz, only; b) full-tensor gravity-gradient data; and c) joint inversion of vertical gravity field and full-tensor gravity-gradient data. Figure 6.19 shows a comparison of the inversion results produced for different gravity data sets with the true model. From the figure, we can clearly see that with the known value of sediment-basement density contrast, all three inversions provide good estimation of the depth-to-basement. The shallower part of my inversion result 49 shows better fitting with the true model than the deeper part, which is reasonable since with the increase of the depth, the sensitivity of the gravity data decreases. One can see that full-tensor inversion produces better result over the inversion of the vertical gravity component. I did not observe any significant improvement of the joint inversion result comparing to the inversion of full-tensor gravity gradiometry data. 6.1.5 Inversion for depth-to-basement and density contrast As we have seen in the previous examples, the simultaneous inversion for depth-to- basement and sediment density is a highly nonunique problem. In this inversion, I set the boundaries for the density contrast from 0 to 1 g/cm3, which corresponds well to the known density contrast values in the Ensenada Bay basement model (Gallardo et al., 2003). Constraint in the density contrast was implemented by doing the inversion in a logarithmic weighted-model space. One needs to notice that the Frechlet derivatives with respect to depth and density contrast are interdependent. In order to avoid the numerical singularity problem, I need to start my inversion from a nonzero depth-to-basement and density. As in the case of the inversion with the known density contrast, I ran the inversions for three different data sets: a) vertical component of the gravity field, gz, only; b) full-tensor gravity-gradient data; and c) joint inversion of the vertical gravity field and the full-tensor gravity-gradient data. For the inversion of each data set, I considered different initial density contrasts to test the stability and robustness of the inversion algorithm. I used seven different initial values of the density contrast ranging from 0.1 g/cm3 to 0.7 g/cm3 with a spacing of 0.1 g/cm3. Table 6.1 presents a list of recovered density contrast values from the gz inversion using different initial values. For all the different selections of initial density contrast, the final density contrast recovered from inversion converges to a reasonable range of 0.29 g/cm3 to 0.33 g/cm3. The initial density contrast surface is a flat surface at a depth of 1 m above the reference sediment-basement interface in order to avoid singularity. Figure 6.20 presents a comparison of the true model and inversion results along a profile at y = 0. The initial value of the density contrast for the results shown in these 50 figures is 0.1 g/cm3. One can clearly see that the depth-to-basement is well recovered from the inversion of full tensor gravity gradiometry data and the joitn inversion. However, I do not observe any improvement of joint inversion result compared to the inversion of full tensor gravity gradiometry data due to the low resolution of vertical gravity component. From Figure 6.20, we can see the estimated density contrast surface smoothed out some variations in the true model. Inversion of the full-tensor gradient data produced a model which is very close to the true model. With the unknown density contrast value, the inversion using vertical gravity component to recover both the depth-to-basement and density contrast value becomes a serious nonunique problem. The gravity gradiometry data themself have higher resolution than the vertical gravity component and the inversion of full tensor gravity gradiometry components can further reduce the model uncertainty. From this figure, we can also also find that the inversion of full tensor gravity gradiometry data produces a better result than the joint inversion of vertical gravity component and full tensor gravity gradiometry data. The recovered density contrast value is 0.325 g/cm3 for the gz inversion, 0.385 g/cm3 for the full-tensor inversion, and 0.373 g/cm3 for the joint inversion of gz and the full-tensor gravity-gradient data. All of these recovered densities are very close to the true value of 0.4 g/cm3. 6.2 Modeling and inversion of magnetic data for the sediment-basement interface model In this section, I present two model studies for the modeling of the magnetic data for the sediment-basement interface. In my models, I assume that the horizontal reference plane P is located at a depth of 1000 m under the earth's surface. In the first model, the variable sediment-basement interface, r , is located above the horizontal plane P , with a maximum elevation of -800 m. Model 2 is a reconstructed model of the Ensenada Bay basement from Gallardo-Delgado et al., (2002). For this model, the left portion of the sediment-basement interface is above the reference plane P , while the right portion is below plane P . The actual interface extends from -1400 m 51 to -600 m in the z direction. For these models, all scalar components of the magnetic field and the total magnetic intensity were computed. The forward modeling results based on the Cauchy-type integral approach are compared with the results obtained using traditional prismatic discretization. 6.2.1 Model 1 In Model 1, the actual sediment-basement interface is a surface above the reference plane P (z - -1000). The observation stations are located at the earth's surface at z - 0. Figure 6.21 provides an illustration of this model with the sediment-basement interface. The black dots indicate the observation stations on the earth's surface with the spacing of 100 m in x and y directions. I assumed that the inducing magnetic field was tilted from the vertical direction and had all three nonzero scalar components as follows: Hj - (30000, 20000, 40000) nT. (6.5) I also assumed that the sediments were nonmagnetized and the magnetic suscepti­bility of the basement was 0.01 in an SI unit. All three scalar components, x, y, and z, of the magnetic field components and the total magnetic intensity were computed for this model using Cauchy-type integral representations and the conventional method. Since the total magnetic intensity data are commonly used in geophysical exploration, I analyzed only the TMI data for all models in this dissertation. Figure 6.22 shows a comparison of the TMI data computed using the Cauchy-type integral (solid blue curve) and the conventional method (red stars) at y - 0. From this figure, we can see that the fields computed using these two methods practically coincide. The fitting between the fields computed from these two different methods in this model is actually better than the previous gravity models due to the finer discretization used here. The computation based on the Cauchy-type integral was very fast since only the sediment-basement interface was discretized instead of the 3D discretization required for the conventional method. 52 6.2.2 Model 2: Ensenada Bay model In this section, I will present the results of inversion of the magnetic data for the Ensenada Bay sediment-basement interface model. The geological settings for the Ensenada Bay area have been discussed in the previous sections. I reconstructed a model based on the vertical geological sections as shown in Figure 6.18. For this model, I assumed a vertical magnetization. The magnetic susceptibility of the sediment was 0 and the magnetic susceptibility of the basement was set to 0.01 in the SI unit. I computer-simulated the total magnetic intensity data with 5% random noise added. The simulated data were used as synthetic observed data for inversion to recover the depth-to-basement interface. In the inversion, the tolerance was set to 5%, which was the noise level of the synthetic data. The inversion converged to this tolerance level after just a few iterations. Figure 6.23 presents the data-fitting for the inversion. Figure 6.24 presents a comparison of the true model and the inversion result at y - 0. We can see that the depth-to-basement is recovered well, especially in the left part, where the actual sediment-basement interface is shallower than the reference plane P. In the right part, the recovered depth-to-basement is slightly underestimated. This can be explained by the decrease in the sensitivity of the data to the basement as the depth increased. 6.3 Modeling and inversion of MT data for the sediment-basement interface model In this section, I will demonstrate my inversion algorithm using several realis­tic synthetic models of the sediment-basement interface. I first assumed that the sediment conductivity was known and the inversion was for the depth-to-basement estimate only. I will also consider the inversion for both the depth-to-basement and the conductivity of the sediments. We should note, however, that in practical applications, one should apply a conventional 3D inversion of the MT data first in order to determine the volumetric distribution of the conductivity in the subsurface. The inverse model produced by conventional MT inversion can be used to create the initial model for the depth-to-basement estimate using the developed novel algorithm. 53 The model represents a sediment-basement interface with asymmetric shape (Fig­ure 6.25) and with a maximum depth of 600 m. The conductivity of the basement is 0.001 S/m, while the conductivity of the sediments is 0.05 S/m. Figure 6.26 shows a vertical cross section of the conductivity distribution for this model. We will discuss below the results of the inversion of the principal MT impedance components. The synthetic MT data were computer-simulated at 425 MT stations located on a rectangular grid as shown in Figure 6.25 using the IE method with 5% random noise added. 6.3.1 Inversion for the depth-to-basement estimate Now, I consider the inversion for the depth-to-basement estimation uisng MT data with known sediment conductivity information. I also tried different starting models for the inversion. 6.3.1.1 Flat surface as starting model For this model, I first assumed that the sediment conductivity was given and the inversion was run for the depth-to-basement only. I used the synthetic observed MT impedance data with frequencies uniformly distributed from 0.1 Hz to 100 Hz in logarithmic space (4 frequencies per decade). The initial model was a horizontal plane at a depth of 300 m. The inversion process was terminated after 30 iterations, when the misfit between the observed and predicted data reached the noise level, as shown in Figure 6.27. Figure 6.28 shows a vertical section of the inversion result with the yellow circles representing the recovered interface and the black stars showing the true sediment-basement interface. As one can see, the inversion did a good job everywhere with the exception of the very bottom of the interface, where the maximum depth recovered from inversion was 568 m, while the true maximum depth of the sedimentary basin was 600 m. Figures 6.29 shows a comparison between the observed and predicted data at a frequency of 1 Hz. From the figure, one can clearly see that the observed data are 54 very close to the predicted data. The data fitting for other frequencies was also very good. 6.3.1.2 Conventional M T inversion result as starting model In a practical case, the conventional MT inversion result can be used to construct an initial model in order to speed up the inversion and reduce the uncertainty. In the next step, the depth-to-basement was updated in the inversion by using the new developed algorithm. Figure 6.30 shows a vertical section of the conventional MT inversion result at y=0. From this figure, one can see that the geometry of the sedimentary basin is well recovered. However, the conductivity distribution is very diffusive and it is hard to determine the sediment basement interface. I manually picked a conductivity-contrast surface (dashed black line with circles) from the conventional MT inversion. The depth-to-basement inversion with the selected surface as initial model is shown in Figure 6.31. One can see that both the shape of the sedimentary basin and its depth are well recovered. The maximum depth recovered from this inversion is 596 m, which is very close to the true maximum depth of 600 m. 6.3.2 Inversion for both the depth-to-basement and the sediment conductivity In this section, I will invert the MT field data to estimate both the depth-to-basement and sediment conductivity. Similar to the previous section, I first use a flat surface as the starting model for the depth-to-basement. Following this, I consider the use of conventional MT inversion result as the starting model. 6.3.2.1 Flat surface as starting model For this model, the inversion was also applied for both the sediment-basement interface and the sediment conductivity. The frequencies are uniformly distributed from 0.1 Hz to 100 Hz in logarithmic space (4 frequencies per decade), and the data were contaminated by 5% random noise. I used a flat surface at 300 m depth and a sediment conductivity of 0.1 S/m as the initial model. 55 Figure 6.32 shows a vertical section of the inversion results with the yellow circles representing the recovered model and the black stars indicating the true sediment-basement interface. The maximum depth recovered from the inversion was 573 m, while the actual maximum depth of the sedimentary basin was 600 m. The inverted sediment conductivity converged to a value of 0.0437 S/m, which was very close to the true value of 0.05 S/m. Figure 6.33 shows the convergence plot for this inversion, in which the normalized misfit reached the noise level after 75 iterations. 6.3.2.2 Conventional M T inversion result as starting model Similar to the previous inversion to recover the depth-to-basement with known conductivity, I ran the inversion to recover both the depth-to-basement and the sediment conductivity with the conventional MT inversion result as an initial model. The initial average sediment conductivity was estimated to be 0.03 S/m from the conventional MT inversion. Figure 6.34 shows the recovered depth-to-basement with the conventional MT inversion result as an initial model. One can see that both the geometry and the depth of the sedimentary basin are well recovered. The recovered sediment conductivity is 0.0454 S/m, which is very close to the true value of 0.05 S/m, for this inversion. One also needs to note that the normalized misfit converged to the noise level after around 35 iterations, which is much less than the previous inversion with a flat surface as initial model. By comparing Figure 6.32 and Figure 6.34, one can see that obvious difference between the true model and recovered model at all locations on Figure 6.32 although the maximum depth is well recovered. The recovered model shows a better fitting with the true model in most of the locations in Figure 6.34 comparing to Figure 6.32. 6.4 Joint inversion of MT and potential-field data for depth-to-basement estimation In the previous section, we have observed the serious nonuniqueness problem of inverting the depth-to-basement and the physical properties of the sediment using 56 a single geophysical data set, especially using gravity data to invert the depth-to-basement and the density contrast between sediment and basement. Potentially, the nonuniqueness problem can be reduced by the joint inversion of gravity field and MT field data. In this section, I will present the synthetic model studies for the joint inversion of MT and gravity field data to recover the depth-to-basement and physical properties of the sediments. In the first subsection, I assumed that the density and conductivity of the sediments were given and that the depth-to-basement would be estimated from both the gravity and MT data. Following this, I assumed that the depth-to-basement, sediment conductivity, and the density contrast were unknown. I expected these parameters to be recovered simultaneously by the joint inversion of the gravity and the MT field. 6.4.1 Inversion for depth-to-basement only In this subsection, I consider the same sediment-basement interface model shown in Figures 6.12 and 6.25 for the gravity and MT problems. The conductivity of the sediment is 0.05 S/m and the density contrast value is 0.4 g/cm3. To be realistic, I have contaminated the model parameters by 5% random noise to produce the geological nosise. Figure 6.35 shows the vertical section of density contrast value and conductivity distribution, at y-0, with geological noise. For this inversion, I assumed that the density contrast and the conductivity values were all given and only the depth-to-basement would be recovered using the joint inversion of the gravity and MT field data. The starting model is a flat surface located at 300 m depth. I have tried other reasonable starting models and all of them produced almost the same result. Figure 6.36 shows the inverted depth-to-basement compared with the true model at y-0. The maximum depth recovered from the inversion is 560 m, which is very close to the true maximum depth of 600 m. One can see that the recovered depth-to-basement is much closer to the true model in comparison with the inversion from the gravity data in Figure 6.13. However, we do not see any improvement in the inversion result compared to the result from the inversion of the MT data shown in Figure 6.28. As a result, for such a simple model with given density contrast 57 and conductivity, the inversion of the MT data only was capable of producing a high-resolution image of the sedimentary basin. 6.4.2 Inversion for depth-to-basement, density, and conductivity Here I consider the same exact model as the one in the previous subsection for the gravity and MT problem. The initial model used for this inversion is a flat sediment-basement interface located at 300 m depth. The initial sediment conductivity was set to be 0.1 S/m, and the initial density contrast was set to be 0.5 g/cm3. It took around 60 iterations for the inversion to converge to the noise level of 5%. Figure 6.37 shows a comparison between the true model and the recovered model using joint inversion. One can clearly see that the shape of the sedimentary basin and the depth are well recovered by the joint inversion even though the values of the density and conductivity are not well known. The final recovered sediment conductivity and the density contrast are 0.056 S/m (true value of 0.05 S/m) and 0.403 g/cm3 (true value of 0.4 g/cm3). Compared to the previous subsections where either the gravity or MT data were inverted separately, the recovered sediment conductivity and density contrast values are very close to their true values from the joint inversion. z(m) 58 y(m) x(m) Figure 6.1. Synthetic sediment-basement interface model with a symmetric shape for the gravity problem. The black dots on the surface indicate the gravity stations. 59 Figure 6.2. Representation of the density contrast distribution for a symmetric model with a constant density contrast, using prismatic cells at y _ 0. mGal 60 gz at y=0 x(m) Figure 6.3. Comparison of forward modeling results obtained using Cauchy-type integral (dotted line) and traditional volume integral (solid line) methods for a symmetric model with a constant density contrast. Depth(m) 61 x(m) Figure 6.4. Inversion result for a symmetric model with a constant density contrast at y= 0 using the vertical gravity component gz. 62 Figure 6.5. The density contrast profile changes linearly with depth for a sedimen-tary- basement interface model with a symmetric shape. 63 Figure 6 .6 . Representation of the density contrast distribution for a symmetric model with density contrast changes linearly with depth, using prismatic cells at y _ 0. mGal 64 gz at y=0 x(m) Figure 6.7. Comparison of forward modeling results obtained using Cauchy-type integral (dotted line) and traditional volume integral (solid line) methods for a symmetric model whose density contrast value changes linearly with depth. 65 Figure 6 .8 . The density contrast profile changes exponentially with depth for a sedimetary-basement interface model with a symmetric shape. 66 Figure 6.9. Representation of the density contrast distribution, using prismatic cells, for a symmetric model whose density contrast changes exponentially with depth. mGal 67 gz at y=0 x(m) Figure 6.10. Comparison of forward modeling results obtained using Cauchy-type integral (dotted line) and traditional-volume integral (solid line) methods for a symmetric model whose density contrast value changes exponentially with depth. Depth(m) 68 x(m) Figure 6.11. Inversion result for a symmetric model whose density contrast value changes exponentially with depth at y-0 using the vertical gravity component gz. z(m) 69 y<m) x(m) Figure 6.12. Synthetic sediment-basement interface model with asymmetric shape for the gravity problem. The black dots on the surface indicate the gravity stations. Depth(m) 70 x(m) Figure 6.13. Inverted depth-to-basement for an asymmetric model using the gz component of the gravity field with known density contrast. Depth(m) 71 x(m) Figure 6.14. Inverted depth-to-basement for an asymmetric model using full tensor gravity gradiometry data with a known density contrast. ,ui)i»d3C| 72 x(m) Figure 6.15. Inverted depth-to-basement for an asymmetric model using the component of the gravity field with an unknown density contrast. ,ui)i»d3C| 73 x(m) Figure 6.16. Inverted depth-to-basement for an asymmetric model using full tensor gravity gradiometry data with an unknown density contrast. Height (km) 74 Figure 6.17. 2D cross section showing the basement, the sedimentary basin, and the bathymetry topography of Ensenada Bay (after Gallardo-Delgado, 2003). 75 y(m> x(m) Figure 6.18. Reconstructed model for the Ensenada Bay basement. The surface in the figure is the actual boundary for sediment and basement, and the referenced sediment-basement interface is a horizontal plane located at z=-1000 m. The black dots are the observation stations for recording the gravity response. 76 (a) (b) -500 I----------------1---------------- 1---------------- r- E. -1000 -1500 ---------------- '---------------- '---------------- '---------------- '---------------- '---------------- -6000 -4000 -2000 0 2000 4000 6000 (c) Figure 6.19. Comparison of inversion results (with given density contrast value) for the true sediment-basement model of Ensenada Bay along a profile at y=0. Panels a, b, and c present the inversion results for three different data sets: a) vertical component of the gravity field, , only; b) full-tensor gravity-gradient data; and c) joint inversion of vertical gravity field and full-tensor gravity-gradient data. In each panel, the solid blue line indicates the true model whereas the red dots show the inversion result. 77 Figure 6.20. Comparison of inversion results (with unknown density contrast value) for the true sediment-basement model of Ensenada Bay along a profile at y=0. Panels a, b, and c present the inversion results for three different data sets: a) vertical component of the gravity field, gz, only; b) full-tensor gravity-gradient data; and c) joint inversion of vertical gravity field and full-tensor gravity-gradient data. In each panel, the solid blue line indicates the true model whereas the red dots show the inversion result. 78 Figure 6.21. Model 1. A realistic model of the sediment-basement interface with the positive anomaly. The lower surface in the figure is the actual boundary between the sediments and basement with a maximum value at -800 m. The reference plane for the sediment-basement interface is located at -1000 m. The black dots are the observation stations for recording the magnetic response. 79 Figure 6.22. Comparison of TMI data computed using the Cauchy-type integral approach and the conventional method for Model 1 with titled magnetization. y(m) y(m) 80 Observed data -6000 4000 -2000 0 2000 4000 6000 nT x(m) Predicted data x(m) Figure 6.23. Comparison of the observed and predicted TMI data for the inversion of the magnetic data for the Ensenada Bay model. The upper panel shows the observed data whereas the lower panel shows the predicted data. 81 Figure 6.24. Comparison between the true sediment-basement interface and the interface recovered from inversion for the Ensenada Bay model along the profile at y = 0. The solid blue curve shows the true model whereas the red stars indicate the inversion result. 82 Figure 6.25. The sediment-basement interface with asymmetric shape. The loca­tions of the MT stations are shown by the red dots. 83 0.05 0.045 0.04 0.035 0.03 0.025 0.02 0.015 0.01 0.005 -2500 -2000 -1500 -1000 -500 0 500 1000 1500 2000 2500o (S/m) x(m) Figure 6.26. A vertical cross section at y = 0 of the sediment-basement interface with asymmetric shape. The black line indicates the actual sediment-basement interface, whereas the prismatic approximation of the interface is shown by the dark red color, reflecting the conductivity of the sediments of 0.05 S/m on the corresponding color scale. Normalized misfit 84 Iteration Number Figure 6.27. Inversion with given conductivity of the sediments: a convergence plot of the normalized misfit functional. (U l )Z 85 Figure 6.28. Inversion with given conductivity of the sediments: a vertical section of the inversion result at y = 0 with the yellow circles representing the recovered model, and the black stars showing the true sediment-basement interface. y(m) y(m) 86 Figure 6.29. Inversion with given conductivity of the sediments: a comparison between the observed and predicted impedance data at a frequency of 1 Hz. 87 Figure 6.30. A vertical section of the conductivity distribution at y=0 from the conventional MT inversion result. The dashed black line with circle indicates the sediment-basement interface that I manually picked from the conventional inversion. 88 Figure 6.31. Inversion with given conductivity of the sediments: a vertical section of the inversion result at y = 0 with the yellow circles representing the recovered model and the black stars showing the true sediment-basement interface. A conventional MT inversion result was used to construct the initial model. 89 Figure 6.32. Inversion with unknown conductivity of the sediments: a vertical section of the inversion result at y = 0 with the yellow circles representing the recovered model and the black stars showing the true sediment-basement interface. Normalized misfit 90 Iteration Number Figure 6.33. Inversion with unknown conductivity of the sediments: a convergence plot of the normalized misfit functional. 91 Figure 6.34. Inversion with unknown conductivity of the sediments and the con­ventional MT inversion used as a starting model: a vertical section of the inversion result at y = 0 with the yellow circles representing the recovered model and the black stars showing the true sediment-basement interface. 92 Figure 6.35. Vertical sections, at y=0, of the sedimentary basin model used for the joint inversion of gravity and MT data. The upper panel shows the density distribution while the lower panel indicates the conductivity distribution. The model parameters are contaminated by 5% random noise. 93 Figure 6.36. Recovered depth-to-basement at y=0 for synthetic Model 2 with known conductivity and density using the joint inversion of gravity and MT field data. The yellow circles represent the recovered model, and the black stars show the true sediment-basement interface. 94 Figure 6.37. Recovered depth-to-basement at y=0 for synthetic Model 2 with unknown conductivity and density using the joint inversion of gravity and MT field data. The yellow circles represent the recovered model, and the black stars show the true sediment-basement interface. 95 Table 6.1. List of recovered density contrast values from the inversion of gz with different intial values of density contrast. The relative error between the recovered density contrast from inversion and the true value is calculated. Intial density (g/cm3) Recovered density (g/cm3) Relative error(%) 0.1000 0.3253 18.67 0.2000 0.3260 18.50 0.3000 0.2900 27.50 0.4000 0.3100 22.50 0.5000 0.3050 23.77 0.6000 0.3228 19.30 0.7000 0.2965 25.87 CHAPTER 7 CASE STUDIES In this chapter, I will show several case studies based on the gravity data released by the USGS in the Big Bear Lake area and the geological models. The isostatic Bouguer data are made available in this area to study the basin structures. In the first section, I will discuss the inversion of the gravity data in this area to recover the depth-to-basement with given information for the density contrast profile with depth. As I mentioned before, MT data can also be used for this application to produce a higher resolution of sedimentary structures and the physical properties such as the sediment conductivity. However, MT data are not made available in this area. I simulated the synthetic MT data based on recovered basin structures from the gravity data. 7.1 Inversion of gravity data at the Big Bear Lake Area In the previous chapters, I have demonstrated the effectiveness of the 3D Cauchy-type integral in the modeling and inversion of gravity and gravity-gradiometry data by several synthetic model studies. In this section, I will describe the application of the developed modeling and inversion algorithm to the field gravity data in Big Bear Lake for the estimation of depth-to-basement. 7.1.1 USGS gravity survey in the Big Bear Lake area Gravity surveys are widely used for basin study. The depth-to-basement can be well estimated based on isostatic Bouguer gravity data since a gravity anomaly is caused primarily by the density contrast between the sediments and the basement. Many gravity measurements were made in the 1960s and 1970s by various groups in order to produce gravity maps covering California at a scale of 1:250,000 for the California Division of Mines and Geology (Roberts et al., 2002). The USGS also 97 conducted a new gravity survey in the Big Bear Lake area. The new survey data were merged with the previous gravity survey to produce a new gravity grid (Roberts et al., 2002). We should note that in this paper, I have gridded and used for the inversion the data from the new USGS survey only. The USGS applied the conventional prism inversion method to the combined new gravity data to recover the depth-to-basement. In their inversion, the subsurface was discretized to a grid of prisms whose horizontal size was 2000 m by 2000 m. The density distribution along each column of prisms was assumed to be known from the well-log data, and the thickness of the prisms was determined by fitting them to the isostatic Bouguer gravity anomaly. The USGS inversion was well constrained by the well-log data and bedrock locations. In addition, at several locations, the thickness of the prisms was assumed to be known and stayed unchanged during the inversion (Roberts et al., 2002). Due to a data ownership issue, the USGS released only the new data they collected; the well-log data were not made available. 7.1.2 Geological background o f the Big Bear Lake area The Big Bear Lake area is located in the southeast part of California. The area is characterized by a deep sedimentary basin surrounded by uplifted bedrock. The USGS produced a basin model from the surface geology, well-logs, and potential field data. Figure 7.1 shows that the whole basin area can be divided into three parts from the northeast to the southwest: Deadman Lake Basin, Surprise Spring Basin, and Joshua Tree Basin. The average depth and density variations between sediment and bedrock may be slightly different. Figure 7.2 presents a digital elevation model of the area. From the surface geology, we can observe three fault belts trending from the northwest to the southeast and one fault belt trending from the west to the east. Different basins in this area are separated by these four main fault structures. 98 7.1.3 Processing o f the USGS data Figure 7.3 presents the released USGS data with the locations of the gravity stations shown by the black dots. As one can see, the original gravity data were collected in an irregular grid. It is well known that gridded data have a significant advantage over scattered data for inversion in terms of robustness (Cordell, 1992), because having regular gridded data helps produce a robust inversion result. There are different gridding methods available. The traditional mathematical gridding approach can produce significant artifacts, especially in areas with a few observation stations. I used a gridding approach based on the equivalent-source concept (Cordell, 1992). According to this concept, in the first step, I determined an equivalent layer with some surface density distribution recovered based on the inversion of the data collected in an irregular grid. In the next step, I computed the gravity data at the regular grid using the equivalent layer as the source. Note that the gridded gravity data can be used directly for inversion if we assume that the isostatic Bouguer anomaly is caused purely by a deficiency in the density of the sediments. By making this assumption, we assume that the density of the bedrock is the same as in the reference density model of the earth's crust. However, in a real case, the density of the bedrock may be different from the reference model. Therefore, the isostatic Bouguer gravity anomaly can be written as a sum of the bedrock component and the sediment component: g = gb + gs. (7.1) The bedrock gravity component, gb, can be estimated initially based on the gravity data observed on the bedrock (Roberts et al., 2002). Figure 7.4 shows the gridded bedrock component. One can see that in the southern part, there is a strong negative anomaly for the bedrock gravity component. The gridded bedrock component of the gravity anomaly was subtracted from the gravity grid in Figure 7.3 to obtain the gravity anomaly caused by the sediment only. Figure 7.5 shows the gravity anomaly obtained after removal of the bedrock component. This grid represents the final data that I used for inversion. One needs to note that the approximation of the bedrock component of the gravity anomaly by interpolating the anomaly observed on outcrops is not a rigorous approach 99 due to the presence of a nearby sedimentary basin with low density. I used an iterative method to remove the bedrock gravity component. In my approach, the bedrock component of the gravity field was initially computed by simple extrapolation from the gravity observations on the bedrock outcrop. Inside the inversion, it is corrected based on the inverted basin depth. The corrections are terminated when there is no significant change in the bedrock component of the gravity field (Roberts et al., 2002). 7.1.4 Inversion o f the USGS gravity data One needs to know the density variation with the depth in order to get an accurate model of the depth-to-basement. As I mentioned above, this information can be obtained from well-log data. The density models of the Deadman Lake and Surprise Spring Basins are slightly different from that of the Joshua Tree Basin. The USGS report states that in the northern part, a density contrast of 400, 350, 300, 250, and 200 kg/m3 with bottom depths of 50 m, 100 m, 150 m, and 300 m is a good approximation of the basin density (Roberts et al., 2002). The USGS report also states that this model may not be suitable for the Deadman Lake Basin well since there are very limited well constraints in Deadman Lake Basin (Roberts et al., 2002). In the southern part (Joshua Tree Basin), a constant density contrast value of 550 kg/m3 is suitable (Roberts et al., 2002). The northern part (the Deadman Lake and Surprise Spring Basins) and the southern part (the Joshua Tree Basin) of the survey area were be inverted separately. In order to speed up the inversion and get the most reasonable result, the well-known Bouguer slab formula (Chakravarthi and Sundararajan, 2006) could be applied to generate an initial model: = gBA Po (7 2) 41.89ApQ + agB where gB is the Bouguer gravity anomaly; Ap0 is the density contrast between sediment and basement on the earth's surface, and this density contrast decreases in the vertical direction with the gradient a. However, my inversion algorithm does not depend on the selection of the starting model. The selection of a flat surface as a starting model produces almost the same result as using the Bouguer slab formula as a starting model. 100 In the inversion, I used a grid size of 300 m by 300 m in the x and y directions, which is much finer than the USGS model grid for prismatic inversion (2000 m by 2000 m). 7.1.5 Inversion o f the gravity data in the Deadman Lake and Surprise Spring Basins In order to take the variable density contrast into account, I need to use some analytical function of depth to approximate the density contrast. For the USGS model, I found that it was better to use equation 6.3 to approximate the true density contrast. The optimized values for the parameters in equation 6.3 are given as follows: a = 251.5 kg/m3, b = -0.007, c = 197 kg/m3, d = 5.2656 x 10-6. (7.3) Figure 7.6 presents plots of the USGS staircase density variation model and my approximation by the exponential function. The results of the inversion are shown in Figure 7.7 overlapped with the DEM (digital elevation model) and the fault structure. One can see that the northwest-southeast trending faults correspond well to the edge of the Surprise Spring and Deadman Lake Basins. The east edge of the recovered Deadman Lake Basin fits well with the mountain belt. Figure 7.8 shows an overlap of the inversion result with the USGS basin and bedrock models. In this figure, one can see that the recovered location of the basin is similar to the USGS model. Figure 7.9 shows a comparison of my inversion result with the inversion result provided by the USGS for the Deadman Lake and Surprise Spring Basins. One can see that the recovered basin geometry obtained by my method