Three-Component Forward Modeling for Transient Electromagnetic Method

In general, the time derivative of vertical magnetic ﬁeld 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 unsatisﬁed with the demand of ﬁeld 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 ﬁeld 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.


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 timedomain electromagnetic method has not yet been properly resolved and is still one of the most difficult problems in 2 International Journal of Geophysics the field of computational geophysics [21][22][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][14][15][24][25][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 (y) 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 n, solving by traditional method, the size of the global system matrix should be 2n × 2n, while by the present algorithm, it will be just n × n, 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.

Time-Domain EM Equations. A Cartesian coordinate
system is defined with z vertically downwards and the y 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 x − z plane, μ is the magnetic permeability of free space, ε is the dielectric permittivity, E and H are the electric and magnetic field vectors, respectively, and J m 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 M is the magnetic current, I is the electric current, and t is sampling time.

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 σ 0 (z) and the (1) also govern the primary fields, thus when σ = σ(z) and when subtraction equations obtained for the secondary fields are ∂ 2 e y,s ∂x 2 + ∂ 2 e y,s ∂z 2 − λ 2 e y,s = ce y,p , 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 s is the Laplace variable, generally a complex number, k y is the along-strike wave number in Fourier domain, σ s is the difference between the total 2D conductivity and the 1D background conductivity used to calculate the primary fields, e x,s , e y,s , h x,s , h y,s , and h z,s are the secondary electromagnetic field components, respectively. e x,p , e y,p , h x,p , h y,p , and h z,p denote the components of primary electromagnetic field in turn, complex unit is denoted by i.

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.

Element Interpolation.
For a linear triangular element e 1 , 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 Φ Δ 1 , Φ Δ 2 , and Φ Δ 3 , respectively. In each element, σ, ε, and μ are constant, and the interpolation or expansion functions are descried with the global coordinates x and z [29][30][31]. Each component Φ of electromagnetic field is represented as where N Δ j (x, z) are the linear interpolation function given by in which Δ is the area of the triangle element e 1 , α Δ j , β Δ j , and γ Δ j are constant coefficients to be determined.

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 e 1 is Substituting (11) into (12) yields From Green's theorem [15], and (13) can be written as where ∂D Δ denotes the contour enclosing D Δ , n Δ is the outward unit vector normal to ∂D Δ , and n x and n z are the direction cosines of the angle between the outward normal vector and respective axes x and z. Replacing e y,s by Φ and employing (9) to (15) which can be written in matrix form as where

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 (e 1 , e 2 , e 3 ,and e 4 ) can be derived, specifically K Δ 1 , K Δ 2 , K Δ 3 , and K Δ 4 . 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 N e 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, e y,s = h y,s = 0; 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

Solution to the Linear Equations System.
Note that the equations of e y,s and h y,s are of the same form, so we obtain from(22) where Φ is the secondary electric or magnetic field component to be solved, namely e y,s or h y,s , and b e and b h are vectors resulting from the primary fields e y,p and h y,p at all nodes, respectively. The resulting linear system of equations in (23)      The EM fields in the Laplace and Fourier domains are transformed into the space y and time t 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 y = 0 are taken into consideration, they are evaluated by Ns k=1 V k f x, k y , z, s k dk y , (30) where V k are the coefficients, N s is an even integer whose optimal value depends upon the computer word length. Gaver-Stehfest algorithm requires only the N s real values s k = k ln 2/t of the Laplace variable s 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 k y . 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.

Layered Earth Model.
To test the algorithm, the threelayered 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 2 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 3 Figure 4(a) has two main abnormal peaks, corresponding to the left and right edges of the horizontal tabular orebody. The profile of component y 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]. 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 z 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 International Journal of Geophysics 9 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.

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 2 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 x-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 x 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.

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, x and z components are sensitive to tectonic information. By the peaks of curves, the x component can reflect information about the plane position and geological occurrence of the abnormity object, while the y component is of great advantage to reveal the variation of geological strata along strike direction.