3-D dc resistivity modelling based on spectral element method with unstructured tetrahedral grids

SUMMARY

In this paper, we propose a spectral element method (SEM) based on unstructured tetrahedral grids for direct current (dc) resistivity modelling. Unlike the tensor-product of 1-D Gauss–Lobatto–Legendre (GLL) quadrature in conventional SEM, we use Proriol–Koornwinder–Dubiner (PKD) polynomials to form the high-order basis polynomials on tetrahedral grids. The final basis functions are established by using Vandermonde matrix. Compared to traditional SEM, our method truly takes into account the high precision of spectral method and the flexibility of finite element method with unstructured grids for modelling the complex underground structures. After addressing the theory on the construction of basis functions and interpolation and integration nodes, we validate our algorithm using the analytical solutions for a layered earth model and the results from other methods for multiple geoelectrical models. We further investigate a dual-track scheme for improving the accuracy of our SEM by increasing the order of interpolation polynomials or by refining the grids.

Electrical properties, Electromagnetic theory, Numerical approximations and analysis, Numerical modelling

1 INTRODUCTION

The direct current (dc) resistivity method is one of the classic exploration methods in geophysics. Due to simple equipment, low cost, diverse observation methods, it is widely used in mining, groundwater and geothermal exploration and in various engineering or environmental studies (Rucker et al. 2009; Abu Rajab & El-Naqa 2013; Weiss et al. 2016).

However, except for some simple models that have analytic solutions, for most complex models the numerical simulation is feasible for data interpretation. Fortunately, the theory of dc resistivity method is relatively simple and the numerical simulation for this method has early entered the 3-D era (Günther et al. 2006; Rücker et al. 2006; Beskardes & Weiss 2018). Since a fast and reliable forward modelling builds the foundation for inversions, much attention has been paid to 3-D modelling. There are mainly three mainstream methods, including integral equation (IE), finite difference (FD) and finite element (FE) method. Each of them have their own pros and cons, but all of them must solve two fundamental problems associated with the accuracy and speed. Among these methods, the IE method only handles the abnormal bodies by the regular grids, it is fast and accurate, however, it cannot simulate complicated models (Xu et al. 1988; Zhdanov et al. 2000; Hvoždara 2012). The FD method can be easily implemented, but it has low accuracy and is limited to structured grids (Spitzer 1995; Hou et al. 2006; Liu & Yin 2014). In comparison, the FE method is more flexible due to its capability of using unstructured grids and thus can handle complex models and irregular interfaces. It has now been developed into an efficient and effective modelling tool (Zhou & Greenhalgh 2001; Li & Spitzer 2002; Yang et al. 2017).

In FE method, there generally exist two strategies for improving the model accuracy: to refine grids (h-type) and to raise the order of interpolation functions (p-type). Between the two strategies, the grid refinement has been more widely used, especially with the application of the adaptive unstructured FE method (Ren & Tang 2010; Ren & Tang 2014; Yin et al. 2016). The other strategy, namely the high-order interpolation generally takes nodes at the equal interval (Silvester 1969). However, with increasing order, the Runge phenomenon (Runge 1901) may occur that can create serious oscillations at the ends of the domain. In that case, a non-linear interpolation may be a good solution and the spectral method (SM) was introduced. Its mathematical basis was set up by Gottlieb & Orszag (1977). Maday & Quarteroni (1981) and Canuto et al. (1988) systematically discussed the polynomial approximation and application of Fourier method. However, as a global method, the SM cannot be used to simulate complex models. To solve this problem, Patera (1984) combined the FE method with the spectral method and put forward the spectral element method (SEM) for fluid dynamics. The classic SEM uses the Gauss–Lobatto–Legendre (GLL) polynomials for basis functions with the numerical integration realized on GLL nodes. In this way, one can obtain a diagonal mass matrix for the modelling and enhancing simultaneously the precision and speed (Komatitsch & Tromp 1999). Komatitsch (2002) used the SEM to model the propagation of the global seismic wave and solved the seismic modelling problem of large scale. After that, the SEM has also been used in EM modelling (Lee et al. 2006; Yin et al. 2017). Since GLL polynomials can only be applied to rectangular or hexahedron, to fit irregular interfaces, one needs to use very small deformed hexahedrons, and as a result, the strict condition for stability must be enforced and the modelling accuracy is reduced. On the other hand, a triangle or tetrahedron can well model the complex and irregular interfaces with relatively homogeneous grids (Wang et al. 2013). This stimulates the research in this paper to develop a spectral element method based on tetrahedral grids (TSEM).

We will deal with three key issues—the basis functions, the interpolation nodes and the numerical integration. In classic SEM, the rectangular or hexahedral elements have independent edge directions, so that one can obtain the nodes for multidimensional interpolations and integrations via tensor product extension of 1-D GLL polynomials (Karniadakids & Sherwin 2005). However, as the triangular and tetrahedral elements are all calculated in simplex region, the edge directions are coupled, it is till now impossible to find similar GLL basis functions. Thus, we try to use Proriol–Koornwinder–Dubiner (PKD) polynomials (Proriol 1957; Koornwinder 1975; Dubiner 1991) that are the compound polynomials built for the coupled directions. PKD polynomials cannot be used directly as the basis functions for they do not satisfy |${\varphi _i}{\rm{(}}{x_j}{\rm{)}} = {\delta _{ij}}$|⁠, thus a Vandermonde matrix is needed to connect PKD polynomials and basis functions. Note that, in this case the expansion with nth polynomials combines all elements of (n–1)th polynomials, so that the highest order in the basis functions is n. This kind of basis functions are called the hierarchical basis functions (Sherwin 1997). When the order of the basis functions is 1, the SEM becomes the traditional FE method.

The SEM has obvious advantage in solving the time-domain problem, because one can build a diagonal mass matrix when the same interpolation and integration points are used, and thus allowing explicit time discretization and reducing the computational cost (Cohen et al. 1994; Archer & Whalen 2005; Liu et al. 2014). Liu et al. (2017) proposed mass-lumping technique based on Lp-norm optimization for triangular elements and obtained a set of nodes that can be used for both interpolation and integration. Geevers et al. (2018) designed the nodes for tetrahedral elements along the same line, the maximum order can be 4. However, to ensure the consistency of the nodes for both interpolation and integration, additional points need to be added to reach a compromise between them (Helenbrook 2009). It implies that the number of interpolation and integration points will increase, while neither of the two point sets are optimal.

In this paper, we study the 3-D forward modelling problem for dc resistivity method using TSEM, we choose different nodes for interpolation and integration, so that we have good interpolation nodes to guarantee the property of the interpolation while we can have a high accuracy by using Gaussian-type integrals (Pasquetti & Rapetti 2006). The Fekete nodes (Fekete 1923) are the optimal set of interpolation nodes for triangles. At the edges of the triangles, they are GLL points that can be used in construction of hybrid grids because they can well match those of a rectangular element. (Taylor & Wingate 1999; Komatitsch et al. 2001; Pasquetti & Rapetti 2004). However, it is not the case for a tetrahedral element. There are many ways to evaluate the property of an interpolation, such as maximization of Vandermonde determinant (Taylor et al. 2001), minimization of Lebesgue constant (Smith 2006), minimization of electrostatic potential (Hesthaven & Teng 2000), etc. In this paper, we take the explicit construction method proposed by Warburton (2006) to calculate the interpolation nodes. Pasquetti & Rapetti (2010) confirmed that the nodes obtained in this way have good interpolation performance. For the numerical integration, we use the symmetric integration nodes proposed by Zhang et al. (2009). We will focus our attention on establishing the SEM basis functions and constructing the interpolation nodes. After validating the accuracy of our algorithm to 3-D dc resistivity modelling, we will analyse the characteristic of dc resistivity for typical abnormal bodies and the influence of the topography.

2 METHODOLOGY

2.1 Governing equations

We first establish a Cartesian coordinate system (x, y, z) in the physical domain |$\Omega $| with x, y axes locating at the earth surface, while the z-axis is positive downwards. In such a coordinate system, the forward problem for dc resistivity method can be written as (Ren & Tang 2010)

$$\begin{eqnarray*} \nabla \cdot \left[ {\sigma {\rm{(}}{{\bf r}}{\rm{)}}\nabla U{\rm{(}}{{\bf r}}{\rm{)}}} \right] = - I\delta {\rm{(}}{{\bf r}} - {{{\bf r}}_0}{\rm{),}} \end{eqnarray*}$$

(1)

where |${{{\bf r}}_0} = {\rm{(}}{x_0}{\rm{,}}{y_0}{\rm{,}}{z_0}{\rm{)}}$| and |${{\bf r}} = {\rm{(}}x{\rm{,}}y{\rm{,}}z{\rm{)}}$| denote, respectively, the source and receiving points, I is the transmitting current,|$\sigma ({{\bf r}})$| is the conductivity, U is the electrical potential and δ is the Dirac delta function.

To solve the electrical potential uniquely, we need to add the boundary conditions to eq. (1). At the interface |${\Gamma _0}$| between the earth and the air, there is no current flowing into the air, thus we assume the following Neumann boundary condition:

$$\begin{eqnarray*} {{\bf n}} \cdot \nabla U{\rm{(}}{{\bf r}}{\rm{) }}\left| {_{{\Gamma _0}}} \right. = 0\,\,{\rm{,}} \end{eqnarray*}$$

(2)

where n denotes the outward normal to the boundary. At the outside boundaries |${\Gamma _\infty }$|⁠, the potential satisfies the following mixed boundary condition (Dey & Marrison 1979), that is

