DEVELOPMENT OF AN EFFICIENT ALGORITHM FOR 3D MAGNETOTELLURIC DATA INVERSION

Abstract

Magnetotelluric (MT) method is a passive source method used to delineate the subsurface conductivity structure of earth. Natural electromagnetic waves in the frequency range 10􀀀5 Hz - 104 Hz are used as source elds. The horizontal electric and magnetic eld 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 eld components are linearly related to the two horizontal magnetic eld components through appropriate transfer function (Cagniard [23],Tikhonov [142]). The objective of the present study is to understand the mathematical, physical and numerical aspect of 3D MT inversion leading to an e cient 3D inversion software, 3DINV FD, for magnetotelluric data. The estimation of model parameters from the physical elds, measured on earth surface, is termed as an inverse modeling. In magnetitelluric method, the earth is parameterized in terms of electrical resistivity which is of special signi cance as it carries information about the lithology, pore

uid, temperature and chemical variations. As the EM eld is a non-linear function of subsurface resistivity distribution, the inverse problem is also non linear in nature. In the present work, the inverse problem is quasi-linearized and then solved iteratively. The inverse problem is solved using Gauss-Newton with conjugate gradient method. For each inversion iteration, a new forward problem, yielding the response of current resistivity model and several pseudo forward problems, for Jacobian matrix computations, are solved. Therefore, an e cient forward modeling algorithm is a prerequisite for an inversion algorithm. The mapping of model parameters to measured elds is known as a forward iii modeling. For generation of MT response, a boundary value problem is solved analytically or numerically. However, for the problems involving complex geometries one has to seek numerical solutions. Due to its simple mathematics and easy implementation, staggered grid nite di erence method (FDM) has been chosen over other numerical techniques for solving the MT boundary value problem. The FDM results in a matrix equation, which is then solved using Bi-Conjugate Gradient Stabilized (BICGSTAB) with DILU preconditioner to compute the MT response. The quasi-linearization of non-linear problem results in a matrix equation which is solved using Conjugate Gradient (CG) method, a semi-iterative matrix solver that dispenses with the necessity of explicit computation of Jacobian matrix. The initial guess is made on the basis of observed anomaly and other a priori information. The inversion algorithm 3DINV FD is the culmination of research that started with the development of a primitive algorithm. The algorithm has been written in FORTRAN 90 language and implemented on an Intel Core i7 3.6 Ghz machine with 4 Gbyte of RAM. 3DINV FD comprises 6887 lines having 44 subroutines and works in double precision arithmetic. The main program has two basic modules - Forward and Inverse. Its special e ciency features which result in cost e ectiveness are - (i) Optimal computational parameters for static divergence correction, (ii) BICGSTAB with DILU preconditioner, which results in fast convergence, (iii)Gaussian noise addition to synthetic data, (iv) Computation of multi frequency response in parallel using OpenMP, (v) Use of logarithm of resistivity to ensure positive values of estimated parameters, (vi)In-built computation of regularization parameter and (vii) CG matrix solver for inverse problem. Besides being e cient, 3DINV FD is versatile on account of its iv features like (i) Inversion of eld/synthetic data, (ii) Error free/erroneous synthetic data and (iii) Inversion of pro ling/sounding data. For establishing the validity of forward modeling algorithm, the published results of various models have been reproduced. The validity of the inversion algorithm 3DINV FD is established by inverting the synthetic data generated from di erent models. To ensure the stability of the algorithm the inversion is performed after adding the Gaussian noise to the synthetic data. Furthermore, to demonstrate the robustness of the algorithm, the data generated from ModEM algorithm (Kelbert et al. [61]) has been inverted successfully. The synthetic experiments designed to understand the e ect of number of sites and their distribution on the inversion, suggest that accurate resolution of the anomalous body data should be acquired along straight pro les whenever possible. And the a priori information about the target body should be taken into account for optimal site selection. The results of various experiments and inversion of synthetic have established the veracity of the algorithm and also amply displayed the capabilities of the inversion algorithm. Also discussed, is the possible scope of future work in various directions for its upgradation.