Electromagnetic inversion using quasi-linear approximation

Three-dimensional electromagnetic inversion continues to be a challenging problem in electrical exploration. We have recently developed a new approach to the solution of this problem based on quasi-linear approximation of a forward modeling operator. It generates a linear equation with respect to the modified conductivity tensor, which is proportional to the reflectivity tensor and the complex anomalous conductivity. We solved this linear equation by using the regularized conjugate gradient method. After determining a modified conductivity tensor, we used the electrical reflectivity tensor to evaluate the anomalous conductivity. Thus, the developed inversion scheme reduces the original nonlinear inverse problem to a set of linear inverse problems. The developed algorithm has been realized in computer code and tested on synthetic 3-D EM data. The case histories include interpretation of a 3-D magnetotelluric survey conducted in Hokkaido, Japan, and the 3-D inversion of the tensor controlled-source audio magnetotelluric data over the Sulphur Springs thermal area, Valles Caldera, New Mexico, U.S.A.

GEOPHYSICS. VOL 65, NO 5 (SEPTEMBER*OCTOBER P. 1301-1513.18 FIGS. Electromagnetic inversion using quasi-linear approximation Michael S. Zhdanov*, Sheng Fang*, and Gabor Hursan* ABSTRACT Three-dimensional electromagnctic inversion contin­ues to be a challenging problem in electrical exploration. We have recently developed a new approach to the so­lution of this problem based on quasi-linear approxima­tion of a forward modeling operator. It generates a linear equation with respect to the modified conductivity ten­sor, which is proportional to the reflectivity tensor and the complex anomalous conductivity. We solved this lin­ear equation by using the regularized conjugate gradient method. After determining a modified conductivity ten­sor, we used the electrical reflectivity tensor to evaluate the anomalous conductivity. Thus, the developed inver­sion scheme reduces the original nonlinear inverse prob­lem to a set of linear inverse problems. The developed algorithm has been realized in computer code and tested on synthetic 3-D EM data. The case histories include in­terpretation of a 3-D magnetotelluric survey conducted in Hokkaido, Japan, and the 3-D inversion of the ten­sor controlled-source audio magnetotelluric data over the Sulphur Springs thermal area. Valles Caldera, New Mexico. U.S.A. INTRODUCTION During the last decade, considerable advances were made in developing a multidimensional electromagnetic (EM) in­terpretation technique. Several papers were published during this period on 3-D inversion of EM data (Eaton, 1989: Mad­den and Mackie, 1989; Smith and Booker, 1991: Lee and Xie, 1993; Oristaglio et al., 1993; Pellerin el al., 1993; Nekut. 1994; Zhdanov and Keller. 1994; Torres-Verdin and Habashy. 1995; Xie and Lee. 1995; Newman and Alumbaugh. 1996). However, the development of effective interpretation schemes for 3-D inhomogeneous geological structures is still one of the most challenging problems in EM geophysics. In our recent publications (Zhdanov and Fang. 1996a. b). we developed a novel approach to 3-D EM forward modeling and inversion based on linearization of the integral equations for scattered EM fields. In this paper, we discuss a new devel­opment of this method using a regularized conjugate gradient inversion scheme, and we present the results of case history studies. For completeness, we also include an overview of the main principles of quasi-linear (QL) inversion. The QL approximation is based on the assumption that the anomalous field E" is linearly related to the normal (back­ground) field Eb in the inhomogeneous domain: E" = XEb. where X is an electrical reflectivity tensor w'hich is a slowly varying function of position and frequency. In the inversion, we introduce a modified material property tensor m. propor­tional to the reflectivity tensor and the anomalous conductivity Ao. In this case, the QL approximation generates a linear equa­tion with respect to the modified material property tensor m. After determining m. we use the electrical reflectivity tensor X to evaluate the anomalous conductivity Arr. The developed algorithm has been realized in the QLINV3D computer code and tested on synthetic 3-D EM data. Syn­thetic inversion examples, with and without random noise, have demonstrated that the algorithm for inverting 3-D EM data is fast and stable. The case histories presented here include in­terpretation of a 3-D magnetotelluric (MT) survey conducted by the New Energy and Industrial Technology Development Organization (NEDO) in the Minamikayabe area in the south­ern part of Hokkaido. Japan, and 3-D inversion of the tensor controlled-source audio magnetotelluric (CSAMT) data col­lected by Phil Wannamaker over the Sulphur Springs thermal area, Valles Caldera. New Mexico, U.S.A. The interpretation of real data demonstrates that the 3-D quasi-linear EM inversion can be a practical tool of geophysical EM data analysis. QUASI-LINEAR INVERSION Consider a 3-D geoelectric model with the normal (back­ground) conductivity oh and local inhomogeneity D with an ar­bitrarily varying conductivity a = ah+An. The model is excited Manuscript received by the Editor June 12,1998; revised manuscript received January 24,2000. ‘University of Utah. Department of Geology and Geophysics. Salt Lake City, Utah 84112-0111. E-mail: m/hdanovt®mines.utah.edu. {Formerly University of Utah. Department of Geology and Geophysics, Salt Lake City, Utah 84112-0111; presently Baker Allas. 10201 Westheimer WA-IA. Houston. Texas. E-mail: [e-mail]. © 2000 Society of Exploration Geophysicists. All rights reserved. 1501 1502 Z h d a n o v e t al. by an electromagnetic field generated by an arbitrary source. We consider a quasi-stationary model of the field, so we ne­glect the displacement currents as is usual in low frequency geophysical applications (Zhdanov and Keller, 1994). The elec­tromagnetic fields in this model can be presented as a sum of the background and anomalous fields: E = Eb + Ea, H = Hb + Ha, (1) where the background field is a field generated by the given sources in the model with the background distribution of con­ductivity ab, and the anomalous field is produced by the anoma­lous conductivity distribution Act. It is well known that the anomalous field in the frequency domain can be presented as an integral over the excess currents in the inhomogeneous domain D (Weidelt, 1975): Fa(Fj) = ///„ CF(Fj 1 F)'^(F)rEb(r) + ^ where F" stands for Ea or Ha observed on the surface of the earth, and GF(rj | r) stands for the electric or magnetic Green's tensor defined for an unbounded conductive medium with the normal conductivity ab. We apply the quasi-linear approach (Zhdanov and Fang, 1996a, b), : ' :i'' Ea(r) % A(r)Eb(r), to equation (2), which yields ^ * /// ' r)ACT(r)[! + ^r)]EbWdv' (4) where X is an electrical reflectivity tensor, and I is the unit tensor. ,K We introduce a new tensor function r . ni(r) = Arr(r)[I + i(r)], ' (5) which we call a modified material property tensor. Equation (5) can be treated as a modification of Ohm's law: jD = ActE = Aa(Eb + Ea) = Act [I + A(r)]Eb. (6) The modified material property tensor can be considered as the ratio of the anomalous, induced currents, jD, to the background electric field. Finally, equation (4) takes the form Fa(rj) « Jjj GF(rj | r)m(r)Eb(r)dv = GF(ih), (7) where GF is the corresponding Green's linear operator. The last equation is a linear one with respect to ift(r). In our algorithm, we adopt the regularized conjugate gradient method for determining the modified material property tensor (Zhdanov, 1993). The reflectivity tensor X can be determined from the follow­ing linear equation inside the inhomogeneous domain D, as long as we know fn : Ea(rj)~///DeE(rj|r)lh(r)Eb(r)dv ■■ :: ^ A(rj)Eb(rj). (8) After determining m and X, it is possible to evaluate the anomalous conductivity distribution Act from equation (5). This inversion scheme reduces the original nonlinear inverse problem to three linear inverse problems: the first one (the quasi-Born inversion) for the parameter rii, another one for the parameter a, and the third one (a correction of the result of the quasi-Born inversion) for the conductivity Act . ft is based on a QL approximation; that is why we call this approach aQI. inversion. Note that we can rewrite equation (7) using matrix notations: F = GFm. (9) Here m is the vector-column of the modified conductivity ten­sor m, F is the vector-column of the data, and the matrix GF is the matrix of the linear Green's operator defined by formula (7). Similarly, we can get from (8) , : XEb = GEm, . ■ (10) where X is a block-diagonal matrix of reflectivity tensors, Eb is the vector-column of the background field, and GE denotes a matrix of linear operators defined by expression (7) with the electric Green's tensor. Equation (5) in matrix notations takes the form .... . in = Act11 • (11) where I is the unit matrix. The solution of the inverse problem is reduced to the inver­sion of linear system (9) with respect to m and then to the in­version of equation (10) with respect to X. After that, we obtain Act using equation (11). We apply the Tikhonov regularization theory to solve this system of linear equations (Tikhonov and Arsenin, 1977). THE REGULARIZED CONJUGATE GRADIENT METHOD FOR SOLVING LINEAR INVERSE PROBLEM EQUATIONS Note that the solution of equation (7) and, correspond­ingly, equation (9) can be nonunique in the general case since there may exist m(r) distributions that produce a zero exter­nal field, similar to the existence of the current distributions producing zero external field (Habashy et al, 1994; Svetov and Gubatenko, 1985). However, according to the uniqueness the­orem by Gusarov (1981) the excess conductivity Acr(r) can be found uniquely in the class of piecewise-analytic functions. Therefore, to make our solution unique, we have to develop a unified approach for simultaneously finding m, X, and Act. We use the method of constrained inversion developed by Zhdanov and Chernyak (1987). A similar approach to 2-D in­verse scaltering problems was discussed also by Kleinman and van den Berg (1993). It is based on the introduction of the parametric functional P"(m) = </>(m) + a S(m), (12) where the misfit functional is specified as : / 0(m) = ||GFm - F|[2 + ||m - Act(I + A)||2, (13) (3)r and the stabilizer is selected to be equal to S(m) = ||m - mapr||2. Q u a s i-L in e a r (14) Tlie a priori model iiiapr is some reference model, selected on the basis of all available geological and geophysical informa­tion about the area under investigation. The scalar multiplier a is a regularization parameter. As one can see from equation (13). the misfit functional takes care not only of the solution of equation (9) but. simultane­ously, also fulfilling of the modified Ohm's law (5). Now. the inversion is reduced to the solution of the mini­mization problem for the parametric functional: Pu(m) = min. (15) One can calculate the first variation of the parametric func­tional under the assumption that A a and X are, temporarily, constants: ^"(m) = 2Re{<$m'[GF*(GFm - F) + (m - Aa[I + X)) + or(m - mapr)]|. where the asterisk denotes transposed complex conjugate matrices. Let us select 3m as Am = -kY'(m), 0 < k < oo. (16) where l"(m) = GF*(Glm-F) + (m-A(7[I + X])-fo'(in-mapr), (17) and k is the length of the iteration step. This selection makes SPa(m) = -2X"Re|l"'(m)l"(m)) < 0. That means that the parametric functional is reduced if we apply perturbation (16) to the model parameters. In order to accelerate the convergence, we apply the regu­larized conjugate gradient method described by the following iteration process: mn-i ~ u>n - m„ &nl (ni„). where the "directions" of ascent r (mn) are selected according to the algorithm described below. In the first step, we use the "direction" of regularized steepest ascent determined by equation (17): I i"(m,i) = l°(nio). Therefore, the first iteration is computed by the formula ni| = m,, + <5mo = mu - *i,[GF*(GFmapr - F)]. (18) where Hlapr = A<7apr[I + Aapr|. and Aorapr and Xa|)r are the parameters of the reference model. The second iteration is where in: = ni| +<Sni| = ni| - /til"(mi). I"(mi) = !•' (m I) +^|I"(mr,) (19) and 1" (mi) = GF*(GFm, - F) + (m, - Anapr|I + X,]) + a(m, - mapr). (20) On the (n + 1 )th step mn+i = m„ -t-<5m„ = mn - k„i"(mn), (21) and f"(mn) = r'(mn) + /JBJ"(mn-i), where l"(m„) = GF'(GFmn - F) + (mn - Arrn_,[I a„]) + a(mn-mapr). (22) The matrix of the reflectivity tensor X„ is determined from the matrix m„ using equation (10): EM In v e rs io n 1503 X„E - G m„. (23) The anomalous conductivity on the (n - l)th iteration can be found using equation (II). Note that the anomalous con­ductivity has to satisfy the condition A<t„_i > -<rb. (24) because the electrical conductivity has to be positive. Notice also that equation (11) should hold for any frequency. In reality, of course, it holds only approximately. Therefore the conduc­tivity AcTn_i can bc found by using the least-squares method of solving equation (II): A<7n_| = Ref ^Tm^a + Xn-.r > ft» ' xT^a-t-x^ra+*„_!)] L in . = a„_|. for a„_i >-fTb. -l and Aa„_| = -^<rb for an_| < -wh. (25) (26) The coefficient k„ can be determined from the condition /,"(nin. |) = P'(m„ - A„i"(mn)) = f(<.'0) = min. Solution of this minimization problem gives the following best estimate for the length of the step: •l"(mn)la(mI1) lc,,(mn)(GF*GF + ufl)l‘'(mn) (27) The conjugate gradient method requires that the vectors i''(m„) will be mutually conjugate. This requirement is fulfilled if the coefficients /),, are determined by the formula (Tarantola. 1987): /».= [| 1" (mn) | ~ ||l-(ni„_i)|1504 Z h d a n o v e t al. Using equations (17), (18), and (27), we can obtain m itera­tively. The regularization parameter a describes the trade-off be­tween the best fitting and reasonable stabilization. In a case where a is selected as too small, the minimization of the para­metric functional />c'(iti) is equivalent to the minimization of the misfit functional </>(ni). Therefore we have no regulariza­tion, which can result in an unstable incorrect solution. When a is loo large, the minimization ol the parametric functional P"(m) is equivalent to the minimization of the stabilizing func­tional iS(m), which will force the solution to be closer to the a priori model. Ultimately, we would expect the final model to be exactly like the a priori model, while the observed data are totally ignored in the inversion. Thus, the critical question in the regularized solution of the inverse problem is the selection of the optimal regularization parameter a. The basic princi­ples used for determining the regularization parameter a are discussed by Tikhonov and Arsenin (1977) and Zhdanov and Keller (1994). The traditional numerical method to determine the param­eter a is based on the following simple idea. Consider for ex­ample the progression of numbers ak = anq k. k =0,1,2,q > 0. (28) For any number ak we can find an element m„k, minimizing P'-1*- (in), and calculate the misfit |GKniUk - F||2. The optimal value of the parameter a is the number a(o, for which we have C/in lo,k0 = s, (29) where S is the level of noise in observed data. Equality (29) is called the misfit condition. lifipii ! IOUr ■ Border of the measured area 19 equidistant profiles 19 equidistant points on each profiles 200 m 100 Ohm-m 200 m 200 m 200 m 540 m Fig. 1. 3-D model of a synclinal structure in the near-surface conductive layer excited by a plane wave, and its subdivision into cells used for the inversion.Q u a s i-L in e a r EM In v e rs io n 1505 In our code, we use the following algorithm of the regularized conjugate gradient method: I"-i (m„) = !"■+> (mn+,) + A,+J"" (mj. where a„ are the subsequent values of the regularization pa­rameter. This method is called the adaptive regularization method (Tikhonov and Arsenin, 1977). In order to avoid di­vergence, we begin an iteration from a big value of a (e.g„ cto = 10000.), then reduce an (ur„ =u„ j/2) on each subsequent iteration, and continuously iterate until the misfit condition (29) is reached. MODEL STUDY OF 3-D QUASI-LINEAR INVERSION METHOD The method was carefully tested on a number of 3-D mod­els describing different typical geoelectrical situations. We will present here just one example of the model study. The model represents a synclinal structure in the near­surface conductive layer with 10 ohm-m resistivity. The hor­izontal and vertical cross-sections of the model are shown in Figure 1. The model was excited by a plane EM wave. The the­oretical observed field was calculated al the nodes of a square grid on the surface using the integral equation (IE) forward modeling code SYSEM (Xiong, 1992). TTie distance between the observation points is 30 m in the X and Y directions. Following the traditional approach used in practical MT ob­servations, we calculated synthetic observed apparent resis­tivities and phases based on two off-diagonal elements of the magnetotelluric tensor at each observation point. The quanti­ties pyx and <t>,x are assigned to the nominal TE mode, whereas pIV and 0,v are assigned to the nominal TM mode. Note that this 2-D nomenclature is artificial and approximate in nature for 3-D structures. However, it is widely used in practical MT observations, and usually only the quantities />„, and <f>n are available for inversion. That is why we use the same approach in the model study. Magnetotelluric apparent resistivities and phases were cal­culated for the frequencies f = 0.5, 1, 2, 5, 10, 20, 50, 100. 200. and 500 Hz. Thus the synthetic MT data set which we use for inversion corresponds to 361 observation points and 10 fre­quencies. As an example, the synthetic observed magnetotel­luric sounding curves at a point (x = 0. y = 90 m) are shown in Figure 2 by solid lines. We defined an inverse area as a rectangular body with sides of 400 m in both horizontal directions and of 450 m in a vertical direction. It was divided into rectangular cells with constant but unknown resistivities. The boundaries of the cells used for inversion are shown by dashed lines in Figure 1. We applied QL inversion simultaneously toTE and TM data with 3% random noise added. We used the adaptive scheme described above for selecting the optimum regularization pa­rameter a. The initial value of a{) was reduced in the process of conjugate gradient minimization of the parametric functional (15) by a factor of two. The behavior of the misfit functional is a good indicator to measure the efficiency of the inversion pro­cess. Figure 3 shows that the normalized misfit can be reduced to as low as 3% in this case, which corresponds to the level of noise in the data. The corresponding quasi-optimal value of the regularization parameter for the final iteration was aq„ = 0.05. The CPU time for numerical calculations using QL inversion was 120 minutes Let us look at the most important result of the inversion: the resistivity distribution. Figure 4 shows the 3-D inversion image of the predicted model by inversion carried out using TE and TM data simultaneously. One can clearly see in I his figure the configuration of the conductive synclinal structure. Figure 2 shows the comparison between synthetic observed data and predicted data for the inversion result presented in Figure 4. One can see that predicted data based on QL approx­imation fit very well the initial (synthetic observed) data be­cause the inversion was based on the QL method. At the same time, there is a small discrepancy between predicted data based on the fE method and the initial (synthetic observed) data. This discrepancy reflects actually the level of QL approximation ac­curacy in the forward modeling solution. It demonstrates also the practical limitations of QL inversion: this inversion can be used for geoelectrical models with relatively small conductivity contrasts between different geoelectrical structures. According to Zhdanov and Fang (1996a). QL approximate forward mod­eling is valid for a conductivity contrast up to 100. Thus, wc should conclude that QL inversion works well only within the limits of QL approximation. Fig. 2. Typical magnetotelluric sounding curves at x =0, y = 90 m above the initial and predicted model for the synthetic data inversion. We plot three curves for each mode: (i) initial synthetic observed data computed using the IE method (solid lines), (2) predicted data computed for inverse model using QL approximation (dashed lines), and (3) predicted data com­puted for inverse model using the IE method (lines formed by crosses). TE mode phase 1 10 100 Frequency [Hz] Observed Predicted QL Predicted full IE TM mode phase -901--------------------------------- -100 -110 -120 --------- Observed ---------Predicted QL ♦ + Predicted lull IE -160 -170 -1801---------------------------- 1 10 100 Frequency [Hz] -1501506 Zhdanov et al. In summary, numerical modeling and inversion results demonstrate that QL inversion can produce reasonable im­ages of 3-D geoelectrical structures. We will discuss below the two case histories of interpreting MT and CSAMT data using QL inversion. FlG. 3. Behavior of the misfit functional during the inversion for the synthetic data contaminated by 3% noise. The New Energy and Industrial Technology Development Organization (NEDO) has been conducting a "Geothermal Development Promotion Survey" in the Minamikayabe area located in the southern part of Hokkaido, Japan (Takasugi et al., 1992). This area is particularly interesting because of its geothermal potential. In 1988, MT sounding and AMT sound­ing were conducted in this area. Figure 5 shows the location of the MT measurement sites. The site spacing is 100 m along and across observational profiles (Figure 5). The frequency range of MT data is from 1 to 130 Hz. Because the survey area at Mi­namikayabe is located very close to Uchiura Bay on the north­eastern side, low frequency data are strongly affected by the "coast effect" (Takasugi et al., 1992). Since we arc interested in the shallow geological structure in this area, it is reasonable to restrict ourselves to frequencies which are greater than 1 Hz so that the coast effect can be neglected. We use MT data for inversion at seven different frequencies between 1 and 100 Hz (96.48, 24.12,6. 3, and 1.5 Hz). In order to reasonably choose the background model for the inversion, we first apply one­dimensional inversion to all observed data and approximate the background model as a three-layer model (the resistivi­ties are 480. 6, and 150 ohm-m; the thicknesses are 30 and 450 m). We use both nominal TE mode and TM mode data for a total of 161 MT soundings in the inversion. In this area, we expect to find a geothermal high conductivity zone in the intrusive resis­tive rock. This assumption is reflected in our three-layer back­ground model. The subdivisions of the inverse area are shown in Figure 6. Il was divided into 16 x 16 cells in a horizontal KAYABE CASE STUDY IN THE MINAM1KAYABE AREA, JAPAN 41:57:29 41:56:13 * 140:53:17 140:53:39 140:54:2 140:54:24 140:54:47 Longitude (East) 0 500 1000 METERS Ohm-m FlG. 4. Volume image of the inversion result for the synthetic model. Fig. 5. The distribution of MT sounding sites in the Mi­namikayabe area. Stars denote the positions of MT sites.Q u a s i-L in e a r EM In v e rs io n 1507 plane (each cell has is 100 tn x 100 m size) and seven layers in a vertical direction. This is consistent with the spacing of 100 m between the observational points. After 320 iterations of the QL inversion scheme we reached the misfit of 8%. The CPU time for this computing was about 200 minutes. Figures 7 and 8 show the inverse results. Figure 7 presents the general 3-D images of the model, whereas Fig­ure 8 shows vertical slices along three profiles. The 3-D inver­sion is consistent with the 2-D interpretation results included in Takasugi et al. (1992); that is, there exists a low resistivity area (middle layer) in the intrusive resistive rock formation. As a general tendency, the resistivity varies with depth approx­imately as high-low-high, and also a resistive zone exists within the low resistivity layer, which is in agreement with the result of Takasugi et al. (1992). However, 3-D inversion produces much more detailed information about the internal structure of the inverse area. Horizontal discretization of the inverse area with the locations of the MT stations g MT-profile 13 MT-profile 12 MT-profile 11 MT-profile 10 MT-profile 9 MT-profile 8 MT-profile 7 MT-profile 6 MT-profile 5 MT-profile 4 MT-profile 3 MT-profile 2 MT-profile 1 -500 0 500 Y [m] Vertical discretization of the inverse area 100 200 ^300 £ N 400 500 600 700 Fig. 6. The subdivision of the inverse area for Japanese MT data (16 x 16 x 7). Top panel is plan view; bottom panel is ver­tical cross-section. The MT stations are marked by stars on the horizontal cross-section (top panel). -500 0 Y[m] 500 I UJ z 500 -500 1- ‘ ! ' t | niooo 316 100 32 10 l: Ohm-m Fig. 7. 3-D resistivity images of the inverse results obtained by QL inversion of the MT data. Predicted model below Profile 2 Fig. 8. The vertical slices of the inverse results for the Japanese MT data along Profiles 2.5. and 9. 1 1000 100 L U , Ohm-m i Ohm-m -500 0 500 Predicted model below Profile 9 1000 100 10 i Ohm-m I I -500 0 5001508 Z h d a n o v e t al. The comparisons of original data with the predicted data are presented in Figures 9 and 10. They show good agreement between the observed and predicted apparent resistivity and phase pseudosections. 3-D INVERSION OF TENSOR CSAMT DATA COLLECTED OVER THE SULPHUR SPRINGS THERMAL AREA, VALLES CALDERA. NEW MEXICO, U.S.A. The Valles Caldera. New Mexico, U.S.A. is one of the Quaternary rhyolite systems of North America. To examine its volcanic history an extensive CSAMT survey was carried out over the Sulphur Springs geothermal area in the caldera (Wannamaker. 1997a, b). Tlie detailed geology of the measured area and the receiver system is shown in Figure 11. Two inde­pendent electric bipoles were used as a transmitter at 13-km distance from the measured area. Frequency soundings were carried out along four profiles at 45 receiver points marked in Figure 11 (Nl, N2, SI, and El). The frequency range was 4-4096 Hz. The tensor CSAMT methodology was based, according to Wannamaker (1997a, b), on measuring five EM field compo­nents (£,. Ey. Hx. Hy, H:) for two independent polarizations of source current bipoles. Plane-wave impedances and H. ten­sor elements were obtained at each frequency. For example, the off-diagonal impedance elements were computed by the formulas: ~ Ex\Hx2 , Z,v = „ „-----„ „ . (30) HxlHyl - HrtHyl and Zyx = Ey2 Hy I - Ev\H y\ny2 HX I Hy2 - H x2 Hy 1 (31) where subscripts 1 and 2 refer to the source bipole number. Apparent resistivities (ptv and pyx) and impedancc phases (0rv and <pyx) were calculated in the usual manner. The quantities pxy and <f>xy were assigned to the nominal TE mode, whereas py, and 4>yx were assigned to the nominal TM mode. Wannamaker (1997a, b) noted that while this 2-D terminology was only ap­proximate, it still was consistent with the strike directions of the major geological structures beneath the observational lines (Figure 11). For our 3-D quasi-linear inversion, we were able to use only these quantities, apparent resistivities (ptv and pvt) and impedance phases (</>xy and (p, v), because the original EM field components data were not available. At the same time, numer­ical analysis conducted by Wannamaker (1997a, b) and our own numerical results have demonstrated that nonplane-wave effects associated with the resistive crystalline basement are Observed p^1 along profile 5 -200 0 200 Y position [m] Observed p]|* along profile 9 Observed <DyI along profile 5 -600 -400 -200 0 200 400 500 Predicted 4?1 along profile 5 -600 -400 -200 0 200 400 600 <*•«'*• Y position [m] Fig. 9. Comparison of apparent resistivity and phase (TM mode) pseudosections between the observed and predicted data along profile 5. Predicted along profile 9 Y position [m] Fig. 10. Comparison of apparent resistivity and phase (TM mode) pseudosections between the observed and predicted data along profile 9. -400 -200 0 200 400 600Ohm-n? Y position [m] Observed ^ along profile 9I Q u a s i-L in e a r EM In v e rs io n observed in the CSAMT data at 8 Hz and lower frequencies. In other words, theoretical CSAMT sounding curves, computed using formulas (30) and (31) practically coincide with the the­oretical MT soundings curves for the geoelectrical structures typical in the caldera. To demonstrate this fact, we present the results of forward modeling for three I -D models which can be considered as typical in the caldera. Models A and B are three-layer models with 12-8-500 ohm-m resistivities and 250-750 m thicknesses, and model C represents a live-layer model with 458-165-7-500-20 ohm-m resistivities and 32-200-750-100000 m thicknesses. We calcu­lated theoretical MT and CSAMT apparent resistivity and phase curves versus frequency for these models. Figure 12 presents the ratio of these curves. One can see that down to / = 8 Hz, the ratios of the corresponding MT and CSAMT ap­parent resistivities and phases arc equal to one. So. we can state that the observed CSAMT data can be treated as plane-wave data. That is why it was reasonable to interpret this data on the 1509 Qvsm-i ■;Qvsm. Fault ;WC23-4 Border of the inverted areal lfAfA,A*AtVA*A*A'± [m] E1 Ir CSAMT profile (dots mark electrode positions, pluses mark selected sites) VC-2B Drill collar location Surficial acid - sulfate alteration VMXXA*:-. a i * jl ■ a i li Surficial deposits undivided Caldera-fill deposits Valles rhyolite Qvrc - Redondo Creek Qvsm - San Antonio Ml. Fault, square on down-thrown side, dashed where inferred Limits of detailed mapping of Goff and Gardner (1980) Caldera structural margin (ring-fracture system) Bandelier Tuff Fig. 11. Detailed geological map of the Sulphur Springs area in the Valles Caldera (after Wannamaker, 1997a) with the receiver profiles and the border of the area used for the inversion shown bv thick black lines.1510 Z h d a n o v e t al. basis of conventional MT formulas in a similar way, as was done for the MT survey in the previous case history. We used for in­version the CSAMT data along three profiles Nl, N2. and SI (Figure 11). Observed frequency sounding TE- and TM-mode CSAMT apparent resistivity and phase pseudosections along the central line (Nl) are shown in Figures 13 and 14. The subdivision of the inverse area in the rectangular cells is shown in Figure 15. We divided the area into 20 by 5 cells in a horizontal plane (see lop panel in Figure 15). Each cell is 500 m long in the X direction and 200 m long in the Y direction, which corresponds to the minimum distance of 200 m between the receivers along the observation lines and the distance of 700 m between the profiles. Note also that the length of the electric receiver bipoles is 250 m. Taking into account the spac­ing of the observational system, it seems unreasonable to select the horizontal size of the inversion grid cells less than 200 m. We used seven different depth levels of the volume grid. The vertical sizes of the substructures increase with depth (bottom panel in Figure 15). After 500 iterations of the QL inversion scheme which re­quired 150 minutes of CPU time, we reached the misfit of 7%. Figure 16 shows the results of inversion for three verti­cal slices passing through the observational lines. Figures 13 and 14 represent also the predicted pseudosections for the TE- and TM-mode apparent resistivity and phase along the central line of the survey. Note that the result of the 3-D inversion is consistent with the previous2-D inversion results (Wannamaker. 1997a). Figure 17 shows 2-D interpretation result obtained by Phil Wannamaker along the central line. One can see strong correlation between these two results. We can clearly locate the typical geoelec­trical layers which correspond to the known geological struc­tures of the area. For example, we observe a resistive Bandelier Tuff in the western part of the cross-section and a low resis­tive Paleozoic sedimentary rocks below the tuff. Also, we can clearly see in the results of 3-D inversion the conductive struc­ture at the depth of about 600-1000 m which can be associated with the Sulphur Creek fault in the western part of the cross­section (Figure 17). Note that the 3-D inversion results produce much more detailed information about deep geological structures of the caldera than the 2-D inversion. For example, similar to 2-D results presented in Wannamaker (1997a), we observe the low resistivities of the Bandelier Tuff southeast of Sulfur Creek under all three profiles. However, the downdrop in the Paleo­zoic rocks along the Redondo Border fault is expressed more strongly in the N1 and N2 profiles, and it is not as clearly seen in the SI profile (Figure 16). In contrast to the 2-D interpretation, one can clearly see the high-resistive zone in the upper 500 m in the central part (near horizontal coordinates .t = 0 and y = 0) of the 3-D inverse image (Figure 16. middle panel; Figure 18. bottom panel). It may possibly correspond to a vapor zone near the well VC-2B (Figure 11) as implied from acid-sulphate Frequency [Hz] Fig. 12. Ratios between MTand CSAMT apparent resistivities and phases versus frequency for three synthetic models. Predicted ■long profile N1 -500 0 500 V position (m) Fig. 13. Observed and predicted TE-mode apparent resistivity and phase pseudosections along the central profile (Nl) of the CSAMT survey. D100 To 1500Ohm-m Observed pjj* along profile Nl Predicted pjj* along profile N1 >1000 -600 0 500 1000 Y position [m] 1500 degree Observed <X>*y along profile Nl -1000 -500 0 500 1000Q u a s i-L in e a r EM In v e rs io n 1511 alteration and subhydrostatic well pressures (Wannamaker, 1997a). Note that this zone was not detected by 2-D interpre­tation of CSAMT data. CONCLUSION We have developed a rapid 3-D electromagnetic inversion algorithm based on the QL approximation of the forward mod­eling. The main advantage of the method is that we reduce the original nonlinear inverse problem to a set of linear inverse problems to obtain a rapid 3-D EM inversion. The QL inverse problem is solved by a regularized conjugate gradient method which ensures stability and rapid convergence. The inversion of the real data indicates that 3-D quasi-linear EM inversion is fast and stable. Tlie inverse results provide reasonable re­covery of the models and real geological features. The main limitation of QL inversion is that it can work well only within the limits of QL approximation. Practically, it means that QL inversion can provide an accurate estimate of the geoelectrical model with a conductivity contrast up to 100. ACKNOWLEDGMENTS Financial support for this work was provided by the National Science Foundation under grant EAR-9614136. The authors acknowledge the support of the University of Utah Consortium of Electromagnetic Modeling and Inversion (CEMI), which includes Advanced Power Technologies, BHP, INCO, Japan National Oil Corporation, Mindeco, MIM Exploration, Naval Research Laboratory, Newmont, RioTinto, Shell International Horizontal discretization of the inverse area with the locations of the data profiles 1000 500 £ 0 -500 -1000 -1000 0 1000 Y[m] 2000 NW------ SE Vertical discretization of the inverse area 200 - 400 N 600 800 1000 -1000 Y[m] 1000 2000 Fig. 15. Tlie subdivision of the inverse area for CSAMT data (5 x 20 x 7). The locations of CSAMT stations used for the inversion are marked by stars on the horizontal cross-section (top panel). These data are obtained by interpolation lo Ihe uniform grid. Predicted model below Profile S1 Observed oyT along profit* N1 -4098 £.1024 J 256 I M -1500 -1000 -500 0 500 Predicted along profit* N1 0 500 Y poaltlon [m] Fig. 14. Observed and predicted TM-mode apparent resistivity and phase pseudosections along the central profile (N l) of the CSAMT survey. 200 E, 400 -C l & H 600 1 1 ft 1 800 1000 m 1100 " 10 3 -1000 0 1000 2000otmwn Predicted model below Profile N1 © ° 800 100 32 10 3 -1000 0 1000 2000 °hm-m Predicted model below Profile N2 S. 600 <D ° 800 1000 100 32 10 3 -1000 0 1000 Y position [m] 2000 Ohm-m Fig. 16. Vertical cross-seclions of the predicted model ob­tained by QL inversion of the CSAMT data. B|l» 10 l500Ohm-m Observed pJ* along proflla N1 -500 0 500 1000 Pr*dlct*d pJ* along proflla Nl -1000 -500 0 500 1000 Y poaltlon [m] 1500Ohm-m1512 Z h d a n o v e t al. r T i VflCR-SULFSPR-Nl 0 j km Fig. 17. Result of the 2-D interpretation along the central line of the survey Nl (after Wannamaker, 1997b). Ohm-m Ohm-m Fig. 18. 3-D resistivity images of the inverse results obtained by QL inversion of the CSAMT data. Exploration and Production. Schlumberger-Doll Research, Western Atlas, Western Mining, Unocal Geothermal Corpo­ration, and Zonge Engineering. We are also thankful to Dr. Takasugi, GERD. Japan, for providing the 3-D MT data. We would like to thank Dr. P. Wannamaker for providing us the tensor CSAMT data and useful discussions. We thank the anonymous reviewers for their useful sugges­tions and comments. REFERENCES Eaton. P. A.. 1989,3-D electromagnctic inversion using integral equa­tions: Geophys. Prosp., 37,407-426. Gusarov. A. L.. 1981, About the uniqueness of the solution of the mag- netotelluric inverse problem, in Mathematical models in geophysical problems: Moscow University Press, 31-80. Habashy. T. M., Oristaglio, M. L., and de Hoop, A. T„ 1994, Simulta­neous nonlinear reconstruction of two-dimensional permittivity and conductivity: Radio Sci.. 29, 1101-1118. Kleinman, R. E., and van den Berg, P. M.. 1993, An extended range- modified gradient technique for profile inversion: Radio Sci.. 28. 877-884. Lee. K. H., and Xie, G.. 1993, A new approach to imaging with low frequency EM fields: Geophysics. 58.780-796. Madden. T. R.. and Mackie, R. L.. 1989, Three dimensional magne­totelluric modeling and inversion: Proc. IEEE, 77,318-332. Nekut. A.. 1994. EM ray trace tomography: Geophysics. 59.371-377. Newman. G. A., and Alumbaugh, D. L., 1996,'3-D EM model­ing and inversion on massively parallel computers: Sandia Report SAND96-0852. Oristaglio, M. L.. Wang, T.. Hohmann. G.. and Tripp, A.. 1993, Resis­tivity imaging of transient EM data by conjugate-gradient method: 63rd Ann. Internat. Mtg.. Soc. Expl. Geophys., Expanded Abstracts. 347-350. Pellerin. L.. Johnston. J.. and Hohmann. G„ 1993, 3-D inversion of EM data: 63rd Ann. Internal. Mtg., Soc. Expl, Geophys., Expanded Abstracts. 360-363. Smith. J. T., and Booker, J. R., 1991. Rapid inversion of two- and three­dimensional magnetotelluric data: J. Geophys. Res.. 96. no. B3.3905- 3922. Svetov, B. S., and Gubatenko, V. P.. 1985. About die equivalence of die system of the extraneous electric and magnetic currents: Radiotech- niqucs and Electronics, no. 4.31-40. Takasugi S., Keisaku T., Noriaki K., and Shigeki M„ 1992, High spa­tial resolution of the resistivity structure revealed by a dense net­work MT measurement-A case study in the Minamikayabe area. Hokkaido. Japan: J. Geomag. Geoelectr.. 44,289-308. Tarantola, A.. 1987, Inverse problem theory: Elsevier. Tikhonov. A. N„ and Arsenin. V. Y., 1977, Solution of ill-posed prob­lems: W. H. Winston and Sons. Torres-Verdin, C„ and Habashy, T. M.. 1995, An overview of the ex­tended Bom approximation as a nonlinear scattering approach and its application to cross-well resistivity imaging: Progress in Electro­magnetic Research Symposium Proc., 323. Wannamaker, P., 1997a, Tensor CSAMT survey over the Sulphur Springs thermal area. Valles Caldera. New Mexico. U.S.A.. Part I: Implications for structure of the western caldera: Geophysics, 62, 451-465.Q u a s i-L in e a r EM In v e rs io n 1513 ---------, 1997b. Tensor CSAMT survey over the Sulphur Springs ther­mal area, Valles Caldera, New Mexico, U.S.A., Part II: Implications for CSAMT methodology: Geophysics, 62,466-476. Weidelt, P., 1975, EM induction in three-dimensional structures: J. Geo­phys.. 41,85-109. Xie. G., and Lee, K. H.. 1995, Nonlinear inversion of 3-D electromag­netic data: Progress in Electromagnetic Research Symposium Proc., 323. Xiong, Z„ 1992, EM modeling of three-dimensional structures by the method of system iteration using integral equations: Geophysics, 57, 1556-1561. Zhdanov, M. S.. 1993. Regularization in inversion theory, tutorial: Col­orado School of Mines. Zhdanov, M. S, and Chernyak, V. V., 1987. An automated method of solving the two-dimensional inverse problem of electromagnetic induction within the earth: Transaction (Doklady) of the USSR Academy of Sciences, Earth Science Sections, 296, no. 4-6, 59­63. Zhdanov, M.S., and Fang, S., 1996a, Quasi-linear approximation in 3-D EM modeling: Geophysics, 61,646-665. ---------, 1996b, 3-D quasi-linear electromagnetic inversion: Radio Sci., 31,741-754. Zhdanov, M. S., and Keller, G. V., 1994, The geoelectrical methods in geophysical exploration: Elsevier.