Numerical Modelling of TEM Fields by the Edge-Based Finite Element Method Using Tetrahedral Element for Central Rectangular Loops Modelado Numérico de Campos TEM Mediante el Método de Elementos Finitos Basados en Aristas Utilizando Elemento Tetraédrico para Bucles

In this study, we have implemented an edge-based finite element method for the numerical modeling of the transient electromagnetic method. We took the Helmholtz equation of the electric field as the governing equation for the edge-based finite element analysis. The modeling domain was discretized using linear tetrahedral mesh supported by Whitney-type vector basis functions. We inferred the equations by applying the Galerkin method. The system of equation was solved using a corrected version of the Bi Conjugate Gradient Stabilized Method (BiCGStab) algorithm to reduce the computational time. We obtained numerical solution for electric field in the Laplace domain; then the field was transformed into the time domain using the Gaver-Stehfest algorithm. Following this, the impulse response of the magnetic field was obtained through the Faraday law of electromagnetic induction as it is considerably more stable and computationally more efficient than inversion using the Fourier Transform. 3D geoelectric models were used to investigate the convergence of the edge-based finite element method with the analytic solution. The results are in good agreement with the analytical solution value for two resistivity contrasts in the 3D geoelectric brick model. We also compared the results of tetrahedral elements with the brick element in the 3D horizontal sheet and 3D conductive brick model. The results indicated that these two elements show very similar errors, but tetrahedral reflects fewer relative errors. For the low resistivity geoelectric model, numerical checks against the analytical solution, integral-equation method, and finite-difference time-domain solutions showed that the solutions would provide accurate results.


INTRODUCCION
The transient electromagnetic method (TEM) is widely applied in many fields such as mineral exploration (Zhdanov 2013), hydrocarbon exploration (Li and Constable 2010), tunnel prediction (Xue et al. 2017), etc. Loop-source Electromagnetic Method (EM) for the frequency and time domains play an important role in a wide range of geophysical investigations such as oil-sands exploration (Batayneh 2001;Yang and Oldenburg 2012;Zhdanov et al. 2013). Since 1980, TEM has been developed and refined intensively (Nedlec 1980). This makes the method relatively young compared to the frequency domain and magnetotelluric methods (Christiansen et al. 2006). The TEM data can be acquired in any one of the configurations, namely Central loop, In-loop, offset loop, and coincident loop methods (Taylor et al. 1992;Zhang and Xiao 2001). However, in practice, the central loop and coincident loop systems have been used more frequently (Nabighian 1988;Nabighian and James 1991;Karmis et al. 2003).
Integral-Equation Method (IEM), Finite-Difference Time-Domain method(FDTD), Finite-Volume, and Finite-Element (FE) Method have widely been applied for forward modeling of electromagnetic data. All of these numerical methods are based on either the background-anomalous field or the primary-secondary method. This method decomposes the total field into a background field and an anomalous field, which are typically evaluated by analytic and numerical methods, respectively. Compared with the IE and FDTD methods, complicated models such as irregular geologic bodies and surface topography can be described accurately by Finite Element Method (FEM) (Li and Hu 2017). In the finite-element modeling of electromagnetic fields, both edge and nodal elements play an important role. The difficulty in the modeling of field strengths that is caused by discontinuities in the properties of the media that can be overcome using edges of the elements. Edges of the elements can be used for computing electromagnetic fields in both homogeneous and inhomogeneous domains. The edge-based FEM adopt the vector basis functions defined on the center of the edges. The computation domain is usually discretized using either structured or unstructured meshes which employs rectangular and tetrahedral elements, respectively (Cai et al. 2014;Cai et al. 2017). Compared with the brick element, the tetrahedral element is preferred to represent complex geometry (Jin 2002). The edge-based FEM has some advantages such as imposing tangential boundary conditions easily (Whitney 1957;Jin 2002;Epov et al. 2007;Nam et al. 2007). Currently, the edge-based FEM is mainly used in numerical simulations of the TEM and magnetotelluric methods (Mitsuhata and Uchida 2004;Nam et al. 2007;Marinenko et al. 2009). Bossavit and Vertie (1982) discussed the construction of edge-based elements with tetrahedral and rectangular bricks for threedimensional eddy-current problems. Xiao et al. (2017) applied an edge-based finite element for modeling of magnetotelluric data. Xiao et al. (2017) applied an edge-based finite element for modeling of magnetotelluric data.
In this paper, the edge-based FEM with tetrahedral elements was employed to conduct the numerical modelling of the TEM for the central-loop configuration. We divided the computing domain into homogenous tetrahedral elements and used the Whitney-type vector basis function. As the divergence of Whitney function is zero and its curl is nonzero, it seems ideal for representing the vector field in the source-free region. In this way, we can obtain the system of the large symmetrical sparse equations which must be solved.
All the numerical methods for solving Ax = b can be divided into two broad categories, direct methods, and iterative methods. The iteration method is an attempt to solve a problem by finding successive approximations to the solution starting from an initial guess. Iterative methods are computationally more attractive than the direct methods particularly for large and sparse systems (Freund et al. 1992, Dreyer 2009). There are a wide variety of iterative techniques, such as Conjugate Gradient (CG), minimum residual, generalized minimum residual, Bi-conjugate gradient, quasi minimal residual, conjugate gradient squared, Bi-conjugate Gradient Stabilized (BiCGStab) and Chebyshev iterations (Cools 2019). We used the Corrected Bi Conjugate Gradient Stabilized Method (CBiCGStab) for the solutions of equations with unsymmetric coefficient matrices. The resulting CBiCGStab algorithm would maintain the favorable properties of the original method while not increasing computational costs (Yang and Brent 2002). After solving the equation with the CBiCGstab method, we derived the Laplace transform of the solution for TEM soundings and used the Gaver-Stehfest algorithm to invert it numerically to the time domain (Stehfest 1970, Davies andMartin 1978). Abate and Valko (2004) showed that the Gaver-Stehfest method would be quite accurate for functions which were not discontinuous or highly oscillatory and that this method would be much faster than the Discrete Fourier Transform method (DFT) since the vertical component of the impulse response of the magnetic field and Electromotive Force (EMF) could be obtained from the electric field at the edges of elements in the time domain with Gaver-Stehfest method.

