Three-Component Forward Modeling for Transient Electromagnetic Method

Abstract

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.

1. Introduction

The transient electromagnetic method has shown great potential in hydrological and hazardous waste site characterization [1, 2], mineral exploration [3], and general geological mapping, and geophysical reconnaissance. However, the behavior of TEM fields is not yet fully understood [4]. 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 [7]. The first published theoretical finite element derivation for the 2.5D electromagnetic problem was by Coggon [8]. Stoyer and Greenfield [9] calculated the frequency-domain response of a 2D earth to a vertical magnetic dipole source with the finite difference method. Lee [10] and Lee and Morrison [11] use FEM to obtain the fields induced by a magnetic dipole source. Leppin [12] presented an FD numerical scheme by which 2.5D TEM scattering problems can be solved directly in the time domain. Everett and Edwards [13] and Unsworth et al. [14] evaluated the time-domain and frequency-domain responses, respectively. Mitsuhata [15] 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 [16], Li and Key [17], Li and Constable [18], Abubakar et al. [19], and Kong et al. [20]. 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 2𝑛×2𝑛, 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 [15], 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. Theory

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∇×𝐄(𝑟,𝑡)=−𝜇𝜕𝐇(𝑟,𝑡)𝜕𝑡−𝐉_𝑚,∇×𝐇(𝑟,𝑡)=𝜀𝜕𝐄(𝑟,𝑡)𝜕𝑡+𝜎𝐄(𝑟,𝑡), (1) 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:𝐉_𝑚=𝜇⋅𝜕𝐌/𝜕𝑡⋅𝛿(𝑥)𝛿(𝑦)𝛿(𝑧), (2) where𝐌=𝐦[1−𝐼(𝑡)]𝐼(𝑡)=0𝑡≤0,1𝑡>0. (3)𝐌 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𝜕²𝑒_𝑦,𝑠𝜕𝑥²+𝜕²𝑒_𝑦,𝑠𝜕𝑧²−𝜆²𝑒_𝑦,𝑠=𝑐𝑒_𝑦,𝑝,𝜕²ℎ_𝑦,𝑠𝜕𝑥²+𝜕²ℎ_𝑦,𝑠𝜕𝑧²−𝜆²ℎ_𝑦,𝑠=𝑐ℎ_𝑦,𝑝, (4)(5) 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ℎ_𝑧,𝑠=−𝑖𝑎₂𝜕ℎ_𝑦,𝑠𝜕𝑧−𝑎₁𝜕𝑒_𝑦,𝑠,ℎ𝜕𝑥_𝑥,𝑠=−𝑖𝑎₂𝜕ℎ_𝑦,𝑠𝜕𝑥+𝑎₁𝜕𝑒_𝑦,𝑠,𝑒𝜕𝑧_𝑥,𝑠=𝑎₃𝜕ℎ_𝑦,𝑠𝜕𝑧−𝑖𝑎₂𝜕𝑒_𝑦,𝑠,𝑒𝜕𝑥_𝑧,𝑠=−𝑎₃𝜕ℎ_𝑦,𝑠𝜕𝑥−𝑖𝑎₂𝜕𝑒_𝑦,𝑠,𝑎𝜕𝑧₁=𝜀𝑠+𝜎𝜆²,𝑎₂=𝑘_𝑦𝜆²,𝑎₃=−𝜇𝑠𝜆², (6)(7) where𝜆²=𝑘²_𝑦−𝑘²,𝑘²=−𝜇𝑠(𝜀𝑠+𝜎),𝑐=𝜇𝑠𝜎_𝑠, (8) 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.

(a)

(b)

Figure 1 