$$\begin{eqnarray*} {{\bf n}} \cdot \nabla U{\rm{(}}{{\bf r}}{\rm{)}} + \frac{{{\rm{cos(}}{{\bf r}} - {{{\bf r}}_0}{\rm{,}}{{\bf n}}{\rm{)}}}}{{\left| {{{\bf r}} - {{{\bf r}}_0}} \right|}}U{\rm{(}}{{\bf r}}{\rm{) }}\left| {_{{\Gamma _\infty }}} \right. = 0\,\,{\rm{.}} \end{eqnarray*}$$

(3)

We use the weighted residual method to solve the forward problem in eqs (1)–(3). For this purpose, we multiply the both sides of eq. (1) by the weighting function W in the Hilbert space and take into account the boundary conditions (Li & Spitzer 2002) and get

$$\begin{eqnarray*} && \iiint_\Omega {\sigma \nabla W \cdot \nabla U\,\,{\rm{d}}\Omega } + \iint_{{\Gamma _\infty }} {\sigma \frac{{{\rm{cos(}}{{\bf r}} - {{{\bf r}}_0}{\rm{,}}{{\bf n}}{\rm{)}}}}{{\left| {{{\bf r}} - {{{\bf r}}_0}} \right|}}WU\,\,{\rm{d}}\Gamma } \nonumber \\ && \qquad = \iiint_\Omega {WI\delta {\rm{(}}{{\bf r}} - {{{\bf r}}_0}{\rm{)}}\,\,{\rm{d}}\Omega }. \end{eqnarray*}$$

(4)

Further, we divide the solution domain into tetrahedral elements and expand the electrical potential according to the basis functions |${\varphi _i}(i = 1{\rm{,}} \ldots {\rm{,}}\,\,Nt)$| with Nt denoting the number of interpolation nodes, that is

$$\begin{eqnarray*} U = {{{\bf U}}^{\rm{T}}}{{\boldsymbol \varphi }} = \sum\limits_{i = 1}^{Nt} {{U_i}{\varphi _i}} . \end{eqnarray*}$$

(5)

In the following, we take an element e as example to do the element analysis. By using Galerkin method, we have

$$\begin{eqnarray*} {{{\bf K }}^e}{{{\bf U}}^e} = {{{\bf B}}^e}, \end{eqnarray*}$$

(6)

where U^e is the vector for the electrical potential, |${{{\bf K }}^e}$| is the coefficient matrix and |${{{\bf B}}^e}$| is the source term. They can be expressed as

$$\begin{eqnarray*} K_{ij}^e = \iiint_{{\Omega _{^e}}} {\sigma \nabla {\varphi _i} \cdot \nabla {\varphi _j}{\rm{d}}\Omega } + \iint_{{\Gamma _e}} {\sigma \frac{{{\rm{cos(}}{{\bf r}} - {{{\bf r}}_0}{\rm{,}}{{\bf n}}{\rm{)}}}}{{\left| {{{\bf r}} - {{{\bf r}}_0}} \right|}}{\varphi _i}{\varphi _j}\,\,{\rm{d}}\Gamma }{\rm{,}} \end{eqnarray*}$$

(7)

$$\begin{eqnarray*} B_i^e = \iiint_{{\Omega _{^e}}} {{\varphi _i}I\delta {\rm{(}}{{\bf r}} - {{{\bf r}}_0}{\rm{)}}\,\,{\rm{d}}\Omega }. \end{eqnarray*}$$

(8)

In eq. (7), Ω_e and Γ_e define the volume and surface of element e, respectively. The second term exists only when |${\Gamma _{\rm{e}}} \in {\Gamma _\infty }$|⁠.

After we finish processing all elements, we can assemble them together to get the matrix equation for the whole solution domain:

$$\begin{eqnarray*} {{\bf K U}} = {{\bf B}}{\rm{,}} \end{eqnarray*}$$

(9)

where the dimension of the electrical potential vector is equal to total number of interpolation nodes in the model domain. Solving eq. (9), we can get the electrical potentials and calculate the apparent resistivity.

2.2 Construction of basis functions

In this section, we first construct the basis functions |${\varphi _i}{\rm{(}}{{{\bf r}}_j}{\rm{)}}$|⁠, which are actually Lagrange polynomials and can be written as (Liu et al. 2017)

$$\begin{eqnarray*} {\varphi _i}{\rm{(}}{{{\bf r}}_j}{\rm{)}} &=& c_i^1{\psi _1}{\rm{(}}{{{\bf r}}_j}{\rm{)}} + c_i^2{\psi _2}{\rm{(}}{{{\bf r}}_j}{\rm{)}} + \cdots + c_i^{Nt}{\psi _{Nt}}{\rm{(}}{{{\bf r}}_j}{\rm{)}} \nonumber \\ &=& \sum\limits_{k = 1}^{Nt} {c_i^k{\psi _k}{\rm{(}}{{{\bf r}}_j}{\rm{)}}} = {\delta _{ij}}, \end{eqnarray*}$$

(10)