2-1. Helmholtz Equation
In TEM, the current flowing in a transmitter loop sets up a magnetic field and induces an eddy current to flow in any good electrical conductor in the ground. These eddy currents set up a secondary magnetic field which can be detected by a receiver loop as time-dependent decaying voltage.
We took the Helmholtz equation of the electric field as the governing equation for the finite element analysis. The Helmholtz equation often arises in the study of electromagnetic problems which represents a time-independent form of the wave equation. The total fields are decomposed into the background electromagnetic fields (̂,̂) and anomaly electromagnetic fields (̂,̂) in the Laplace domain to avoid the source singularity problem (Cai et al. 2015): Without regard to displacement current, when no electric or magnetic abnormal body exists, we can calculate the Helmholtz equation in the Laplace domain using Maxwell's equations and Eq 1, as explained in (Mitsuhata and Uchida 2004;Marinenko et al. 2009;Cai et al. 2015): where = − and stands for the Laplace domain variable, 0 is the magnetic permeability, stands for angular frequency, ̂ and ̂a re the background and anomaly electric fields respectively in the Laplace domain, = + is the conductivity of the geoelectric model, denotes the background conductivity without the electric or magnetic abnormal body and stands for the abnormal conductivity.

2-2. Galerkin Method
Taking equation (2) as the governing equation of the finite element analysis, we infer the finite element equations by applying the Galerkin method (Jin 2002;Li et al. 2017). A popular and useful variation of the method of weighted residuals is the Galerkin method. (Li et al. 2017). Suppose a residual vector r for equation (2) is Then, based on the Galerkin method, the following is obtained in each element (Cai et al. 2015): In equation 4, N stands for the vector consisting of vector basis functions in different directions and is the volume differential. Figure 1 is an illustration of a tetrahedral element with node and edge indexing and in this figure, the relationship between node and edge is shown. From Fig 1 it can be seen that a direction is assigned to each edge and it points from one node to another ( Cai et al. 2015). In node-based FEM, according to the node numbers defined in figure 1, the anomalous electric field ̂ in a tetrahedral element can be determined based on the amount of the anomalous electric field (̂) in each node ( Jin 2002;Cai et al. 2017):