Example of the subdivision of a two-dimensional domain.

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Φ^Δ(𝑥,𝑧)=³_𝑗=1𝑁^Δ_𝑗(𝑥,𝑧)Φ^Δ_𝑗, (9) where 𝑁^Δ_𝑗(𝑥,𝑧) are the linear interpolation function given by𝑁^Δ_𝑗1(𝑥,𝑧)=2Δ(𝛼^Δ_𝑗+𝛽^Δ_𝑗𝑥+𝛾^Δ_𝑗𝑧),𝑗=1,2,3, (10) 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𝜕𝑟=²𝑒_𝑦,𝑠𝜕𝑥²+𝜕²𝑒_𝑦,𝑠𝜕𝑧²−𝜆²𝑒_𝑦,𝑠−𝑐𝑒_𝑦,𝑝, (11) and thus the weighted residual for element 𝑒₁ is𝑅^Δ_𝑖=_𝐷Δ𝑁^Δ_𝑖𝑟𝑑𝑥𝑑𝑧. (12) Substituting (11) into (12) yields𝑅^Δ_𝑖=_𝐷Δ𝑁^Δ_𝑖(𝜕²𝑒_𝑦,𝑠𝜕𝑥²+𝜕²𝑒_𝑦,𝑠𝜕𝑧²−𝜆²𝑒_𝑦,𝑠−𝑐𝑒_𝑦,𝑝)𝑑𝑥𝑑𝑧. (13) From Green’s theorem [15],_𝐷𝜙𝜕𝜑𝜕𝑥𝑑𝑥𝑑𝑧=−_𝐷𝜕𝜙𝜕𝑥𝜑𝑑𝑥𝑑𝑧+_𝜕𝐷𝜙𝜑𝑛_𝑥𝑑ℓ, (14) and (13) can be written as𝑅^Δ_𝑖=−_𝐷Δ(𝜕𝑁^Δ_𝑖𝜕𝑥𝜕𝑒_𝑦,𝑠+𝜕𝑥𝜕𝑁^Δ_𝑖𝜕𝑧𝜕𝑒_𝑦,𝑠𝜕𝑧+𝜆²𝑁^Δ_𝑖𝑒_𝑦,𝑠−)𝑑𝑥𝑑𝑧_𝐷Δ𝑐𝑁^Δ_𝑖𝑒_𝑦,𝑝𝑑𝑥𝑑𝑧+_𝜕𝐷Δ𝑁^Δ_𝑖𝐃⋅𝐧^Δ𝑑ℓ, (15) where𝐃⋅𝐧^Δ=𝜕𝑒_𝑦,𝑠𝑛𝜕𝑥^𝑥+𝜕𝑒_𝑦,𝑠𝑛𝜕𝑧^𝑧, (16)𝜕𝐷_Δ 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 𝑅^Δ_𝑖=−³_𝑗=1𝐷Δ(𝜕𝑁^Δ_𝑖𝜕𝑥𝜕𝑁^Δ_𝑗+𝜕𝑥𝜕𝑁^Δ_𝑖𝜕𝑧𝜕𝑁^Δ_𝑗𝜕𝑧+𝜆²𝑁^Δ_𝑖𝑁^Δ_𝑗)Φ^Δ_𝑗−𝑑𝑥𝑑𝑧_𝐷Δ𝑐𝑁^Δ_𝑖𝑒_𝑦,𝑝𝑑𝑥𝑑𝑧+_𝜕𝐷Δ𝑁^Δ_𝑖𝐃⋅𝐧^Δ𝑑ℓ (17) which can be written in matrix form as{𝑅^Δ}=−[𝐾^Δ]{Φ^Δ}−{𝑏^Δ}+{𝑔^Δ}, (18) where𝐾^Δ_𝑖𝑗=_𝐷Δ(𝜕𝑁^Δ_𝑖𝜕𝑥𝜕𝑁^Δ_𝑗+𝜕𝑥𝜕𝑁^Δ_𝑖𝜕𝑧𝜕𝑁^Δ_𝑗𝜕𝑧+𝜆²𝑁^Δ_𝑖𝑁^Δ_𝑗𝑏)𝑑𝑥𝑑𝑧,𝑖,𝑗=1,2,3,_Δ𝑖=_𝐷Δ𝑐𝑁^Δ_𝑖𝑒_𝑦,𝑝𝑔𝑑𝑥𝑑𝑧,_Δ𝑖=_𝜕𝐷Δ𝑁^Δ_𝑖𝐃⋅𝐧^Δ𝑑ℓ. (19)

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𝐊^■=𝐊^Δ₁+𝐊^Δ₂+𝐊^Δ₃+𝐊^Δ₄, (20) similarly,𝐛^■=𝐛^Δ₁+𝐛^Δ₂+𝐛^Δ₃+𝐛^Δ₄,𝐠^■=𝐠^Δ₁+𝐠^Δ₂+𝐠^Δ₃+𝐠^Δ₄, (21)(22) therefore the element equation of the rectangle element is{𝑅^■}=−[𝐾^■]{Φ^■}−{𝑏^■}+{𝑔^■}. (23)