where |${\psi _i}{\rm{(}}{{{\bf r}}_j}{\rm{)}}$|are the orthogonal functions, |${\delta _{ij}} = \Big\{\ _{{0,}\,\,{i\,\,\ne\,\,j}}^{{1,}\,\,{i\,\,=\,\,j}}.\,{c_i^k}$| are expansion coefficients that are connected with the general Vandermonde matrix V as following:

$$\begin{eqnarray*} {{\bf Vc}} = {{\bf I}}, \end{eqnarray*}$$

(11)

where I is a unit matrix and

$$\begin{eqnarray*} {{\bf V}} = \left[ {\begin{array}{@{}*{4}{c}@{}} {{\psi _1}{\rm{(}}{{{\bf r}}_{\rm{1}}}{\rm{)}}} & \quad{{\psi _2}{\rm{(}}{{{\bf r}}_{\rm{1}}}{\rm{)}}} & \quad \cdots & \quad{{\psi _{Nt}}{\rm{(}}{{{\bf r}}_{\rm{1}}}{\rm{)}}}\\ {{\psi _1}{\rm{(}}{{{\bf r}}_{\rm{2}}}{\rm{)}}} & \quad{{\psi _2}{\rm{(}}{{{\bf r}}_{\rm{2}}}{\rm{)}}} & \quad \cdots & \quad{{\psi _{Nt}}{\rm{(}}{{{\bf r}}_{\rm{2}}}{\rm{)}}}\\ \vdots & \quad \vdots & \quad \ddots & \quad \vdots \\ {{\psi _1}{\rm{(}}{{{\bf r}}_{Nt}}{\rm{)}}} & \quad{{\psi _2}{\rm{(}}{{{\bf r}}_{Nt}}{\rm{)}}} & \quad \cdots & \quad{{\psi _{Nt}}{\rm{(}}{{{\bf r}}_{Nt}}{\rm{)}}} \end{array}} \right], \end{eqnarray*}$$

(12)

$$\begin{eqnarray*} {{\bf c}} = \left[ {\begin{array}{@{}*{4}{c}@{}} {c_1^1} & \quad {c_2^1} & \quad \cdots & \quad {c_{Nt}^1}\\ {c_1^2} & \quad {c_2^2} & \quad \cdots & \quad {c_{Nt}^2}\\ \vdots & \quad \vdots & \quad \ddots & \quad \vdots \\ {c_1^{Nt}} & \quad {c_2^{Nt}} & \quad \cdots & \quad {c_{Nt}^{Nt}} \end{array}} \right]. \end{eqnarray*}$$

(13)

Altogether, we can obtain the following relation

$$\begin{eqnarray*} {{\boldsymbol \varphi }} &=& \left[ {\begin{array}{@{}*{1}{c}@{}} {{\varphi _1}{\rm{(}}{{\bf r}}{\rm{)}}}\\ {{\varphi _2}{\rm{(}}{{\bf r}}{\rm{)}}}\\ \vdots \\ {{\varphi _{Nt}}{\rm{(}}{{\bf r}}{\rm{)}}} \end{array}} \right] = \left[ {\begin{array}{@{}*{4}{c}@{}} {c_1^1} & \quad {c_1^2} & \quad \cdots & \quad {c_1^{Nt}}\\ {c_2^1} & \quad {c_2^2} & \quad \cdots & \quad {c_2^{Nt}}\\ \vdots & \quad \vdots & \quad \ddots & \quad \vdots \\ {c_{Nt}^1} & \quad {c_{Nt}^2} & \quad \cdots & \quad {c_{Nt}^{Nt}} \end{array}} \right]\left[ {\begin{array}{@{}*{1}{c}@{}} {{\psi _1}{\rm{(}}{{\bf r}}{\rm{)}}}\\ {{\psi _2}{\rm{(}}{{\bf r}}{\rm{)}}}\\ \vdots \\ {{\psi _{Nt}}{\rm{(}}{{\bf r}}{\rm{)}}} \end{array}} \right] = {{{\bf c}}^{\rm{T}}}{{\boldsymbol \psi }} \nonumber \\ &=& {{\rm{(}}{{{\bf V}}^{{\rm{ - 1}}}}{\rm{)}}^{\rm{T}}}{{\boldsymbol \psi }} \end{eqnarray*}$$

(14)

To obtain eq. (14), we need to find the polynomials |${\psi _i}{\rm{(}}{{{\bf r}}_j}{\rm{)}}$| with good orthogonality and to ensure the numerical stability and the reversibility of the Vandermonde matrix. The polynomials are defined in the reference domain of right-angled tetrahedron with vertex |${\rm{(}} - {\rm{1,}} - {\rm{1,}} - {\rm{1),}}\,\,\,\,{\rm{(1,}} - {\rm{1,}} - {\rm{1),}}\,\,\,\,{\rm{(}} - {\rm{1,1,}} - {\rm{1),}}{\rm{(}} - {\rm{1,}} - {\rm{1,1)}}$| that can be expressed as |${\Omega _s} = \{ { {{\rm{(}}\xi {\rm{,}}\eta {\rm{,}}\zeta {\rm{)}}} |{\rm{(}} \xi {\rm{,}}\eta {\rm{,}}\zeta {\rm{)}}\geq\ -1{\rm{,}}\,\,\,\,\,\,\xi + \eta + \zeta \le - 1} \}$|⁠. For this purpose, we define an n-dimensional (the order of polynomials) general canonical function (Hesthaven & Warburton 2008):

$$\begin{eqnarray*} {\psi _N} = {\xi ^i}{\eta ^j}{\zeta ^k}, \end{eqnarray*}$$

(15)

where |$\xi {\rm{,}}\,\,\,\,\eta {\rm{,}}\,\,\,\,\zeta $| are direction parameters, i, j, k are the orders that satisfy |$i{\rm{,}}j{\rm{,}}k \ge 0{\rm{,}}\,\,\,\,i + j + k \le n$|⁠. |$Nt = \frac{{{\rm{(}}n + 1{\rm{)(}}n + 2{\rm{)(}}n + 3{\rm{)}}}}{6}$| is the number of interpolation points in tetrahedral, |$1 \le N \le Nt$| is a composite index connected to the order n. Corresponding to the permutation of i, j, k, we have |$\sum\limits_{N = 1}^{Nt} \Leftrightarrow \sum\limits_{i = 0}^n {\sum\limits_{j = 0}^{n - i} {\sum\limits_{k = 0}^{n - i - j} {} } } $| and

$$\begin{eqnarray*} N &=& 1 + \frac{{{\rm{(}}11 + 12n + 3{n^2}{\rm{)}}i}}{6} + \frac{{{\rm{(}}2n + 3{\rm{)}}j}}{2} + k - \frac{{{\rm{(}}2 + n{\rm{)}}{i^2}}}{2} \nonumber \\ &&- ij - \frac{{{j^2}}}{2} + \frac{{{i^3}}}{6}. \end{eqnarray*}$$

(16)

Then, we can use Gram–Schmidt orthogonalization to get PKD polynomials

$$\begin{eqnarray*} {\psi _N}{\rm{(}}\xi {\rm{,}}\eta {\rm{,}}\zeta {\rm{)}} &=& \sqrt {\displaystyle\frac{8}{{{\gamma _{ijk}}}}} {P_i}^{{\rm{(}}0{\rm{,}}0{\rm{)}}}\left( { - 2\displaystyle\frac{{1 + \xi }}{{\eta + \zeta }} - 1} \right)P_j^{{\rm{(}}2i + 1{\rm{,}}0{\rm{)}}} \left( {2\displaystyle\frac{{1 + \eta }}{{1 - \zeta }} - 1} \right)\nonumber \\ && \times {{\left( {2 - 2\displaystyle\frac{{1 + \eta }}{{1 - \zeta }}} \right)}^i}P_k^{{\rm{(}}2i + 2j + 2{\rm{,}}0{\rm{)}}}{\rm{(}}\zeta {\rm{)(}}1 - \zeta {{\rm{)}}^{i + j}} \nonumber \\ &=& \!\sqrt {\!\displaystyle\frac{8}{{{\gamma _{ijk}}}}} {P_i}^{{\rm{(}}0{\rm{,}}0{\rm{)}}}{\rm{(}}r{\rm{)}}P_j^{{\rm{(}}2i \!+\! 1{\rm{,}}0{\rm{)}}}{\rm{(}}s{\rm{)(1}} \!-\! s{{\rm{)}}^i}P_k^{{\rm{(}}2i + 2j + 2{\rm{,}}0{\rm{)}}}{\rm{(}}t{\rm{)(1}} \!-\! t{{\rm{)}}^{i + j}}{\rm{,}} \nonumber \\ \end{eqnarray*}$$

(17)

$$\begin{eqnarray*} r = - 2\frac{{1 + \xi }}{{\eta + \zeta }} - 1{\rm{,}}\,\,\,\,s = 2\frac{{1 + \eta }}{{1 - \zeta }} - {\rm{1,}}\,\,\,\,t = \zeta , \end{eqnarray*}$$

(18)

where |${\gamma _{ijk}} = \frac{2}{{2i + 1}}\frac{{{2^{2i + 2}}}}{{2i + 2j + 2}}\frac{{{2^{2i + 2j + 3}}}}{{2i + 2j + 2k + 3}}$|⁠,|$P_i^{{\rm{(}}\alpha {\rm{,}}\beta {\rm{)}}}{\rm{(}}x{\rm{)}}$| are orthogonal Jacobian polynomials of ith order with the orthogonal weight |$w{\rm{(}}x{\rm{)}} = {{\rm{(}}1 - x{\rm{)}}^\alpha }{{\rm{(}}1 + x{\rm{)}}^\beta }$|⁠. As we know, the Jacobian polynomial is a solution to the singular Sturm–Liouville eigenvalue problem that can be solved by recurrence method. When |$\alpha = \beta = 0$|⁠, |$P_i^{{\rm{(}}\alpha {\rm{,}}\beta {\rm{)}}}{\rm{(}}x{\rm{)}}$| becomes the Legendre polynomial, when |$\alpha = \beta = \pm1/2$|⁠, it becomes the Chebyshev polynomial.

2.3 Coordinate transformation

From the above discussion, it can be seen that the physical domain and reference domain are not identical, thus we need to build a mapping relation between the physical coordinate (x, y, z) and reference coordinate (ξ, η, ζ) using the affine transformation (cf. Fig. 1a). Besides, to obtain the coefficient matrix, we need to utilize orthogonality and calculate the partial derivative of the basis functions. To simplify the calculation, we introduce new variables (r, s, t) for the polynomials in eq. (18), where (r, s, t) belong to a standard hexahedron domain |${\Omega _H} = \{ { {{\rm{(}}r{\rm{,}}s{\rm{,}}t{\rm{)}}} | - 1 \le r{\rm{,}}s{\rm{,}}t \le 1} \}$|⁠. The mapping between the reference tetrahedral and the standard hexahedron is called collapsed transformation (cf. Fig. 1b). Thus, to accomplish the calculation related to the basis functions, we need to establish the reversible mapping relations between three coordinate systems (physical, reference and standard hexahedron), namely |${\rm{(}}x{\rm{,}}y{\rm{,}}z{\rm{)}} \Leftrightarrow {\rm{(}}\xi {\rm{,}}\eta {\rm{,}}\zeta {\rm{)}} \Leftrightarrow {\rm{(}}r{\rm{,}}s{\rm{,}}t{\rm{)}}$| (Karniadakis & Sherwin 2005). Fig. 1 shows the detailed mapping relations between these coordinate systems.

Figure 1.

(a) Affine transformation between physical coordinate |${\rm{(}}x{\rm{,}}y{\rm{,}}z{\rm{)}}$| and reference coordinate |${\rm{(}}\xi {\rm{,}}\eta {\rm{,}}\zeta {\rm{)}}$|⁠; (b) Collapsed transformation between reference coordinate |${\rm{(}}\xi {\rm{,}}\eta {\rm{,}}\zeta {\rm{)}}$| and standard hexahedral coordinate |${\rm{(}}r{\rm{,}}s{\rm{,}}t{\rm{)}}$|⁠.(refer to Karniadakis & Sherwin 2005)

To realize the affine transformation, we need only to design a rotationally invariant barycentric coordinates |${\rm{(}}{\nu _1}{\rm{,}}{\nu _2}{\rm{,}}{\nu _3}{\rm{,}}{\nu _4}{\rm{)}}$| or so called volume coordinates to build the relations between the two coordinate systems in Fig. 1(a). The inner nodes satisfy |${\nu _1} + {\nu _2} + {\nu _3} + {\nu _4} = 1{\rm{,}}\,\,\,\,{\rm{0}} \le {\nu _{\rm{i}}} \le 1$|⁠. Assume a tetrahedron in the physical domain with the vertexes|${\rm{(}}{x_1}{\rm{,}}{y_1}{\rm{,}}{z_1}{\rm{),}}\,\,{\rm{(}}{x_2}{\rm{,}}{y_2}{\rm{,}}{z_2}{\rm{),}}\,\,{\rm{(}}{x_3}{\rm{,}}{y_3}{\rm{,}}{z_3}{\rm{),}}\,\,{\rm{(}}{x_4}{\rm{,}}{y_4}{\rm{,}}{z_4}{\rm{)}}$|⁠, then any point inside tetrahedron can be expressed as

$$\begin{eqnarray*} \left\{ {\begin{array}{@{}*{1}{c}@{}} {x = {\nu _3}{x_1} + {\nu _4}{x_2} + {\nu _2}{x_3} + {\nu _1}{x_4}}\\ {y = {\nu _3}{y_1} + {\nu _4}{y_2} + {\nu _2}{y_3} + {\nu _1}{y_4}}\\ {z = {\nu _3}{z_1} + {\nu _4}{z_2} + {\nu _2}{z_3} + {\nu _1}{z_4}} \end{array}} \right.. \end{eqnarray*}$$

(19)

Note that the relations in eq. (19) still exists in the reference domain. Thus, for the right-angled tetrahedron in the reference coordinate system with the known coordinates of four vertexes, namely |${\rm{(}} - {\rm{1,}} - {\rm{1,}} - {\rm{1),}}\,\,\,\,{\rm{(1,}} - {\rm{1,}} - {\rm{1),}}\,\,\,\,{\rm{(}} - {\rm{1,1,}} - {\rm{1),}}\,\,\,\,{\rm{(}} - {\rm{1,}} - {\rm{1,1)}}$|⁠, we can obtain from eq. (19) the coordinates in the reference domain, that is

$$\begin{eqnarray*} \left\{ {\begin{array}{@{}*{1}{c}@{}} {\xi = - {\nu _3} + {\nu _4} - {\nu _2} - {\nu _1}}\\ {\eta = - {\nu _3} - {\nu _4} + {\nu _2} - {\nu _1}}\\ {\zeta = - {\nu _3} - {\nu _4} - {\nu _2} + {\nu _1}} \end{array}} \right.. \end{eqnarray*}$$

(20)

From eq. (20), we have

$$\begin{eqnarray*} {\nu _1} =&& \frac{{\zeta + 1}}{2},\quad {\nu _2} = \frac{{\eta + 1}}{2},\quad {\nu _3} = - \frac{{\xi + \eta + \zeta + 1}}{2},\nonumber \\ && {\nu _4} = \frac{{\xi + 1}}{2}. \end{eqnarray*}$$

(21)

Insert eq. (21) into (19), we have

$$\begin{eqnarray*} \left\{ {\begin{array}{@{}*{1}{c}@{}} {x = - \displaystyle\frac{{\xi + \eta + \zeta + 1}}{2}{x_1} + \displaystyle\frac{{\xi + 1}}{2}{x_2} + \displaystyle\frac{{\eta + 1}}{2}{x_3} + \displaystyle\frac{{\zeta + 1}}{2}{x_4}}\\ {y = - \displaystyle\frac{{\xi + \eta + \zeta + 1}}{2}{y_1} + \displaystyle\frac{{\xi + 1}}{2}{y_2} + \displaystyle\frac{{\eta + 1}}{2}{y_3} + \displaystyle\frac{{\zeta + 1}}{2}{y_4}}\\ {z = - \displaystyle\frac{{\xi + \eta + \zeta + 1}}{2}{z_1} + \displaystyle\frac{{\xi + 1}}{2}{z_2} + \displaystyle\frac{{\eta + 1}}{2}{z_3} + \displaystyle\frac{{\zeta + 1}}{2}{z_4}} \end{array}} \right.. \nonumber\\ \end{eqnarray*}$$

(22)

Using eq. (22) we can easily make transformation between two coordinate systems, so that the relations between two domains in Fig. 1 (a) are established. The coordinate transformation in Fig. 1(b) is obtained via eq. (18).

2.4 Construction of interpolation nodes

As reasonable nodes are crucial for the stability of numerical solutions, we will in the following use the specific spectral interpolation node sets to construct the above discussed Vandermonde matrix and to locate the nodal positions in each element. However, a simple location mapping from GLL nodes to a right-angled tetrahedral element is not appropriate, because a lot of interpolation points defined in this way will become asymmetrical clustering to the vertices, resulting in ill-conditioned operators (Hesthaven & Warburton 2008). Karniadakis & Sherwin (2005) have shown that the Lebesgue constant defines the farthest distance between the constructed nodes and the optimal nodes. In other words, one needs to find a set of points that possess greater closeness to GLL. In this paper, we will search a set of nodes within the tetrahedron that are most possibly close to GLL nodes, we will follow the concept of Warburton (2006) to constructed the warp & blend nodes. The procedure is as follow:

• 1-D deformation function

We first define 1-D deformation function |$\omega $|⁠. For this purpose, we assume the n + 1 equidistant nodes within [−1, 1] as the initial position, then we use the Lagrange polynomials to map the equidistant nodes |$x_i^{{\rm{Equi}}}$| to GLL nodes |$x_i^{{\rm{GLL}}}$| and define the deformation function

$$\begin{eqnarray*} \omega {\rm{(}}x{\rm{)}} = \sum\limits_{i = 1}^{n + 1} {{\rm{(}}x_i^{{\rm{GLL}}} - x_i^{{\rm{Equi}}}{\rm{)}}} L_i^{{\rm{Equi}}}, \end{eqnarray*}$$

(23)

where |$L_i^{{\rm{Equi}}} = \mathop \prod \limits_{\scriptstyle j = 1\atop \scriptstyle j \ne i}^{n + 1} ( {\frac{{x - x_i^{{\rm{Equi}}}}}{{x_i^{{\rm{Equi}}} - x_j^{{\rm{Equi}}}}}} )$|⁠.

• 2-D warping and blending function

To obtain the construction functions for 2-D interpolation nodes, we apply the 1-D deformation function in the warp transformation for the three edge directions of a triangle and combine with blend transformation. To construct a set of symmetric nodes, we consider an equilateral triangle with the coordinates of the vertexes |${\rm{(}} - {\rm{1,}}\frac{{ - {\rm{1}}}}{{\sqrt {\rm{3}} }}{\rm{),(1,}}\frac{{ - {\rm{1}}}}{{\sqrt {\rm{3}} }}{\rm{),(0,}}\frac{{\rm{2}}}{{\sqrt {\rm{3}} }}{\rm{)}}$|⁠. Similar to 1-D case, we try to find the equidistributed nodes and denote them by the barycentric coordinates |${\rm{\{ }}{\nu _i}{\rm{\} ,}}\,\,i = 1{\rm{,}}\,\,2{\rm{,}}\,\,3$|⁠. Then, we use the warp and blend transformation to obtain the warp and blend nodes (Warburton 2006).

We first define 2-D warping functions |${{{\bf w}}_i}\,\,{\rm{(}}i = 1{\rm{,}}\,\,2{\rm{,}}\,\,3{\rm{)}}$| along three edge directions as (see Fig. 2)

$$\begin{eqnarray*} \begin{array}{@{}*{1}{l}@{}} {{{{\bf w}}_1}{\rm{(}}{\nu _{\rm{1}}}{\rm{,}}{\nu _{\rm{2}}}{\rm{,}}{\nu _{\rm{3}}}{\rm{)}} = \omega {\rm{(}}{\nu _{\rm{3}}} - {\nu _{\rm{2}}}{\rm{)}}\left( {\begin{array}{@{}*{1}{c}@{}} 1\\ 0 \end{array}} \right)}\\ \\ {{{{\bf w}}_2}{\rm{(}}{\nu _{\rm{1}}}{\rm{,}}{\nu _{\rm{2}}}{\rm{,}}{\nu _{\rm{3}}}{\rm{)}} = \omega {\rm{(}}{\nu _{\rm{1}}} - {\nu _{\rm{3}}}{\rm{)}}\frac{{\rm{1}}}{{\rm{2}}}\left( {\begin{array}{@{}*{1}{c}@{}} { - 1}\\ {\sqrt 3 } \end{array}} \right)}\\ \\ {{{{\bf w}}_3}{\rm{(}}{\nu _{\rm{1}}}{\rm{,}}{\nu _{\rm{2}}}{\rm{,}}{\nu _{\rm{3}}}{\rm{)}} = \omega {\rm{(}}{\nu _{\rm{2}}} - {\nu _{\rm{1}}}{\rm{)}}\frac{{\rm{1}}}{{\rm{2}}}\left( {\begin{array}{@{}*{1}{c}@{}} { - 1}\\ { - \sqrt 3 } \end{array}} \right)} \end{array}, \end{eqnarray*}$$

(24)

and extend the edge warp into the triangle by blending in the edge normal direction and obtain the blending functions |${b_i}{\rm{(}}i = 1{\rm{,}}\,\,2{\rm{,}}\,\,3{\rm{)}}$|⁠, that is

$$\begin{eqnarray*} {b_i}{\rm{(}}{\nu _{\rm{1}}}{\rm{,}}{\nu _{\rm{2}}}{\rm{,}}{\nu _{\rm{3}}}{\rm{)}} = \prod\limits_{\scriptstyle j = 1\atop \scriptstyle j \ne i}^3 {\left( {\frac{{2{\nu _j}}}{{2{\nu _j} + {\nu _i}}}} \right)} . \end{eqnarray*}$$

(25)

Altogether, we can obtain the construction functions for the 2-D interpolation nodes as

$$\begin{eqnarray*} {{\boldsymbol \lambda }}{\rm{(}}{\nu _{\rm{1}}}{\rm{,}}{\nu _{\rm{2}}}{\rm{,}}{\nu _{\rm{3}}}{\rm{)}} =&& {\rm{(}}1 + {{\rm{(}}\alpha {\nu _{\rm{1}}}{\rm{)}}^{\rm{2}}}{\rm{)}}{b_1}{{{\bf w}}_1} + {\rm{(}}1 + {{\rm{(}}\alpha {\nu _{\rm{2}}}{\rm{)}}^2}{\rm{)}}{b_2}{{{\bf w}}_2}\,\, + {\rm{(}}1\nonumber \\ &&+ {{\rm{(}}\alpha {\nu _{\rm{3}}}{\rm{)}}^2}{\rm{)}}{b_3}{{{\bf w}}_3}{\rm{.}} \end{eqnarray*}$$

(26)

From eq. (26), the nodes located at edges coincide the GLL nodes, while the nodes inside the triangle are stretched towards the edges. Besides, we have introduced a parameter |$\alpha $| in eq. (26) to strengthen the nodes further toward the edges. |$\alpha $| can be obtained by optimizing the Lebesgue constant (Warburton 2006). For the nodes at the vertexes, |${\nu _{\rm{1}}} = {\nu _{\rm{3}}} = 0$| or |${\nu _{\rm{1}}} = {\nu _{\rm{2}}} = 0$|⁠, we do not need to process, the interpolation nodes stay at the location of the vertexes.

Figure 2.

Directions of 2-D warping functions.

In the practical construction of the interpolation nodes, we first determine the equidistributed points within an equilateral triangle and use eqs (24) and (25) to calculate w_i and b_i(i = 1,2,3), and then use (26) to calculate the 2-D interpolation nodes.

• 3-D warping and blending function

For 3-D case, we can follow a similar process within an equilateral tetrahedron with the barycentric point as the origin, then the coordinates for the four vertexes are |$( { - {\rm{1,}}\frac{{ - {\rm{1}}}}{{\sqrt {\rm{3}} }}{\rm{,}}\frac{{ - {\rm{1}}}}{{\sqrt {\rm{6}} }}} ){\rm{,}}\,\,( {{\rm{1,}}\frac{{ - {\rm{1}}}}{{\sqrt {\rm{3}} }}{\rm{,}}\frac{{ - {\rm{1}}}}{{\sqrt {\rm{6}} }}} ){\rm{,}}\,\,( {{\rm{0,}}\frac{{\rm{2}}}{{\sqrt {\rm{3}} }}{\rm{,}}\frac{{ - {\rm{1}}}}{{\sqrt {\rm{6}} }}} ){\rm{,}}\,\,( {{\rm{0,0,}}\frac{{\rm{3}}}{{\sqrt {\rm{6}} }}} )$|⁠. By converting the equidistributed nodes into barycentric coordinates |${\rm{\{ }}{\nu _i}{\rm{\} ,}}\,\,i = 1{\rm{,}}\,\,2{\rm{,}}\,\,3{\rm{,}}\,\,4$| and using 3-D warping functions |${{{\bf w}}_i}\,\,{\rm{(}}i = 1{\rm{,}}\,\,2{\rm{,}}\,\,3{\rm{,}}\,\,4{\rm{)}}$| and blending functions |${b_i}{\rm{(}}i = 1{\rm{,}}\,\,2{\rm{,}}\,\,3{\rm{,}}\,\,4{\rm{)}}$|⁠, we can obtain the construction function for the 3-D interpolation nodes as

$$\begin{eqnarray*} {{\bf G}}{\rm{(}}{\nu _{\rm{1}}}{\rm{,}}{\nu _{\rm{2}}}{\rm{,}}{\nu _{\rm{3}}}{\rm{,}}{\nu _{\rm{4}}}{\rm{)}} =&& {\rm{(}}1 + {{{\rm{(}}\alpha {\nu _{\rm{1}}}{\rm{)}}}^2}{\rm{)}}{b_1}{{{\bf w}}_1} + {\rm{(}}1 + {{{\rm{(}}\alpha {\nu _{\rm{2}}}{\rm{)}}}^2}{b_2}{{{\bf w}}_2}\nonumber \\ &&+ \, {\rm{(}}1 + {{{\rm{(}}\alpha {\nu _{\rm{3}}}{\rm{)}}}^2}{\rm{)}}{b_3}{{{\bf w}}_3} + {\rm{(}}1 + {{{\rm{(}}\alpha {\nu _{\rm{4}}}{\rm{)}}}^2}{\rm{)}}{b_4}{{{\bf w}}_4} , \end{eqnarray*}$$

(27)

$$\begin{eqnarray*} {b_i}{\rm{(}}{\nu _{\rm{1}}}{\rm{,}}{\nu _{\rm{2}}}{\rm{,}}{\nu _{\rm{3}}}{\rm{,}}{\nu _{\rm{4}}}{\rm{)}} = \prod\limits_{\scriptstyle j = 1\atop \scriptstyle j \ne i}^4 {\left( {\frac{{2{\nu _j}}}{{2{\nu _j} + {\nu _i}}}} \right)} . \end{eqnarray*}$$

(28)

Here, the warping functions for each face of the tetrahedron need to use the 2-D construction functions. Taking into account that each of the four faces have different orientations, we can write the warping functions as

$$\begin{eqnarray*} \begin{array}{@{}*{1}{l}@{}} {{{{\bf w}}_1}{\rm{(}}{\nu _{\rm{1}}}{\rm{,}}{\nu _{\rm{2}}}{\rm{,}}{\nu _{\rm{3}}}{\rm{,}}{\nu _{\rm{4}}}{\rm{)}} = {\lambda _1}{\rm{(}}{\nu _{\rm{2}}}{\rm{,}}{\nu _{\rm{3}}}{\rm{,}}{\nu _{\rm{4}}}{\rm{)}}{{{\bf t}}_{1{\rm{,}}1}} + {\lambda _2}{\rm{(}}{\nu _{\rm{2}}}{\rm{,}}{\nu _{\rm{3}}}{\rm{,}}{\nu _{\rm{4}}}{\rm{)}}{{{\bf t}}_{1{\rm{,}}2}}}\\ {{{{\bf w}}_2}{\rm{(}}{\nu _{\rm{1}}}{\rm{,}}{\nu _{\rm{2}}}{\rm{,}}{\nu _{\rm{3}}}{\rm{,}}{\nu _{\rm{4}}}{\rm{)}} = {\lambda _1}{\rm{(}}{\nu _{\rm{1}}}{\rm{,}}{\nu _{\rm{3}}}{\rm{,}}{\nu _{\rm{4}}}{\rm{)}}{{{\bf t}}_{2{\rm{,}}1}} + {\lambda _2}{\rm{(}}{\nu _{\rm{1}}}{\rm{,}}{\nu _{\rm{3}}}{\rm{,}}{\nu _{\rm{4}}}{\rm{)}}{{{\bf t}}_{2{\rm{,}}2}}}\\ {{{{\bf w}}_3}{\rm{(}}{\nu _{\rm{1}}}{\rm{,}}{\nu _{\rm{2}}}{\rm{,}}{\nu _{\rm{3}}}{\rm{,}}{\nu _{\rm{4}}}{\rm{)}} = {\lambda _1}{\rm{(}}{\nu _{\rm{1}}}{\rm{,}}{\nu _{\rm{4}}}{\rm{,}}{\nu _{\rm{2}}}{\rm{)}}{{{\bf t}}_{3{\rm{,}}1}} + {\lambda _2}{\rm{(}}{\nu _{\rm{1}}}{\rm{,}}{\nu _{\rm{4}}}{\rm{,}}{\nu _{\rm{2}}}{\rm{)}}{{{\bf t}}_{3{\rm{,}}2}}}\\ {{{{\bf w}}_4}{\rm{(}}{\nu _{\rm{1}}}{\rm{,}}{\nu _{\rm{2}}}{\rm{,}}{\nu _{\rm{3}}}{\rm{,}}{\nu _{\rm{4}}}{\rm{)}} = {\lambda _1}{\rm{(}}{\nu _{\rm{1}}}{\rm{,}}{\nu _{\rm{3}}}{\rm{,}}{\nu _{\rm{2}}}{\rm{)}}{{{\bf t}}_{4{\rm{,}}1}} + {\lambda _2}{\rm{(}}{\nu _{\rm{1}}}{\rm{,}}{\nu _{\rm{3}}}{\rm{,}}{\nu _{\rm{2}}}{\rm{)}}{{{\bf t}}_{4{\rm{,}}2}}} \end{array}, \end{eqnarray*}$$

(29)

where |${\lambda _1}$| and |${\lambda _2}$| are the two components of the 2-D construction function |${{\bf \lambda }}$| in eq. (26), |${{{\bf t}}_{i,j}}$| are two vectors forming orthogonal axes in the plane of each face (see Fig. 3).

Figure 3.

Construction of 3-D interpolation nodes. |${{{\bf t}}_{i,j}}$| denote the vectors forming orthogonal axes in the plane of each face.

From the above discussion, we can summarize the process of determining the interpolation nodes: a series of equidistributed points are stretched towards certain directions at certain distances, so that they are put into the optimal locations closing to the vertexes and edges. Once these nodes are obtained in an equilateral tetrahedron, we can map them into reference domain of a right-angled tetrahedron and obtain the interpolation nodes. Fig. 4 shows the distributions of the nodes for 2nd, 5th and 8th order.

Figure 4.

Distribution of interpolation nodes in a right-angled tetrahedron for the order of n = 2, 5 and 8. The red points are nodes at vertexes, the yellow ones are for edge nodes, the green ones are for surface nodes, while the blue points are for interior nodes.(refer to Warburton 2006)

For the nth interpolation, we have Nt interpolation nodes, including vertex nodes, edge nodes, face nodes, and interior nodes. When n = 1, we have only vertex nodes, with four nodes for each tetrahedron. When |$n \ge 2$|⁠, there exist edge nodes, the number of nodes at each edge is n–1 (exclude the vertex nodes). When |$n \ge 3$|⁠, there exist face nodes, the number of the face nodes is (n + 1)(n + 2)/2–3n (exclude the edge and vertex nodes). When|$n \ge 4$|⁠, there exist interior nodes, the number of nodes is (n – 1)(n – 2)(n– 3)/6. Table 1 gives the number of totalnodes and interior nodes.

Table 1.

Number of total nodes Nt and interior nodes Ni for nth-order interpolation.

Order/n	1	2	3	4	5	6	7	8
Nt	4	10	20	35	56	84	120	165
Ni	0	0	0	1	4	10	20	35

Finally, we can get the complete basis functions |${\varphi _i}$| and interpolation nodes for our modelling. Fig. 5 shows 1-D basis functions of 4th order. It consists of five functions and their intersections are GLL interpolation nodes. Fig. 6 shows 2-D basis functions of the 4th order. It consists of 15 functions, with the locations of the maximum value corresponding to the interpolation nodes. For the 3-D case, there are 35 functions in total that are connected with variables |${\rm{(}}\xi {\rm{,}}\eta {\rm{,}}\zeta {\rm{)}}$|⁠. This results in the difficulty in displaying them, so we do not show them here.

Figure 5.

1-D basis functions of 4th order.

Figure 6.

2-D basis functions of 4th order.(refer to Hesthaven & Warburton 2008)

2.5 Numerical integration

After obtaining the basis functions and the interpolation nodes, we can calculate the entries of the coefficient matrix in eqs (6)–(8) for each element. For that purpose, we need to numerically calculate the triple integration of the gradient of the basis functions. Refer to eq. (7), we have

$$\begin{eqnarray*} \begin{array}{@{}*{1}{l}@{}} {\iiint\limits_{{\Omega _e}} {\sigma \nabla {\varphi _i} \cdot \nabla {\varphi _j}}\,\,{\rm{d}}\Omega }\\ { = \iiint\limits_{{\Omega _s}} {\sigma \nabla {\varphi _i} \cdot \nabla {\varphi _j}}\left| {{{{\bf J}}_{3D}}} \right|\,\,{\rm{d}}\xi {\rm{d}}\eta {\rm{d}}\zeta }\\ { = \sigma \sum\limits_{k{\rm{,}}l = 1}^{Nt} {{{({{{\bf V}}^{ - 1}})}_{ki}}{{({{{\bf V}}^{ - 1}})}_{lj}}\iiint\limits_{{\Omega _s}} {\nabla {\psi _k} \cdot \nabla {\psi _l}\left| {{{{\bf J}}_{3D}}} \right|{\rm{d}}\xi {\rm{d}}\eta {\rm{d}}\zeta }} }\\ { \approx \sigma \sum\limits_{k{\rm{,}}l = 1}^{Nt} {{{({{{\bf V}}^{ - 1}})}_{ki}}{{({{{\bf V}}^{ - 1}})}_{lj}}\sum\limits_{m = 1}^{Nf} {{w_m}{{\left[ {\left. {\nabla {\psi _k} \cdot \nabla {\psi _l}} \right]} \right.}_m}} } {{\left| {{{{\bf J}}_{3D}}} \right|}_m}}\\ { = \sigma \sum\limits_{k{\rm{,}}l = 1}^{Nt} {{{({{{\bf V}}^{ - 1}})}_{ki}}{{({{{\bf V}}^{ - 1}})}_{lj}}} \sum\limits_{m = 1}^{Nf} {{w_m}\left| {{{{\bf J}}_{3D}}{\rm{(}}{\xi _m}{\rm{,}}{\eta _m}{\rm{,}}{\zeta _m}{\rm{)}}} \right| } }\\ \times {\left( {\left[ {\displaystyle\frac{{\partial {\psi _k}{\rm{(}}{\xi _m}{\rm{,}}{\eta _m}{\rm{,}}{\zeta _m}{\rm{)}}}}{{\partial \xi }}\displaystyle\frac{{\partial \xi }}{{\partial x}} + \displaystyle\frac{{\partial {\psi _k}{\rm{(}}{\xi _m}{\rm{,}}{\eta _m}{\rm{,}}{\zeta _m}{\rm{)}}}}{{\partial \eta }}\displaystyle\frac{{\partial \eta }}{{\partial x}} + \displaystyle\frac{{\partial {\psi _k}{\rm{(}}{\xi _m}{\rm{,}}{\eta _m}{\rm{,}}{\zeta _m}{\rm{)}}}}{{\partial \zeta }}\displaystyle\frac{{\partial \zeta }}{{\partial x}}} \right]} \right. }\\ \,\,\times {\left[ {\displaystyle\frac{{\partial {\psi _l}{\rm{(}}{\xi _m}{\rm{,}}{\eta _m}{\rm{,}}{\zeta _m}{\rm{)}}}}{{\partial \xi }}\displaystyle\frac{{\partial \xi }}{{\partial x}} + \displaystyle\frac{{\partial {\psi _l}{\rm{(}}{\xi _m}{\rm{,}}{\eta _m}{\rm{,}}{\zeta _m}{\rm{)}}}}{{\partial \eta }}\displaystyle\frac{{\partial \eta }}{{\partial x}} + \displaystyle\frac{{\partial {\psi _l}{\rm{(}}{\xi _m}{\rm{,}}{\eta _m}{\rm{,}}{\zeta _m}{\rm{)}}}}{{\partial \varsigma }}\displaystyle\frac{{\partial \zeta }}{{\partial x}}} \right] }\\ \,\, + {\left[ {\displaystyle\frac{{\partial {\psi _k}{\rm{(}}{\xi _m}{\rm{,}}{\eta _m}{\rm{,}}{\zeta _m}{\rm{)}}}}{{\partial \xi }}\displaystyle\frac{{\partial \xi }}{{\partial y}} + \displaystyle\frac{{\partial {\psi _k}{\rm{(}}{\xi _m}{\rm{,}}{\eta _m}{\rm{,}}{\zeta _m}{\rm{)}}}}{{\partial \eta }}\displaystyle\frac{{\partial \eta }}{{\partial y}} + \displaystyle\frac{{\partial {\psi _k}{\rm{(}}{\xi _m}{\rm{,}}{\eta _m}{\rm{,}}{\zeta _m}{\rm{)}}}}{{\partial \zeta }}\displaystyle\frac{{\partial \zeta }}{{\partial y}}} \right] }\\ \,\,\times {\left[ {\displaystyle\frac{{\partial {\psi _l}{\rm{(}}{\xi _m}{\rm{,}}{\eta _m}{\rm{,}}{\zeta _m}{\rm{)}}}}{{\partial \xi }}\displaystyle\frac{{\partial \xi }}{{\partial y}} + \displaystyle\frac{{\partial {\psi _l}{\rm{(}}{\xi _m}{\rm{,}}{\eta _m}{\rm{,}}{\zeta _m}{\rm{)}}}}{{\partial \eta }}\displaystyle\frac{{\partial \eta }}{{\partial y}} + \displaystyle\frac{{\partial {\psi _l}{\rm{(}}{\xi _m}{\rm{,}}{\eta _m}{\rm{,}}{\zeta _m}{\rm{)}}}}{{\partial \zeta }}\displaystyle\frac{{\partial \zeta }}{{\partial y}}} \right] }\\ \,\,+ {\left[ {\displaystyle\frac{{\partial {\psi _k}{\rm{(}}{\xi _m}{\rm{,}}{\eta _m}{\rm{,}}{\zeta _m}{\rm{)}}}}{{\partial \xi }}\displaystyle\frac{{\partial \xi }}{{\partial z}} + \displaystyle\frac{{\partial {\psi _k}{\rm{(}}{\xi _m}{\rm{,}}{\eta _m}{\rm{,}}{\zeta _m}{\rm{)}}}}{{\partial \eta }}\displaystyle\frac{{\partial \eta }}{{\partial z}} + \displaystyle\frac{{\partial {\psi _k}{\rm{(}}{\xi _m}{\rm{,}}{\eta _m}{\rm{,}}{\zeta _m}{\rm{)}}}}{{\partial \zeta }}\displaystyle\frac{{\partial \zeta }}{{\partial z}}} \right] }\\ \,\,\times {\left. {\left[ {\displaystyle\frac{{\partial {\psi _l}{\rm{(}}{\xi _m}{\rm{,}}{\eta _m}{\rm{,}}{\zeta _m}{\rm{)}}}}{{\partial \xi }}\displaystyle\frac{{\partial \xi }}{{\partial z}} + \displaystyle\frac{{\partial {\psi _l}{\rm{(}}{\xi _m}{\rm{,}}{\eta _m}{\rm{,}}{\zeta _m}{\rm{)}}}}{{\partial \eta }}\displaystyle\frac{{\partial \eta }}{{\partial z}} + \displaystyle\frac{{\partial {\psi _l}{\rm{(}}{\xi _m}{\rm{,}}{\eta _m}{\rm{,}}{\zeta _m}{\rm{)}}}}{{\partial \zeta }}\displaystyle\frac{{\partial \zeta }}{{\partial z}}} \right]} \right)}, \end{array}\nonumber \\ \end{eqnarray*}$$

(30)

where Nf is the number of the numerical integration points,|${w_m}$| is the weight for the mth point. The interpolation bases ψ_i (⁠|$1 \le i \le Nt$|⁠) are given in eq. (17). We use the Gauss-type integration points given by Zhang et al. (2009) for the tetrahedrons. Table 2 gives the number of integration nodes for tetrahedral elements.

Table 2.

Number of integration nodes for tetrahedral elements.

Order/n	1	2	3	4	5	6	7	8
Nf	1	4	8	14	14	24	35	46

For the straight-sided tetrahedrons used in this paper, the Jacobian determinant |$| {{{{\bf J}}_{\rm{3D}}}} |$| in eq. (30) is a constant, which can be written as

$$\begin{eqnarray*} \left| {{{{\bf J}}_{\rm{3D}}}} \right| = \left| {\begin{array}{@{}*{3}{c}@{}} {\displaystyle\frac{{\partial x}}{{\partial \xi }}}&{\displaystyle\frac{{\partial x}}{{\partial \eta }}}&{\displaystyle\frac{{\partial x}}{{\partial \zeta }}}\\ {\displaystyle\frac{{\partial y}}{{\partial \xi }}}&{\displaystyle\frac{{\partial y}}{{\partial \eta }}}&{\displaystyle\frac{{\partial y}}{{\partial \zeta }}}\\ {\displaystyle\frac{{\partial z}}{{\partial \xi }}}&{\displaystyle\frac{{\partial z}}{{\partial \eta }}}&{\displaystyle\frac{{\partial z}}{{\partial \zeta }}} \end{array}} \right| = \left| {\begin{array}{@{}*{3}{c}@{}} {\displaystyle\frac{{{x_2} - {x_1}}}{2}}&{\displaystyle\frac{{{x_3} - {x_1}}}{2}}&{\displaystyle\frac{{{x_4} - {x_1}}}{2}}\\ {\displaystyle\frac{{{y_2} - {y_1}}}{2}}&{\displaystyle\frac{{{y_3} - {y_1}}}{2}}&{\displaystyle\frac{{{y_4} - {y_1}}}{2}}\\ {\displaystyle\frac{{{z_2} - {z_1}}}{2}}&{\displaystyle\frac{{{z_3} - {z_1}}}{2}}&{\displaystyle\frac{{{z_4} - {z_1}}}{2}} \end{array}} \right|. \end{eqnarray*}$$

(31)

It is a special case of transfinite blending functions in the general case of curved edges. Through collapse transform, the partial derivatives of the basis functions with respect to the reference coordinates can be written as

$$\begin{eqnarray*} \left\{ {\begin{array}{@{}*{1}{l}@{}} {\displaystyle\frac{{\partial {\psi _i}}}{{\partial \xi }} = \displaystyle\frac{{\partial {\psi _i}}}{{\partial r}}\displaystyle\frac{{\partial r}}{{\partial \xi }}}\\ {\displaystyle\frac{{\partial {\psi _i}}}{{\partial \eta }} = \displaystyle\frac{{\partial {\psi _i}}}{{\partial r}}\displaystyle\frac{{\partial r}}{{\partial \eta }} + \displaystyle\frac{{\partial {\psi _i}}}{{\partial s}}\displaystyle\frac{{\partial s}}{{\partial \eta }}}\\ {\displaystyle\frac{{\partial {\psi _i}}}{{\partial \zeta }} = \displaystyle\frac{{\partial {\psi _i}}}{{\partial r}}\displaystyle\frac{{\partial r}}{{\partial \zeta }} + \displaystyle\frac{{\partial {\psi _i}}}{{\partial s}}\displaystyle\frac{{\partial s}}{{\partial \zeta }} + \displaystyle\frac{{\partial {\psi _i}}}{{\partial t}}} \end{array}} \right., \end{eqnarray*}$$

(32)

while the terms |$\frac{{\partial {\rm{(}}\xi {\rm{,}}\eta {\rm{,}}\zeta {\rm{)}}}}{{\partial {\rm{(}}x{\rm{,}}y{\rm{,}}z{\rm{)}}}}$| in eq. (30) can be calculated through |$| {{{{\bf J}}_{\rm{3D}}}} |$|⁠. Moreover, the second term in eq. (7) can be written as

$$\begin{eqnarray*} \begin{array}{@{}*{1}{l}@{}} {\iint\limits_{{\Gamma _e}} {\sigma \displaystyle\frac{{{\rm{cos(}}{{\bf r}} - {{{\bf r}}_0}{\rm{,}}{{\bf n}}{\rm{)}}}}{{\left| {{{\bf r}} - {{{\bf r}}_0}} \right|}}{\varphi _i}{\varphi _j}{\rm{d}}\Gamma }}\\ { = \sigma \displaystyle\frac{{{\rm{cos(}}{{\bf r}} - {{{\bf r}}_0}{\rm{,}}{{\bf n}}{\rm{)}}}}{{\left| {{{\bf r}} - {{{\bf r}}_0}} \right|}}\sum\limits_{k,l = 1}^{Ntt} {{{({{{\bf V}}^{ - 1}})}_{ki}}{{({{{\bf V}}^{ - 1}})}_{lj}}\iint\limits_{{\Gamma _s}} {{\psi _k}{\psi _l}\left| {{{{\bf J}}_{\rm{2D}}}} \right|{\rm{d}}\xi {\rm{d}}\eta }} }\\ { = \sigma \displaystyle\frac{{{\rm{cos(}}{{\bf r}} - {{{\bf r}}_0}{\rm{,}}{{\bf n}}{\rm{)}}}}{{\left| {{{\bf r}} - {{{\bf r}}_0}} \right|}}\sum\limits_{k,l = 1}^{Ntt} {{{({{{\bf V}}^{ - 1}})}_{ki}}{{({{{\bf V}}^{ - 1}})}_{lj}}\left| {{{{\bf J}}_{\rm{2D}}}} \right|{W_{kl}}} } \end{array}, \end{eqnarray*}$$

(33)

where the PKD polynomial are given by

$$\begin{eqnarray*} {\psi _N} =&& \sqrt {\frac{2}{{{\gamma _{i^{\prime}j^{\prime}}}}}} P_{i^{\prime}}^{{\rm{(}}0{\rm{,}}0{\rm{)}}}{\rm{(}}r{\rm{)}}P_{j^{\prime}}^{{\rm{(}}2i^{\prime} + 1{\rm{,}}0{\rm{)}}}{\rm{(}}s{\rm{)(}}1 - s{{\rm{)}}^{i^{\prime}}},\nonumber \\ r =&& 2\frac{{1 + \xi }}{{1 - \eta }} - 1{\rm{,}}\,\,\,\,s = \eta \end{eqnarray*}$$

(34)

$$\begin{eqnarray*} N =&& j^{\prime} + {\rm{(}}n + 1{\rm{)}}i^{\prime} + 1 - \frac{{i^{\prime}}}{2}{\rm{(}}i^{\prime} - 1{\rm{)}},\nonumber \\ &&{\rm{(}}i^{\prime},j^{\prime}{\rm{)}} \ge 0,i^{\prime} + j^{\prime} \le Ntt. \end{eqnarray*}$$

(35)

where |${\gamma _{i^{\prime}j^{\prime}}} = \frac{2}{{2i^{\prime} + 1}}\frac{{{2^{2i^{\prime} + 2}}}}{{2i^{\prime} + 2j^{\prime} + 2}}$|⁠. In eqs (33) and (35), |$Ntt = \frac{{{\rm{(}}n + 1{\rm{)(}}n + 2{\rm{)}}}}{2}$| is the number of interpolation points, |${\Gamma _s}$| is the reference domain mapped from |${\Gamma _e}$|⁠, |$| {{{{\bf J}}_{2D}}} |$| is the 2-D Jacobian determinant that can be calculated similarly to eq. (31), while W_kl can be calculated by

$$\begin{eqnarray*} {W_{kl}} =&& \iint\limits_{{\Gamma _s}} {{\psi _k}{\psi _l}{\rm{d}}\xi {\rm{d}}\eta }\nonumber \\ && { = \int_{ - 1}^{{1}}{{\int_{ - 1}^{{1}}{{{\psi _k}}}}}{\psi _l}\left| {{{{\bf J}}_H}} \right|{\rm d}r{\rm d}s}\nonumber \\ && = \int_{ - 1}^{{1}}\int_{ - 1}^{{1}}{{\displaystyle\frac{2}{{\sqrt {{\gamma _{\ell m}}{\gamma _{pq}}} }}P_\ell ^{{\rm{(}}0,0{\rm{)}}}{\rm{(}}r{\rm{)}}P_m^{{\rm{(}}2\ell + 1{\rm{)}}}{\rm{(}}s{\rm{)(}}1 - s{{\rm{)}}^\ell }P_p^{{\rm{(}}0,0{\rm{)}}}{\rm{(}}r{\rm{)}} }}\nonumber \\ && \quad \times {{P_q^{{\rm{(}}2p + 1{\rm{)}}}{\rm{(}}s{\rm{)(}}1 - s{{\rm{)}}^p}}}\left| {{{{\bf J}}_H}} \right|{\rm d}r{\rm d}s\nonumber \\ && { = \left\{ {\begin{array}{@{}*{2}{l}@{}} {{\rm{1}},}&{\ell = p\,\,{\rm{and}}\,\,m = q;} \\ {0,}&{{\rm{otherwise}}} \end{array}} \right.} . \end{eqnarray*}$$

(36)

In eq. (36), |${\bf J}_{H}\,=\,| _{\frac{\partial \xi}{\partial s}\,\,\,\,\frac{\partial \eta}{\partial s}}^{\frac{\partial \xi}{\partial r}\,\,\,\,\frac{\partial \eta}{\partial r}}|\,=\,\frac{1-s}{2},\,\,{\ell}$|⁠, m, p, q and k, l satisfy the similar relations in eq. (35). In this paper, we use the open source code TetGen (Hang 2015) to discretize the model domain into tetrahedral elements. After obtaining the elements in eqs (7) and (8), we assemble them together to form the equations system of eq. (9) and solved it by the open source code MUMPS (Amestoy et al. 2000).

3 NUMERICAL EXPERIMENTS

3.1 Accuracy verification

To verify the accuracy of our TSEM, we design a 2-layer earth model. The first layer has a resistivity of 1 |$\Omega \cdot {\rm{m}}$| and a thickness of 10 m, the basement has a resistivity of 100 |$\Omega \cdot {\rm{m}}$|⁠. The point current source of 1A is located at the origin. We calculate the potentials at interval of 1 m within an area of |$100\,{\rm{m}} \times 100\,{\rm{m}}$|⁠. To demonstrate the advantages of our high-order algorithm that can improve the modelling accuracy by raising the order of basis functions in the interpolation, we assume a very coarse initial mesh. We only refine the mesh around the source for we use total field method here. Fig. 7 shows the mesh for the research area at the earth surface.

Figure 7.

Mesh for the research area at the earth surface.

Fig. 8 shows the results in contours using our TSEM for different orders. From the figure, it is seen that when n = 1, the contour becomes deformed. With increasing order, the contours become more and more smooth and circular, implying that the TSEM delivers more and more accurate results. Moreover, we compare our modelling results with the analytical solutions given by Wait (1990). From Fig. 9, it is seen that with the increasing order the relative errors decrease. The maximum error for n = 1 is 451.5 per cent, while it is 2.3 per cent for n = 4. This demonstrates that our TSEM method can model the electrical potential effectively. With higher order of TSEM, the modelling accuracy can be very much improved.

Figure 8

Electrical potentials for a point source at the surface of a 2-layer earth. The TSEM order is changing. (a) n=1 (b) n=2 (c) n=3 (d) n=4. Note that with increasing order, the contours become more and more smooth and circular that is close to the true distribution of potential for a point source

Figure 9.

Errors E in percentage related to the 1-D analytical solutions for the model in Fig. 7 for different order of TSEM. (a) n = 1 (b) n = 2 (c) n = 3 (d) n = 4. With increasing order, the relative errors reduce sharply.

Meanwhile, we have chosen a profile along the y-axis to our model results and the relative errors against the analytical solutions are given in Fig. 10. From the figure, we can draw the similar conclusion that with increasing order in our TSEM the accuracy is vastly improved.

Figure 10

Modelling results of our TSEM in comparison to the analytical solutions along the profile z = 0 m, x = 0 m for the layered model given in Fig. 7. (a) Electrical potential; (b) relative errors.

3.2 Efficiency analysis

To illustrate another advantage of efficiency of our TSEM, we take the model designed by Ren & Tang (2010) and compare the results from our TSEM with their adaptive finite element method (AFEM) results. Refer to Fig. 11, the model domain has a dimension of |$2000\,\,{\rm{m}} \times 2000\,\,{\rm{m}} \times 1200\,\,{\rm{m}}$|⁠. We assume that there exists a vertical contact that divides the half-space into two parts, with resistivities of 10 and 100 |$\Omega \cdot {\rm{m}}$|⁠, respectively. A cubic body with an edge length of 2 m and a resistivity of 10 |$\Omega \cdot {\rm{m}}$| is buried in the resistive section. The top depth of the body is 2 m. We choose the Schlumberger array for our modelling, the potential electrodes MN are, respectively, located at (0, −3.4 m, 0) and (0, −2.6 m, 0). The current electrodes AB are symmetric about the center of MN, AB/2 ranges from 0.5 to 100 m. The current is 1A. After obtaining the potential using TSEM, we use the definition given by Telford et al. (1990) to calculate the apparent resistivity.

Figure 11.

A cube buried in a half-space with a vertical contact (refer to Ren & Tang 2010).

To make the evaluation on our TSEM, we take the semi-analytical results from boundary integral (BI) method as the standard solution (Hvozdara & Kaikkonen 1994) and compare our results with those calculated by the AFEM with the open-source code at https://software.seg.org/2010/0002/index.html. Fig. 12(a) shows the comparison between the apparent resistivities calculated by FE method from adaptively refined grids (Ap_0, Ap_1, Ap_2, Ap_3) and those calculated from our TSEM at the coarse grids Ap_0 for different orders (n = 1, n = 2, n = 3, n = 4), while Fig. 12(b) shows the relative errors against the BI solutions. One sees that the FE method improves the modelling accuracy by refining the grids, while our TSEM can do that by enhancing the order. The precision of our results at the second order for the coarse Ap_0 grids has surpassed that of the FE results for Ap_3 (refine grids three times). The high order TSEM can achieve even higher accuracy.

Figure 12.

(a) Comparison of apparent resistivity between our TSEM and the AFEM by Ren & Tang (2010) for the model in Fig. 11; (b) relative errors.

Fig. 13 shows the comparison of apparent resistivity for different grids/orders of TSEM for the model in Fig. 11 and the relative errors against the semi-analytical solutions. It can be observed that in contrast to traditional FE method that can only improve the modelling accuracy by refining grids, our TSEM can improve the accuracy by both refining grids and enhancing the order. In fact, from Fig. 13(b), one sees that the accuracy for n = 4 at coarse grids Ap_0 is comparable to that for n = 3 at fine grids Ap_1.

Figure 13.

(a) Comparison of apparent resistivity for different grids/orders of TSEM for the model in Fig. 11; (b) relative errors. In the figure, the original modelling for n = 1, 2, 3, 4 takes the Ap_0 grids, while the new modelling for n_refine = 1, 2, 3, 4 takes the refined Ap_1 grids.

Table 3 gives comparison of our TSEM with AFEM results in terms of grids, DOFs and relative errors against the semi-analytical solutions from Hvozdara & Kaikkonen (1994) for the model in Fig. 11. It is seen that the DOFs for AFEM at Ap_3 (three-time refinements) is 114 310, for 4th TSEM, however, the DOFs is only 46 413, roughly 40 per cent of the DOFs for AFEM, and the accuracy of our TSEM is 17–18 times higher than AFEM. In addition, the DOFs of AFEM increase very fast even when the adaptive technique is used.

Table 3.

Comparison of TSEM with AFEM in terms of grids and DOFs.

Method	AFEM				TSEM
	Ap_0	Ap_1	Ap_2	Ap_3	n = 1	n = 2	n = 3	n = 4
No. of grids	3215	46 503	209 189	682 514	3215	3215	3215	3215
DOFs	927	8140	35 193	114 310	927	6341	20 188	46 413
Max error (per cent)	27.54	8.86	6.32	3.97	27.54	2.03	0.71	0.22
Average error (per cent)	12.34	5.12	3.60	2.07	12.34	0.92	0.21	0.12

Summarizing the above discussion, we find that the TSEM of this paper has the advantage of improving the modelling accuracy by enhancing the order and refining the grids. In practice, we can take a dual strategy for high-accuracy numerical modelling, namely we first enhance the order on coarse grids to improve the accuracy until further enhancement of the order does not obviously improve the accuracy. At this time, we can refine the grids and further improve the accuracy.

3.3 Responses to typical models

In this section, we design a topographic earth model with two spheres embedded. Fig. 14 gives the model parameters. The 2-D hill and valley have an expansion of 50 m, the peak of the hill and the depth of the valley are both 10 m, the peak and valley bottom are, respectively, located at x = −25 m, x = 25 m. Two spheres with a radius of 5 m are buried 5 m deep under the hill and the valley, respectively. The resistivities are |${\rho _1} = 1\,\Omega \cdot {\rm{m}}$| to the left and |${\rho _2} = 10\,000\,\Omega \cdot {\rm{m}}$| to the right sphere. The earth half-space has a resistivity of 100 |$\Omega \cdot {\rm{m}}$|⁠. The survey line is along the x-axis. Figs 14(a)–(c) display the tetrahedral grids for the topographic earth with abnormal bodies. Fig. 14(d) displays the distribution of the interpolation nodes at the cross section of a quarter of sphere. We use the 3rd order of TSEM to calculate the potential and the apparent resistivities for gradient array and pole–dipole array (forward and reverse measurement, Sheriff 2002).

Figure 14.

A topographic earth with two spherical bodies embedded. The white dashed box denotes an area to be enlarged in the next figure. (a) Tetrahedral grids; (b) enlarged model in profile; (c) grids of the resistive sphere; (d) nodes for 3rd order of TSEM.

Fig. 15 shows the apparent resistivities for both configurations. From Fig. 15(a), one sees that the apparent resistivity for pure topographic earth has a mirror relation with the topography. With the introduction of the abnormal bodies, the anomalies for both spherical bodies are superposed by the topographic response. For the forward and reverse pole–dipole array, the responses |$\rho _{_{\rm{S}}}^{\rm{A}}$| and |$\rho _{_{\rm{S}}}^{\rm{B}}$| for the pure topography demonstrate the standard form for a hill (the left-hand side) and valley (the right-hand side) topography, with the cross points located over the peak of the hill or valley bottom. The addition of two spheres results in the deformation of the anomalies. Whereas the anomaly of the right-hand side resistive body is strengthened by the valley response, the left conductive body is overwhelmed by the response of the hill topography.

Figure 15.

Apparent resistivities for (a) gradient array and (b) forward and reverse pole–dipole array for the model in Fig. 14. The TSEM of 3rd order has been used in the modelling.

From the above modelling results we observed that the topographic effect for dc resistivity method is very serious. To remove this effect, we can use the topographic correction given by Fox et al. (1980), that is

$$\begin{eqnarray*} \rho _s^G = {\rho _s}/{\rm{(}}\rho _s^T/\rho {\rm{)}}, \end{eqnarray*}$$

(37)

where |${\rho _s}$| is the measured apparent resistivity, |$\rho _s^T$| is the apparent resistivity for the pure topography of resistivity |$\rho $|⁠, |$\rho _s^G$| is the apparent resistivity after topographic correction. Fig. 16 shows the apparent resistivities for two configurations after topographic correction to the data in Fig. 15. It is seen that the anomalies for the conductive and resistive spheres are clearly observed. The central locations of the abnormal bodies can also be distinguished.

Figure 16.

Apparent resistivities after topographic corrections. The anomalies for the conductive and resistive spherical bodies are clearly observed. The central locations of the abnormal bodies can also be distinguished.

4 CONCLUSIONS

By incorporating the SEM with unstructured grids, we have successfully developed an algorithm for dc resistivity modelling. In contrast to traditional SEM based on structured grids, our TSEM truly combines the merits of quick convergence of the spectral method and the unstructured finite element method for complex models. The modelling accuracy can be improved via a dual-track procedure of increasing the order of the interpolation and refining grids. Numerical experiments demonstrated that we can obtain accurate model results iteratively by increasing the TSEM order first until no further improvement is achieved, and then we refine the grids and repeat the process, so that we can quickly achieve the expected accuracy. By revisiting the typical topographic earth models, we confirmed again the seriousness of topographic effect and importance of topographic correction. However, we only test our TSEM based on node interpolation on dc resistivity method. Our future research will be focusing on EM modelling based on vector interpolation.

ACKNOWLEDGEMENTS

The paper is financially supported by the National Natural Science Foundation of China (41774125, 41530320), the Strategic Priority Research Program of the Chinese Academy of Sciences (XDA14020102), the Key National Research Project of China (2017YFC0601900, 2016YFC0303100) and the S&T Program of Beijing (Z181100005718001).