2-3. Edge-Based FEM
(5) Where is the number of nodes in each element and is the interpolation function and is obtained by is the volume of the element, , , and can be determined according to the nodes coordinate (Cai et al. 2017).
In the edge-based FEM, we have to determine vector basis function instead of scaler interpolation function in FEM. The edge number ( ) and the associated nodes 1 and 2 for edge are defined in Fig 1. The continuity of the tangential field across all element edges must be guaranteed; thus must have a tangential component along the th edge and none along the other edges. They must be ideal for representing the vector field in source-free region; thus another property of this function is that each satisfies the divergence condition .
= 0 within the element. Therefore, for basis vector function, we examine the vector function This vector basis function is called Whitney-form function. It is not difficult to see that where is the unit vector pointing from node 1 to node 2 and denotes the length of the edge connecting node 1 and node 2 .
1 2 has a constant tangential component along edge ( 1, 2 ) and has no tangential component along the other five edges. For the vector basis function to be in the same direction as the unit vector, we define the vector basis function for the edge i as follows: where . = 1.
Thus possesses all the necessary properties to be a vector basis function for the edge ith associated with nodes (i 1 ,i 2 ) (Newman et al. 1986;Harrington and Mautz 1991). The anomalous electric field is defined at the center of each edge and is denoted as , and the field inside the tetrahedral element can be represented by the following equation (Xiong 2011;Tiaaojie et al. 2019): For using Galerkin method, by substituting Eq 3 into Eq 4, we have and by substituting Eq 10 into Eq 11, one can find the original differential equation as follows: In these relations, ̂ is the vector that represents the unknown anomalous electric field at each edge in the Laplace domain and ̂ is the vector consisting of the background electric field in the Laplace domain on the edge. Using the expression given for in the Eq (9), we evaluate matrices and (Xiao et al. 2018).
For evaluation of the right-hand side b, the electric field components excited by rectangular Loop in the homogeneous half-space are required to calculate.

