Integral electric current method in 3-D electromagnetic modeling for large conductivity contrast

Abstract-We introduce a new approach to 3-D electromagnetic (EM) modeling for models with large conductivity contrast. It is based on the equations for integral current within the cells of the discretization grid, instead of the electric field or electric current themselves, which are used in the conventional integral-equation method. We obtain these integral currents by integrating the current density over each cell. The integral currents can be found accurately for the bodies with any conductivity. As a result, the method can be applied, in principle, for the models with high conductivity contrast. At the same time, knowing the integral currents inside the anomalous domain allows us to compute the EM field components in the receivers using the standard integral representations of the Maxwell's equations. We call this technique an integral-electric-current method. The method is carefully tested by comparison with an analytical solution for a model of a sphere with large conductivity embedded in the homogenous whole space.

1282 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 45, NO. 5, MAY 2007 Integral Electric Current Method in 3-D Conductivity Contrast Michael S. Zhdanov, Vladimir I. Dmitriev, and Alexander V. Gribenko Abstract-We introduce a new approach to 3-D electromagnetic (EM) modeling for models with large conductivity contrast. It is based on the equations for integral current within the cells of the discretization grid, instead of the electric field or electric current themselves, which are used in the conventional integral-equation method. We obtain these integral currents by integrating the current density over each cell. The integral currents can be found accurately for the bodies with any conductivity. As a result, the method can be applied, in principle, for the models with high- conductivity contrast. At the same time, knowing the integral currents inside the anomalous domain allows us to compute the EM field components in the receivers using the standard integral representations of the Maxwell's equations. We call this tech­nique an integral-electric-current method. The method is carefully tested by comparison with an analytical solution for a model of a sphere with large conductivity embedded in the homogenous whole space. Index Terms-Electromagnetic (EM) modeling, high conductiv­ity contrast, integral equations. 1. Introduction NE OF THE difficult problems in electromagnetic (EM) modeling is accurate numerical solution for models with large conductivity contrast. This problem appears, for example, in modeling EM data for mineral exploration when we have a conductive target embedded in relatively resistive host rocks. The study of the topography effect on EM data requires the solution of a similar problem, because the contrast in con­ductivity between the conductive earth and nonconductive air can be as large as 10s -1010 times. Well-logging is another area where one should take into account a strong contrast between the cased borehole, for example, and surrounding rock formations. Manuscript received November 28, 2005; revised November 21, 2006. This work was supported by the University of Utah Consortium for Electromagnetic Modeling and Inversion (CEMI), which includes Baker Atlas Logging Services, BGP China National Petroleum Corporation, BHP Billiton World Explo­ration, Inc., ENI S.p.A., ExxonMobil Upstream Research Company, INCO Exploration, Norsk Hydro, Petrobras, Rio Tinto-Kennecott, Schlumberger Oilfield Services, Shell International Exploration and Production, Inc., Statoil, Sumitomo Metal Mining Company, and Zonge Engineering and Research Organization. M. S. Zhdanov and A. V. Gribenko are with the University of Utah, Salt Lake City, UT 84112 USA. V. I. Dmitriev is with the Department of Computational Mathematics and Cybernetics, Moscow State University, Moscow 119992, Russia. Digital Object Identifier 10.1109/TGRS.2007.893562 In this paper, we introduce a new approach to the solution of this problem based on the integral-equation (IE) method. The basic principles of the IE method were outlined in the pioneer papers [8], [9], [12], [17], [18], and [22], Over decades, these methods were further developed and improved in a large number of publications (see, for example, [2], [6], [7], [II], [14], [19], [21], [23], [24], and [26]). However, most existing IE methods fail for large conductivity contrast, because they use the boxcar basis functions to approximate the electric field within the conductive body [10], [11], [16], [19], The development of accurate EM modeling methods for the models with large conductivity contrast is considered one of the most challenging problems in EM geophysics. The conventional IE algorithms are usually written for the electric field or electric current components within the domain with anomalous conductivity. This domain is divided in the number of cells, which are selected to be so small that the field components vary slowly within the cell. If the conductivity of the body and/or frequency are high, it is difficult to satisfy this condition. The EM field varies extremely fast within a good conductor, which may result in errors of numerical modeling. In order to overcome this difficulty, Newman and Hohman used a special grouping of the boxcar basis functions to form current loops within the conductor [16], Farquharson and Oldenburg implemented the more sophisticated edge element basis func­tions to avoid the inaccuracy of the conventional boxcar basis function approach [10], In this paper, we consider a novel approach for solving this problem. We develop a new form of the IE method, which is based on the equations for integral current within the cells, instead of the electric field or electric current themselves. We obtain these integral currents by integrating the current density over each cell. The integral currents can be found accurately for a body with any conductivity. We do not use anymore the requirements that the field varies slowly inside the cell, because we deal with the integral of this field. As a result, the method can be applied, in principle, for mod­els with arbitrary conductivity contrast. At the same time, knowing the integral currents inside the anomalous domain allows us to compute the EM field components in the re­ceivers using the standard integral representations of Maxwell's equations. We call this technique an integral-electric-current (1EC) method. We will present below the detailed description of the 1EC method and will illustrate it by numerical modeling results. 0196-2892/S25.00 © 2007 IEEE Authorized licensed use limited to: The University of Utah. Downloaded on May 17,2010 at 22:44:50 UTC from IEEE Xplore. Restrictions apply. ZHDANOV et al.: IEC METHOD IN 3-D EM MODELING FOR LARGE CONDUCTIVITY CONTRAST 1283 II. Formulation of the IE Method the centers of the cells D„ of the anomalous domain D We consider, first, the basic IEs of 3-D EM forward model­ing, written for the total electric E and magnetic H fields E(r') H(r') G^r^r) ■ [A<j(r)E(r)] dv + E (V) ■■GE [Ac?(r)E(r)] + E (V) Gh (r1r) ■ [Ac?(r)E(r)] dv + H (r') -Gh [Aff(r)E(r)] +Hb(r') (1) (2) „ __ r rpl rp2 rp.N rpl rp.2 rp.N rpl rp2 rp.N]T eD - iPx J ^x J ' ' ' ^x J ^y J ^y J • • • ^y J ^z J ^z J • • • ^ z J -D r rpbA. rpb.2 fP:b,N rpbA. rp.b,'. x ■, -c'Ve ■, ■ ■ ■ -c'Ve J ^y ' ^y pb,N pb, 1 p6,2 pb.N]T •C/y z J •C/z i - ■ ■ ^z J • These vectors have the order 3N. The 3N x 3N matrix Gd is formed by the volume integrals over the elementary cells Dn of the components of the corre­sponding electric Green's tensor Ge, acting inside domain D G r,pn ^ad (6) where G£;(rj|r) and Gjy(rj|r) are the electric and magnetic Green's tensors defined for an unbounded conductive medium with the complex background conductivity ab = er - iue; Ge and Gh are corresponding Green's linear operators; and Eb, IT1' are the background electric and magnetic fields; domain D corresponds to the volume with the anomalous conductivity distribution a(r) = + Acf(r), r € D. Equation (1) written for the points r' located inside domain D, r' € D, gives us an IE with respect to electric field E(r). The main problem is to solve this IE. The conventional approach to discretization of the integral (1) is based on dividing domain D into N elementary cells, Dn, formed by some rectangular grid in the domain D = U *=1Dn, and assuming that Ac?(r) has the constant value Aan within the cell. Note that the coefficients Aan can be represented as the components of a vector cr of the order N a = [A<7i, Ac?2, • • • , AC?Ar]T where superscript "T" denotes transposition. We also assume that each cell Dn is so small that the electric field is approximately constant within the cell, E(r) « E(r„), where rn is a center point of rectangular cell Dn. Under this condition, (1) takes the form where (jjjf't-n.ir,, r}t!r. a, j3 = x, y, z; p, n = 1,2,... N. Note that (3) or equivalent matrix (4) provides an ad­equate approximation of the original IE, if the following conditions hold. 1) The linear size h of elementary cell Dn is much smaller than the wave length Ab of the EM field in the background medium h <C Ab- (7) 2) h is much smaller than the wave length Aa of the EM field in a medium with anomalous conductivity h An. (8) The first condition (7) usually holds for typical geophysical EM modeling problems. The second condition may fail in the case of high anomalous conductivity, which is the subject of this paper. E(rp) £/ n=l J 1 GE(rp\r)dv ■ A-TnE(r„) ■ p = 1,2,... N. (3) Thus, inside the anomalous domain D, the discrete analog of (1) can be written as [27] en = Gi,nv.i, ■'D (4) where 5 is a (3N x 3N) diagonal matrix of anomalous conductivities a = diag(A<7i,..., AaN,a-[,..., A aN, ai,..., i\oN) (5) eo and are the vectors of the total and background electric fields formed by the x, y, and z components of these fields at III. IEs for Integral Currents Our goal is to construct a discrete analog of integral (1), which would provide an accurate approximation only under condition (7). We will consider, first, (1) and (2), written for calculation of the EM field in the receivers located outside domain D. Let us denote by R the minimal distance from the receivers to domain D with the anomalous conductivity: R = minreD lr/ - r|, where r' is the observation point. We assume that the linear size h of elementary cell Dn is much smaller than the distance to the observation point r' (9) and that it is also much smaller than the wave length Ab in the background medium [condition (7)]. Under these conditions the Green's tensor G^r'lr) slowly varies inside cell Dn, if point r moves within this cell, and the observation point r' is far Authorized licensed use limited to: The University of Utah. Downloaded on May 17,2010 at 22:44:50 UTC from IEEE Xplore. Restrictions apply. 1284 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 45, NO. 5, MAY 2007 away from Dn. It is possible, therefore, to write the integral representations (1) and (2) in the form N n=1 N H(r> * £ n=1 GE(r'[ Dn j(r )dv j(r)dv ■ E (: (10) D „ Hb(r') (11) where r„ is the center point of rectangular cell Dn, and j(r) is the excess electric current Note that integral j(r) = Acr(r)E(r). j(r')<fo' = I„ (12) (13) D„ is the IEC within elementary cell Dn. Substituting (13) into (10) and (11), we obtain b/„/\ E(r')=£Gi3(r'|r„)-I„+Eb( n=l N H(r') £ G a (V r„ ) • l„ • H''( b/„/\ (14) (15) domain D, r' € D, requires holding condition (8), which may fail in the case of high anomalous conductivity. That is why we need to modify (16). At the same time, we should emphasize that in the case of the domain (16) and (1), we do not need to impose condition (9), which was used only in developing formulas (14) and (15) for calculation of the EM field in the receivers located outside domain D from the IEC given inside domain D. In order to obtain an equation with respect to integral cur­rents, we integrate both sides of (16) over elementary cell Dp and assume that anomalous conductivity is constant within the cell Dp, Act = Actp Actp£ D„ G£;(r'|r)dt/ Dp •j(r)efe + Ib (17) where Ib is the integral current in the cell Dp due to background field Eb Ip = I I I j" (r')dv' = j j j ACT(r')Eb(r')dt/. (18) Dp Dp Let us introduce the notation G/.-lr' rh/r' = G(r). (19) D„ Thus, the EM field components can be calculated in the re­ceivers, if we know the integral currents within the cells of the grid. Note also that we obtained formulas (14) and (15) without imposing any restriction on the behavior of the electric field or electric current within the elementary cell of the grid. In other words, we do not use condition (8), and therefore the conductivity of the anomaly can be arbitrarily high. This result indicates that we can develop the IE method for the models with high-conductivity contrast, if the corresponding equations are written not for the electric field E(r„) but for the integral currents I„. In order to obtain a system of linear equations with respect to electric currents, let us multiply both sides of (1) by Act|V), assuming that Ao^r') ^ 0 within domain D. As a result, we have j(r') = A5(r') JJJ G/.-(r' r) • j(r),/r • jb(V) (16) D where jb(r') = ACT(r')Eb(r')! r' C 1) is the induced current due to background field Eb. We should note that integral (16) with respect to the electric currents, in principle, equivalent to the original integral (1) with respect to the electric field. These equations were analyzed in many publications (e.g., [1], [5], [18], and [22]). The discretiza­tion of (16), similar to (1) written for the points r' located inside Note that integral G_E.P(r) represents smoothing of the Green's tensor, and it is a relatively slow varying function. Therefore, we can take this expression outside of the integral over Dn in formula (17) Ip ACTp 'y G l-\p( 1'/, I I/, ' I|: (20) where in the case of a slow varying field G_E,P(r) and small cells Dn, one can assign the points r„ to the centers of the cells. Thus, we have arrived at a system of linear equations with respect to the integral currents within each cell, which, using matrix notation, can be written in the form ctGdId rb lD (21) where ct is a (3N x 3N) diagonal matrix of anomalous conduc­tivities, Id and ijp are the vectors of the total and background electric field intensities formed by the x, y, and z components of these fields at the centers of the cells Dn of the anomalous domain D *-D III. II.... If. Tl. Tl.... If. Tl. It... I;V1 T L x > x > x > y > y > y ' fb _ f jb,l jb.2 rb.N rb.l rb. 2 rb.N rb.l rb. 2 r&.jV] ^ D L X ■ X * * * * t / ?/ / ?/ ** * * ■*?/ / ?; / ?; ** * * ^ *-x -*y These vectors have the order 3N. Note that the background IEC can be found by the corresponding numerical integration according to (18). Authorized licensed use limited to: The University of Utah. Downloaded on May 17,2010 at 22:44:50 UTC from IEEE Xplore. Restrictions apply. ZHDANOV et al.: IF.C METHOD IN 3-D F.M MODELING FOR LARGE CONDUCTIVITY CONTRAST 1285 The 3Ar x 3Ar matrix Gj-> is formed by the volume integrals over the elementary cells D.„ of the components of the corre­sponding electric Green's tensor G /?, acting inside domain D. Due to the reciprocity principle [25], elements of this matrix, G'J™;, can be written as iOpn Uaf3 ‘D'n G'Rrt73(r|r)))ffc G/?,3a'(r»|r)cfo BI where B MiBMJ M l{!n M, diag and where = G'^, a, ,3 = x, y, z; p, n = 1,2,... jV where Gp^a are the elements of the corresponding matrix (6) of the conventional integral (3) for the electric field. In other words, matrix Gp is a transposed matrix of the original linear system (4), Gp, for the vector of electric field ep (which justifies the notation we use for Gp). Thus, forward EM modeling based on the IE method is reduced to the solution of the matrix (21) for the unknown vector If) of IEC components inside domain D. The equation is a 3Ar x 3Ar linear system (22) (23) and 1 is identity tensor. We have reduced the EM forward modeling problem to the solution of the matrix (22) with respect to the IEC. The use of the integral current instead of the conventional current density constitutes the key new idea of our method. The integral currents can be found accurately for a body with any conduc­tivity. We do not use anymore the traditional for IE method requirements that the electric field (or electric current) varies slowly inside the cell, because we deal with the integral of this field (or of this current) over the entire cell. As a result, the method can be applied, in principle, for models with arbitrary conductivity contrast. Matrix B is a 3Ar x 3 Ar dense matrix. We use the contraction IE method to precondition matrix (22) [13], [27] Fig. 1. Panel (a) shows a model of a prismatic conductive body with a resistivity of 0.1 £2 - m embedded in the two-layered background. Panel (b) presents the vertical distribution of the electric field within a conductive prism computed using three different discretizations in the vertical direction: 5 cells (stars), 15 cells (circles), and 25 cells (crosses). with M2 defined as the diagonal matrix M2 = diag (u^1, <G1 ' • • • 'a iv ' ail ' a21 ' • • • 2(J5j i 2 \/ n'h < -l -l -l -i\ tt.y , ttj , , . . . , tt.y J 1,2,..., jV. (26) (24) where Mi is the 3Ar x 3Ar diagonal matrix of the square root of the background conductivity, similar to matrix (5) M. 2 (25) One can use different types of iterative methods for the solution of this problem. Detailed analysis of the different solvers is given in [13]. After determining the integral current, we can find the components of the EM field by substituting this current in (14) and (15). IV. Numerical Analysis of thl Ellctric Currlnt Distribution Insidl thl: Conductive: Body Consider a model of a prismatic conductive body with a resistivity of 0.1 fi • m embedded in a two-layered background (Fig. 1). The resistivity contrast between the second layer and the prism is 104. The incident field is a vertically propagated plane EM wave at the 25-Hz frequency, containing both the TM and TE modes. We investigated the effect of different vertical discretizations of the prismatic body on the electric current calculations. Three different discretizations were used in the vertical direction: 5, 15, and 25 cells. The discretizations in Authorized licensed use limited to: The University of Utah. Downloaded on May 17,2010 at 22:44:50 UTC from IEEE Xplore. Restrictions apply. 1286 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 45, NO. 5, MAY 2007 Fig. 2. y component of the electric field Ey obtained using three discretiza­tions in the vertical direction, 5, 15, and 25 cells, respectively. eE.u = 171(25) _ ZT'(^) ^y_____hy 771(25) h^5) - i45) tt{25) fix Fig. 4. Model of two conductive bodies embedded within different layers of two-layered horizontally homogeneous background medium. Fig. 1, panel (b), presents the vertical distribution of the elec­tric field within the conductive prism computed using different vertical disretizations. To produce these plots, we selected a central vertical column of the cells within the prism for each discretization. The electric field, E(rn), was calculated in the center of each elementary cell from this column according to the following formula based on expressions (13) and (12): E(r„) = I„/(A anDn) (27) Fig. 3. Two left panels show the real part (top) and imaginary part (bottom) of the x component of the magnetic field Hx obtained using three discretizations in the vertical direction, 5, 15, and 25 cells, respectively. Two right panels show the relative errors in real (top) and imaginary (bottom) parts of the component. the x and y directions remained the same: 10 and 20 cells, respectively. The horizontal components of the TE mode elec­tric and magnetic fields obtained using all three discretizations are shown in the left panels of Figs. 2 and 3, respectively. The right panels in these figures present the relative errors, eev and ehx > in the real (top) and imaginary (bottom) parts of the corresponding components computed as the difference between the field for the finest discretization (25) and the field, obtained for the coarsest discretization (5), normalized by the field at the finest discretization (25) One can see that these errors do not exceed 1.5%. using the IEC, In, computed for this cell with the IEC method (where Dn and Aan are the volume and the anomalous conductivity of the corresponding elementary cell, respec­tively). Fig. 1, panel (b), shows that the electric field computed for the finest discretization (25 cells in the vertical direction) describes well the skin effect within the conductive body, while the field on the coursed discretization of five vertical cells is practically insensitive to the skin effect. At the same time, the difference between the observed EM field components at the surface is within just 1.5% (Figs. 2 and 3). This remarkable property of the IEC solution is related to the main principle of the IEC method, which is based on computing the IEC, In, within every cell. In this case, the electric field computed according to formula (27) should also describe the averaged electric field within the cell, which corresponds well to the plots shown in Fig. 1, panel b. One can see that the plots of the horizontal electric field components for the coarser discretization describe the average values of the same plot for the finer discretization. The plots of the vertical component of the electric field behaves a little bit differently, because the vertical field is 103 times smaller than the horizontal fields. Nevertheless, the plots for 15-cell and 25-cell discretizations practically coincide, which is a clear manifestation that we reached the optimal level of discretization at 25 cells in the vertical direction. The solution will not change if we will use y Authorized licensed use limited to: The University of Utah. Downloaded on May 17,2010 at 22:44:50 UTC from IEEE Xplore. Restrictions apply. ZHDANOV et al: IEC METHOD IN 3-D EM MODELING FOR LARGE CONDUCTIVITY CONTRAST 1287 Fig. 5. Plots of the real and imaginary parts of the magnetic field Hx, computed using the conventional IE method and the IEC method. The crosses show the data obtained by the IE, while the solid lines present the results of numerical modeling with the IEC method with the cell size of 10 x 10 x 10 m3 for the following resistivity contrasts, c = pn/pi, between the homogeneous background layer and the conductive upper body: (a) c = 102, (b) c = 104, (c) c = 105, and (d) c = 106. Circles show the results obtained using the conventional IE method for fine discretization grid with the cell size of 10 x 10 x 1 m3. the finer discretization. Thus, another important property of the new IEC method is that it does not require a very fine discretization to produce an accurate result, because it does not operate with the discretized electric field but with the integral currents, instead. V. Comparison Between the IEC and the Conventional IE Method In this section, we compare the conventional IE method with the new IEC algorithm. We consider a model shown in Fig. 4. The model consists of two conductive bodies embedded within different layers of two-layered horizontally homogeneous back­ground with the resistivity of the first and second layers equal to 100 and 1000 Q • m, respectively, and with a thickness of the first layer equal to 100 m. In our numerical experiments, the resistivity of the lower body stays constant at 0.1 Q • m, while the resistivity of the upper body changes: p\ = 1, 0.01, 0.001, and 0.0001 Q • m. The incident field is an E-polarized (TE mode) vertically propagated plane EM wave at a frequency of 25 Hz. We use the same discretization of the conductive bodies for both the conventional IE and the IEC methods with the cell size of 10 x 10 x 10 m3. Fig. 5 presents the real and imaginary parts of the mag­netic field component Hx computed using the conventional IE method (crosses) and the IEC method (solid lines), for the resistivity contrast between the homogeneous background layer and the conductive upper body, c = pn/pi, equal to 102, 104, 105, and, 106, respectively. One can see that for this model, the two methods generate practically the same result for low- resistivity contrast. However, they produce different results with the increase of the resistivity contrast, as one would expect. We have repeated the calculations using the conventional IE method with more fine discretization with the cell size of 10 x 10 x 1 m3. These results are shown by circles in panel (c) of Fig. 5 for the highest conductivity contrast c = 106, where two methods have diverged. One can see that in this case the IE method produces the result which is closer to one generated by the IEC method. Note that the limitations of the computer mem­ory did not allow us to run the modeling for the conductivity contrast c = 106 for the smaller cell's size than 10 x 10 x 1 m3 using the conventional IE method. At the same time, for the conductivity contrasts up to c = 105, the results obtained by the conventional IE method with the fine discretization and the results of IEC modeling on the relatively coarse grid are practically the same. This example shows that, the conventional method requires more cells than the IEC method to get the same accuracy. VI. Comparison Between the IEC Method and Analytical Solution for a Conductive Sphere In order to check the accuracy of the new IEC method, we apply this technique to model a response of the conductive sphere excited by the vertically propagated plane EM wave. This problem represents one of a few EM problems which allow for an analytical solution. The mathematical solution of this problem has been described in several publications (see, for example, [3], [4], [15], and [20]). This problem is usually solved by means of the Debye potentials. We compare this analytical solution with numerical modeling using the IEC method. Authorized licensed use limited to: The University of Utah. Downloaded on May 17,2010 at 22:44:50 UTC from IEEE Xplore. Restrictions apply. 1288 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 45, NO. 5, MAY 2007 Fig. 6. (a) Model of a conductive sphere with a radius of 50 m embedded in the homogeneous whole space with a background resistivity of 1000 Q • m; and (b) approximation of the sphere with a model formed by cubic cells with a side of 6.25 m. Fig. 6(a) shows a model of the conductive sphere with a radius of 50 m embedded in the homogeneous whole space with a background (normal) resistivity of pn = 1000 ft • m. In the model study, we use different resistivities of the sphere: Pd = 100,10,1, 0.1, and 0.01 Q • m. The incident field is an E-polarized (TE mode) vertically propagated plane EM wave at a frequency of 25 Hz. The origin of the Cartesian coordinate system is located in the center of the sphere. The receiver profile runs from -410 to 410 m in the x direction at an elevation of 350 m above the center of the sphere. The receivers are located every 20 m. To calculate the sphere response by the IEC method, we approximated the sphere with a model formed by cells with a side of 6.25 m [see Fig. 6(b)]. Using both the analytical solution and the IEC method, we computed an apparent magnetotelluric resistivity for a sphere model according to the formula o \HXJ Note that, according to the method of Debye potentials [4], the EM field components are represented in the form of series. Therefore, the result may depend on the number of the terms kept in these series in calculations. However, these series con­verge extremely fast. Fig. 7 represents the plot of the maximum value of the apparent resistivity versus the number of terms used in the series in analytical calculations for the model with maximum conductivity contrast (le + 5). One can see that the result practically does not change after adding the third term. Fig. 8 shows the plots of the real and imaginary parts of the apparent resistivity, Repa and Impa, for the different resis­tivity contrasts between the homogeneous background and the conductive sphere, c = pn/pd, equal to 10, 102, 103, 104, and 105, respectively. The solid lines correspond to the analytical solution, while the dashed lines present the numerical IEC results. One can see that the difference between the analytical and the numerical IEC solutions does not exceed 0.15% at the extremum value of the apparent resistivity for the highest conductivity contrast of 105. This result demonstrates that the developed new method of integral current equations produces 1000.3811 1000.3811 Number of terms Fig. 7. Plot of the maximum value of the apparent resistivity versus the number of terms used in the series in analytical calculations for the model with maximum conductivity contrast (le + 5 ). an accurate result even for the models with high-conductivity contrast. VII. Conclusion For a long time the main limitation of the IE method was modeling the EM field for models with high-conductivity contrast. In this paper, we have developed a new approach to the construction of the IE method. It is based on using IECs, calculated over the elementary cells of the discretization grid, instead of the electric field itself within the cells, as is commonly used in the conventional IE method. As a result, the method is capable of modeling the EM response in geoelectrical structures with high contrast of conductivity. The method was carefully tested. We compared the numer­ical modeling results with the exact analytical solution for a model of a conductive sphere. Future work will be directed to application of the new method for examining the complex mod­els of geological targets with the large conductivity contrast, typical for mineral exploration. Authorized licensed use limited to: The University of Utah. Downloaded on May 17,2010 at 22:44:50 UTC from IEEE Xplore. Restrictions apply. ZHDANOV et al: IEC METHOD IN 3-D EM MODELING FOR LARGE CONDUCTIVITY CONTRAST 1289 Real pari of apparent resistivity 710­708 Imaginary part of apparent resistivity C=1e+1 E ; •v ; ; ^ ; 1 : C=1e+1 - IEC O 704 --- IEC -- analytical --- analytical -300 -200 -100 (b) 0 X,m -500 -400 -300 -200 -100 0 X,m 200 300 400 500 (c) X,m Fig. 8. Plots of the real and imaginary parts of tlie apparent resistivity K.cpa and lmpa computed using the analytical solution and the IEC method. The solid lines show the data obtained by tlie analytical solution, while the dashed lines present tlie results of numerical modeling with the IEC method for the following resistivity contrasts c = pn/Pd between the homogeneous background and the conductive sphere: (a) c = 10, (b) c = 102, (c) c = 103, (d) c = 104, and (e) e = l'o5. Acknowledgment The authors would like to thank S. K. Lee for developing a Matlab code for an analytical solution for a conductive sphere. References [1J A. Abubakar and P. M. van der Berg, "Iterative forward and inverse algo­rithms based on domain integral equations for three-dimensional electric and magnetic objects," J. Comput. Pins., vol. 195, no. 1, pp. 236-262, Mar. 2004. [2] D. B. Avdeev, A. V. Kuvshinov, O. V. Pankratov, and G. A. Newman, "High-performance three-dimensional electromagnetic modeling using modified Neumann series. Wide-band numerical solution and examples," J. Geomagn. Geoelectr., vol. 49, no. 11/12, pp. 1519-1539, 1997. [3] C. A. Balanis, Advanced Engineering Electromagnetics. Hoboken, NJ: Wiley, 1989. [4] M. N. Berdichevsky and M. S. Zhdanov, Advanced Theory of Deep Geomagnetic Sounding. Amsterdam, The Netherlands: Elsevier, 1984. p. 408.' [5] W. C. Chew, Waves and Fields in Inhomogeneous Media. New York: Van Nostrand, 1990. p. 608. [6] T. J. Cui and W. C. Chew, "Fast algorithm for electromagnetic scatter­ing bv buried 3-D dielectric objects of large size," IEEE Trans. Geosci. Remote Sens., vol. 37, pt. 2, no.5, pp. 2597-2608, Sep. 1999. [7J T. J. Cui, W. C. Chew, A. A. Aydiner, and Y. H. Zhang, "Fast-forward solvers for tlie low-frequency detection of buried dielectric objects," IEEE Trans. Geosci. Remote Sens., vol. 41, pt. 1, no. 9, pp. 2026-2036, Sep. 2003. [8] V. I. Dmitriev, "Electromagnetic fields in inhomogeneous media," in Proc. Comput. Center, 1969, pp. 54-67, (in Russian). Authorized licensed use limited to: The University of Utah. Downloaded on May 17,2010 at 22:44:50 UTC from IEEE Xplore. Restrictions apply. 1290 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING. VOL. 45. NO. 5. MAY 2007 [91 V. I. Dmitriev and N. I. Nesmeyanova, "Integral equation method in three­dimensional problems of low-frequency electrodynamics," in Computat. Math. Modeling, vol. 3. New York: Plenum, 1992, pp. 313-317. [10| C. G. Harquharson and D. W. Oldenburg, "An integral equation so­lution to the geophysical electromagnetic forward modeling prob­lem," in Three-Dimensional Electromagnetics, M. S. Zhdanov and P. H. Wannamaker, iids. Amsterdam, The Netherlands: Blsevier, 2002, pp. 3-19. |111 C. G. Harquharson, K. Duckworth, and D. W. Oldenburg, Comparison of Integral Equation and Physical Scale Modelling of the Electromag­netic Response of Models with Large Conductivity Contrasts. Banff, AB: Scie. Meeting Can. Geophys. Union, 2003. [12| G. W. Hohmann, "Three-dimensional induced polarization and HM modeling," Geophysics, vol. 40, no. 2, pp. 309-324, 1975. [13) G. Hursan and M. S. Zhdanov, "Contraction integral equation method in three-dimensional electromagnetic modeling," Radio Sci., vol. 37, no. 6, p. 1089, 2002. DO I: 10.1029/2001RS002513. |14| Q. H. Liu, Z. Q. Zhang, and X. M. Xu, "The hybrid extended born approx­imation and CG-!!!!T method for electromagnetic induction problems," IEEE Trans. Geosci. Remote Sens., vol. 39, no. 2, pp. 347-355, !!eb. 2001. |15| H. W. March, "The field of a magnetic dipole in the presence of a con­ducting sphere," Geophysics, vol. 18, no. 3, pp. 671-684, Jul. 1953. |16| G. A. Newman and G. W. Hohmann, "Transient electromagnetic re­sponses of high-contrast prism in a layered earth," Geophysics, vol. 53, no. 5, pp. 691-706, May 1988. |17| A. P. Raiche, "An integral equation approach to three-dimensional mod­elling," Geophys. J. R. Astron. Soc., vol. 36, no. 2, pp. 363-376, 1974. |18| L. A. Tabarovsky, Application of Integral Equation Method to Geoelec­trical Problems. Novosibirsk, Russia: Nauka, 1975. (in Russian). |19| P. H. Wannamaker, "Advances in three-dimensional magnetotelluric modeling using integral equations," Geophysics, vol. 56, no. 11, pp. 1716-1728, Nov. 1991. |20| S. H. Ward and G. W. Hohmann, "Hlectromagnetic theory for geophysical applications," in Electromagnetic Methods in Applied Geophysics, vol. 1, M. N. Nabighian, Hd. Tulsa, OK: Soc. Hxploration Geophysicists, 1987. [21] J. T. Weaver, Mathematical Methods for Geo-Electromagnetic Induction. Taunton, U.K.: Res. Studies Press, 1994. p. 316. [221 P. Weidelt, "HM induction in three-dimensional structures," J. Geophys., vol. 41, no. 1,pp. 85-109, 1975. [231 Z. Xiong, "Hlectromagnetic modelling of three-dimensional structures by the method of system iterations using integral equations," Geophysics, vol. 57, no. 12, pp. 1556-1561, Dec. 1992. [24| Z. Xiong and A. Kirsch, "Three-dimensional earth conductivity inver­sion," J. Comput. Appl. Math., vol. 42, no. 1, pp. 109-121, Sep. 1992. [25[ M. S. Zhdanov, Integral Transforms in Geophysics. Berlin, Germany: Springer-Verlag, 1988. p. 367. [26| M. S. Zhdanov and S. Hang, "Quasi linear series in 3D HM modeling," Radio Sci., vol. 32, no. 6, pp. 2167-2188, 1997. [271 M. S. Zhdanov, Geophysical Inverse Theory and Regularization Prob­lems. Amsterdam, The Netherlands: Hlsevier, 2002. p. 628. Vladimir I. Dmitriev was bom in Baku, Azerbaijan, on October 20, 1932. He received the B.Sc. degree in physics from Moscow State University (MSU), Moscow, Russia, and the Ph.D. degree in physics and mathematics, in 1955 and 1967, respectively. He is a Director of the Mathematical Physics laboratory and a Deputy Dean of the Department of Computational Mathematics and Cybernetics at MSU. He has made considerable contributions in the development of the theory of electromagnetic (HM) methods applied to exploration geophysics and geoelectrics. He participated in the foundation of the methods used for solving inverse problems in HM geophysics and magnetotelluric data interpretation, and played a leading role in developing numerical modeling and solution methods for forward and inverse problems in electrodynamics. He has authored more than 300 scientific papers, including 15 monographs. Dr. Dmitriev was honored with the following awards: USSR State Award (1986), MSU Lomonosov Award (1983), and USSR Minister's Cabinet Award (1986). He is a member of the Russian Academy of Natural Sciences. Alexander V. Gribenko received the M.S. degree in geophysics from Novosibirsk State University, Novosibirsk, Russia, in 1998, the M.S. degree in geosciences from Texas Tech University, Lubbock, in 2001, and the Ph.D. degree in geophysics from University of Utah, Salt Lake City, in 2005. He previously worked as a Research Associate with the Institute of Geophysics of the Siberian Branch of the Russian Academy of Sciences from 1996 to 1999. He worked as a Research and Teaching Assistant with Texas Tech University from 1999 to 2001, and with the University of Utah from 2001 to 2005. His current research interests include 3-D modeling and inversion of marine controlled source electromagnetics (CSHM) and induction logging data. Michael S. Zhdanov received the Ph.D. degree in physics and mathematics from Moscow State Uni­versity, Moscow, Russia, in 1970. He spent more than 20 years as a Professor with the Moscow Academy of Oil and Gas and as Head of the Department of Deep Hlectromagnetic Study, and Deputy Director of IZMIRAN, and later as Director of the Geoelectromagnetic Research In­stitute, Russian Academy of Sciences in Moscow. In 1990, he received an honorary diploma of the Gauss Professorship from the Gottingen Academy of Sciences, Germany. In 1991, He was elected a full member of the Russian Academy of Natural Sciences. He became an Honorary Professor with the China National Center of Geological Hxploration Technology in 1997 and a Hellow of the Hlectromagnetics Academy, in 2002. He joined the University of Utah as Hull Professor in 1993 and has been Director of the Consortium for Hlectromagnetic Modeling and Inversion (CHMI) since 1995. Authorized licensed use limited to: The University of Utah. Downloaded on May 17,2010 at 22:44:50 UTC from IEEE Xplore. Restrictions apply.