The system of equations is then obtained by expanding (23) and then assembling it for all elements, giving{𝑅}=^𝑁𝑒_𝑒=1{𝑅^■}=^𝑁𝑒_𝑒=1(−[𝐾^■]{Φ^■}−{𝑏^■}+{𝑔^■})={0}, (24) which can also be written as[𝐾]{Φ}=−{𝑏}+{𝑔}, (25) 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,𝑒_𝑦,𝑠=ℎ_𝑦,𝑠=0; (26) 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[𝐾]{Φ}=−{𝑏}. (27)

3.5. Solution to the Linear Equations System

Note that the equations of 𝑒_𝑦,𝑠 and ℎ_𝑦,𝑠 are of the same form, so we obtain from(22)𝐊Φ=−𝐒, (28) where 𝐒=𝐛_𝑒𝐛_ℎ, (29)Φ 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 [32] 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 𝑂(𝑁^1.5) operations for matrix factorization and 𝑂(𝑁log𝑁) 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 [34] is used to approximately invert the Laplace transform. Since the fields at 𝑦=0 are taken into consideration, they are evaluated by 𝐹(𝑥,𝑦=0,𝑧,𝑡)≃ln2𝜋𝑡^∞₀^𝑁𝑠_𝑘=1𝑉_𝑘𝑓(𝑥,𝑘_𝑦,𝑧,𝑠_𝑘)𝑑𝑘_𝑦, (30) 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 𝑠_𝑘=𝑘ln2/𝑡 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. Results

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 m×50 m, the equivalent area of the receiving loop is 43 m² 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.

Figure 2 

Comparison between FEM results with analytical results for the normalized vertical EMF on a 3-layer earth.

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 m×100 m square loop source on the surface. The equivalent area of the receiving loop is 43 m² and the transmitted current is 1.0 A.

Figure 3 

Horizontal tabular orebody model.

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 [35].

(a)

(b)

(c)

Figure 4 

TEM profiles of horizontal tabular orebody for (a) the horizontal EMF in 𝑥 direction, (b) the horizontal EMF in 𝑦 direction, and (c) the vertical EMF.

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 m² and the transmitted current of the source is 10.0 A.

Figure 5 

Compound orebody model of two vertical slabs.

The horizontal and vertical components of 𝜕𝐁/𝜕t 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.

(a)

(b)

(c)

Figure 6 

TEM profiles of compound orebody for (a) the horizontal EMF in 𝑥 direction, (b) the horizontal EMF in 𝑦 direction, and (c) the vertical EMF.

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 m×50 m, equivalent area 2500 m² of the receiving coils, and the current amplitude of the source is 10.0 A.

Figure 7 

Inclined tabular orebody model.

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 [36]. 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)

(b)

(c)

Figure 8 

TEM profiles of inclined tabular orebody for (a) the horizontal EMF in 𝑥 direction, (b) the horizontal EMF in 𝑦 direction, and (c) the vertical EMF.

5. Conclusions

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.

Acknowledgments

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.