ELECTROMAGNETIC FIELD EXCITED BY RECTANGULAR LOOP
For solving equation 12, the electric field components in the homogeneous half-space need to be calculated. A rectangular coordinate system with its origin at the center of the rectangular loop is adopted (Fig 2). The rectangular transmitting loop belongs to the TE exciting mode, and the electric field contains the horizontal component only. The z-axis is downward vertically. The small vertical magnetic dipole coordinate is ( , , 0) and the rectangular loop area is integrated to obtain and at the point ( ′ , ′ , ). and are electric field component along x and y axis respectively.
In the equations 18 and 19, stands for the transmitter current, s is the Laplace variable and is equal to − , stands for angular frequency, is the conductivity of the homogeneous half-space, is the wavenumber, stands for the space between the receiver and the transmitter, ( ′ , ′ , ) is the coordinate of the receiver point, 0 is the Bessel function of the first kind, zero-order and is the transmission coefficient, which can be found from the following (Li et al. 2011 The receiver is on the surface and at the center of the transmitter, then = 0. The EMF can be obtained by the solution of the magnetic impulse response from the Maxwell equation on the ground. Therefore, the vertical component of the magnetic field response excited by the rectangular loop in the Laplace domain is: where the following relations are used We use equations 18 and 19 to evaluate ̂ and , then we should solve the system of equations ̂= to calculate the component of horizontal electric field in each edge and element.

4-1. Iteration Algorithm
All elements should be integrated as a whole to form a Global Matrix A and right-hand side b after they are analyzed. Then, we can obtain the system of equations .̂= where ̂ a is the vector which represents the unknown electric field anomaly field in the Laplace domain at each edge. For the solution of equations, we use the BiCGstab method which developed by Vost (1992). This method is an extension of the CG method and has widely been used in the solution of differential equations. To reduce computing time, Jacques et al. (2000) propose a modified parallel version of the BiCGStab method (MBiCGStab). They applied the developed solver with preconditioner to linear electromagnetic systems. Motivated by the same ideas, Yang and Brent (2002) proposed the Corrected BiCGStab method(CBiCGstab). This method is reorganized without changing the numerical stability, all inner products of a single iteration step are independent and subsequently, the communication time required for the inner product can be overlapped efficiently with a computation time of vector updates. Therefore, the time of computations on computers can significantly be reduced. The resulting CBiCGStab algorithm maintains the favorable properties of the original method while not increasing computational costs. The efficiency of this method is demonstrated by the numerical experimental results obtained from a massively parallel distributed memory system (Yang and Brent 2002). The convergence of the proposed CBiCGstab is almost the same as the MBiCGStab and original BiCGStab version. Experimental results of the speed up for these algorithms are given in Figure 3.These results are based on timing measurements of a fixed number of iterations.

4-2. Inverse Laplace Transform (Gaver-Stehfest Algorithm)
Using the component of horizontal electric field that is obtained from CBiCGstab method, we can determine impulse response of magnetic field from Maxwell equation in Laplace domain. This field must be transformed to the time domain with inverse Laplace transform. More generally, for the TEM problem, there are several inverse Laplace transform methods such as Fourier transform and Gaver-stehfest algorithm. The Fourier series technique for the Laplace inversion is based on choosing the contour of integration in the inversion integral, then converting the inversion integral into the Fourier transform, and finally approximating the transform by a Fourier series (Abate and Valko 2004). The Gaver-Stehfest algorithm was first proposed by Stehfest (1970) for numerical inverse Laplace transform. Abate and Valko (2004) compared the Fourier transform with Gaver-Stehfest algorithm for functions with oscillatory behavior in TEM (functions of type − ). also،they were compared these two methods with analytical inverse Laplace transform.
In the Gaver-Stehfest and Fourier methods, N is the expansion terms of the series. In the Gaver-Stehfest method For N=28, time of calculations was acceptable (t=0.641s) and the amount of error was minimum (standard deviation was 1.47 − 11). Tests of the Fourier method showed that the accuracy of computations to be lower than the Gaver-Stehfest method. In the Fourier method, the time of calculation was about 50 minutes and the amount of standard deviation was 7.63 − 03. The tests performed showed that in the investigated case (function of type − ), the faster and more accurate procedure would be the Gaver-Stehfest algorithm series (Abate and Valko 2004); thus, we used the Gaver-Stehfest algorithm to transform the electric field to the time domain from the Laplace domain. ̂( ) is the electromagnetic field in the Laplace domain and E(t) is the electromagnetic field in the time domain. The Gaver-Stehfest method uses the summation (Harrington and Mautz 1991;Nabighian and James 1991;krougly et al. 2013): where N is the expansion terms in the series. coefficients only depend on the number of expansion terms N and they are: Since the 2 in equation 26 is a natural number, N must be even. We achieved the best results with N = 28 which this number is called optimum value. It was not possible to get better results by increasing N beyond the optimum, as we increased the number of terms N in the computation, we quickly discovered that the numerical inversion would become unstable and our function would be subject to numerical errors.
We have coded the algorithms using Matlab and used DPLOT for graphs. To test background electromagnetic fields which have obtained from algorithms, we used Electromagnetic Analysis Software (http://www.hgg.au.dk).

EXAMPLES
In this study, for edge-based FEM, we used a tetrahedral element instead of the brick element. The major disadvantage of the brick element is that it is restricted to a limited class of geometries. For comparing time and accuracy between tetrahedral and brick elements, we chose a model that could be solved with both and would often be used in electromagnetic interpretations. This model had already been investigated by Li et al. (2011) with a brick element. The 3D thin horizontal sheet and rectangular body are common configurations for conductive zones (such as the mineralization of massive pyrite and sulfide body in the ground).
The transient electromagnetic field diffuses outward and downward from the exciting source. A particular advantage of TEM Systems is the fact that the measurement is taken when the transmitted fields are switched off. To further verify the correctness of the algorithm, it is compared with the analytical solutions, IE and finite difference time domain method(FDTD) in the conductive 3D geoelectric brick model.

5-1. Example A: 3D Low Resistivity Horizontal Sheet
The parameters of the geoelectric model are shown in Figure 4. The air resistivity is set as 10 10 Ω , the homogeneous half-space resistivity is 10 2 Ω , and the resistivity of the horizontal low resistivity sheet is 1 Ω . The side of the transmitter loop is 100 m long, the equivalent of the receiving coil is 1 2 and the transmitting current is 1 A. In this model, in addition to using tetrahedral elements, we used brick elements to compare the results. It should be noted that to compare two kinds of the element, we should try to preserve the same condition for both. It means that the number of the unknown edges in tetrahedral and brick meshes need to be almost equal. Three Matrix size can be adopted as shown in Tables 1 and 2 for the mesh size in brick and tetrahedral meshes. The number of nodes and edges depends on the maximum length of the edge in meshes. We increased the nodes and edges in the three schemes to compare the results. The number of edges in each scheme at both meshes are approximately equal. Due to different geometric shapes in brick and tetrahedral elements, we were not able to provide a completely identical number for edges in both meshes. The calculated results for the horizontal sheet in different mesh schemes for tetrahedral and brick meshes are shown in Figs 5 and 6. It can be seen that scheme 1 gives poor results and scheme 3 is the best scheme in brick and tetrahedral meshes (Figs 5 and 6).
The EMF curves of brick and tetrahedral meshes had the same general behavior. In the tetrahedral mesh (Figs 5 and 6), prior to 3 ms, the EMF curves calculated by the three mesh schemes are basically uniform. From among 4 ms, the curve of tetrahedral meshes have less error than the brick meshes. For both tetrahedral and brick meshes, the relative errors are shown in figure 7. For tetrahedral mesh, it is observed from the curves of the relative error (Fig. 7), that the relative errors of the various mesh schemes were all within 2% in 0.8 ms, and when time is 3.4 ms, the relative errors for the three mesh schemes are 12.5%, 7%, and 4.5% or so, respectively. In 13 ms, the relative errors for these three mesh schemes are 27%, 13%, and 6% or so, respectively. For brick element meshes, the relative errors of the various mesh schemes are all within 4% in 0.8 ms, and when time is 3.4 ms, the relative errors for the three mesh analytic schemes are 18%, 10%, and 6% or so, respectively.
At time 13 ms, the relative errors for the three mesh analytic schemes are 34%, 17.5%, and 7.5% or so, respectively (Fig 7). As the figure shows, early-time relative errors in brick and tetrahedral meshes are less than the late-time relative errors. The phenomenon can be interpreted by the 'smoke ring' theory of transient electromagnetic field dispersion, in the horizontal direction, the electromagnetic field has less scope for dispersion in the early stage, and the induction current is limited to the horizontal low resistivity sheet.
With increasing time, the induction current disperses out of the horizontal low resistivity sheet. The comparison between brick and tetrahedral meshes shows that relative errors in scheme 1,2 and 3 for tetrahedral, are less than the brick element. It is shown in Fig. 8 that the number of edges is important in amount of relative errors. The time in these curves is 8 ms. It is observed from figure 8 that with an increase in the number of edges in the model, the relative error between brick and tetrahedral meshes would decrease. Given that the speed of the edge-based FEM is slower than FEM, sometimes the number of edges should be reduced and in this situation, the tetrahedral mesh has much less error than the brick mesh.
For tetrahedral mesh, the multi-measuring point EMF curves are shown in Fig 9. The measuring line was 1500 m long perpendicular to the strike with 10 measuring points in total. The EMF curve at the different moments is one canyon flanked by two mountains. The forming phenomenon is related to the number of magnetic lines cut by the upright low resistivity sheet. In the early stage, the electromagnetic field is mainly centralized to the surface and the EMF received in the measuring points is mainly related to the number of magnetic lines cutting the top of the sheet and showed background field. These curves have been drawn for scheme 3 in the tetrahedral mesh (Fig 9).

5-2. Example B: 3D Vertical Geoelectric brick model with Different Resistivity Contrast
In our design, the resistivity of the homogeneous half-space was set as 20 Ω , and 100 Ω respectively, and low resistivity sheet resistivity was set as 2 Ω (Fig10). The transmitting side was 100 m long, and the transmitting current was 1A. The equivalent area for the receiving coil was 1 2 . Erath was divided into 2892 tetrahedral elements that consisted of 17358 edges and 697 nodes. The obtained results in two resistivity contrasts are shown in Fig 11. The edge-based FE solution fitted right in with the analytical solution value for two resistivity contrast models. For the model with a homogeneous half-space resistivity of 20 Ω , the relative error was less than 1.18% at 0.2 ms and was less than 2.28% at 1.29 ms. Fig 11. The EMF for two resistivity contrast in tetrahedral mesh for example B As for models with homogeneous half-space resistivity of 100 Ω , the relative error reaches 17% at 1.5 ms. It is seen in Fig. 12 that with the increase in the homogeneous half-space resistivity, the EMF would decrease, while the EMF error solved by the edge-based FEM with tetrahedral mesh would become bigger and bigger, so the edge-based finite element would be suitable for low resistivity contrast. In the initial stages, currents penetrate in the vicinity of the surface and curve show background field. The curves are sharper than a 3D horizontal sheet (Fig 9).

5-3. Example C: 3D Low Resistivity Geoelectric brick Model
To further verify the correctness of the algorithm, the algorithm in this paper was contrasted with the Integral Equation and FDTD methods ( Newmann et al. 1986;Wang anf Hohmann 1993). Model C is shown in Fig 14.   Fig 14. The 3D conductive brick model (example C) The side of the transmitting loop was 100 m long, and the observing value was the vertical component of the EMF. Erath was divided into 682 brick and 1487 tetrahedral elements. The numbers of nodes in brick and tetrahedral meshing were 880 and 389 respectively. The brick mesh consisted of 8184 edges and there were 8922 edges in tetrahedral meshes. To further verify the correctness of the edge-based FEM with tetragedral mesh, it is compared with the analytical solutions, edge-based FEM with brick mesh (Li et al. 2011), IE and finite difference time domain method(FDTD) in the conductive 3D geoelectric brick model, the behavior of the four curves stayed the same over the time. The edge-based FE solution fitted well with the other methods from 0 to 1.5 ms. From 1.5 to 7 ms, the edge-based FE solution was bigger than the analytical solution; however, the IEM and FDTD method solutions were less than the analytical solution. (Fig 15). The IE and FDTD methods reflect the same error rates in early and late time however, from 1.5 ms to 7 ms, tetrahedral meshes reflect fewer errors than brick elements. Between the four curves, edge-based FEM with tetrahedral element has the best agreement with analytical solution. For this model (example C), the multi-measuring point EMF curves for tetrahedral mesh are shown in Fig 16(edge-based-FE solution). The measuring line was 160 m long perpendicular to the strike with 16 measuring points in total. The distance between two maximums is proportional to the width of the body and these points are seen on the boundaries of the model.

CONCLUSIONS
The 3D numerical modeling is performed for the central loop configuration of the TEM using the edgebased finite element method. Although the brick element is a simple mesh and an easy finite-element analysis, the mesh is not used to precisely describe the complicated geoelectric model and surface topography. Thus, in this study, the linear tetrahedron element was used to replace the brick element. We validated the method using several 3D geoelectric models with different mesh schemes and different resistivity contrasts. The numerical studies showed that using the tetrahedral element can produce an accurate result for TEM modeling. The edge-based finite element solutions of these models show a good agreement with the analytical and integral equation solutions. In the 3D horizontal sheet, the numerical results also demonstrated that the tetrahedral mesh commits fewer errors than the brick mesh. To reduce the computational time, the number of unknown parameters and edges need to be reduced and this will cause more errors in solutions. To overcome this problem, choosing a tetrahedral mesh is a suitable option. By taking the low resistivity 3D brick geoelectric model as an example, the results of the EMF obtained by the edge-based finite element method, integral equation method, and finite difference time domain method are consistent with each other. In this paper, the brick and sheet geoelectric models are applied for comparing tetrahedral and brick meshes. Thus in the subsequent studies, complicated and combined models can also be used. This method can also be applied and tested for separate transmitterreceiver loops in the transient electromagnetic method.