3D SIMULATION OF MAGNETOTELLURIC DATA USING FINITE DIFFERENCE EIGENMODE METHOD

Abstract

Magnetotelluric method is used to delineate the subsurface conductivity structure of earth using natural electromagnetic waves in the frequency range 10" Hz - 10 Hz as source field. These natural fields are generated mainly by thunderstorm activity (>1 Hz) and the interaction of solar wind with the earth's magnetosphere (<1 Hz) (Kaufman and Keller, 1981). The horizontal electric and magnetic field components are measured at the earth's surface and analyzed to infer electrical resistivity distribution in the earth's interior. The two orthogonal horizontal electric field components are linearly related to the two horizontal magnetic field components through appropriate transfer function (Tikhonov, 1950 and Cagniard, 1953). The depth of penetration of electromagnetic (EM) wave depends upon its frequency and conductivity distribution of medium. The EM fields are studied using Maxwell's equations, coupled in electric (E) and magnetic field (B) vectors. These equations are transformed into vector Helmholtz equation for decoupled E-field or B-field. The vector Helmholtz equation is used to solve for the response ofa given earth model. Typical model parameters are geometry ofthe target and spatial distribution ofconductivity. The estimation ofmodel parameters from the physical fields, measured on earth surface, is termed as an inverse problem, while the mapping of model parameters to measured fields is known as a forward problem. For a good inversion algorithm, an efficient forward modeling code is needed. This work deals with the development of an efficient 3Dforward modeling algorithm. The popular numerical modeling schemes can be broadly classified into Integral Equation Methods (IEM) and Differential Equation Methods (DEM) (Finite Difference Method (FDM), Finite Element Method (FEM)). While IEM can be efficiently used only in for computing the responses of confined targets buried in a layered earth, the DEMs are capable of modeling arbitrary complex distributions of conductivity. The coefficient matrix in case of IEM is full but small in size, while in case of DEM it is large but grossly sparse. In DEMs use of staggered grid is popular, particularly in 3D case, because its use analytically incorporates the divergence equation of magnetic field. FDM with staggered grid is used in the present study. Instead of using FDM to solve the complete Boundary Value Problem (BVP) with sources, we have first studied fundamental nature of the eigenvalue problem obtained in case of source free BVP. Eigenvalues and eigenvectors, collectively known as eigenmodes, exhibit the basic characteristics of the response to a given physical property distribution in the model. After estimating the eigenmodes for a given geometry and physical property distribution, the EM response for a given source frequency can be obtained through superposition of the eigenvectors. In geophysical applications, similar approach was implemented by Druskin et al. (1994, 1999)and Stuntebeck(2003). In the eigenmode method, the responses for additional frequencies can be obtained in negligible time. In contrast, in case of traditional use of FDM to generate multifrequency responses, one has to re-run the code for each frequency. During evaluation of superposition coefficients, the eigenvalues appear in the denominator, implying that the smaller eigenvalues contribute more significantly to the field. Therefore, one need compute only a subset of the smallest eigenvalues and corresponding eigenvectors for a given degree of accuracy of field values. For evaluation ofthis subset, the iterative methods serve better than the direct methods, particularly in case of 3Dproblems where the matrix size is extremely large. IV The most widely used methods for evaluating a subset of eigenmodes are Krylov subspace projection methods. In these methods, only product ofthe matrix with a vector is needed and, therefore, only non-zero elements of the sparse coefficient matrix need be stored. Lanczos and Arnoldi methods are two popular Krylov subspace methods. Former is used for symmetric matrices while the latter is used for non-symmetric matrices. Before launching the development of3D code, we gained experience ofeigenmode method by developing ID and 2D forward modeling codes. The FDM coefficient matrix is symmetric. In case of 3D, the symmetric coefficient matrix is of large size, which is reduced in size by using the current divergence equation, to eliminate the z-component of electric field from expressions. This step transforms the symmetric coefficient matrix to a nonsymmetric one, albeit of much smaller size. So, Lanczos method is used to obtain the eigenmodes in ID and 2D case while Arnoldi method is used in case of 3D. The eigenmode evaluation subprogram of our algorithm is adapted from the routines of ARPACK (1997) software which is based on Implicit Restarted Lanczos/Arnoldi Method (IRLM/IRAM) given by Sorensen et al. (1992). ARPACK works in different modes such as 'regular', 'shift and invert' etc. The regular mode is efficient in obtaining largest magnitude eigenvalues while invert mode is efficient in obtaining smallest magnitude ones. Since we are interested in the smallest eigenvalues, shift and invert mode is used. Further, to circumvent the problem of loss of Lanczos vector orthogonality in case of degenerate eigenvalues, their complete reorthogonalization has been employed. The development of3D algorithm was carried out on a PC. As a result, we had to introduce several I/O detours and had to work under severe limitations imposed on the size of the grid. Therefore, we designed several appropriate experiments using the coarse grid to validate the 3D algorithm. The organization of seven chapters in the thesis is presented next. In chapter one, the literature review is presented. In chapter two, the theory for 3D Magnetotellurics using electric field vector Helmholtz equation, obtained from Maxwell's equations, is discussed. Different types of boundary conditions such as domain and interface boundary conditions are described. The eigenmode problem is formulated and the eigenmodes are used to obtain the EM response for multi-frequency case. Thederivations of response functions, i.e. apparent resistivity and phase corresponding to both 2D-TE and 2D-TM modes are discussed. In chapter three implementation of FDM on staggered grid is described. The domain is discretized into a grid comprising cuboids. We have followed the convention that electric field components are defined on midpoints of edges while magnetic field components are defined at the centers of surfaces. The derivation of matrix equation from the governing partial differential equation and boundary conditions is presented next. The coefficient matrix obtained is symmetric and about one third of its eigenvalues are zero. These spurious zero eigenvalues do not contribute to field synthesis. This knowledge is made use of in reducing the coefficient matrix size to the number of non-zero eigenvalues by eliminating the vertical component of electric field and working only with the horizontal components. This step transformed the symmetric coefficient matrix to a non-symmetric one. A brief review of ARPACK subprograms adapted to determine the eigenmodes is presented. In chapter four, the details of various stages of development of algorithm, MT3DEA, are discussed. Starting with symmetric matrix eigenmode evaluation using Singular Value Decomposition (SVD), the Lanczos and Arnoldi methods in 'regular' and in 'shift and invert' modes are presented. In invert mode a matrix equation need be solved so efficient matrix solvers based on Conjugate Gradient Method and various preconditioners used are described next. Finally, the algorithm is presented along with flow charts of important subprograms. In chapter five, the synthetic experiments designed to test and validate the algorithms MT_2D_EA and MT3DEA are discussed. First we performed different tests such as grid convergence and no contrast case to check the consistency and accuracy. Then, we compared the results of2D version ofour algorithm with published results. We studied two 2D models (simple and complex) taken from COMMEMI (Zhdanov et al., 1997) and obtained good match with the average values given in the paper. The RMS errors for simple and complex models are 0.01 and 0.06 respectively. Next we studied the impact on the field values of using different percentages of eigenmodes. We observed that for obtaining accurate field values, 5% eigenmodes were sufficient for the conductive block model whereas for the resistive block 20% eigenmodes were needed for same accuracy. In the multi-frequency experiment, we studied Weaver (1976) model. We used two grids for time periods Is and 10s and generated the responses (true) for these grids. Next we generated the response at Is using 10s eigenmodes and vice versa and found excellent fit with the corresponding true responses. In case of 3D, additional experiment conducted was to verify that 3D apparent resistivity values converge to corresponding 2D values as the strike length in one direction is extended. We compared our 3D response with the vu published apparent resistivity values of the model described as 3D-2 in COMMEMI report and found a good fit. In chapter six, we have used MT_3D_EA algorithm to generate 3D models whose responses are commensurate with the MT field data acquired by Israil et al. (2008) in Garhwal Himalaya. Tyagi (2007) and Israil et al. (2008) analyzed this data using WingLink software, and proposed the first 2D geoelectric model. We used this model as base model for our study. Using MT2DEA algorithm, we generated responses of this model at two time periods and found excellent match with the corresponding WINGLINK responses. Due to limited computer resources, we could not run the complex model using MTJDEA algorithm. So, for 3D study we designed a simplified 3D model retaining the dominant feature of conducting block. We generated the 3D responses for 4 strike length values (20 km, 50 km, 70 km and 100 km) of the conducting body. At 100 km strike length the 3D response of the model matches well with the 2D response. Finally, we experimented with the strike-length and the depth of burial of the block and generated equivalent 3D models that would explain the conducting anomaly in the observed data. The 3D geometry of the conductive block, buried under the Roorkee-Gangotri profile near MCT, can be taken as 70 km strike, 20-26 km width and 4 km depth and its resistivity is estimated as 8 Q-m. However, the detailed 3D study suggests that the conductive block can be approximated as a 2D one. In chapter seven, we have discussed the strategies for further improvement of our algorithm.