Three-Component Forward Modeling for Transient Electromagnetic Method
In general, the time derivative of vertical magnetic field is considered only in the data interpretation of transient electromagnetic (TEM) method. However, to survey in the complex geology structures, this conventional technique has begun gradually to be unsatisfied with the demand of field exploration. To improve the integrated interpretation precision of TEM, it is necessary to study the three-component forward modeling and inversion. In this paper, a three-component forward algorithm for 2.5D TEM based on the independent electric and magnetic field has been developed. The main advantage of the new scheme is that it can reduce the size of the global system matrix to the utmost extent, that is to say, the present is only one fourth of the conventional algorithm. In order to illustrate the feasibility and usefulness of the present algorithm, several typical geoelectric models of the TEM responses produced by loop sources at air-earth interface are presented. The results of the numerical experiments show that the computation speed of the present scheme is increased obviously and three-component interpretation can get the most out of the collected data, from which we can easily analyze or interpret the space characteristic of the abnormity object more comprehensively.
The transient electromagnetic method has shown great potential in hydrological and hazardous waste site characterization [1, 2], mineral exploration , and general geological mapping, and geophysical reconnaissance. However, the behavior of TEM fields is not yet fully understood . The need for further theoretical insight is reflected by the increasing demands placed on transient electromagnetic methods for petroleum, mineral, and geothermal exploration. Forward modeling is one of the most common and effective methods that help us understand the physical significance of the electromagnetic responses [5, 6]. Computer solutions for this method have been mainly confined to the vertical component of magnetic field time derivative. Until now, a limited number of solutions have appeared in the literatures which are relevant to the TEM three-component responses of a 3D source over 2D earth, which is the so-called 2.5D.
At present the 2.5D model represents the only way of interpreting controlled source electromagnetic data in terms of a complex earth, due to the prohibitive amount of computer time and storage required for a complex 3D model . The first published theoretical finite element derivation for the 2.5D electromagnetic problem was by Coggon . Stoyer and Greenfield  calculated the frequency-domain response of a 2D earth to a vertical magnetic dipole source with the finite difference method. Lee  and Lee and Morrison  use FEM to obtain the fields induced by a magnetic dipole source. Leppin  presented an FD numerical scheme by which 2.5D TEM scattering problems can be solved directly in the time domain. Everett and Edwards  and Unsworth et al.  evaluated the time-domain and frequency-domain responses, respectively. Mitsuhata  described a CSEM modeling method to obtain the responses for 2D models including topography. The new development of using FEM in EM applications can be found in the works of Key and Weiss , Li and Key , Li and Constable , Abubakar et al. , and Kong et al. . Unfortunately, because of the complexity of the problem itself, 2.5D forward modeling for time-domain electromagnetic method has not yet been properly resolved and is still one of the most difficult problems in the field of computational geophysics [21–23]. Moreover, most of the forward schemes for controlled-source EM methods have been carried out by solving the boundary value problem of coupled electromagnetic fields [13–15, 24–26].
In this paper, a three-component forward algorithm for 2.5D transient electromagnetic method based on the independent electric and magnetic field has been developed. By dint of the operation rule of curl and divergence, the present modeling method transforms the governing Maxwell equations into Laplace domain and then expands them in component format. Applying the Fourier transform to the electric and magnetic components in the strike direction, provided conductivity is invariant in this direction, I can get 2.5D TEM boundary value problem in the Laplace and along-strike spatial Fourier domain [27, 28]. Compared with the conventional algorithm, the main advantage of the new scheme is that it can reduce the size of the global system matrix to the utmost extent. For example, if the total number of nodes created by finite element method is solving by traditional method, the size of the global system matrix should be while by the present algorithm, it will be just that is to say, the present is only one fourth of the former. Noting that the left sides of electric and magnetic equations are the same form and only the right sides are different, when the simple Dirichlet condition is assigned to the boundary , forward modeling for many source locations will not increase the computing burden of solving matrix evidently. In order to illustrate the feasibility and usefulness of the present algorithm, several typical geoelectric models of the transient electromagnetic responses produced by loop sources at air-earth interface are presented.
2.1. Time-Domain EM Equations
A Cartesian coordinate system is defined with vertically downwards and the direction is set to be the strike direction. The total electric field and magnetic fields generated by a specified magnetic source satisfy the Maxwell equations where is the electrical conductivity which is assumed to vary only in plane, is the magnetic permeability of free space, is the dielectric permittivity, and are the electric and magnetic field vectors, respectively, and is the impressed magnetic source. The magnetic dipole source placed at the origin in the Cartesian coordinates is described at the surface as follows: where is the magnetic current, is the electric current, and is sampling time.
2.2. Laplace and Fourier Domain-Independent EM Equations
I transform the governing Maxwell (1) into Laplace domain and express them in component form. Applying the Fourier transform to the electric and magnetic components in the strike direction, we obtain 2.5D total electromagnetic fields in the Laplace and along-strike spatial Fourier domain [27, 28]. It is clear that the primary fields can be calculated for a simple, one-dimensional conductivity structure and the (1) also govern the primary fields, thus when and when subtraction equations obtained for the secondary fields are assuming that the media in the vicinity of the source is homogenous and omitting some irrelevant items, the other components of electric filed and magnetic field not along the strike could be written as where and is the Laplace variable, generally a complex number, is the along-strike wave number in Fourier domain, is the difference between the total 2D conductivity and the 1D background conductivity used to calculate the primary fields, , , , , and are the secondary electromagnetic field components, respectively. , , , , and denote the components of primary electromagnetic field in turn, complex unit is denoted by .
3. Finite Element Analysis
3.1. Domain Discretization
In this paper, the rectangle element subdivided into four sympeda triangle elements (Figure 1(a)) has been adopted, which does not increase the total number of FEM nodes and is also feasible for the meshing for inclined physical interface.
3.2. Element Interpolation
For a linear triangular element , there are three nodes located at the vertices of the triangle (Figure 1(b)). Assume that the nodes are numbered counterclockwise by number 1, 2, and 3, with the corresponding values of unknown function denoted by , and , respectively. In each element, , , and are constant, and the interpolation or expansion functions are descried with the global coordinates and [29–31]. Each component of electromagnetic field is represented as where are the linear interpolation function given by in which is the area of the triangle element , , , and are constant coefficients to be determined.
3.3. Formulation via the Galerkin Method
To be general, I use (4) as an example to formulate the system of equation using the Galerkin method. The residual associated with (4) is and thus the weighted residual for element is Substituting (11) into (12) yields From Green’s theorem , and (13) can be written as where denotes the contour enclosing , is the outward unit vector normal to , and and are the direction cosines of the angle between the outward normal vector and respective axes and Replacing by and employing (9) to (15), I derive the element equation which can be written in matrix form as where
3.4. Assembly to Form the System of Equations
Considering that node 5 (Figure 1(a)) has no relation to all other nodes but to the four corner points of the rectangle element, we can eliminate it in the system of linear algebraic equations before solving the equations system, which will reduce the rank of the overall matrix and decrease the computational complexity [29, 30].
According to (19), the matrices of the four triangle elements and can be derived, specifically , , and . The summation of them is the element matrix of the rectangular element similarly, therefore the element equation of the rectangle element is
The system of equations is then obtained by expanding (23) and then assembling it for all elements, giving which can also be written as where is the total number of the FEM elements. This global system matrix does not include the central points of each rectangle element and the rank of it is just the total number of nodes of rectangular meshes. Considering the truncated boundary condition and in order to avoid undesirable effects from the boundaries, the boundaries are placed far away from the area of interest, so that the secondary field can be ignored, that is to say, at the same time, to eliminate the reflection from the boundaries, the boundaries are expanded by increasing the node-spacing gradually [5, 6, 15]. Therefore, (25) can be simplified as
3.5. Solution to the Linear Equations System
Note that the equations of and are of the same form, so we obtain from(22) where is the secondary electric or magnetic field component to be solved, namely or , and and are vectors resulting from the primary fields and at all nodes, respectively. The resulting linear system of equations in (23) is then solved using a multifrontal LU decomposition method developed by Davis and Duff  that is solved for all source excitations simultaneously, in which the solutions for many source locations can be achieved by factorizing the stiffness matrix only once. For a problem with total nodes, the optimal multifrontal LU method takes operations for matrix factorization and operations for any additional source calculation [19, 33]. However, because the delaying time and Laplace variable information are contained in the stiffness matrix , each new delaying time or Laplace variable requires an additional stiffness matrix factorization. Thus, the total run time is dependent on the number of cells in the model times the number of delaying sampling times and Laplace variables.
3.6. Electromagnetic Fields in Domain
Once the EM components and in Laplace and along-strike spatial Fourier domains are computed numerically, the other secondary field components can be evaluated from (6). These complementary components are calculated from the space derivatives of and which are evaluated from the space derivatives of the interpolation function and in each element.
The EM fields in the Laplace and Fourier domains are transformed into the space and time domains by discrete inverse Fourier and inverse Laplace transformations, respectively. In this study, the Gaver-Stehfest algorithm  is used to approximately invert the Laplace transform. Since the fields at are taken into consideration, they are evaluated by where are the coefficients, is an even integer whose optimal value depends upon the computer word length. Gaver-Stehfest algorithm requires only the real values of the Laplace variable and eliminates the need for complex arithmetic anywhere in the algorithm. The integration in (30) is carried out by cubic spline interpolation in the logarithmic wave number domain . Finally, the primary fields of the central loop source are evaluated directly from Fourier sine or cosine transform expressions and added to the secondary fields from (30) to yield the total electromagnetic fields.
4.1. Layered Earth Model
To test the algorithm, the three-layered earth model computed consists of layer 1 resistivity = 100 Ω·m, layer 2 resistivity = 10 Ω·m, layer 3 resistivity = 1000 Ω·m, and layer 1 and layer 2 thickness = 100 m, and 50 m, respectively. The source is a central loop of dimensions 50 0 m, the equivalent area of the receiving loop is 43 m2 and the transmitted current is 1.0 A.
As a check for our forward algorithm, let us compare the finite-element numerical solution to the analytical solution for the vertical electromotive force (EMF) of a 3-layered earth model shown in Figure 2 at 17 sampling times, ranging from 0.001 to 10 ms. The two solutions demonstrate an excellent agreement with each other. The FEM solution shows a slightly smaller decay rate at early times (before 0.01 ms) than does the analytical solution.
4.2. Horizontal Tabular Model
As illustrated in Figure 3, the 10 Ω·m body is 320 m wide, 10 m in depth extent, and is 50 m deep. It is embedded in a 1000 Ω·m homogeneous half-space and excited by a 100 00 m square loop source on the surface. The equivalent area of the receiving loop is 43 m2 and the transmitted current is 1.0 A.
The transient electromagnetic responses by the finite element method along the positive -axis are shown in Figure 4. The horizontal and vertical electromotive forces all correspond to the magnetic field measured with horizontal coils. Results are plotted as a function of the location of sounding and are shown for 17 values of , the number of delaying sampling times. It is obvious that the profile of horizontal EMF of component showed in Figure 4(a) has two main abnormal peaks, corresponding to the left and right edges of the horizontal tabular orebody. The profile of component displayed in Figure 4(b) reflects the information for long direction extension of the orebody along strike direction. Figure 4(c) shows the profile of vertical components and is the most reliable information for TEM data interpretation at present. At the early times, for the coupling of the overburden layer and loops that is much stronger than that of the orebody and loops, the signals observed mainly reflect the effect of the overburden layer. With the delay of the sampling time, the shape of curves is going to release the information of the deposit gradually. However, at the late stage the definition of the abnormality is decreasing. Compared to the oblique tabular orebody, the horizontal one can be coupled with the loops, mostly, so that the amplitude of the abnormality is distinct and usually has the shape of one peak with flat top, which is consistent with the results in physical simulation .
4.3. Compound Conductive Block Model
In this model, two 10 Ω·m orebodies are located in a 1000 Ω·m uniform half-space (Figure 5). The two bodies are 10 m thick, 100 m in depth extent, 170 m away from each other, and are buried at a depth of 20 m. A square loop 200 m on a side is laid on the surface, the equivalent area of the receiving loop is 10000 m2 and the transmitted current of the source is 10.0 A.
The horizontal and vertical components of at 17 delay times along the positive -axis are shown in Figure 6. Profiles of component , as shown in Figure 6(a), have two zero-value points just locating at the same position as the two projection points on the earth surface which are projected from the vertical central axes of the two vertical orebodies, that is to say, the component reflects the horizontal position of the two vertical slabs exactly. In Figure 6(b), the extension of the two vertical tabular orebodies along strike direction is showed clearly, while Figure 6(c) represents the profiles of component in TEM responses. It can be concluded from the above that the distance of the two vertical tabular orebodies and the sampling delay time determine the spacial characteristic of the abnormality. In other words, only when the separation distance is large enough can each of them be distinguished by the number of peaks in the abnormality curves. If the distance is too small, the response will be similar as that of a thick tabular orebody. At early stage each body can be discriminated but with the postponement of sampling time, the response will behave like that of a single sheet tabular orebody.
4.4. Inclined Tabular Model
The model is a dip tabular body embedded in a homogeneous half-space with resistivity 1000 Ω·m (Figure 7). The body with dip of 45 degree is 7 m thick, 50 m in vertical depth extent, and is buried at a depth of 10 m. It has a resistivity of 10 Ω·m and is excited by a central loop on the surface, in which the square transmitting loop of dimensions 50 0 m, equivalent area 2500 m2 of the receiving coils, and the current amplitude of the source is 10.0 A.
Figure 8 shows the profile variations of the horizontal and vertical EMF along the -axis for sampling time from 0.001 ms to 10 ms. It can be seen clearly in Figure 8(a) that the curves are unsymmetric according to the origin points and what is more different from those in Figures 4 and 6 is that the sequence of the peak values is inversed. Specifically, the plus peak values come first and then the minus ones. Therefore, we can say that the component can help us to identify the occurrence information . From Figure 8(b), the information of the dipping tabular orebody along strike direction can be obtained. With the postponing of the sampling time, the minimum value of the curve is moving right and up. In Figure 8(c), the responses have the characteristic like that of an equiaxed body at the early and middle times. The curves are different from time to time and the maximum value is moving toward the dipping side with the right hand of the curve becoming gentler. Then the curves are becoming somewhat like having two peaks. As the time goes by, the minimum value is shifting to the dipping side and the double peak is vanishing gradually. At the late stage, the coupling of the circulation current downward and the orebody is getting weaker and weaker so that the curves are showing the same characteristic as early stage.
A three-component finite-element forward modeling algorithm for 2.5D transient electromagnetic method based on the independent electric and magnetic fields has been developed. The solution possesses a number of appealing features. First, compared with the algorithm for the coupling electric and magnetic fields, due to the form of equations controlling the electric and magnetic fields along strike direction are alike, the new scheme presented in this paper can reduce the size of the global system matrix to the utmost extent, thus the present is only one fourth of the former. Once the factorization of the global matrix is performed, solutions for magnetic fields can be obtained swiftly by simple back-substitution, which simplify the whole solving process greatly. Second, interpretation by three-component data can get the most out of the collected data, from which we can analyze or interpret the space characteristic of the abnormity object more comprehensively. Specifically, and components are sensitive to tectonic information. By the peaks of curves, the component can reflect information about the plane position and geological occurrence of the abnormity object, while the component is of great advantage to reveal the variation of geological strata along strike direction.
The author thanks Professor Ruan B. Y. of Guilin University of Technology for sharing his many helpful discussions. This work was supported by the National High Technology Research and Development Program of China Grant 2007AA06Z134, Higher School Doctor Subject Special Scientific Research Foundation Grant 20070533104, and National Natural Science Foundation of China Grant 40974077. He also wishes to thank the editor and the referees for their constructive comments.
L. Pellerin and D. L. Alumbaugh, “Tools for electromagnetic investigation of the shallow subsurface,” Leading Edge, vol. 16, pp. 1631–1640, 1997.View at: Google Scholar
M. N. Nabighian and J. C. Macnae, “Time domain electromagnetic prospecting methods,” in Electromagnetic Methods in Applied Geophysics, M. N. Nabighian, Ed., vol. 2, pp. 427–520, Society for Exploration Geophysics, Tulsa, Okla, USA, 1991.View at: Google Scholar
G. A. Newman, G. W. Hohmann, and W. L. Anderson, “Transient electromagnetic response of a Three-Dimensional body in a layered earth,” Geophysics, vol. 51, pp. 1608–1627, 1986.View at: Google Scholar
J.-S. Shen and W.-B. Sun, “A study on the 2.5-D electromagnetic modeling by the finite element method and the wave number selection,” Computing Techniques for Geophysical and Geochemical Exploration, vol. 30, no. 2, pp. 135–144, 2008.View at: Google Scholar
M. N. Nabighian, Electromagnetic Methods in Applied Geophysics, 1, Theory, Society of Exploration Geophysicists, Tulsa, Okla, USA, 1988.
J. H. Coggon, “Electromagnetic and electrical modeling by the finite element method,” Geophysics, vol. 36, pp. 132–155, 1971.View at: Google Scholar
C. H. Stoyer and R. J. Greenfield, “Numerical solutions of the response of a two-dimensional earth to an oscillating magnetic dipole source,” Geophysics, vol. 41, pp. 519–530, 1976.View at: Google Scholar
K. H. Lee, Electromagnetic scattering by a two-dimensional inhomogeneity due to an oscillating magnetic dipole source, Ph.D. thesis, University of California, Berkeley, Calif, USA, 1978.
K. H. Lee and H. F. Morrison, “A numerical solution for the eleccromagnetic scattering by a two- dimensional inhomogeneity,” Geophysics, vol. 50, pp. 1163–1165, 1985.View at: Google Scholar
M. Leppin, “Electromagnetic modeling of 3-D sources over 2-D inhomogeneities in the time domain,” Geophysics, vol. 57, pp. 994–1003, 1992.View at: Google Scholar
M. E. Everett and R. N. Edwards, “Transient marine electromagnetics: the 2.5-D forward problem,” Geophysical Journal International, vol. 113, pp. 545–561, 1993.View at: Google Scholar
M. J. Unsworth, B. J. Travis, and A. D. Chave, “Electromagnetic induction by a finite electric dipole source over a 2-D earth,” Geophysics, vol. 58, pp. 198–214, 1993.View at: Google Scholar
Y. Mitsuhata, “2-D electromagnetic modeling by finite-element method with a dipole source and topography,” Geophysics, vol. 65, pp. 465–475, 2000.View at: Google Scholar
G. Q. Xue, X. Li, and Q. Y. Di, “The progress of TEM in theory and application,” Progress in Geophysics, vol. 22, pp. 1195–1200, 2007.View at: Google Scholar
G. Y. Lv, “Current status and development trend of transient EM method,” Computing Techniques for Geophysical and Geochemical Exploration, vol. 29, pp. 111–115, 2007.View at: Google Scholar
G. Q. Xue, X Li, and Q. Y. Di, “Research progress in TEM forward modeling and inversion calculation,” Progress in Geophysics, vol. 23, pp. 1165–1172, 2008.View at: Google Scholar
Y. L. Meng and Y. Z. Luo, “Finite element method for solving electromagnetic field with harmonic dipolar source on two-dimensional geoelectric model,” in Annual of the Chinese Geophysical Society, Seismological Press, Beijing, China, 1994.View at: Google Scholar
H. J. Wang, Methods study on 2.5-D forward and inversion for transient electromagnetic with central loop, Ph.D. thesis, China University of Geosciences, Beijing, China, 2001.
B. Xiong and Y.-Z. Luo, “Finite element modeling of 2.5-D TEM with block homogeneous conductivity,” Chinese Journal of Geophysics, vol. 49, pp. 590–597, 2006.View at: Google Scholar
K. J. Xu, Study on 2.5D complex resistivity electromagnetic forward and inversion, Ph.D. thesis, Jilin University, Jilin, China, 2007.
B. Xiong, “A new formula based on the independent electric and magnetic fields for 2.5-D forward of TEM,” in Proceedings of the 19th International Workshop on Electromagnetic Induction in the Earth, pp. 536–539, 2008.View at: Google Scholar
Y. Z. Luo and G. Q. Zhang, Application of Electronic Computer in Electrical Prospecting, Wuhan College of Geology Press, 1987.
S. Z. Xu, FEM in Geophysics, Science Press, Beijing, China, 1994.
J. M. Jin, J. G. Wang, and D. B. Ge, The Finite Element Method in Electromagnetics, Xidian University Press, Xi'an, China, 1997.
T. A. Davis and I. S. Duff, “An unsymmetric-pattern multifrontal method for sparse LU factorization,” SIAM Journal on Matrix Analysis and Applications, vol. 18, pp. 140–158, 1997.View at: Google Scholar
T. A. Davis, Direct Methods for Sparse Linear Systems, SIAM, Philadelphia, Pa, USA, 2007.
H. Stehfest, “Numerical inversion of Laplace transforms,” Communications of the ACM, vol. 13, pp. 47–49, 1970.View at: Google Scholar
Z. L. Niu, The Theory of Time-Domain Electromagnetic Methods, Central South University of Technology Press, Hunan, China, 1992.
J. T. Liu, H. M. Gu, and X. Y. Hu, “Three-components interpretation for transient electromagnetic method,” Yangtze River, vol. 39, pp. 114–116, 2008.View at: Google Scholar